Airy, Bessel and hypergeometric functions

[Waldek Hebisch]

This is mainly to let you know why bug fixes take so much
time. There is old bug, namely our numeric Airy functions
fails for positive real arguments (and gives quite wrong
results for some complex arguments). Problem is that
our current formula for Airy Ai uses Bessel functions and
square root. For some arguments (including positive reals)
we get wrong branch of square root and conseqently
completely wrong result. Mathematically, fix
looks easy. Namely there is formula for Airy Ai:

Ai(x) = (3)^(-2/3)*(1/Gamma(2/3))*H01(2/3, (1/9)*x^3) -
(3)^(-1/3)*(1/Gamma(1/3)* x*H01(4/3, (1/9)*x^3)

where H01 is hypergeometric function. H01 is analytic
for all complex x, so no problem with branches. But
it has two problem, one fundamental. First, for
positve reals H01 grows quite fast while Ai decays fast.
So we are subracting two large numbers to get tiny
difference, which leads to significant loss of
accuracy (this is fundamental problem). This one can
be corrected by expressing Ai in terms of Bessel K
function (this only for limited set of arguments, as
outherwise we would get back problems with branches).
However, our Bessel K uses formula that is similar
to expression for Ai above, and in particular has
the sane problem with loss of accuracy that we are
trying to avoid...

Other problem is that for some arguments our H01 is
quite inaccurate. Now, H01 is closely related to Bessel
J function:

H01(p, x) = (sqrt(-x))^(-(p-1))*Gamma(p)J_{p-1}(2 sqrt(-x))

For efficiency it makes some sense to have separate
special method for Bessel J and for H01, but one
can use common method just doing little adjustment
implied by equation above. Certainly, there is no
reson to have significant difference in accuracy.
For positive reals formula for H01 means that we take
square root of negative number, that is evaluate
Bessel J at imaginary arguments. Alternativly, one
can use Bessel I function. Again, for complex arguments
Bessel J and Bessel I are closely related and one
can use common method. But our (Boot) code for
Bessel J and Bessel K differ quite a lot. In particular
Bessel I raises error for most arguments. This means
we have not only problem with Airy functions, but
also we have problem with Bessel functions (and H01).

I looked at methods used to compute Bessel functions.
There are some:
- direct use of power series (good only for small x)
- asymptotic expansions (needs large x or p)
- Chebyshev series (moderate arguments)
- recurences

There are several asymptotic expansions, each with its
validity conditions and with different errors. It seems
that most expansions are given without error bounds, so
it is not clear where they are really useful. Recurences
have problem with stability, some are claimed stable,
but this is probably valid only for real p and real x.
My calculations suggest that in "interesting" ranges
of p and x recurences are unstable. What looked
like easy bugfix expanded to a research program...

Anyway, it took some time and I probably will postopone
fixing this one past current release. But there are
few other "easy looking" potential fixes...

BTW: Obvious thing to do is to look at other implementations.
ATM I have not found "full range" implementation, that
is one giving good accuracy for all arguments where
function may be computed (there good looking implementations
for real p). According to various papers Mathematica
is using multiple procision which gives good accuracy,
but is so slow for some arguments that range is
limited by compute time. AFAICS the same it true of
mpmath.

--
Waldek Hebisch

[Kurt Pagani]

Maybe you are already aware of

COMPUTATION OF BESSEL AND AIRY FUNCTIONS AND OF RELATED GAUSSIAN QUADRATURE
FORMULAE by WALTER GAUTSCHI
https://www.cs.purdue.edu/homes/wxg/selected_works/section_02/169.pdf

The code is in Fortran though.
https://www.cs.purdue.edu/archives/2001/wxg/codes/

IMO it would be easier, however, if we took it from GSL:
https://www.gnu.org/software/gsl/doc/html/specfunc.html
Bill wrote an interface: https://github.com/billpage/gsla/blob/master/gsl.lisp

[riccard...@gmail.com]

Here is a link to ARB Airy implementation:

https://github.com/fredrik-johansson/arb/blob/master/acb_hypgeom/airy.c

Riccardo

[Tobias Neumann]

I'm sure there are tons of these implementations available, but here's another one

(one that I also use in production):

https://www.netlib.org/amos/ https://dl.acm.org/doi/10.1145/7921.214331

The paper gives a nice visual presentation of which method is used in which region.

It's quite complicated indeed.

What is your policy on integrating such functions in FriCAS? Only native SPAD

implementations, or are interfaces to C/Fortran OK? Consequently, is arbitrary precision a requirement

or is double precision enough?

Personally I would not mind having specialized C/Fortran codes like arb, flint etc. that are well packaged

and self-contained included or optional in FriCAS. And with that, generally, some C/Fortran

interfacing, even at the cost of supporting less Lisp implementations (I think in practice only ecl and sbcl work anyway?)

Tobias

[Waldek Hebisch]

On Thu, Apr 08, 2021 at 07:17:42PM +0200, Kurt Pagani wrote:
> Maybe you are already aware of
>
> COMPUTATION OF BESSEL AND AIRY FUNCTIONS AND OF RELATED GAUSSIAN QUADRATURE
> FORMULAE by WALTER GAUTSCHI
> https://www.cs.purdue.edu/homes/wxg/selected_works/section_02/169.pdf

Thanks for reference.

> The code is in Fortran though.
> https://www.cs.purdue.edu/archives/2001/wxg/codes/
>
> IMO it would be easier, however, if we took it from GSL:
> https://www.gnu.org/software/gsl/doc/html/specfunc.html
> Bill wrote an interface: https://github.com/billpage/gsla/blob/master/gsl.lisp

Both ilutrate what I wrote: Gautschi has 0<p<1 and real x. GSL
has real arguments. As I wrote, there are good looking
implementations for real p (and complex x). A have a toy
multiple precision implementation, it is fast enough
if arguments are resonably small. OTOH our Bessel J and
H01 also seem to havve resonable accuracy when arguments
are small. The problem is large arguments. Also, it
would be silly to use multiple precision implementation
to get double accuracy, but do not provide real multiple
precision implementation. But with multiple precision
results problem of large arguments gets worse...

Just to be clear: in case of Airy we can try to hide
bug, that is provide low accuracy version. That
would be easy, but the accuracy problem would come back
later. The point is to sort out the mess so that
it works. And ATM it seem that nobody else did this
in really satisfactory way, so appearently we should
be first to do this.

--
Waldek Hebisch

[Kurt Pagani]

Oh, I must have had some other Fortran code in mind (Gautschi has only $p\in
(0,1)$, indeed).

https://jblevins.org/mirror/amiller/

I guess it was the one below, but there are some others (CTRL-F Bessel ...)
https://jblevins.org/mirror/amiller/cbsslj.f90
cbsslj.f90 Complex Bessel function J_{\nu}(z) where both the argument, z, and
the order, \nu, are complex.

Anyway, it makes hardly sense to convert F90 to fricas, but it may be helpful
for comparing purposes.

[Waldek Hebisch]

On Thu, Apr 08, 2021 at 11:05:45AM -0700, Tobias Neumann wrote:
> I'm sure there are tons of these implementations available, but here's
> another one
> (one that I also use in production):
>
> https://www.netlib.org/amos/ https://dl.acm.org/doi/10.1145/7921.214331
> The paper gives a nice visual presentation of which method is used in which
> region.
> It's quite complicated indeed.

It is one of "good loking" that I mentioned. It should be
fine when you are satisfied with real p (it suppors complex x).
But we need complex p.

> What is your policy on integrating such functions in FriCAS? Only native
> SPAD
> implementations, or are interfaces to C/Fortran OK? Consequently, is
> arbitrary precision a requirement
> or is double precision enough?

My policy is "never say never". It is pragmatic choice: how badly
we need given function, what is effort to reimplement compared to
resue etc. Having said that, up to now looking at various
factors in most cases my conclusion was to provide implementation
in Spad. Factors are:
- I am mostly looking at long time perspecitive, adding a
single function in moderate improvement and can wait
some time
- beside double precision version we would like to have multiple
precision one (not as absolute requirement, but strong
preference)
- eventually we would like to cover a lot of functions and
full range of arguments. For this it is natural to
modify/tweak existing code and Spad it preferable there
- one can use symbolic version to help create numeric one. One
can use multiple precision version to help develop double
precision one
- we would like not only code but also explanation why it
works

ATM we have a few functions that actual implementation
is multiple-precision only, double variant is just wrapper for
multiple precision. For elliptic functions we have separate
double precision version using popular Carlson method. And we
have multiple-precison version based on Gauss-Landen transforms.
There are still holes and bugs, but overall IMO it one of best
implementations available. Let me add there there is some
not entirely trivial math there to avoid errors. I have
somewhat rough (unfinished) paper describing the issue.
One reason that it is unfinished is that there are few things
to work out: there is still one corner case in ellipticPi that
should be handled (and currently remain unhadled) and few
uniplemented functions like Weirestrass Sigma and Jacobi Theta.
When finished we should have very good implementation with
explanation how this is done.

For Airy, Bessel and hypegometric F01 there is extra constraint:
we have existing code and new version should improve it
as opposed to regression. Since existing version is double
precision one, it is natural to improve double precision
version before providing multiple precision one. Actually,
the core isssue is similar in both cases: to have good
implementation we need to handle large complex arguments.

> Personally I would not mind having specialized C/Fortran codes like arb,
> flint etc. that are well packaged
> and self-contained included or optional in FriCAS. And with that,
> generally, some C/Fortran
> interfacing, even at the cost of supporting less Lisp implementations (I
> think in practice only ecl and sbcl work anyway?)

There is Closure CL, AFAIK it works fine. I do not know how
many FriCAS users buid it on top of Closure CL, but it is
quite good Lisp implementation. One higlight is precise
garbage collector: with sbcl (and probably ecl too) it can
happen that there is free memory but garbage collector can
not reclaim it and consequently one get "out of memory"
error. In practice when one get "out of memory" erorr
it is hard to tell exactly what is the reason. Theoretically
Closure CL should behave better in this case.

AFAICS Clisp works, but is least preferable due to performance.
Old GCL mostly works, but there is a bug not present when
using older implementations. Newest (git checkout) GCL
currently does not work.

--
Waldek Hebisch

[Tobias Neumann]

There is Closure CL, AFAIK it works fine. I do not know how
many FriCAS users buid it on top of Closure CL, but it is
quite good Lisp implementation. One higlight is precise
garbage collector: with sbcl (and probably ecl too) it can
happen that there is free memory but garbage collector can
not reclaim it and consequently one get "out of memory"
error. In practice when one get "out of memory" erorr
it is hard to tell exactly what is the reason. Theoretically
Closure CL should behave better in this case.

The precise garbage collection sounds good indeed.

I recently tried FriCAS compiled with Clozure CL to compare SBCL/ECL/CCL performance,

but for me it started failing with some very basic things:

 test := 1.1

   >> System error:

   The value NIL is not of the expected type CCL::UVECTOR.

[Dima Pasechnik]

On Tue, Apr 13, 2021 at 6:36 PM Tobias Neumann <tobias....@gmail.com> wrote:
>>
>> There is Closure CL, AFAIK it works fine. I do not know how
>> many FriCAS users buid it on top of Closure CL, but it is
>> quite good Lisp implementation. One higlight is precise
>> garbage collector: with sbcl (and probably ecl too) it can

ECL is using the well-known Boehm–Demers–Weiser GC - which is very
widely used, and is rather robust and fast, it's multithreaded etc -
if used correctly.

In SageMath (which embeds ECL to run Maxima, FriCAS, etc) we rarely notice
any ECL's GC problems.

>> happen that there is free memory but garbage collector can
>> not reclaim it and consequently one get "out of memory"
>> error. In practice when one get "out of memory" erorr
>> it is hard to tell exactly what is the reason. Theoretically
>> Closure CL should behave better in this case.
>
>
> The precise garbage collection sounds good indeed.
> I recently tried FriCAS compiled with Clozure CL to compare SBCL/ECL/CCL performance,
> but for me it started failing with some very basic things:
>
> test := 1.1
> >> System error:
> The value NIL is not of the expected type CCL::UVECTOR.
>

> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/7576ed80-58ec-4d1c-aab9-7ad713eb1300n%40googlegroups.com.

[Richard Fateman]

Here's a paper that discusses Bessel function computation in the context of arbitrary precision.

Miller's algorithm  (backward recurrence) is one method, but there are other methods for

large argument.

https://people.eecs.berkeley.edu/~fateman/papers/hermite.pdf

Have fun.

[Waldek Hebisch]

Which version of Clozure CL did you try? ATM my testing version
is 11.5.


--
Waldek Hebisch

[Waldek Hebisch]

On Tue, Apr 13, 2021 at 01:58:32PM -0700, Richard Fateman wrote:
> Here's a paper that discusses Bessel function computation in the context of
> arbitrary precision.
> Miller's algorithm (backward recurrence) is one method, but there are
> other methods for
> large argument.
>
> https://people.eecs.berkeley.edu/~fateman/papers/hermite.pdf
>
> Have fun.

Earlier I wrote:
> > >
> > > BTW: Obvious thing to do is to look at other implementations.
> > > ATM I have not found "full range" implementation, that
> > > is one giving good accuracy for all arguments where
> > > function may be computed (there good looking implementations
> > > for real p). According to various papers Mathematica
> > > is using multiple procision which gives good accuracy,
> > > but is so slow for some arguments that range is
> > > limited by compute time. AFAICS the same it true of
> > > mpmath.

Thanks for reference, but it seems that your paper have
similar limitation as other works: it concentrates on
real order (or maybe even integer orders).

--
Waldek Hebisch

[Tobias Neumann]

> test := 1.1
> >> System error:
> The value NIL is not of the expected type CCL::UVECTOR.

Which version of Clozure CL did you try? ATM my testing version
is 11.5.

I am using 1.12 (it's the newest release)

"Clozure Common Lisp Version 1.12 (v1.12-37-g27c4bc11) LinuxX8664" 

https://github.com/Clozure/ccl/releases/tag/v1.11.5

Is this the one you are using? Seems that release is already a bit older, from 2017. I can also try that

[Waldek Hebisch]

I somehow missed that Clozure CL folks declared their github repository
to be a release. AFAICS they changed how binum routines work, so
we need to adapt GMP interface. For me build without GMP works,
but bignum tests take longer...


--
Waldek Hebisch

[Waldek Hebisch]

On Tue, Apr 13, 2021 at 03:49:16PM -0700, Tobias Neumann wrote:
>
> >

I have now pushed change so that Clozure CL 12 should work with GMP.

--
Waldek Hebisch