MINICERN/mathlib/gen/c/r2sp.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_SINGLE)
  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 CDIGAM factor
  23 C
  24 C  where AA, BB, and Z are defined below
  25 C
  26       IMPLICIT REAL(A-H,O-Z)
  27       COMPLEX X,ETA,ZL,P,AA,BB,Z,C309R2,CDIGAM
  28       COMPLEX DD,G,F,AI,BI,T
  29       DOUBLE PRECISION AR,BR,GR,GI,DR,DI,TR,TI,UR,UI,FI,FI1,DEN
  30
  31       PARAMETER(TBBB = 3D0/2D0)
  32
  33       ABSC(AA)=ABS(REAL(AA))+ABS(AIMAG(AA))
  34       NINTC(AA)=NINT(REAL(AA))
  35
  36       AA=ZL+1-ETA*P
  37       BB=2*ZL+2
  38       Z=2*P*X
  39       IF(REAL(BB) .LE. 0 .AND. ABS(BB-NINTC(BB)) .LT.
  40      1 SQRT(SQRT(ACC8))**3 .AND. REAL(BB)+LIMIT .GE. TBBB) THEN
  41        NITS=-1
  42        RETURN
  43       END IF
  44       IF(LIMIT .LE. 0) THEN
  45        C309R2=0
  46        ERR=0
  47        NITS=1
  48        RETURN
  49       END IF
  50       TA=1
  51       RK=1
  52       IF(KIND .LE. 0 .AND. ABSC(Z)*ABSC(AA) .GT. ABSC(BB)*1.0) THEN
  53        DR=1
  54        DI=0
  55        GR=1
  56        GI=0
  57        AR=REAL(AA)
  58        BR=REAL(BB)
  59        FI=0
  60        DO 20 I = 2,LIMIT
  61        FI1=FI+1
  62        TR=BR*FI1
  63        TI=AIMAG(BB)*FI1
  64        DEN=1/(TR*TR+TI*TI)
  65        UR=(AR*TR+AIMAG(AA)*TI)*DEN
  66        UI=(AIMAG(AA)*TR-AR*TI)*DEN
  67        TR=UR*GR-UI*GI
  68        TI=UR*GI+UI*GR
  69        GR=REAL(Z)*TR-AIMAG(Z)*TI
  70        GI=REAL(Z)*TI+AIMAG(Z)*TR
  71        DR=DR+GR
  72        DI=DI+GI
  73        ERR=ABS(GR)+ABS(GI)
  74        IF(ERR .GT. FPMAX) GO TO 60
  75        RK=ABS(DR)+ABS(DI)
  76        TA=MAX(TA,RK)
  77        IF(ERR .LT. RK*EPS .OR. I .GE. 4 .AND. ERR .LT. ACC16) GO TO 30
  78        FI=FI1
  79        AR=AR+1
  80    20  BR=BR+1
  81 C
  82    30  C309R2=DR+(0,1)*DI
  83        ERR=ACC16*TA/RK
  84       ELSE
  85 C
  86 C*    If REAL*16 arithmetic is not available, (or already using it!),
  87 C*    then use KIND > 0
  88 C
  89        G=1
  90        F=1
  91        IF(KIND .GE. 2) F=CDIGAM(AA)-CDIGAM(BB)-CDIGAM(G)
  92        DD=F
  93        DO 40 I = 2,LIMIT
  94        AI=AA+(I-2)
  95        BI=BB+(I-2)
  96        R=I-1
  97        G=G*Z*AI/(BI*R)
  98 C
  99 C                       multiply by (psi(a+r)-psi(b+r)-psi(1+r))
 100 C
 101        IF(KIND .EQ. 2) F=F+1/AI-1/BI-1/R
 102        T=G*F
 103        DD=DD+T
 104        ERR=ABSC(T)
 105        IF(ERR .GT. FPMAX) GO TO 60
 106        RK=ABSC(DD)
 107        TA=MAX(TA,RK)
 108        IF(ERR .LT. RK*EPS .OR. ERR .LT. ACC8 .AND. I .GE. 4) GO TO 50
 109    40  CONTINUE
 110
 111    50  ERR=ACC8*TA/RK
 112        C309R2=DD
 113       END IF
 114    60 NITS=I
 115       RETURN
 116       END
 117 #endif