Hi Eike,
Eike Rathke schrieb:
Hi Regina,
On Sunday, 2009-06-28 20:20:17 +0200, Regina Henschel wrote:
Searching around I come across an iterative solution from Peter
Deuflhard [1] (German). My tests as Basic macro results in an accuracy
of at least 12 digits. But it needs approximately order+x*1.1+30
iterations.
[1] http://www.mathematik.uni-dortmund.de/ieem/BzMU/BzMU2007/Deuflhard.pdf
So the question is, how much iterations are acceptable?
Depends of course on the cost per iteration. Hard to say. How much time
is a user willing to wait for an accurate result? Which still doesn't
say anything about a series of values to be calculated. But, since
a non-accurate solution wouldn't help much, I guess having to wait
a little isn't too much.. so I just throw a number and say 10000. Does
that help anything?
With 10000 iterations you get convergence for x<6000 - might be x<10000
- depending on order N. When I calculate the asymptotic formula with
MuPad in high precision, then I get only 2 to 3 (estimated in the mean)
accurate digits even for x>200*N^2. From point of accuracy this
asymptotic formula is unusable.
For example BESSELJ(3000;200) results in -7,79235815417491E-003 with the
current implementation with asymptotic formula, but the true value is
-1.186524260848996e-2. Gnumeric results -1.186524260848841E-002, my
Basic version of Deuflhard's algorithm results -1,18652426084901E-002.
But calculating this needs 3147 iterations. For order 200 you would need
x>520000 to result in at least 1 digit accuracy with the asymptotic
formula.
Second problem: The current implementation restricts the order to
integer values. In mathematic the order need not to be integer. ODF1.2
also allows non-integer order. Are there plans to implement a solution
in OOo with non-integer orders? I personally have no idea how to do it.
No plans. Is there any practical benefit in calculating non-integer
orders?
I don't know. The Bessel-functions are used in technical tasks. I also
don't know, which ranges are really necessary for x and for order N.
Does any spreadsheet application do it?
Excel rounds down to integer. The comment in ODF spec 16.15.11 says that
Gnumeric allow fractional N, but my Gnumeric 1.6.3 rounds down too. In
the comment Kspread is said to allow fractional N, but I haven't got
that program.
So the question is, shall I implement Deuflhard's algorithm nevertheless?
Perhaps warn the user in the help, that large input values results in
long lasting calculations?
What is "large" and what is "long lasting"?
I would have to implement it in C++ to see how long it really lasts. The
Basic macro lasts several seconds for x>10000. The effort is linear in x.
Or return an error, if neither asymptotic
formula results at least 2 digits (How tell the user, that the result is
not accurate?)
A few functions set errNoConvergence if a function can't determine
a result.
That would be better than returning a result, which has no correct
digit, as the current solution does.
nor Deuflhard's algorithm calculates in reasonable time
(where to cut?)?
"reasonable time" depends on user perspective and intention. I don't
think there is a general reasonable timeout.
So I would prefer to implement an iterative solution and drop the
inaccurate asymptotic part. The algorithm should converge, but I
hesitate to use a while-loop without any cut. I need to test, for which
x I can get a result in 1 second. Of cause the help should warn the user
about high x-values. The asymptotic formula can be mentioned in the help
in Wiki. It is not long and the user can write it directly as
spreadsheet formula, if he needs huge x-values. That way he is aware of
the accuracy problem.
Btw, newer g++ compilers (>= 4.3) have some bessel functions built-in,
as part of the TR1 extensions, see
http://en.wikipedia.org/wiki/C%2B%2B_Technical_Report_1, it appears that
MSVC does not have them though.
The MSVC Express edition has not any of that functions.
Also the Boost library provides bessel
functions, I don't know off-hand which library version introduced them
and if it is in the one OOo uses. Maybe time to finally upgrade ;-)
You are right to remember me to Boost, I have overlooked it. I had a
look at Boost now. The Bessel functions (and a lot of others) are
included in library "Math Toolkit" since version 1.35.0. I see version
1.34.0 in OOo. There is a description in
http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/bessel.html.
Having an actual version of the Boost library would solve a lot of other
accuracy problems too. I would prefer upgrading Boost library. When can
the upgrade be done? Using it then would mean to rewrite most of
interpr3, interpr6 and analysis part of scaddin, which might take a
while. I guess rewrite can be done step by step. In the long run we need
no longer to struggle with accuracy of special functions and get free
time for fixing other issues.
Nevertheless, I can provide a patch with the Deuflhard algorithm likely
in time for OOo3.2 when give me a "go".
BTW, I still need a decision on the UI name of the left tailed
F-distribution to finish that work. I have FDISTL in mind.
kind regards
Regina
---------------------------------------------------------------------
To unsubscribe, e-mail: [email protected]
For additional commands, e-mail: [email protected]
