[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / gapnc64.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:05  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
#if !defined(CERNLIB_DOUBLE)
      FUNCTION RGAPNC(A,X)
#endif
#if defined(CERNLIB_DOUBLE)
      FUNCTION DGAPNC(A,X)
#include "gen/imp64.inc"
#endif
C     Calculates the incomplete gamma function P(A,X) as defined by
C     GSTAR(A,X) in Ref. 1. Based on
C     1. W. Gautschi, ALGORITHM 542 Incomplete Gamma Functions,
C        ACM Trans. Math. Software 5 (1979) 482-489
C     2. W. Gautschi, A computational procedure for incomplete gamma
C        functions, ACM Trans. Math. Software 5 (1979) 466-481

      CHARACTER NAME*(*)
      CHARACTER*80 ERRTXT
#if !defined(CERNLIB_DOUBLE)
      PARAMETER (NAME = 'RGAPNC')
#endif
#if defined(CERNLIB_DOUBLE)
      PARAMETER (NAME = 'RGAPNC/DGAPNC')
#endif
      PARAMETER (EPS = 5D-14)
      PARAMETER (ALH = -0.69314 71805 59945 31D0)
      PARAMETER (Z1 = 1, HALF = Z1/2, QUAR = Z1/4)
      PARAMETER (C1 = 3*Z1/2, KMAX = 600, EPS1 = EPS/100)

      DIMENSION C(25)

      DATA C( 1) / 0.57721 56649 01532 86D0/
      DATA C( 2) /-0.65587 80715 20253 88D0/
      DATA C( 3) /-0.04200 26350 34095 24D0/
      DATA C( 4) / 0.16653 86113 82291 49D0/
      DATA C( 5) /-0.04219 77345 55544 34D0/
      DATA C( 6) /-0.00962 19715 27876 97D0/
      DATA C( 7) / 0.00721 89432 46663 10D0/
      DATA C( 8) /-0.00116 51675 91859 07D0/
      DATA C( 9) /-0.00021 52416 74114 95D0/
      DATA C(10) / 0.00012 80502 82388 12D0/
      DATA C(11) /-0.00002 01348 54780 79D0/
      DATA C(12) /-0.00000 12504 93482 14D0/
      DATA C(13) / 0.00000 11330 27231 98D0/
      DATA C(14) /-0.00000 02056 33841 70D0/
      DATA C(15) / 0.00000 00061 16095 10D0/
      DATA C(16) / 0.00000 00050 02007 64D0/
      DATA C(17) /-0.00000 00011 81274 57D0/
      DATA C(18) / 0.00000 00001 04342 67D0/
      DATA C(19) / 0.00000 00000 07782 26D0/
      DATA C(20) /-0.00000 00000 03696 81D0/
      DATA C(21) / 0.00000 00000 00510 04D0/
      DATA C(22) /-0.00000 00000 00020 58D0/
      DATA C(23) /-0.00000 00000 00005 35D0/
      DATA C(24) / 0.00000 00000 00001 23D0/
      DATA C(25) /-0.00000 00000 00000 12D0/

#if !defined(CERNLIB_DOUBLE)
      GLGAMA(V)=ALGAMA(V)
#endif
#if defined(CERNLIB_DOUBLE)
      GLGAMA(V)=DLGAMA(V)
#endif

      HST=0
      IF(X .LT. 0) THEN
       WRITE(ERRTXT,101) X
       CALL MTLPRT(NAME,'C334.1',ERRTXT)
       GO TO 99
      ELSEIF(X .EQ. 0) THEN
       IF(A .LT. 0) THEN
        WRITE(ERRTXT,102) A
        CALL MTLPRT(NAME,'C334.2',ERRTXT)
       ELSEIF(A .EQ. 0) THEN
        HST=1
       ENDIF
       GO TO 99
      ELSE
       ALX=LOG(X)
      ENDIF
      IF(X .LT. QUAR) THEN
       ALFA=ALH/ALX
      ELSE
       ALFA=X+QUAR
      ENDIF
      MA=HALF-A
      AEPS=A+MA

      IF(MA .GT. 0) THEN
       IF(AEPS .NE. 0) THEN
        SG=(-1)**(MA-1)*SIGN(Z1,A)*SIGN(Z1,AEPS)
        ALGP1=GLGAMA(1+AEPS)-LOG(ABS(AEPS))
        IF(MA .NE. 1) ALGP1=ALGP1+GLGAMA(1-AEPS)-GLGAMA(MA-AEPS)
       ELSE
        SG=0
        ALGP1=0
       ENDIF
      ELSE
       SG=SIGN(Z1,A)
       ALGP1=GLGAMA(1+A)
      ENDIF
      IF(A .GT. ALFA) THEN
       TERM=1
       SUM=1
       DO 1 K = 1,KMAX
       TERM=X*TERM/(A+K)
       SUM=SUM+TERM
       IF(ABS(TERM) .LE. EPS*SUM) GO TO 2
    1  CONTINUE
       GO TO 98
    2  HST=EXP(A*ALX-X+LOG(SUM)-ALGP1)
      ELSEIF(X .GT. C1) THEN
       P=0
       S=1-A
       Q=(X+S)*(X-1-A)
       R=4*(X+S)
       TERM=1
       SUM=1
       RHO=0
       DO 3 K = 2,KMAX
       P=P+S
       Q=Q+R
       R=R+8
       S=S+2
       T=P*(1+RHO)
       RHO=T/(Q-T)
       TERM=RHO*TERM
       SUM=SUM+TERM
       IF(ABS(TERM) .LE. EPS*SUM) GO TO 4
    3  CONTINUE
       GO TO 98
    4  IF(A .LT. 0) THEN
        HST=1
        IF(MA .LT. 0 .OR. AEPS .NE. 0)
     1   HST=1-SG*EXP(LOG(ABS(A)*SUM/(X+1-A))-X+A*ALX-ALGP1)
       ELSEIF(A .EQ. 0) THEN
        HST=1
       ELSE
        HST=1-EXP(A*ALX-X+LOG(A*SUM/(X+1-A))-ALGP1)
       ENDIF
      ELSE
       AE=A
       IF(A .LT. HALF) THEN
        IF(A .LT. -HALF) AE=AEPS
        SUM=C(25)
        DO 12 K = 24,1,-1
   12   SUM=AE*SUM+C(K)
        GA=-SUM/(1+AE*SUM)
        Y=AE*ALX
        IF(ABS(Y) .GE. 1) THEN
         U=GA-(EXP(Y)-1)/AE
        ELSE
         SUM=1
         TERM=1
         DO 7 K = 2,KMAX
         TERM=Y*TERM/K
         SUM=SUM+TERM
         IF(ABS(TERM) .LE. EPS1*SUM) GO TO 8
    7    CONTINUE
         GO TO 98
    8    U=GA-SUM*ALX
        ENDIF
       ELSE
        U=EXP(GLGAMA(A))-X**A/A
       ENDIF
       P=AE*X
       Q=AE+1
       R=AE+3
       TERM=1
       SUM=1
       DO 9 K = 2,KMAX
       P=P+X
       Q=Q+R
       R=R+2
       TERM=-P*TERM/Q
       SUM=SUM+TERM
       IF(ABS(TERM) .LE. EPS1*SUM) GO TO 10
    9  CONTINUE
       GO TO 98
   10  H=U+X**(AE+1)*SUM/(AE+1)
       IF(A .LT. -HALF) THEN
        H=H*EXP(X-AE*ALX)
        DO 13 J = 1,MA
   13   H=(1-X*H)/(J-AE)
        HST=1
        IF(MA .LT. 0 .OR. AEPS .NE. 0)
     1   HST=1-SG*EXP(LOG(ABS(A)*H)-X+A*ALX-ALGP1)
       ELSEIF(A .EQ. 0) THEN
        HST=1
       ELSE
        HST=1-A*H*EXP(-ALGP1)
       ENDIF
      ENDIF
#if defined(CERNLIB_DOUBLE)
   99 DGAPNC=HST
#endif
#if !defined(CERNLIB_DOUBLE)
   99 RGAPNC=HST
#endif
      RETURN

   98 WRITE(ERRTXT,103) A,X
      CALL MTLPRT(NAME,'C334.3',ERRTXT)
      GO TO 99
  101 FORMAT('ILLEGAL ARGUMENT X = ',1P,D15.8,' < 0')
  102 FORMAT('ILLEGAL ARGUMENTS  A = ',1P,D15.8,' < 0, X = 0')
  103 FORMAT('PROBLEMS WITH CONVERGENCE, A = ',1P,D15.8,' X = ',D15.8)
      END
Commit	Line	Data
fe4da5cc	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-SGEXP(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(AALX-X+LOG(ASUM/(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=HEXP(X-AEALX)
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-SGEXP(LOG(ABS(A)H)-X+A*ALX-ALGP1)
195	ELSEIF(A .EQ. 0) THEN
196	HST=1
197	ELSE
198	HST=1-AHEXP(-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