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 On Thursday, April 8, 2021 at 1:17:40 PM UTC-4 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 > > 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 > > > > On 08.04.2021 18:13, Waldek Hebisch wrote: > > 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. > > > -- 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/41782bd5-6af1-41c9-b77f-7f22fc7ac89fn%40googlegroups.com.
