Imprecision of Airy Ai

[Grégory Vanuxem]

Hello,

In FriCAS, airyAi has a problem of precision. Here is the output using
DoubleFloat, MathLink (Wolfram engine) with Julia, and Maxima:

(15) -> airyAi 3.7

(15) 0.0017455720028082336
Type: DoubleFloat
(16) -> airyAi jWSReal "3.7"

(16) 0.0017455720006099784
Type: JuliaWSReal
(17) -> airyBi 3.7

(17) 47.5607474995892
Type: DoubleFloat
(18) -> airyBi jWSReal "3.7"

(18) 47.560747499589446
Type: JuliaWSReal

In Maxima:
(%i1) airy_ai(3.7);
(%o1) 0.0017455720006099781
(%i2) airy_bi(3.7);
(%o2) 47.56074749958946

So as you can see it is not the coercion to DoubleFloat that is to
blame. It seems at first glance that the airyAi computation occurs on
a ~32 bits floating point number. I have not checked. I have another
package that operates also on 64 bits floats that comforts this but
with different function naming schemes, and I plan to uniformize this
right now.

(32) -> )expose JF64SF2
JuliaFloat64SpecialFunctions2 is now explicitly exposed in frame
frame0
(32) -> airyai(jf64 3.7)

(32) 0.001745572000609978
Type: JuliaFloat64

So my question is, do I fill an issue for this particular special
function, and populate later this issue with other function issues
after I encounter others like this one? Or do I open separate issues?
It is also probable that even if I fill an issue for different
functions tested, I encounter some other later. Frankly I would prefer
not to encounter the airyAy imprecision problem for other functions.
And yes, I know machine floating point numbers are inherently
imprecise to represent real numbers.

I have a test file for the other package on 64 bits floats mentioned
above*, I can modify and use it to have a better view if necessary on
a set of special functions.

- Greg

* I use it to test yet another package, sorry, that operates on high
precision floating point numbers.

(33) -> )expose JFSF2
JuliaFloatSpecialFunctions2 is now explicitly exposed in frame
frame0
(33) -> airyai(jfloat "3.7")

(33)
0.001745572000609979136759728232445420447338438146589128698021602847417405690
988777
Type: JuliaFloat
(34) -> airyAi jWSReal "3.7"

(34) 0.0017455720006099784
Type: JuliaWSReal
(35) -> airyAi 3.7

(35) 0.0017455720028082336
Type: DoubleFloat

[Waldek Hebisch]

On Fri, Jun 07, 2024 at 06:10:51AM +0200, Grégory Vanuxem wrote:
> Hello,
>
> In FriCAS, airyAi has a problem of precision.

This is known problem. In fact, this is problem with Bessel
functions. I wrote about this in the past. I did some
work on a new better implementation, but it is still not finished
(not ready to be included).

> (17) -> airyBi 3.7
>
> (17) 47.5607474995892
> Type: DoubleFloat
> (18) -> airyBi jWSReal "3.7"
>
> (18) 47.560747499589446
> Type: JuliaWSReal
>
> In Maxima:
> (%i1) airy_ai(3.7);
> (%o1) 0.0017455720006099781
> (%i2) airy_bi(3.7);
> (%o2) 47.56074749958946
>
> So as you can see it is not the coercion to DoubleFloat that is to
> blame. It seems at first glance that the airyAi computation occurs on
> a ~32 bits floating point number.

No, simply formula used to compute airyAi is badly conditioned
leading to significant loss of accuracy.

Note airyBi above: there is some discrepancy in results, it
is normal when different methods are in use. More precisely,
_computing_ a function using DoubleFloat leads to loss of
accuracy and last few digits can be wrong. In case of airyAi
currently there is substantial loss of accuracy and there is
known method giving better accuracy, so we should get better
method. For airyBi I consider small loss of accuracy as
unavidable.

> I have not checked. I have another
> package that operates also on 64 bits floats that comforts this but
> with different function naming schemes, and I plan to uniformize this
> right now.
>
> (32) -> )expose JF64SF2
> JuliaFloat64SpecialFunctions2 is now explicitly exposed in frame
> frame0
> (32) -> airyai(jf64 3.7)
>
> (32) 0.001745572000609978
> Type: JuliaFloat64
>
> So my question is, do I fill an issue for this particular special
> function, and populate later this issue with other function issues
> after I encounter others like this one? Or do I open separate issues?
> It is also probable that even if I fill an issue for different
> functions tested, I encounter some other later. Frankly I would prefer
> not to encounter the airyAy imprecision problem for other functions.

As I wrote this is proble with Bessel functions. There is also
problem with complex polygamma function (I do not know how bad
it is). Basically all functions computed in 'sfsfun.boot' have
smaller or bigger problems. But they are better than nothing,
for example most of the time you will get resonable plot (which
is main use of them).

More generaly, our DoubleFlat functions normally do computations
using DoubleFlat (there are few fakes that delegate real work
to multiple precision routines) and that leads to natural loss
of accuracy.

Concerning Githib issues, I am not sure if it is worth opening
issue about floating point inaccuracy: this is huge topic, was
discussed in the mailing list and is unlikely to gather substantial
useful information.

--
Waldek Hebisch

[Grégory Vanuxem]

Hello,

Ok, maybe a quick note about "substantial loss of accuracy" for some
DoubleFloat functions, if it does not exist, can be added somewhere.

> More generaly, our DoubleFlat functions normally do computations
> using DoubleFlat (there are few fakes that delegate real work
> to multiple precision routines) and that leads to natural loss
> of accuracy.
>
> Concerning Githib issues, I am not sure if it is worth opening
> issue about floating point inaccuracy: this is huge topic, was
> discussed in the mailing list and is unlikely to gather substantial
> useful information.

Yes.

- Greg