MINICERN/mathlib/gen/c/cdigam64.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 CDIGAM(Z)
  12 #include "gen/defc64.inc"
  13      +  CDIGAM
  14 #endif
  15 #if defined(CERNLIB_DOUBLE)
  16       FUNCTION WDIGAM(Z)
  17 #include "gen/imp64.inc"
  18 #include "gen/defc64.inc"
  19      +   WDIGAM
  20 #endif
  21 #include "gen/defc64.inc"
  22      + Z,U,V,H,R,P
  23       CHARACTER NAME*(*)
  24       CHARACTER*80 ERRTXT
  25 #if !defined(CERNLIB_DOUBLE)
  26       PARAMETER (NAME = 'CDIGAM')
  27 #endif
  28 #if defined(CERNLIB_DOUBLE)
  29       PARAMETER (NAME = 'CDIGAM/WDIGAM')
  30 #endif
  31       DIMENSION C(6)
  32
  33       PARAMETER (Z1 = 1, HF = Z1/2)
  34       PARAMETER (PI = 3.14159 26535 89793 24D0)
  35
  36 #include "gen/gcmpfun.inc"
  37 CSEQ,GCMPLX.
  38
  39       DATA C(1) / 8.33333 33333 33333 33D-2/
  40       DATA C(2) /-8.33333 33333 33333 33D-3/
  41       DATA C(3) / 3.96825 39682 53968 25D-3/
  42       DATA C(4) /-4.16666 66666 66666 67D-3/
  43       DATA C(5) / 7.57575 75757 57575 76D-3/
  44       DATA C(6) /-2.10927 96092 79609 28D-2/
  45
  46       U=Z
  47       X=U
  48       A=ABS(X)
  49       IF(GIMAG(U) .EQ. 0 .AND. -A .EQ. INT(X)) THEN
  50        H=0
  51        WRITE(ERRTXT,101) X
  52        CALL MTLPRT(NAME,'C307.1',ERRTXT)
  53       ELSE
  54        IF(X .LT. 0) U=-U
  55        V=U
  56        H=0
  57        IF(A .LT. 15) THEN
  58         H=1/V
  59         DO 1 I = 1,14-INT(A)
  60         V=V+1
  61     1   H=H+1/V
  62         V=V+1
  63        END IF
  64        R=1/V**2
  65        P=R*C(1)
  66        DO 2 I = 6,1,-1
  67     2  P=R*(C(I)+P)
  68        H=LOG(V)-HF/V-P-H
  69        IF(X .LT. 0) THEN
  70         V=PI*U
  71         X=V
  72         A=SIN(X)
  73         X=COS(X)
  74         Y=TANH(GIMAG(V))
  75         H=H+1/U+PI*GCMPLX(X,-A*Y)/GCMPLX(A,X*Y)
  76        END IF
  77       ENDIF
  78 #if defined(CERNLIB_DOUBLE)
  79       WDIGAM=H
  80 #endif
  81 #if !defined(CERNLIB_DOUBLE)
  82       CDIGAM=H
  83 #endif
  84       RETURN
  85   101 FORMAT(1X,'ARGUMENT EQUALS NON-POSITIVE INTEGER = ',1P,E15.1)
  86       END