Inverse erf

ROOT
1

hi,

i can’t find the inverse error function in TMath?
is it not implemented?

thanks,
gerhard

2

Hi Gerhard,

No it is not (yet?) available. Are you proposing to provide an implementation?

Cheers,
Philippe.

3

hi,

yes i can implement one, technically.
but i am not sure which algorithm to use.
on the web i found some approximations using polynoms.
is it acceptable to have a hardcoded table with parameters
in root?
maybe somebody more senior than i am can propose
where to find a good algorithm, i can then take care of
the implementation.
if not, i’ll just pull something out of my fingers …

cheers,
gerhard

4

We will include the inverse error function in TMath in teh coming release 4.00/08. Meanwhile, you can use the code below

Rene

Double_t ErfInverse(Double_t y) {
// returns the inverse error function obtained
// y must be <-1<y<1

Int_t kMaxit = 50;
Double_t kEps = 1e-14;
Double_t kConst = 0.8862269254527579; // sqrt(pi)/2.0

if(TMath::Abs(y) <= kEps) return kConst*y;

// Newton iterations
Double_t erfi, derfi, Y0,Y1,DY0,DY1;
if(TMath::Abs(y) < 1.0) {
erfi = kConstTMath::Abs(y);
Y0 = TMath::Erf(0.9erfi);
derfi = 0.1*erfi;
for (Int_t iter=0; iter<kMaxit; iter++) {
Y1 = 1. - TMath::Erfc(erfi);
DY1 = TMath::Abs(y) - Y1;
if (TMath::Abs(DY1) < kEps) {if (y < 0) return -erfi; else return erfi;}
DY0 = Y1 - Y0;
derfi *= DY1/DY0;
Y0 = Y1;
erfi += derfi;
if(TMath::Abs(derfi/erfi) < kEps) {if (y < 0) return -erfi; else return erfi;}
}
}
return 0; //did not converge
}