search

LSSTDESC / N5K

Non-local No-Nonsense Non-Limber Numerical Knockout
BSD 3-Clause "New" or "Revised" License
7 stars 7 forks source link

Bessel functions of high order #6

Open tilmantroester opened 3 years ago

tilmantroester commented 3 years ago

Both GSL and boost seem to struggle with Bessel functions of high order (ell>200 or so). GSL is a bit more verbose and gives underflow errors, while boost quietly ignores that issue unless one sets the policy to explicitly throw errors on under/overflows. Since the current setup requests ells up to 2000, this might be a common problem.

Has anyone compared GSL or boost against asymptotic forms (e.g., https://dlmf.nist.gov/10.57)?

damonge commented 3 years ago

The implementation of j_ell I've used in the past is the one in CAMB (which also lives in CCL!), and tends to be more or less stable. I haven't looked into the code enough to see if it's implicitly using the asymptotic forms.

tilmantroester commented 3 years ago

The CCL implementation seems to match the other ones, as far as I can tell. The cases I looked at had small arguments where the result is just really small, so the underflow is warranted.

ntessore commented 3 years ago

If you want something like BesselJ[2000, 5] which is about 1e-5000 you might want doi:10.1002/cnm.972, in particular for products of Bessel functions.

ntessore commented 3 years ago

This is actually pretty useful. Here is a rough implementation, perhaps it helps.

$ cc logbessel.c
$ # compute ln[j_10000(5.0)] and compare to Mathematica
$ ./a.out
jnu = +1.000000 exp(-72950.7471329014661023)
     [+1.000000 exp(-72950.7471329014661023)]
EiffL commented 3 years ago

For what it's worth, I've made a small Python/C++ interface for some Boost Bessel funtions, would be trivial to extend it to provide Python access to all gsl/boost Bessel utilities: https://github.com/EiffL/N5K/blob/master/n5k/cxx/bessel_tools.cpp Let me know if that would be useful and I can make a PR just for these tools