MINICERN/mathlib/gen/c/algama64.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_IBM)
  11 #if defined(CERNLIB_DOUBLE)
  12       FUNCTION DLGAMA(X)
  13 C
  14 #include "gen/imp64.inc"
  15 C
  16       CHARACTER*(*) NAME
  17       PARAMETER(NAME='ALGAMA/DLGAMA')
  18 #endif
  19 #if !defined(CERNLIB_DOUBLE)
  20       FUNCTION ALGAMA(X)
  21 C
  22       CHARACTER*(*) NAME
  23       PARAMETER(NAME='ALGAMA')
  24 #endif
  25 C
  26       DIMENSION P1(7),Q1(7),P2(7),Q2(7),P3(7),Q3(7),C(5)
  27
  28       PARAMETER (Z1 = 1, HF = Z1/2, HF1 = 1+HF)
  29       CHARACTER*80 ERRTXT
  30       DATA P1
  31      1/+3.84287 36567 45991D+0, +5.27068 93753 00983D+1,
  32      2 +5.55840 45723 51531D+1, -2.15135 13573 72570D+2,
  33      3 -2.45872 61722 29242D+2, -5.75008 93603 04123D+1,
  34      4 -2.33590 98949 51284D+0/
  35       DATA Q1
  36      1/+1.00000 00000 00000D+0, +3.37330 47907 07074D+1,
  37      2 +1.93877 84034 37713D+2, +3.08829 54973 42428D+2,
  38      3 +1.50068 39064 89095D+2, +2.01068 51344 33395D+1,
  39      4 +4.57174 20282 50299D-1/
  40       DATA P2
  41      1/+4.87402 01396 83863 6D+0, +2.48845 25168 57407 6D+2,
  42      2 +2.17973 66058 89591 5D+3, +3.79751 24011 52511 8D+3,
  43      3 -1.97780 70769 84164 6D+3, -3.69298 34005 59128 2D+3,
  44      4 -5.60177 73537 80387 7D+2/
  45       DATA Q2
  46      1/+1.00000 00000 00000 0D+0, +9.50999 17418 20893 8D+1,
  47      2 +1.56120 45277 92863 5D+3, +7.23400 87928 94807 1D+3,
  48      3 +1.04595 76594 05895 9D+4, +4.16994 15153 20023 1D+3,
  49      4 +2.76785 83623 80410 1D+2/
  50       DATA P3
  51      1/-6.88062 40094 59425D+3, -4.30699 69819 57098D+5,
  52      2 -4.75045 94653 43956D+6, -2.94234 45930 32234D+6,
  53      3 +3.63218 04931 54257D+7, -3.35677 82814 54576D+6,
  54      4 -2.48043 69488 28593D+7/
  55       DATA Q3
  56      1/+1.00000 00000 00000D+0, -1.42168 29839 65146D+3,
  57      2 -1.55528 90280 85353D+5, -3.41525 17108 01107D+6,
  58      3 -2.09696 23255 80444D+7, -3.45441 75093 34395D+7,
  59      4 -9.16055 82863 71317D+6/
  60       DATA C
  61      1/ 1.12249 21356 561D-1,  7.95916 92961 204D-2,
  62      1 -1.70877 94611 020D-3,  9.18938 53320 467D-1,
  63      2  1.34699 05627 879D+0/
  64
  65
  66 #if defined(CERNLIB_DOUBLE)
  67       ENTRY DLOGAM(X)
  68 #endif
  69 #if !defined(CERNLIB_DOUBLE)
  70       ENTRY ALOGAM(X)
  71 #endif
  72
  73       IF(X .LE. 0) THEN
  74        H=0
  75        WRITE(ERRTXT,101) X
  76        CALL MTLPRT(NAME,'C304.1',ERRTXT)
  77       ELSE IF(X .EQ. 1 .OR. X .EQ. 2) THEN
  78        H=0
  79       ELSE IF(X .LE. HF) THEN
  80        Y=X+1
  81        AP=P1(1)
  82        AQ=Q1(1)
  83        DO 2 I = 2,7
  84        AP=P1(I)+Y*AP
  85     2  AQ=Q1(I)+Y*AQ
  86        H=-LOG(X)+X*AP/AQ
  87       ELSE IF(X .LE. HF1) THEN
  88        AP=P1(1)
  89        AQ=Q1(1)
  90        DO 3 I = 2,7
  91        AP=P1(I)+X*AP
  92     3  AQ=Q1(I)+X*AQ
  93        H=(X-1)*AP/AQ
  94       ELSE IF(X .LE. 4) THEN
  95        AP=P2(1)
  96        AQ=Q2(1)
  97        DO 4 I = 2,7
  98        AP=P2(I)+X*AP
  99     4  AQ=Q2(I)+X*AQ
 100        H=(X-2)*AP/AQ
 101       ELSE IF(X .LE. 12) THEN
 102        AP=P3(1)
 103        AQ=Q3(1)
 104        DO 5 I = 2,7
 105        AP=P3(I)+X*AP
 106     5  AQ=Q3(I)+X*AQ
 107        H=AP/AQ
 108       ELSE
 109        Y=1/X**2
 110        H=(X-HF)*LOG(X)-X+C(4)+(C(1)+Y*(C(2)+Y*C(3)))/
 111      1                                        ((C(5)+Y)*X)
 112       ENDIF
 113 #if defined(CERNLIB_DOUBLE)
 114       DLGAMA=H
 115 #endif
 116 #if !defined(CERNLIB_DOUBLE)
 117       ALGAMA=H
 118 #endif
 119       RETURN
 120   101 FORMAT('NON-POSITIVE ARGUMENT  X = ',1P,E15.6)
 121       END
 122 #endif