(*--------------------------------------------------------------------------*)
(*                 Cinv --- Inverse chi-square routine                      *)
(*--------------------------------------------------------------------------*)

FUNCTION Cinv( P, V: REAL; VAR Ifault: INTEGER ) : REAL;

(*--------------------------------------------------------------------------*)
(*                                                                          *)
(*      Function:  Cinv                                                     *)
(*                                                                          *)
(*      Purpose:  Compute inverse chi-square (percentage point)             *)
(*                                                                          *)
(*      Calling Sequence:                                                   *)
(*                                                                          *)
(*         Cp := Cinv( P, V: REAL; VAR Ifault: INTEGER ) : REAL;            *)
(*                                                                          *)
(*            P      --- tail probability to get percentage point for       *)
(*            V      --- degrees of freedom                                 *)
(*            Ifault --- error indicator                                    *)
(*                       = 0:  no error                                     *)
(*                       = 1:  probability too small to accurately compute  *)
(*                             percentage point.                            *)
(*                       = 2:  degrees of freedom given as <= 0.            *)
(*                       = 3:  chi-square probability evaluation failed.    *)
(*                                                                          *)
(*            Cp    --- resulting chi-square percentage point               *)
(*                                                                          *)
(*      Calls:    GammaIn                                                   *)
(*                Ninv                                                      *)
(*                ALGama                                                    *)
(*                                                                          *)
(*      Remarks:                                                            *)
(*                                                                          *)
(*         This is basically AS 91 from Applied Statistics.                 *)
(*         An initial approximation is improved using a Taylor series       *)
(*         expansion.                                                       *)
(*                                                                          *)
(*--------------------------------------------------------------------------*)

CONST
   E       = 1.0E-8;
   Dprec   = 8;
   MaxIter = 100;

VAR
   XX   : REAL;
   C    : REAL;
   Ch   : REAL;
   Q    : REAL;
   P1   : REAL;
   P2   : REAL;
   T    : REAL;
   X    : REAL;
   B    : REAL;
   A    : REAL;
   G    : REAL;
   S1   : REAL;
   S2   : REAL;
   S3   : REAL;
   S4   : REAL;
   S5   : REAL;
   S6   : REAL;
   Cprec: REAL;

   Iter : INTEGER;

LABEL 9000;

BEGIN (* Cinv *)
                                   (* Test arguments for validity *)
   Cinv   := -1.0;
   Ifault := 1;

   IF ( P < E ) OR ( P > ( 1.0 - E ) ) THEN GOTO 9000;

   Ifault := 2;
   IF( V <= 0.0 ) THEN GOTO 9000;

                                   (* Initialize *)
   P      := 1.0 - P;
   XX     := V / 2.0;
   G      := ALGama( XX );
   Ifault := 0;
   C      := XX - 1.0;
                                   (* Starting approx. for small chi-square *)

   IF( V < ( -1.24 * LN( P ) ) ) THEN
      BEGIN

         Ch := Power( P * XX * EXP( G + XX * LnTwo ) , ( 1.0 / XX ) );

         IF Ch < E THEN
            BEGIN
               Cinv := Ch;
               GOTO 9000;
            END;
      END
   ELSE                         (* Starting approx. for v <= .32 *)
      IF ( V <= 0.32 ) THEN
         BEGIN

            Ch := 0.4;
            A  := LN( 1.0 - P );

            REPEAT
               Q  := Ch;
               P1 := 1.0 + Ch * ( 4.67 + Ch );
               P2 := Ch * ( 6.73 + Ch * ( 6.66 + Ch ) );
               T  := -0.5 + ( 4.67 + 2.0 * Ch ) / P1 -
                     ( 6.73 + Ch * ( 13.32 + 3.0 * Ch ) ) / P2;
               Ch := Ch - ( 1.0 - EXP( A + G + 0.5 * Ch + C * LnTwo ) *
                            P2 / P1 ) / T;
            UNTIL( ABS( Q / Ch - 1.0 ) <= 0.01 );

         END
      ELSE
         BEGIN
                                   (* Starting approx. using Wilson and *)
                                   (* Hilferty estimate                 *)
            X  := Ninv( P );
            P1 := 2.0 / ( 9.0 * V );
            Ch := V * PowerI( ( X * SQRT( P1 ) + 1.0 - P1 ) , 3 );

                                   (* Starting approx. for P ---> 1     *)

            IF ( Ch > ( 2.2 * V + 6.0 ) ) THEN
               Ch := -2.0 * ( LN( 1.0 - P ) - C * LN( 0.5 * Ch ) + G );

         END;

                                   (* We have starting approximation.    *)
                                   (* Begin improvement loop.            *)
   REPEAT
                                   (* Get probability of current approx. *)
                                   (* to percentage point.               *)
      Q     := Ch;
      P1    := 0.5 * Ch;
      P2    := P - GammaIn( P1, XX, Dprec, MaxIter, Cprec, Iter, Ifault );

      IF( Ifault <> 0 ) OR ( Iter > MaxIter ) THEN
         Ifault := 3
      ELSE
         BEGIN
                                      (* Calculate seven-term Taylor series *)

            T  := P2 * EXP( XX * LnTwo + G + P1 - C * LN( Ch ) );
            B  := T / Ch;
            A  := 0.5 * T - B * C;

            S1 := ( 210.0 + A * ( 140.0 + A * ( 105.0 + A * ( 84.0 + A *
                  ( 70.0 + 60.0 * A ) ) ) ) ) / 420.0;
            S2 := ( 420.0 + A * ( 735.0 + A * ( 966.0 + A * ( 1141.0 +
                    1278.0 * A ) ) ) ) / 2520.0;
            S3 := ( 210.0 + A * ( 462.0 + A * ( 707.0 + 932.0 * A ) ) )
                  / 2520.0;
            S4 := ( 252.0 + A * ( 672.0 + 1182.0 * A ) + C * ( 294.0 + A *
                  ( 889.0 + 1740.0 * A ) ) ) / 5040.0;
            S5 := ( 84.0 + 264.0 * A + C * ( 175.0 + 606.0 * A ) ) / 2520.0;
            S6 := ( 120.0 + C * ( 346.0 + 127.0 * C ) ) / 5040.0;
            Ch := Ch + T * ( 1.0 + 0.5 * T * S1 - B * C * ( S1 - B *
                  ( S2 - B * ( S3 - B * ( S4 - B * ( S5 - B * S6 ) ) ) ) ) );

         END;

   UNTIL ( ABS( ( Q / Ch ) - 1.0 ) <= E ) OR ( Ifault <> 0 );

   IF Ifault = 0 THEN
      Cinv := Ch
   ELSE
      Cinv := -1.0;

9000:  ;

END   (* Cinv *);