search

Axect / puruspe

PURe RUSt SPEcial library
Apache License 2.0
11 stars 6 forks source link

Error estimates from special functions #5

Open aoanla opened 2 years ago

aoanla commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

Any contribution is always welcome :smile: 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!