MINICERN/mathlib/gen/c/cgamma64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:55  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       FUNCTION CGAMMA(Z)
  12 #include "gen/defc64.inc"
  13      + CGAMMA
  14 #endif
  15 #if defined(CERNLIB_DOUBLE)
  16       FUNCTION WGAMMA(Z)
  17 #include "gen/imp64.inc"
  18 #include "gen/defc64.inc"
  19      +  WGAMMA
  20 #endif
  21 #include "gen/defc64.inc"
  22      +       Z,U,V,F,H,S
  23       CHARACTER NAME*(*)
  24       CHARACTER*80 ERRTXT
  25 #if !defined(CERNLIB_DOUBLE)
  26       PARAMETER (NAME = 'CGAMMA')
  27 #endif
  28 #if defined(CERNLIB_DOUBLE)
  29       PARAMETER (NAME = 'CGAMMA/WGAMMA')
  30 #endif
  31       DIMENSION C(0:15)
  32
  33       PARAMETER (Z1 = 1, HF = Z1/2)
  34
  35 #if defined(CERNLIB_QF2C)
  36 #include "gen/gcmpfun.inc"
  37 #endif
  38
  39       DATA PI /3.14159 26535 89793 24D0/
  40       DATA C1 /2.50662 82746 31000 50D0/
  41
  42       DATA C( 0) / 41.62443 69164 39068D0/
  43       DATA C( 1) /-51.22424 10223 74774D0/
  44       DATA C( 2) / 11.33875 58134 88977D0/
  45       DATA C( 3) / -0.74773 26877 72388D0/
  46       DATA C( 4) /  0.00878 28774 93061D0/
  47       DATA C( 5) / -0.00000 18990 30264D0/
  48       DATA C( 6) /  0.00000 00019 46335D0/
  49       DATA C( 7) / -0.00000 00001 99345D0/
  50       DATA C( 8) /  0.00000 00000 08433D0/
  51       DATA C( 9) /  0.00000 00000 01486D0/
  52       DATA C(10) / -0.00000 00000 00806D0/
  53       DATA C(11) /  0.00000 00000 00293D0/
  54       DATA C(12) / -0.00000 00000 00102D0/
  55       DATA C(13) /  0.00000 00000 00037D0/
  56       DATA C(14) / -0.00000 00000 00014D0/
  57       DATA C(15) /  0.00000 00000 00006D0/
  58
  59 #if !defined(CERNLIB_QF2C)
  60 #include "gen/gcmpfun.inc"
  61 #endif
  62
  63       U=Z
  64       X=U
  65       IF(GIMAG(U) .EQ. 0 .AND. -ABS(X) .EQ. INT(X)) THEN
  66        F=0
  67        H=0
  68        WRITE(ERRTXT,101) X
  69        CALL MTLPRT(NAME,'C305.1',ERRTXT)
  70       ELSE
  71        IF(X .GE. 1) THEN
  72         F=1
  73         V=U
  74        ELSEIF(X .GE. 0) THEN
  75         F=1/U
  76         V=1+U
  77        ELSE
  78         F=1
  79         V=1-U
  80        END IF
  81        H=1
  82        S=C(0)
  83        DO 1 K = 1,15
  84        H=((V-K)/(V+(K-1)))*H
  85     1  S=S+C(K)*H
  86        H=V+(4+HF)
  87        H=C1*EXP((V-HF)*LOG(H)-H)*S
  88        IF(X .LT. 0) H=PI/(SIN(PI*U)*H)
  89       ENDIF
  90 #if !defined(CERNLIB_DOUBLE)
  91       CGAMMA=F*H
  92 #endif
  93 #if defined(CERNLIB_DOUBLE)
  94       WGAMMA=F*H
  95 #endif
  96       RETURN
  97   101 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER = ',1P,E15.1)
  98       END