Error estimates from special functions

[aoanla]

Something I've noticed that's lacking in (all?) of the Rust-native special functions crates is any actual error estimates on their functions.
For scientific use, it's very important to have an estimate of the accuracy of a particular function evaluation - which is why mature scientific libraries in C [like the Gnu Scientific Library ( https://www.gnu.org/software/gsl/ ) or the MPFR ( https://en.wikipedia.org/wiki/GNU_MPFR )] either report their estimated error in their return value, or allow a required precision to be specified.

Is there any possibility of error estimation in the special functions in this crate?

[Axect]

Thanks for the good point.
I also realize that error reporting or guaranteed precision are very important. So far, the purpose of this crate has been to implement special functions in Rust without any dependencies. Thus, I just followed up implementations in "Numerical Recipes 3rd ed" for convenience. In this book, epsilon is set to 64bit double precision (2^(-53) ~ 1e-16), but I did not rigorously prove that the whole algorithm guarantees the precision.
If I can afford it later, I will document the precision for each algorithm and implement it so that the user can provide the desired precision if possible.

[aoanla]

I'm willing to help contribute to this crate in this matter [and also with regards to other implementations of other special functions - the Bessel J,Y,H series, for example] if it would be appreciated. Whilst the Numerical Recipes approaches are generally very good , I would be surprised if they always reach 0.5ulp precision in their range... I will spend a bit of time writing some error estimation on them.
(You may be interested in: https://members.loria.fr/PZimmermann/papers/glibc234-20210907.pdf , a study by the authors of the Gnu MPFR on the precision limits of various popular math libraries... which shows just how inaccurate an implementation can be.)

[Axect]

Any contribution is always welcome 😄
Especially, I'm not an expert on floating point errors, so your contribution would be really appreciated.
Thanks for the good reference too. That is what I was looking for!