axect / puruspe

PURe RUSt SPEcial library

PURUSPE

PURe RUSt SPEcial function library. There are no dependencies.

Implemented Functions

Gamma functions

• ln_gamma : Logarithmic gamma function
• gamma : Gamma function
• gammp : Regularized lower gamma function
• gammq : Regularized upper gamma function
• invgammp : Inverse regularized lower gamma function

Beta function

• beta : Beta function
• betai : Regularized Incomplete beta function
• invbetai : Inverse regularized incomplete beta function

Error functions

• erf : Error function
• erfc : Complementary Error function
• inverf : Inverse error function
• inverfc : Inverse complementary error function

Bessel functions

• Jn : Bessel function of the first kind (integer order)
• Yn : Bessel function of the second kind (integer order)
• In : Modified Bessel function of the first kind (integer order)
• Kn : Modified Bessel function of the second kind (integer order)
• Jnu_Ynu : Bessel function of the first kind and second kind (fractional order)
• Inu_Knu : Modified Bessel function of the first kind and second kind (fractional order)
• besseljy : Bessel function of the first and second kinds (include derivatives)
• besselik : Modified Bessel function of the second kind (include derivatives)

Reference

• Press, William H., and William T. Vetterling. Numerical Recipes. Cambridge: Cambridge Univ. Press, 2007.

puruspe's Issues

gammpapprox produces incorrect results

There is a bug in the implementation of gammpapprox that causes it to return incorrect values (i.e. values that are very far away from the expected values). This also causes invgammp to panic in some scenarios (when a >= ASWITCH).

Example:

extern crate puruspe;
use puruspe::*;

fn main() {
    let a = 100.;
    for x in [1., 10., 50., 75., 99., 100., 101., 110., 1000.] {
        println!("{}\t{} vs {}", x, gammp(a, x), gammp_no_aswitch(a, x));
    }
}

1       0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000032263012916501365 vs 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000039812808185605916
10      0.0000000000000000000000000000000000000000000000000000000004374842621019958 vs 0.000000000000000000000000000000000000000000000000000000000000005398589727602311
50      0.00002593229505674661 vs 0.00000000032000653242666364
75      271.67102316179927 vs 0.003352441497853259
99      38355.05310862807 vs 0.4733043303523406
100     -39439.69225878758 vs 0.5132987982280703
101     -36230.83923948253 vs 0.5528962934789909
110     -12825.391829109325 vs 0.8417213299556574
1000    1 vs 1

Error estimates from special functions

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?