MINICERN/mathlib/gen/c/gammf64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:54  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION DGAMMF(X)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15       CHARACTER*(*) NAME
  16       PARAMETER(NAME='GAMMF/DGAMMF')
  17 #endif
  18 #if !defined(CERNLIB_DOUBLE)
  19       FUNCTION GAMMF(X)
  20 C
  21       CHARACTER*(*) NAME
  22       PARAMETER(NAME='GAMMF')
  23 #endif
  24       CHARACTER*80 ERRTXT
  25 C
  26       DIMENSION C(0:15)
  27
  28       PARAMETER (PI = 3.14159 26535 89793 24D0)
  29
  30       DATA C( 0) /3.65738 77250 83382 44D0/
  31       DATA C( 1) /1.95754 34566 61268 27D0/
  32       DATA C( 2) /0.33829 71138 26160 39D0/
  33       DATA C( 3) /0.04208 95127 65575 49D0/
  34       DATA C( 4) /0.00428 76504 82129 09D0/
  35       DATA C( 5) /0.00036 52121 69294 62D0/
  36       DATA C( 6) /0.00002 74006 42226 42D0/
  37       DATA C( 7) /0.00000 18124 02333 65D0/
  38       DATA C( 8) /0.00000 01096 57758 66D0/
  39       DATA C( 9) /0.00000 00059 87184 05D0/
  40       DATA C(10) /0.00000 00003 07690 81D0/
  41       DATA C(11) /0.00000 00000 14317 93D0/
  42       DATA C(12) /0.00000 00000 00651 09D0/
  43       DATA C(13) /0.00000 00000 00025 96D0/
  44       DATA C(14) /0.00000 00000 00001 11D0/
  45       DATA C(15) /0.00000 00000 00000 04D0/
  46
  47       U=X
  48       IF(U .LE. 0) THEN
  49        IF(U .EQ. INT(X)) THEN
  50         WRITE(ERRTXT,101) U
  51         CALL MTLPRT(NAME,'C303.1',ERRTXT)
  52         H=0
  53         GO TO 9
  54        ELSE
  55         U=1-X
  56        END IF
  57       ENDIF
  58     8 F=1
  59       IF(U .LT. 3) THEN
  60        DO 1 I = 1,INT(4-U)
  61        F=F/U
  62     1  U=U+1
  63       ELSE
  64        DO 2 I = 1,INT(U-3)
  65        U=U-1
  66     2  F=F*U
  67       END IF
  68       H=U+U-7
  69       ALFA=H+H
  70       B1=0
  71       B2=0
  72       DO 3 I = 15,0,-1
  73       B0=C(I)+ALFA*B1-B2
  74       B2=B1
  75     3 B1=B0
  76       H=F*(B0-H*B2)
  77       IF(X .LT. 0) H=PI/(SIN(PI*X)*H)
  78 #if defined(CERNLIB_DOUBLE)
  79     9 DGAMMF=H
  80 #endif
  81 #if !defined(CERNLIB_DOUBLE)
  82     9  GAMMF=H
  83 #endif
  84       RETURN
  85
  86   101 FORMAT('ARGUMENT IS NON-POSITIVE INTEGER = ',1P,E15.1)
  87       END