1CS    REAL FUNCTION ALGAMA(X)
2      DOUBLE PRECISION FUNCTION DLGAMA(X)
3C----------------------------------------------------------------------
4C
5C This routine calculates the LOG(GAMMA) function for a positive real
6C   argument X.  Computation is based on an algorithm outlined in
7C   references 1 and 2.  The program uses rational functions that
8C   theoretically approximate LOG(GAMMA) to at least 18 significant
9C   decimal digits.  The approximation for X > 12 is from reference
10C   3, while approximations for X < 12.0 are similar to those in
11C   reference 1, but are unpublished.  The accuracy achieved depends
12C   on the arithmetic system, the compiler, the intrinsic functions,
13C   and proper selection of the machine-dependent constants.
14C
15C
16C*********************************************************************
17C*********************************************************************
18C
19C Explanation of machine-dependent constants
20C
21C beta   - radix for the floating-point representation
22C maxexp - the smallest positive power of beta that overflows
23C XBIG   - largest argument for which LN(GAMMA(X)) is representable
24C          in the machine, i.e., the solution to the equation
25C                  LN(GAMMA(XBIG)) = beta**maxexp
26C XINF   - largest machine representable floating-point number;
27C          approximately beta**maxexp.
28C EPS    - The smallest positive floating-point number such that
29C          1.0+EPS .GT. 1.0
30C FRTBIG - Rough estimate of the fourth root of XBIG
31C
32C
33C     Approximate values for some important machines are:
34C
35C                           beta      maxexp         XBIG
36C
37C CRAY-1        (S.P.)        2        8191       9.62E+2461
38C Cyber 180/855
39C   under NOS   (S.P.)        2        1070       1.72E+319
40C IEEE (IBM/XT,
41C   SUN, etc.)  (S.P.)        2         128       4.08E+36
42C IEEE (IBM/XT,
43C   SUN, etc.)  (D.P.)        2        1024       2.55D+305
44C IBM 3033      (D.P.)       16          63       4.29D+73
45C VAX D-Format  (D.P.)        2         127       2.05D+36
46C VAX G-Format  (D.P.)        2        1023       1.28D+305
47C
48C
49C                           XINF        EPS        FRTBIG
50C
51C CRAY-1        (S.P.)   5.45E+2465   7.11E-15    3.13E+615
52C Cyber 180/855
53C   under NOS   (S.P.)   1.26E+322    3.55E-15    6.44E+79
54C IEEE (IBM/XT,
55C   SUN, etc.)  (S.P.)   3.40E+38     1.19E-7     1.42E+9
56C IEEE (IBM/XT,
57C   SUN, etc.)  (D.P.)   1.79D+308    2.22D-16    2.25D+76
58C IBM 3033      (D.P.)   7.23D+75     2.22D-16    2.56D+18
59C VAX D-Format  (D.P.)   1.70D+38     1.39D-17    1.20D+9
60C VAX G-Format  (D.P.)   8.98D+307    1.11D-16    1.89D+76
61C
62C**************************************************************
63C**************************************************************
64C
65C Error returns
66C
67C  The program returns the value XINF for X .LE. 0.0 or when
68C     overflow would occur.  The computation is believed to
69C     be free of underflow and overflow.
70C
71C
72C Intrinsic functions required are:
73C
74C      LOG
75C
76C
77C References:
78C
79C  1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for
80C     the Natural Logarithm of the Gamma Function,' Math. Comp. 21,
81C     1967, pp. 198-203.
82C
83C  2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May,
84C     1969.
85C
86C  3) Hart, Et. Al., Computer Approximations, Wiley and sons, New
87C     York, 1968.
88C
89C
90C  Authors: W. J. Cody and L. Stoltz
91C           Argonne National Laboratory
92C
93C  Latest modification: June 16, 1988
94C
95C----------------------------------------------------------------------
96*
97*  Some modifs from Bruno (25 Feb 2005):
98*    - return Nan (in place of XINF) if x <= 0
99*    - return Inf (in place of XINF) if x > XBIG
100*
101*  In fact an error indicator should be returned in these cases to prevent
102*  the user...
103*
104      INTEGER I
105CS    REAL
106      DOUBLE PRECISION
107     1    C,CORR,D1,D2,D4,EPS,FRTBIG,FOUR,HALF,ONE,PNT68,P1,P2,P4,
108     2    Q1,Q2,Q4,RES,SQRTPI,THRHAL,TWELVE,TWO,X,XBIG,XDEN,XINF,
109     3    XM1,XM2,XM4,XNUM,Y,YSQ,ZERO
110      DIMENSION C(7),P1(8),P2(8),P4(8),Q1(8),Q2(8),Q4(8)
111C----------------------------------------------------------------------
112C  Mathematical constants
113C----------------------------------------------------------------------
114CS    DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/,
115CS   1     FOUR,THRHAL,TWO,PNT68/4.0E0,1.5E0,2.0E0,0.6796875E0/,
116CS   2     SQRTPI/0.9189385332046727417803297E0/
117      DATA ONE,HALF,TWELVE,ZERO/1.0D0,0.5D0,12.0D0,0.0D0/,
118     1     FOUR,THRHAL,TWO,PNT68/4.0D0,1.5D0,2.0D0,0.6796875D0/,
119     2     SQRTPI/0.9189385332046727417803297D0/
120C----------------------------------------------------------------------
121C  Machine dependent parameters
122C----------------------------------------------------------------------
123CS    DATA XBIG,XINF,EPS,FRTBIG/4.08E36,3.401E38,1.19E-7,1.42E9/
124      DATA XBIG,XINF,EPS,FRTBIG/2.55D305,1.79D308,2.22D-16,2.25D76/
125C----------------------------------------------------------------------
126C  Numerator and denominator coefficients for rational minimax
127C     approximation over (0.5,1.5).
128C----------------------------------------------------------------------
129CS    DATA D1/-5.772156649015328605195174E-1/
130CS    DATA P1/4.945235359296727046734888E0,2.018112620856775083915565E2,
131CS   1        2.290838373831346393026739E3,1.131967205903380828685045E4,
132CS   2        2.855724635671635335736389E4,3.848496228443793359990269E4,
133CS   3        2.637748787624195437963534E4,7.225813979700288197698961E3/
134CS    DATA Q1/6.748212550303777196073036E1,1.113332393857199323513008E3,
135CS   1        7.738757056935398733233834E3,2.763987074403340708898585E4,
136CS   2        5.499310206226157329794414E4,6.161122180066002127833352E4,
137CS   3        3.635127591501940507276287E4,8.785536302431013170870835E3/
138      DATA D1/-5.772156649015328605195174D-1/
139      DATA P1/4.945235359296727046734888D0,2.018112620856775083915565D2,
140     1        2.290838373831346393026739D3,1.131967205903380828685045D4,
141     2        2.855724635671635335736389D4,3.848496228443793359990269D4,
142     3        2.637748787624195437963534D4,7.225813979700288197698961D3/
143      DATA Q1/6.748212550303777196073036D1,1.113332393857199323513008D3,
144     1        7.738757056935398733233834D3,2.763987074403340708898585D4,
145     2        5.499310206226157329794414D4,6.161122180066002127833352D4,
146     3        3.635127591501940507276287D4,8.785536302431013170870835D3/
147C----------------------------------------------------------------------
148C  Numerator and denominator coefficients for rational minimax
149C     Approximation over (1.5,4.0).
150C----------------------------------------------------------------------
151CS    DATA D2/4.227843350984671393993777E-1/
152CS    DATA P2/4.974607845568932035012064E0,5.424138599891070494101986E2,
153CS   1        1.550693864978364947665077E4,1.847932904445632425417223E5,
154CS   2        1.088204769468828767498470E6,3.338152967987029735917223E6,
155CS   3        5.106661678927352456275255E6,3.074109054850539556250927E6/
156CS    DATA Q2/1.830328399370592604055942E2,7.765049321445005871323047E3,
157CS   1        1.331903827966074194402448E5,1.136705821321969608938755E6,
158CS   2        5.267964117437946917577538E6,1.346701454311101692290052E7,
159CS   3        1.782736530353274213975932E7,9.533095591844353613395747E6/
160      DATA D2/4.227843350984671393993777D-1/
161      DATA P2/4.974607845568932035012064D0,5.424138599891070494101986D2,
162     1        1.550693864978364947665077D4,1.847932904445632425417223D5,
163     2        1.088204769468828767498470D6,3.338152967987029735917223D6,
164     3        5.106661678927352456275255D6,3.074109054850539556250927D6/
165      DATA Q2/1.830328399370592604055942D2,7.765049321445005871323047D3,
166     1        1.331903827966074194402448D5,1.136705821321969608938755D6,
167     2        5.267964117437946917577538D6,1.346701454311101692290052D7,
168     3        1.782736530353274213975932D7,9.533095591844353613395747D6/
169C----------------------------------------------------------------------
170C  Numerator and denominator coefficients for rational minimax
171C     Approximation over (4.0,12.0).
172C----------------------------------------------------------------------
173CS    DATA D4/1.791759469228055000094023E0/
174CS    DATA P4/1.474502166059939948905062E4,2.426813369486704502836312E6,
175CS   1        1.214755574045093227939592E8,2.663432449630976949898078E9,
176CS   2      2.940378956634553899906876E10,1.702665737765398868392998E11,
177CS   3      4.926125793377430887588120E11,5.606251856223951465078242E11/
178CS    DATA Q4/2.690530175870899333379843E3,6.393885654300092398984238E5,
179CS   2        4.135599930241388052042842E7,1.120872109616147941376570E9,
180CS   3      1.488613728678813811542398E10,1.016803586272438228077304E11,
181CS   4      3.417476345507377132798597E11,4.463158187419713286462081E11/
182      DATA D4/1.791759469228055000094023D0/
183      DATA P4/1.474502166059939948905062D4,2.426813369486704502836312D6,
184     1        1.214755574045093227939592D8,2.663432449630976949898078D9,
185     2      2.940378956634553899906876D10,1.702665737765398868392998D11,
186     3      4.926125793377430887588120D11,5.606251856223951465078242D11/
187      DATA Q4/2.690530175870899333379843D3,6.393885654300092398984238D5,
188     2        4.135599930241388052042842D7,1.120872109616147941376570D9,
189     3      1.488613728678813811542398D10,1.016803586272438228077304D11,
190     4      3.417476345507377132798597D11,4.463158187419713286462081D11/
191C----------------------------------------------------------------------
192C  Coefficients for minimax approximation over (12, INF).
193C----------------------------------------------------------------------
194CS    DATA C/-1.910444077728E-03,8.4171387781295E-04,
195CS   1     -5.952379913043012E-04,7.93650793500350248E-04,
196CS   2     -2.777777777777681622553E-03,8.333333333333333331554247E-02,
197CS   3      5.7083835261E-03/
198      DATA C/-1.910444077728D-03,8.4171387781295D-04,
199     1     -5.952379913043012D-04,7.93650793500350248D-04,
200     2     -2.777777777777681622553D-03,8.333333333333333331554247D-02,
201     3      5.7083835261D-03/
202C----------------------------------------------------------------------
203      Y = X
204      IF ((Y .GT. ZERO) .AND. (Y .LE. XBIG)) THEN
205            IF (Y .LE. EPS) THEN
206                  RES = -LOG(Y)
207               ELSE IF (Y .LE. THRHAL) THEN
208C----------------------------------------------------------------------
209C  EPS .LT. X .LE. 1.5
210C----------------------------------------------------------------------
211                  IF (Y .LT. PNT68) THEN
212                        CORR = -LOG(Y)
213                        XM1 = Y
214                     ELSE
215                        CORR = ZERO
216                        XM1 = (Y - HALF) - HALF
217                  END IF
218                  IF ((Y .LE. HALF) .OR. (Y .GE. PNT68)) THEN
219                        XDEN = ONE
220                        XNUM = ZERO
221                        DO 140 I = 1, 8
222                           XNUM = XNUM*XM1 + P1(I)
223                           XDEN = XDEN*XM1 + Q1(I)
224  140                   CONTINUE
225                        RES = CORR + (XM1 * (D1 + XM1*(XNUM/XDEN)))
226                     ELSE
227                        XM2 = (Y - HALF) - HALF
228                        XDEN = ONE
229                        XNUM = ZERO
230                        DO 220 I = 1, 8
231                           XNUM = XNUM*XM2 + P2(I)
232                           XDEN = XDEN*XM2 + Q2(I)
233  220                   CONTINUE
234                        RES = CORR + XM2 * (D2 + XM2*(XNUM/XDEN))
235                  END IF
236               ELSE IF (Y .LE. FOUR) THEN
237C----------------------------------------------------------------------
238C  1.5 .LT. X .LE. 4.0
239C----------------------------------------------------------------------
240                  XM2 = Y - TWO
241                  XDEN = ONE
242                  XNUM = ZERO
243                  DO 240 I = 1, 8
244                     XNUM = XNUM*XM2 + P2(I)
245                     XDEN = XDEN*XM2 + Q2(I)
246  240             CONTINUE
247                  RES = XM2 * (D2 + XM2*(XNUM/XDEN))
248               ELSE IF (Y .LE. TWELVE) THEN
249C----------------------------------------------------------------------
250C  4.0 .LT. X .LE. 12.0
251C----------------------------------------------------------------------
252                  XM4 = Y - FOUR
253                  XDEN = -ONE
254                  XNUM = ZERO
255                  DO 340 I = 1, 8
256                     XNUM = XNUM*XM4 + P4(I)
257                     XDEN = XDEN*XM4 + Q4(I)
258  340             CONTINUE
259                  RES = D4 + XM4*(XNUM/XDEN)
260               ELSE
261C----------------------------------------------------------------------
262C  Evaluate for argument .GE. 12.0,
263C----------------------------------------------------------------------
264                  RES = ZERO
265                  IF (Y .LE. FRTBIG) THEN
266                        RES = C(7)
267                        YSQ = Y * Y
268                        DO 450 I = 1, 6
269                           RES = RES / YSQ + C(I)
270  450                   CONTINUE
271                  END IF
272                  RES = RES/Y
273                  CORR = LOG(Y)
274                  RES = RES + SQRTPI - HALF*CORR
275                  RES = RES + Y*(CORR-ONE)
276            END IF
277         ELSE
278C----------------------------------------------------------------------
279C  Return for bad arguments
280C----------------------------------------------------------------------
281*           modif from Bruno (see comment at the beginning)
282*            RES = XINF
283            if (X .le. 0.d0) then
284              CALL returnananfortran(RES)
285            else
286C this means that X > XBIG and so that log(gamma) overflows
287C bad trick to get Inf
288               RES = 2*XINF
289            endif
290*           end modif
291      END IF
292C----------------------------------------------------------------------
293C  Final adjustments and return
294C----------------------------------------------------------------------
295CS    ALGAMA = RES
296      DLGAMA = RES
297      RETURN
298C ---------- Last line of DLGAMA ----------
299      END
300
301
302