code 41501;
real procedure PHINV(PROB); value PROB; real PROB;
begin real EPS;
    real procedure INVERF(X); value X; real X;
    begin real ABSX, P, BETAX;
        real array A[0 : 23];
        real procedure CHEPOLSER(N, X, A);
        value N, X; integer N; real X; array A;
           begin integer K; real H, R, S, TX;
              TX:= X + X; R:= A[N];
               H:= A[N - 1] + R × TX;
               for K:= N - 2 step -1 until 1 do
               begin S:= R; R:= H;
                   H:= A[K] + R × TX - S
               end;
               CHEPOLSER:= A[0]- R + H × X
           end CHEPOLSER;
           
           ABSX:= ABS(X);
           if ABSX ≤ 0.8 then
           begin
               A[ 0]:=                 0.99288   53766   1894;
               A[ 1]:=                 0.12046   75161   4310;
               A[ 2]:=                 0.01607   81993   4210;
               A[ 3]:=                 0.00268   67044   3716;
               A[ 4]:=                 0.00049   96347   3024;
               A[ 5]:=                 0.00009   88982   1860;
               A[ 6]:=                 0.00002   03918   1276;
               A[ 7]:=                 0.00000   43272   7162;
               A[ 8]:=                 0.00000   09380   8141;
               A[ 9]:=                 0.00000   02067   3472;
               A[10]:=                 0.00000   00461   5970;
               A[11]:=                 0.00000   00104   1668;
               A[12]:=                 0.00000   00023   7150;
               A[13]:=                 0.00000   00005   4393;
               A[14]:=                 0.00000   00001   2555;
               A[15]:=                 0.00000   00000   2914;
               A[16]:=                 0.00000   00000   0680;
               A[17]:=                 0.00000   00000   0159;
               A[18]:=                 0.00000   00000   0037;
               A[19]:=                 0.00000   00000   0009;
               A[20]:=                 0.00000   00000   0002;
               A[21]:=                 0.00000   00000   0001;
            INVERF:= CHEPOLSER(21, X × X / 0.32 - 1, A) × X
        end else
        if 1 - ABSX ≥ 2510-4 then
        begin
            A[ 0]:=                      0.91215   88034   1755;
            A[ 1]:=                     -0.01626   62818   6766;
            A[ 2]:=                      0.00043   35564   7295;
            A[ 3]:=                      0.00021   44385   7007;
            A[ 4]:=                      0.00000   26257   5108;
            A[ 5]:=                     -0.00000   30210   9105;
            A[ 6]:=                     -0.00000   00124   0606;
            A[ 7]:=                      0.00000   00624   0661;
            A[ 8]:=                     -0.00000   00005   4015;
            A[ 9]:=                     -0.00000   00014   2321;
            A[10]:=                      0.00000   00000   3438;
            A[11]:=                      0.00000   00000   3358;
            A[12]:=                     -0.00000   00000   0146;
            A[13]:=                     -0.00000   00000   0081;
            A[14]:=                      0.00000   00000   0005;
            A[15]:=                      0.00000   00000   0002;
            BETAX:= SQRT(- LN((1 + ABSX) × (1 - ABSX)));
            P:= -1.54881 30423 7326 × BETAX +
                2.56549 01231 4782;
            P:= CHEPOLSER(15, P, A);
            INVERF:= if X < 0 then - BETAX × P
                else BETAX × P
        end else
        begin
              A[ 0]:=    0.95667 97090   2049;
              A[ 1]:=   -0.02310 70043   0907;
              A[ 2]:=   -0.00437 42360   9751;
              A[ 3]:=   -0.00057 65034   2265;
              A[ 4]:=   -0.00001 09610   2231;
              A[ 5]:=    0.00002 51085   4703;
              A[ 6]:=    0.00001 05623   3607;
              A[ 7]:=    0.00000 27544   1233;
              A[ 8]:=    0.00000 04324   8450;
              A[ 9]:=   -0.00000 00205   3034;
              A[10]:=   -0.00000 00438   9154;
              A[11]:=   -0.00000 00176   8401;
              A[12]:=   -0.00000 00039   9129;
              A[13]:=   -0.00000 00001   8693;
              A[14]:=    0.00000 00002   7292;
              A[15]:=    0.00000 00001   3282;
              A[16]:=    0.00000 00000   3183;
              A[17]:=    0.00000 00000   0167;
              A[18]:=   -0.00000 00000   0204;
              A[19]:=   -0.00000 00000   0097;
              A[20]:=   -0.00000 00000   0022;
              A[21]:=   -0.00000 00000   0001;
              A[22]:=    0.00000 00000   0001;
              A[23]:=    0.00000 00000   0001;
              BETAX:=   SQRT(- LN((1 +   ABSX) × (1 - ABSX)));
              P:= -0.55945 76313 29832 × BETAX +
                  2.28791 57162 6336;
              P:= CHEPOLSER(23, P, A);
              INVERF:= if X < 0 then - BETAX × P
                  else BETAX × P
          end
     end INVERSE ERROR FUNCTION;
     EPS:= 10-14;
     if PROB < EPS ∨ 1 - PROB < EPS then
     STATAL3 ERROR(“PHINV”, 1, PROB);
    PHINV:= INVERF(2 × PROB - 1) × 1.41421356237310
end PHINV;
eop