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