]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:01:56 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | #if defined(CERNLIB_DOUBLE) | |
11 | FUNCTION C309R2(X,ETA,ZL,P,EPS,LIMIT,KIND,ERR,NITS, | |
12 | 1 FPMAX,ACC8,ACC16) | |
13 | C | |
14 | C *** evaluate the HYPERGEOMETRIC FUNCTION | |
15 | C i | |
16 | C F (AA;BB; Z) = SUM (AA) Z / ( (BB) i! ) | |
17 | C 1 1 i i i | |
18 | C | |
19 | C to accuracy EPS with at most LIMIT terms. | |
20 | C If KIND = 0 : using extended precision but real arithmetic only, | |
21 | C 1 : using normal precision in complex arithmetic, | |
22 | C or 2 : using normal complex arithmetic, but with WDIGAM factor | |
23 | C | |
24 | C where AA, BB, and Z are defined below | |
25 | C | |
26 | IMPLICIT DOUBLE PRECISION(A-H,O-Z) | |
27 | COMPLEX*16 X,ETA,ZL,P,AA,BB,Z,C309R2,WDIGAM | |
28 | COMPLEX*16 DD,G,F,AI,BI,T | |
29 | #if (!defined(CERNLIB_UNIX))&&(!defined(CERNLIB_QMALPH)) | |
30 | REAL*16 AR,BR,GR,GI,DR,DI,TR,TI,UR,UI,FI,FI1,DEN | |
31 | #endif | |
32 | #if defined(CERNLIB_UNIX)||defined(CERNLIB_QMALPH) | |
33 | DOUBLE PRECISION AR,BR,GR,GI,DR,DI,TR,TI,UR,UI,FI,FI1,DEN | |
34 | #endif | |
35 | ||
36 | PARAMETER(TBBB = 3D0/2D0) | |
37 | ||
38 | #if defined(CERNLIB_QF2C) | |
39 | #include "defdr.inc" | |
40 | #endif | |
41 | ||
42 | ABSC(AA)=ABS(DREAL(AA))+ABS(DIMAG(AA)) | |
43 | NINTC(AA)=NINT(DREAL(AA)) | |
44 | ||
45 | C309R2 = 0. | |
46 | ||
47 | AA=ZL+1-ETA*P | |
48 | BB=2*ZL+2 | |
49 | Z=2*P*X | |
50 | IF(DREAL(BB) .LE. 0 .AND. ABS(BB-NINTC(BB)) .LT. | |
51 | 1 SQRT(SQRT(ACC8))**3 .AND. DREAL(BB)+LIMIT .GE. TBBB) THEN | |
52 | NITS=-1 | |
53 | RETURN | |
54 | END IF | |
55 | IF(LIMIT .LE. 0) THEN | |
56 | C309R2=0 | |
57 | ERR=0 | |
58 | NITS=1 | |
59 | RETURN | |
60 | END IF | |
61 | TA=1 | |
62 | RK=1 | |
63 | IF(KIND .LE. 0 .AND. ABSC(Z)*ABSC(AA) .GT. ABSC(BB)*1.0) THEN | |
64 | DR=1 | |
65 | DI=0 | |
66 | GR=1 | |
67 | GI=0 | |
68 | AR=DREAL(AA) | |
69 | BR=DREAL(BB) | |
70 | FI=0 | |
71 | DO 20 I = 2,LIMIT | |
72 | FI1=FI+1 | |
73 | TR=BR*FI1 | |
74 | TI=DIMAG(BB)*FI1 | |
75 | DEN=1/(TR*TR+TI*TI) | |
76 | UR=(AR*TR+DIMAG(AA)*TI)*DEN | |
77 | UI=(DIMAG(AA)*TR-AR*TI)*DEN | |
78 | TR=UR*GR-UI*GI | |
79 | TI=UR*GI+UI*GR | |
80 | GR=DREAL(Z)*TR-DIMAG(Z)*TI | |
81 | GI=DREAL(Z)*TI+DIMAG(Z)*TR | |
82 | DR=DR+GR | |
83 | DI=DI+GI | |
84 | ERR=ABS(GR)+ABS(GI) | |
85 | IF(ERR .GT. FPMAX) GO TO 60 | |
86 | RK=ABS(DR)+ABS(DI) | |
87 | TA=MAX(TA,RK) | |
88 | IF(ERR .LT. RK*EPS .OR. I .GE. 4 .AND. ERR .LT. ACC16) GO TO 30 | |
89 | FI=FI1 | |
90 | AR=AR+1 | |
91 | 20 BR=BR+1 | |
92 | C | |
93 | 30 C309R2=DR+(0,1)*DI | |
94 | ERR=ACC16*TA/RK | |
95 | ELSE | |
96 | C | |
97 | C* If REAL*16 arithmetic is not available, (or already using it!), | |
98 | C* then use KIND > 0 | |
99 | C | |
100 | G=1 | |
101 | F=1 | |
102 | IF(KIND .GE. 2) F=WDIGAM(AA)-WDIGAM(BB)-WDIGAM(G) | |
103 | DD=F | |
104 | DO 40 I = 2,LIMIT | |
105 | AI=AA+(I-2) | |
106 | BI=BB+(I-2) | |
107 | R=I-1 | |
108 | G=G*Z*AI/(BI*R) | |
109 | C | |
110 | C multiply by (psi(a+r)-psi(b+r)-psi(1+r)) | |
111 | C | |
112 | IF(KIND .EQ. 2) F=F+1/AI-1/BI-1/R | |
113 | T=G*F | |
114 | DD=DD+T | |
115 | ERR=ABSC(T) | |
116 | IF(ERR .GT. FPMAX) GO TO 60 | |
117 | RK=ABSC(DD) | |
118 | TA=MAX(TA,RK) | |
119 | IF(ERR .LT. RK*EPS .OR. ERR .LT. ACC8 .AND. I .GE. 4) GO TO 50 | |
120 | 40 CONTINUE | |
121 | ||
122 | 50 ERR=ACC8*TA/RK | |
123 | C309R2=DD | |
124 | END IF | |
125 | 60 NITS=I | |
126 | RETURN | |
127 | END | |
128 | #endif |