{ ******************************************************************
  Inverse of Normal distribution function
  Translated from C code in Cephes library (http://www.moshier.net)
  ****************************************************************** }

unit uinvnorm;

interface

uses
  utypes, uminmax, upolev;

function InvNorm(P : Float) : Float;
{ ------------------------------------------------------------------
  Inverse of Normal distribution function

  Returns the argument, X, for which the area under the Gaussian
  probability density function (integrated from minus infinity to X)
  is equal to P.
  ------------------------------------------------------------------ }

implementation

  function InvNorm(P : Float) : Float;
  const
    P0 : TabCoef = (
        8.779679420055069160496E-3,
      - 7.649544967784380691785E-1,
        2.971493676711545292135E0,
      - 4.144980036933753828858E0,
        2.765359913000830285937E0,
      - 9.570456817794268907847E-1,
        1.659219375097958322098E-1,
      - 1.140013969885358273307E-2,
        0, 0);

    Q0 : TabCoef = (
      - 5.303846964603721860329E0,
        9.908875375256718220854E0,
      - 9.031318655459381388888E0,
        4.496118508523213950686E0,
      - 1.250016921424819972516E0,
        1.823840725000038842075E-1,
      - 1.088633151006419263153E-2,
        0, 0, 0);

    P1 : TabCoef = (
      4.302849750435552180717E0,
      4.360209451837096682600E1,
      9.454613328844768318162E1,
      9.336735653151873871756E1,
      5.305046472191852391737E1,
      1.775851836288460008093E1,
      3.640308340137013109859E0,
      3.691354900171224122390E-1,
      1.403530274998072987187E-2,
      1.377145111380960566197E-4);

    Q1 : TabCoef = (
      2.001425109170530136741E1,
      7.079893963891488254284E1,
      8.033277265194672063478E1,
      5.034715121553662712917E1,
      1.779820137342627204153E1,
      3.845554944954699547539E0,
      3.993627390181238962857E-1,
      1.526870689522191191380E-2,
      1.498700676286675466900E-4,
      0);

    P2 : TabCoef = (
      3.244525725312906932464E0,
      6.856256488128415760904E0,
      3.765479340423144482796E0,
      1.240893301734538935324E0,
      1.740282292791367834724E-1,
      9.082834200993107441750E-3,
      1.617870121822776093899E-4,
      7.377405643054504178605E-7,
      0, 0);

    Q2 : TabCoef = (
      6.021509481727510630722E0,
      3.528463857156936773982E0,
      1.289185315656302878699E0,
      1.874290142615703609510E-1,
      9.867655920899636109122E-3,
      1.760452434084258930442E-4,
      8.028288500688538331773E-7,
      0, 0, 0);

    P3 : TabCoef = (
        2.020331091302772535752E0,
        2.133020661587413053144E0,
        2.114822217898707063183E-1,
      - 6.500909615246067985872E-3,
      - 7.279315200737344309241E-4,
      - 1.275404675610280787619E-5,
      - 6.433966387613344714022E-8,
      - 7.772828380948163386917E-11,
        0, 0);

    Q3 : TabCoef = (
        2.278210997153449199574E0,
        2.345321838870438196534E-1,
      - 6.916708899719964982855E-3,
      - 7.908542088737858288849E-4,
      - 1.387652389480217178984E-5,
      - 7.001476867559193780666E-8,
      - 8.458494263787680376729E-11,
        0, 0, 0);

  var
    X, Y, Z, Y2, X0, X1 : Float;
    Code : Integer;
  begin
    if (P <= 0.0) or (P >= 1.0) then
      begin
        InvNorm := DefaultVal(FDomain, Sgn(P) * MaxNum);
        Exit;
      end;

    Code := 1;
    Y := P;
    if Y > (1.0 - 0.13533528323661269189) then  { 0.135... = exp(-2) }
      begin
        Y := 1.0 - Y;
        Code := 0;
      end;
    if Y > 0.13533528323661269189 then
      begin
        Y := Y - 0.5;
        Y2 := Y * Y;
        X := Y + Y * (Y2 * PolEvl(Y2, P0, 7) / P1Evl(Y2, Q0, 7));
        X := X * Sqrt2Pi;
        InvNorm := X;
        Exit;
      end;

    X := Sqrt(- 2.0 * Ln(Y));
    X0 := X - Ln(X) / X;
    Z := 1.0 / X;
    if X < 8.0 then
      X1 := Z * PolEvl(Z, P1, 9) / P1Evl(Z, Q1, 9)
    else if X < 32.0 then
      X1 := Z * PolEvl(Z, P2, 7) / P1Evl(Z, Q2, 7)
    else
      X1 := Z * PolEvl(Z, P3, 7) / P1Evl(Z, Q3, 7);
    X := X0 - X1;
    if Code <> 0 then
      X := - X;
    InvNorm := X;
  end;

end.