(*-------------------------------------------------------------------------*)
(*          Ninv2 -- Find Percentage Point of Normal Distribution          *)
(*-------------------------------------------------------------------------*)

FUNCTION Ninv2( P : REAL ) : REAL;

(*-------------------------------------------------------------------------*)
(*                                                                         *)
(*       Function:  Ninv2                                                  *)
(*                                                                         *)
(*       Purpose:   Finds percentage point of normal distribution          *)
(*                  (high accuracy)                                        *)
(*                                                                         *)
(*       Calling Sequence:                                                 *)
(*                                                                         *)
(*            Perc  := Ninv2( P : REAL ) : REAL;                           *)
(*                                                                         *)
(*                 P      --- input probability                            *)
(*                 Perc   --- resultant percentage point                   *)
(*                                                                         *)
(*       Calls:   Ninv                                                     *)
(*                                                                         *)
(*       Remark:  This method provides about 12-13 decimal digits of       *)
(*                accuracy.  The method used is to improve the             *)
(*                approximation produced by NINV using a Taylor series     *)
(*                on the approximation error.                              *)
(*                                                                         *)
(*-------------------------------------------------------------------------*)

VAR
   Xp:    REAL;
   P1:    REAL;
   Z:     REAL;
   X3:    REAL;
   X2:    REAL;
   X1:    REAL;
   Phi:   REAL;

BEGIN (* Ninv2 *)

   Xp     := Ninv( P );

   P1     := SigNorm( Xp );

   Phi    := SQRT( 1.0 / ( 2.0 * PI ) ) * EXP( -( Xp * Xp ) / 2.0 );

   Z      := ( P - P1 ) / Phi;

   X3     := ( 2.0 * ( Xp * Xp ) + 1.0 ) * Z / 3.0;
   X2     := ( X3 + Xp ) * Z / 2.0;
   X1     := ( ( X2 + 1.0 ) * Z );

   Ninv2  := Xp + X1;

END   (* Ninv2 *);