Numerical instablity of ErfInverse with very small number

Dear ROOTers,

I calculate the statistical significance for cases where the p-value is very small,

the macro being:

// b = background, d = data
void ss( int b, int d ){

Double_t p = TMath::Poisson( d, b );
Double_t Z = TMath::Sqrt(2)*TMath::ErfInverse( 1 - p );
cout << p << endl;
cout << Z << endl;


In my case the order of magnitude is b = 300, d = 600
, the corresponding p value if 8e-53, the ErfInverse functino returns 0 sigma for such value.
If I test the ErfInverse function for gradually smaller numbers, it is fine until e-17 when it suddenly drops to zero.

It looks like some some numerical unstability or a too rough precision in the estimation of the integral.
Do you confirm ? Is there a way to bypass this problem ?

Thanks in advance,




these are common numerical instabilities that, as far I am concerned, can be faced adequately with one technique:

Rewriting your formula; playing with math.

For example look what ROOT outputs for these expressions:

(const double)0.00000000000000000e+00


that according to the definition of erfc

they are equal.

In a similar fashion you can write erf^-1(1-x) like erfc^-1(x) according to eq. 8 of


root [36] TMath::ErfcInverse(1.e-85)

Check it with Mathematica. I think it works fine.


Hej Leonidas,

Thanks for the trick,

I am left for a particularly extreme case where b = 3000 and d = 6000,
which gives me:

root [12] ss( 3000, 6000, 1 )
Error in TMath::NormQuantile: probability outside (0, 1)
p = 0
Z = -0

Is there a way to go that far in the tails ? Otherwise I may simply quote a lower bound on the significance,




for getting the significance, you should not use this function but

ROOT::Math::normal_quantile_c(p, 1)

which is the inverse of the integral of the normal distribution from +inf to x.
In your case by doing 1.-p if p is less than the double eps ( ~ 2 E-16 ) the result will be always 1.