This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / r2dp.F
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