MINICERN/mathlib/gen/c/gapnc64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:05  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       FUNCTION RGAPNC(A,X)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       FUNCTION DGAPNC(A,X)
  15 #include "gen/imp64.inc"
  16 #endif
  17 C     Calculates the incomplete gamma function P(A,X) as defined by
  18 C     GSTAR(A,X) in Ref. 1. Based on
  19 C     1. W. Gautschi, ALGORITHM 542 Incomplete Gamma Functions,
  20 C        ACM Trans. Math. Software 5 (1979) 482-489
  21 C     2. W. Gautschi, A computational procedure for incomplete gamma
  22 C        functions, ACM Trans. Math. Software 5 (1979) 466-481
  23
  24       CHARACTER NAME*(*)
  25       CHARACTER*80 ERRTXT
  26 #if !defined(CERNLIB_DOUBLE)
  27       PARAMETER (NAME = 'RGAPNC')
  28 #endif
  29 #if defined(CERNLIB_DOUBLE)
  30       PARAMETER (NAME = 'RGAPNC/DGAPNC')
  31 #endif
  32       PARAMETER (EPS = 5D-14)
  33       PARAMETER (ALH = -0.69314 71805 59945 31D0)
  34       PARAMETER (Z1 = 1, HALF = Z1/2, QUAR = Z1/4)
  35       PARAMETER (C1 = 3*Z1/2, KMAX = 600, EPS1 = EPS/100)
  36
  37       DIMENSION C(25)
  38
  39       DATA C( 1) / 0.57721 56649 01532 86D0/
  40       DATA C( 2) /-0.65587 80715 20253 88D0/
  41       DATA C( 3) /-0.04200 26350 34095 24D0/
  42       DATA C( 4) / 0.16653 86113 82291 49D0/
  43       DATA C( 5) /-0.04219 77345 55544 34D0/
  44       DATA C( 6) /-0.00962 19715 27876 97D0/
  45       DATA C( 7) / 0.00721 89432 46663 10D0/
  46       DATA C( 8) /-0.00116 51675 91859 07D0/
  47       DATA C( 9) /-0.00021 52416 74114 95D0/
  48       DATA C(10) / 0.00012 80502 82388 12D0/
  49       DATA C(11) /-0.00002 01348 54780 79D0/
  50       DATA C(12) /-0.00000 12504 93482 14D0/
  51       DATA C(13) / 0.00000 11330 27231 98D0/
  52       DATA C(14) /-0.00000 02056 33841 70D0/
  53       DATA C(15) / 0.00000 00061 16095 10D0/
  54       DATA C(16) / 0.00000 00050 02007 64D0/
  55       DATA C(17) /-0.00000 00011 81274 57D0/
  56       DATA C(18) / 0.00000 00001 04342 67D0/
  57       DATA C(19) / 0.00000 00000 07782 26D0/
  58       DATA C(20) /-0.00000 00000 03696 81D0/
  59       DATA C(21) / 0.00000 00000 00510 04D0/
  60       DATA C(22) /-0.00000 00000 00020 58D0/
  61       DATA C(23) /-0.00000 00000 00005 35D0/
  62       DATA C(24) / 0.00000 00000 00001 23D0/
  63       DATA C(25) /-0.00000 00000 00000 12D0/
  64
  65 #if !defined(CERNLIB_DOUBLE)
  66       GLGAMA(V)=ALGAMA(V)
  67 #endif
  68 #if defined(CERNLIB_DOUBLE)
  69       GLGAMA(V)=DLGAMA(V)
  70 #endif
  71
  72       HST=0
  73       IF(X .LT. 0) THEN
  74        WRITE(ERRTXT,101) X
  75        CALL MTLPRT(NAME,'C334.1',ERRTXT)
  76        GO TO 99
  77       ELSEIF(X .EQ. 0) THEN
  78        IF(A .LT. 0) THEN
  79         WRITE(ERRTXT,102) A
  80         CALL MTLPRT(NAME,'C334.2',ERRTXT)
  81        ELSEIF(A .EQ. 0) THEN
  82         HST=1
  83        ENDIF
  84        GO TO 99
  85       ELSE
  86        ALX=LOG(X)
  87       ENDIF
  88       IF(X .LT. QUAR) THEN
  89        ALFA=ALH/ALX
  90       ELSE
  91        ALFA=X+QUAR
  92       ENDIF
  93       MA=HALF-A
  94       AEPS=A+MA
  95
  96       IF(MA .GT. 0) THEN
  97        IF(AEPS .NE. 0) THEN
  98         SG=(-1)**(MA-1)*SIGN(Z1,A)*SIGN(Z1,AEPS)
  99         ALGP1=GLGAMA(1+AEPS)-LOG(ABS(AEPS))
 100         IF(MA .NE. 1) ALGP1=ALGP1+GLGAMA(1-AEPS)-GLGAMA(MA-AEPS)
 101        ELSE
 102         SG=0
 103         ALGP1=0
 104        ENDIF
 105       ELSE
 106        SG=SIGN(Z1,A)
 107        ALGP1=GLGAMA(1+A)
 108       ENDIF
 109       IF(A .GT. ALFA) THEN
 110        TERM=1
 111        SUM=1
 112        DO 1 K = 1,KMAX
 113        TERM=X*TERM/(A+K)
 114        SUM=SUM+TERM
 115        IF(ABS(TERM) .LE. EPS*SUM) GO TO 2
 116     1  CONTINUE
 117        GO TO 98
 118     2  HST=EXP(A*ALX-X+LOG(SUM)-ALGP1)
 119       ELSEIF(X .GT. C1) THEN
 120        P=0
 121        S=1-A
 122        Q=(X+S)*(X-1-A)
 123        R=4*(X+S)
 124        TERM=1
 125        SUM=1
 126        RHO=0
 127        DO 3 K = 2,KMAX
 128        P=P+S
 129        Q=Q+R
 130        R=R+8
 131        S=S+2
 132        T=P*(1+RHO)
 133        RHO=T/(Q-T)
 134        TERM=RHO*TERM
 135        SUM=SUM+TERM
 136        IF(ABS(TERM) .LE. EPS*SUM) GO TO 4
 137     3  CONTINUE
 138        GO TO 98
 139     4  IF(A .LT. 0) THEN
 140         HST=1
 141         IF(MA .LT. 0 .OR. AEPS .NE. 0)
 142      1   HST=1-SG*EXP(LOG(ABS(A)*SUM/(X+1-A))-X+A*ALX-ALGP1)
 143        ELSEIF(A .EQ. 0) THEN
 144         HST=1
 145        ELSE
 146         HST=1-EXP(A*ALX-X+LOG(A*SUM/(X+1-A))-ALGP1)
 147        ENDIF
 148       ELSE
 149        AE=A
 150        IF(A .LT. HALF) THEN
 151         IF(A .LT. -HALF) AE=AEPS
 152         SUM=C(25)
 153         DO 12 K = 24,1,-1
 154    12   SUM=AE*SUM+C(K)
 155         GA=-SUM/(1+AE*SUM)
 156         Y=AE*ALX
 157         IF(ABS(Y) .GE. 1) THEN
 158          U=GA-(EXP(Y)-1)/AE
 159         ELSE
 160          SUM=1
 161          TERM=1
 162          DO 7 K = 2,KMAX
 163          TERM=Y*TERM/K
 164          SUM=SUM+TERM
 165          IF(ABS(TERM) .LE. EPS1*SUM) GO TO 8
 166     7    CONTINUE
 167          GO TO 98
 168     8    U=GA-SUM*ALX
 169         ENDIF
 170        ELSE
 171         U=EXP(GLGAMA(A))-X**A/A
 172        ENDIF
 173        P=AE*X
 174        Q=AE+1
 175        R=AE+3
 176        TERM=1
 177        SUM=1
 178        DO 9 K = 2,KMAX
 179        P=P+X
 180        Q=Q+R
 181        R=R+2
 182        TERM=-P*TERM/Q
 183        SUM=SUM+TERM
 184        IF(ABS(TERM) .LE. EPS1*SUM) GO TO 10
 185     9  CONTINUE
 186        GO TO 98
 187    10  H=U+X**(AE+1)*SUM/(AE+1)
 188        IF(A .LT. -HALF) THEN
 189         H=H*EXP(X-AE*ALX)
 190         DO 13 J = 1,MA
 191    13   H=(1-X*H)/(J-AE)
 192         HST=1
 193         IF(MA .LT. 0 .OR. AEPS .NE. 0)
 194      1   HST=1-SG*EXP(LOG(ABS(A)*H)-X+A*ALX-ALGP1)
 195        ELSEIF(A .EQ. 0) THEN
 196         HST=1
 197        ELSE
 198         HST=1-A*H*EXP(-ALGP1)
 199        ENDIF
 200       ENDIF
 201 #if defined(CERNLIB_DOUBLE)
 202    99 DGAPNC=HST
 203 #endif
 204 #if !defined(CERNLIB_DOUBLE)
 205    99 RGAPNC=HST
 206 #endif
 207       RETURN
 208
 209    98 WRITE(ERRTXT,103) A,X
 210       CALL MTLPRT(NAME,'C334.3',ERRTXT)
 211       GO TO 99
 212   101 FORMAT('ILLEGAL ARGUMENT X = ',1P,D15.8,' < 0')
 213   102 FORMAT('ILLEGAL ARGUMENTS  A = ',1P,D15.8,' < 0, X = 0')
 214   103 FORMAT('PROBLEMS WITH CONVERGENCE, A = ',1P,D15.8,' X = ',D15.8)
 215       END