Axect / puruspe

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

gammpapprox produces incorrect results #4

Closed main-- closed 2 years ago

main-- commented 2 years ago

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

I'm sorry for late reply (I didn't notice this issue). Thanks to your good point, I can detect errata in source code - Now, I fix it at Ver 0.2.0.

At Ver 0.2.0 :

extern crate puruspe;
use puruspe::*;

fn main() {
    let a = 100.;
    for x in [1., 10., 50., 75., 99., 100., 101., 110., 1000.] {
        println!("{}\t{}", x, gammp(a, x));
    }
}
1       0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000398128081855225
10      0.000000000000000000000000000000000000000000000000000000000000005398589727602037
50      0.00000000032000653242666364
75      0.0033524414978532177
99      0.47330433035233394
100     0.5132987983276017
101     0.5528962934790187
110     0.8417213299556715
1000    1

Thank you again for your interest and good point! :smile: