[javascript] erfinv(x) implementations (and erf(x)/erfc(x) implementations ported from cephes)

erf.js

// erf(x): error function (see: https://en.wikipedia.org/wiki/Error_function)
// by https://github.com/jeremybarnes/cephes/blob/master/cprob/ndtr.c

const P = [
    2.46196981473530512524E-10,
    5.64189564831068821977E-1,
    7.46321056442269912687E0,
    4.86371970985681366614E1,
    1.96520832956077098242E2,
    5.26445194995477358631E2,
    9.34528527171957607540E2,
    1.02755188689515710272E3,
    5.57535335369399327526E2,
];
const Q = [
    1.32281951154744992508E1,
    8.67072140885989742329E1,
    3.54937778887819891062E2,
    9.75708501743205489753E2,
    1.82390916687909736289E3,
    2.24633760818710981792E3,
    1.65666309194161350182E3,
    5.57535340817727675546E2,
];
const R = [
    5.64189583547755073984E-1,
    1.27536670759978104416E0,
    5.01905042251180477414E0,
    6.16021097993053585195E0,
    7.40974269950448939160E0,
    2.97886665372100240670E0,
];
const S = [
    2.26052863220117276590E0,
    9.39603524938001434673E0,
    1.20489539808096656605E1,
    1.70814450747565897222E1,
    9.60896809063285878198E0,
    3.36907645100081516050E0,
];
const T = [
    9.60497373987051638749E0,
    9.00260197203842689217E1,
    2.23200534594684319226E3,
    7.00332514112805075473E3,
    5.55923013010394962768E4,
];
const U = [
    3.35617141647503099647E1,
    5.21357949780152679795E2,
    4.59432382970980127987E3,
    2.26290000613890934246E4,
    4.92673942608635921086E4,
];
function polevl(x, c) {
    return c.reduce((r, c) => r * x + c, 0);
}
function p1evl(x, c) {
    return c.reduce((r, c) => r * x + c, 1);
}

function erf(x) {
    if (Math.abs(x) > 1) return 1 - erfc(x);
    const z = x * x;
    return x * polevl(z, T) / p1evl(z, U);
}
// erfc(x) = 1 - erf(x)
const MAXLOG = Math.log(Number.MAX_VALUE);
function erfc(x0) {
    const x = Math.abs(x0);
    if (x < 1) return 1 - erf(x);
    const z = -x0 * x0;
    if (z < -MAXLOG) return x0 < 0 ? 2 : 0;
    const [p, q] = x < 8 ? [P, Q] : [R, S];
    const y = Math.exp(z) * polevl(x, p) / p1evl(x, q);
    return x0 < 0 ? 2 - y : y;
}
// erfce(x) = exp(x**2) * erfc(x)
function erfce(x) {
    console.assert(x > 1);
    const [p, q] = x < 8 ? [P, Q] : [R, S];
    return polevl(x, p) / p1evl(x, q);
}
// ndtr(a) = (1 + erf(a/sqrt(s))) / 2
const SQRTH = 2 ** 0.5 / 2;
function ndtr(a) {
    const x = a * SQRTH;
    const z = Math.abs(x);
    if (z < 1) return (1 + erf(x)) / 2;
    const y = erfc(z) / 2;
    return x > 0 ? 1 - y : y;
}

// check
const x = Array.from(Array(60), (_, i) => -3 + i * 0.1);
x.forEach(x => console.log(`erf(${x.toFixed(1).padStart(4)}) = ${erf(x)}`));

console.log("-".repeat(80));
const y = x.map(x => erf(x));
for (let i = 0; i < 20; i++) {
    const max = 1 - 0.1 * i;
    const min = max - 0.1;
    const plot = y.map(y => min < y && y <= max ? "*" : " ").join("");
    const r = `${max.toFixed(1).padStart(4)} ~ ${min.toFixed(1).padStart(4)}`;
    console.log(`${r}: ${plot}`);
}
console.log("-".repeat(80));

erfinv.js

// inverse of erf(x)
// see: https://en.wikipedia.org/wiki/Error_function

function* cgen() {
    yield 1;
    const cs = [1];
    for (let k = 1;; k++) {
        let r = 0;
        for (let m = 0; m < k; m++) {
            r += cs[m] * cs[k - 1 - m] / (m + 1) / (2 * m + 1);
        }
        yield r;
        cs[k] = r;
    }
}

function erfinv1(z) {
    // see erf^-1(z) at "Inverse functions"
    let r = 0;
    const c = cgen();
    const pz = z / 2 * (Math.PI ** 0.5);
    for (let k = 0; k < 2000; k++) {
        const i = 2 * k + 1, cv = c.next().value;
        r += cv / i *  (pz ** i);
    }
    return r;
}

function erfinv2(z) {
    // see: erf^-1(z) at "Approximation with elementary functions"
    const a = 0.140012, lnzz = Math.log(1 - z * z), pa2 = 2 / Math.PI / a;
    return Math.sign(z) *
        (((pa2 + lnzz / 2) ** 2 - lnzz / a) ** 0.5 - pa2 - lnzz / 2) ** 0.5;
}

// erfinv(x) = ndtri((x+1)/2) / sqrt(2)
// ndtri(y): inverse function of integral of normal distribution
// by https://github.com/jeremybarnes/cephes/blob/master/cprob/ndtri.c
const P0 = [
    -5.99633501014107895267E1,
    9.80010754185999661536E1,
    -5.66762857469070293439E1,
    1.39312609387279679503E1,
    -1.23916583867381258016E0,
];
const Q0 = [
    1.95448858338141759834E0,
    4.67627912898881538453E0,
    8.63602421390890590575E1,
    -2.25462687854119370527E2,
    2.00260212380060660359E2,
    -8.20372256168333339912E1,
    1.59056225126211695515E1,
    -1.18331621121330003142E0,
];
const P1 = [
    4.05544892305962419923E0,
    3.15251094599893866154E1,
    5.71628192246421288162E1,
    4.40805073893200834700E1,
    1.46849561928858024014E1,
    2.18663306850790267539E0,
    -1.40256079171354495875E-1,
    -3.50424626827848203418E-2,
    -8.57456785154685413611E-4,
];
const Q1 = [
    1.57799883256466749731E1,
    4.53907635128879210584E1,
    4.13172038254672030440E1,
    1.50425385692907503408E1,
    2.50464946208309415979E0,
    -1.42182922854787788574E-1,
    -3.80806407691578277194E-2,
    -9.33259480895457427372E-4,
];
const P2 = [
    3.23774891776946035970E0,
    6.91522889068984211695E0,
    3.93881025292474443415E0,
    1.33303460815807542389E0,
    2.01485389549179081538E-1,
    1.23716634817820021358E-2,
    3.01581553508235416007E-4,
    2.65806974686737550832E-6,
    6.23974539184983293730E-9,
];
const Q2 = [
    6.02427039364742014255E0,
    3.67983563856160859403E0,
    1.37702099489081330271E0,
    2.16236993594496635890E-1,
    1.34204006088543189037E-2,
    3.28014464682127739104E-4,
    2.89247864745380683936E-6,
    6.79019408009981274425E-9,
];
// https://github.com/jeremybarnes/cephes/blob/master/cprob/polevl.c
function polevl(x, c) {
    return c.reduce((r, c) => r * x + c, 0);
}
function p1evl(x, c) {
    return c.reduce((r, c) => r * x + c, 1);
}
const expm2 = Math.exp(-2), s2pi = (2 * Math.PI) ** 0.5, sqrt2i = 2 ** -0.5;

function ndtri(y0) {
    if (y0 <= 0) return -Infinity;
    if (y0 >= 1) return Infinity;
    const ncode = y0 > 1 - expm2;
    const y = ncode ? 1 - y0 : y0;
    if (y > expm2) {
        const y1 = y - 0.5;
        const y2 = y1 * y1;
        const y3 = y2 * polevl(y2, P0) / p1evl(y2, Q0);
        return (y1 + y1 * y3) * s2pi;
    }
    const x = (-2 * Math.log(y)) ** 0.5;
    const z = 1 / x;
    const [p, q] = x < 8 ? [P1, Q1] : [P2, Q2];
    const x0 = x - Math.log(x) * z;
    const x1 = z * polevl(z, p) / p1evl(z, q);
    return ncode ? x0 - x1 : x1 - x0;
}
function erfinv3(z) {
    return ndtri((z + 1) / 2) * sqrt2i;
}


function f32(x) {
    return x.toPrecision(6);
}

console.log("\n[erfinv1]");
console.log(f32(erfinv1(0.5)));
console.log(f32(erfinv1(0.999)));
console.log(f32(erfinv1(0.999999))); // bad result
console.log(f32(erfinv1(-0.5)));
printHist(Array.from(Array(1000), _ => randn(erfinv1)), -4.5, 4.5, 9);

console.log("\n[erfinv2]");
console.log(f32(erfinv2(0.5)));
console.log(f32(erfinv2(0.999)));
console.log(f32(erfinv2(0.999999)));
console.log(f32(erfinv2(-0.5)));
printHist(Array.from(Array(10000), _ => randn(erfinv2)), -4.5, 4.5, 9);

console.log("\n[erfinv3]");
console.log(f32(erfinv3(0.5)));
console.log(f32(erfinv3(0.999)));
console.log(f32(erfinv3(0.999999)));
console.log(f32(erfinv3(-0.5)));
printHist(Array.from(Array(10000), _ => randn(erfinv3)), -4.5, 4.5, 9);

//[values from wolframalpha]
//erfinv(0.5) = 0.476936...
//erfinv(0.999) = 2.32675...
//erfinv(0.999999) = 3.45891...
//erfinv(-0.5) = -0.476936...

// for randn and histgram
function randn(erfinv) {
    return Math.sqrt(2) * erfinv(2 * Math.random() - 1);
}
function hist(vs, min, max, n) {
    const step = (max - min) / n;
    const h = Array.from(Array(n), _ => 0);
    for (const v of vs) {
        const i = (v - min) / step | 0;
        if (0 <= i && i < h.length) h[i]++;
    }
    return h;
}
function printHist(vs, min, max, n) {
    const step = (max - min) / n;
    const h = hist(vs, min, max, n);
    const mins = Array.from(Array(n), (_, i) => min + step * i);
    const l = h.reduce((l, v) => l > v ? l : v, 0);
    console.log("-".repeat(80));
    h.forEach((count, i) => {
        const min = mins[i].toFixed(2).padStart(6);
        const max = (mins[i] + step).toFixed(2).padStart(6);
        const countStr = count.toString().padStart(6);
        const bar = "*".repeat(count / l * 50 >>> 0);
        console.log(`${min} ~ ${max} (${countStr}) ${bar}`);
    });
    console.log("-".repeat(80));
}

result

$ node erfinv.js

[erfinv1]
0.476936
2.32611
2.70058
-0.476936
--------------------------------------------------------------------------------
 -4.50 ~  -3.50 (     1)
 -3.50 ~  -2.50 (     7)
 -2.50 ~  -1.50 (    69) *********
 -1.50 ~  -0.50 (   244) *********************************
 -0.50 ~   0.50 (   359) **************************************************
  0.50 ~   1.50 (   248) **********************************
  1.50 ~   2.50 (    65) *********
  2.50 ~   3.50 (     7)
  3.50 ~   4.50 (     0)
--------------------------------------------------------------------------------

[erfinv2]
0.476919
2.31882
3.45003
-0.476919
--------------------------------------------------------------------------------
 -4.50 ~  -3.50 (     1)
 -3.50 ~  -2.50 (    56)
 -2.50 ~  -1.50 (   583) *******
 -1.50 ~  -0.50 (  2476) ********************************
 -0.50 ~   0.50 (  3848) **************************************************
  0.50 ~   1.50 (  2413) *******************************
  1.50 ~   2.50 (   578) *******
  2.50 ~   3.50 (    43)
  3.50 ~   4.50 (     2)
--------------------------------------------------------------------------------

[erfinv3]
0.476936
2.32675
3.45891
-0.476936
--------------------------------------------------------------------------------
 -4.50 ~  -3.50 (     1)
 -3.50 ~  -2.50 (    55)
 -2.50 ~  -1.50 (   614) ********
 -1.50 ~  -0.50 (  2423) *******************************
 -0.50 ~   0.50 (  3798) **************************************************
  0.50 ~   1.50 (  2448) ********************************
  1.50 ~   2.50 (   597) *******
  2.50 ~   3.50 (    57)
  3.50 ~   4.50 (     7)
--------------------------------------------------------------------------------