GNU Octave: libcruft/amos/gamln.f Source File

Go to the documentation of this file.
00001       FUNCTION GAMLN(Z,IERR)
00002 C***BEGIN PROLOGUE  GAMLN
00003 C***DATE WRITTEN   830501   (YYMMDD)
00004 C***REVISION DATE  830501   (YYMMDD)
00005 C***CATEGORY NO.  B5F
00006 C***KEYWORDS  GAMMA FUNCTION,LOGARITHM OF GAMMA FUNCTION
00007 C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
00008 C***PURPOSE  TO COMPUTE THE LOGARITHM OF THE GAMMA FUNCTION
00009 C***DESCRIPTION
00010 C
00011 C         GAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR
00012 C         Z.GT.0.  THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES
00013 C         GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION
00014 C         G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN.  THE FUNCTION WAS MADE AS
00015 C         PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE
00016 C         10 DIGITS IN A WORD, RLN=AMAX1(-ALOG10(R1MACH(4)),0.5E-18)
00017 C         LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY.
00018 C
00019 C         SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100
00020 C         VALUES IS USED FOR SPEED OF EXECUTION.
00021 C
00022 C     DESCRIPTION OF ARGUMENTS
00023 C
00024 C         INPUT
00025 C           Z      - REAL ARGUMENT, Z.GT.0.0E0
00026 C
00027 C         OUTPUT
00028 C           GAMLN  - NATURAL LOG OF THE GAMMA FUNCTION AT Z
00029 C           IERR   - ERROR FLAG
00030 C                    IERR=0, NORMAL RETURN, COMPUTATION COMPLETED
00031 C                    IERR=1, Z.LE.0.0E0,    NO COMPUTATION
00032 C
00033 C***REFERENCES  COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
00034 C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
00035 C***ROUTINES CALLED  I1MACH,R1MACH
00036 C***END PROLOGUE  GAMLN
00037 C
00038       INTEGER I, I1M, K, MZ, NZ, IERR, I1MACH
00039       REAL CF, CON, FLN, FZ, GLN, RLN, S, TLG, TRM, TST, T1, WDTOL, Z,
00040      * ZDMY, ZINC, ZM, ZMIN, ZP, ZSQ
00041       REAL R1MACH
00042       DIMENSION CF(22), GLN(100)
00043 C           LNGAMMA(N), N=1,100
00044       DATA GLN(1), GLN(2), GLN(3), GLN(4), GLN(5), GLN(6), GLN(7),
00045      1     GLN(8), GLN(9), GLN(10), GLN(11), GLN(12), GLN(13), GLN(14),
00046      2     GLN(15), GLN(16), GLN(17), GLN(18), GLN(19), GLN(20),
00047      3     GLN(21), GLN(22)/
00048      4     0.00000000000000000E+00,     0.00000000000000000E+00,
00049      5     6.93147180559945309E-01,     1.79175946922805500E+00,
00050      6     3.17805383034794562E+00,     4.78749174278204599E+00,
00051      7     6.57925121201010100E+00,     8.52516136106541430E+00,
00052      8     1.06046029027452502E+01,     1.28018274800814696E+01,
00053      9     1.51044125730755153E+01,     1.75023078458738858E+01,
00054      A     1.99872144956618861E+01,     2.25521638531234229E+01,
00055      B     2.51912211827386815E+01,     2.78992713838408916E+01,
00056      C     3.06718601060806728E+01,     3.35050734501368889E+01,
00057      D     3.63954452080330536E+01,     3.93398841871994940E+01,
00058      E     4.23356164607534850E+01,     4.53801388984769080E+01/
00059       DATA GLN(23), GLN(24), GLN(25), GLN(26), GLN(27), GLN(28),
00060      1     GLN(29), GLN(30), GLN(31), GLN(32), GLN(33), GLN(34),
00061      2     GLN(35), GLN(36), GLN(37), GLN(38), GLN(39), GLN(40),
00062      3     GLN(41), GLN(42), GLN(43), GLN(44)/
00063      4     4.84711813518352239E+01,     5.16066755677643736E+01,
00064      5     5.47847293981123192E+01,     5.80036052229805199E+01,
00065      6     6.12617017610020020E+01,     6.45575386270063311E+01,
00066      7     6.78897431371815350E+01,     7.12570389671680090E+01,
00067      8     7.46582363488301644E+01,     7.80922235533153106E+01,
00068      9     8.15579594561150372E+01,     8.50544670175815174E+01,
00069      A     8.85808275421976788E+01,     9.21361756036870925E+01,
00070      B     9.57196945421432025E+01,     9.93306124547874269E+01,
00071      C     1.02968198614513813E+02,     1.06631760260643459E+02,
00072      D     1.10320639714757395E+02,     1.14034211781461703E+02,
00073      E     1.17771881399745072E+02,     1.21533081515438634E+02/
00074       DATA GLN(45), GLN(46), GLN(47), GLN(48), GLN(49), GLN(50),
00075      1     GLN(51), GLN(52), GLN(53), GLN(54), GLN(55), GLN(56),
00076      2     GLN(57), GLN(58), GLN(59), GLN(60), GLN(61), GLN(62),
00077      3     GLN(63), GLN(64), GLN(65), GLN(66)/
00078      4     1.25317271149356895E+02,     1.29123933639127215E+02,
00079      5     1.32952575035616310E+02,     1.36802722637326368E+02,
00080      6     1.40673923648234259E+02,     1.44565743946344886E+02,
00081      7     1.48477766951773032E+02,     1.52409592584497358E+02,
00082      8     1.56360836303078785E+02,     1.60331128216630907E+02,
00083      9     1.64320112263195181E+02,     1.68327445448427652E+02,
00084      A     1.72352797139162802E+02,     1.76395848406997352E+02,
00085      B     1.80456291417543771E+02,     1.84533828861449491E+02,
00086      C     1.88628173423671591E+02,     1.92739047287844902E+02,
00087      D     1.96866181672889994E+02,     2.01009316399281527E+02,
00088      E     2.05168199482641199E+02,     2.09342586752536836E+02/
00089       DATA GLN(67), GLN(68), GLN(69), GLN(70), GLN(71), GLN(72),
00090      1     GLN(73), GLN(74), GLN(75), GLN(76), GLN(77), GLN(78),
00091      2     GLN(79), GLN(80), GLN(81), GLN(82), GLN(83), GLN(84),
00092      3     GLN(85), GLN(86), GLN(87), GLN(88)/
00093      4     2.13532241494563261E+02,     2.17736934113954227E+02,
00094      5     2.21956441819130334E+02,     2.26190548323727593E+02,
00095      6     2.30439043565776952E+02,     2.34701723442818268E+02,
00096      7     2.38978389561834323E+02,     2.43268849002982714E+02,
00097      8     2.47572914096186884E+02,     2.51890402209723194E+02,
00098      9     2.56221135550009525E+02,     2.60564940971863209E+02,
00099      A     2.64921649798552801E+02,     2.69291097651019823E+02,
00100      B     2.73673124285693704E+02,     2.78067573440366143E+02,
00101      C     2.82474292687630396E+02,     2.86893133295426994E+02,
00102      D     2.91323950094270308E+02,     2.95766601350760624E+02,
00103      E     3.00220948647014132E+02,     3.04686856765668715E+02/
00104       DATA GLN(89), GLN(90), GLN(91), GLN(92), GLN(93), GLN(94),
00105      1     GLN(95), GLN(96), GLN(97), GLN(98), GLN(99), GLN(100)/
00106      2     3.09164193580146922E+02,     3.13652829949879062E+02,
00107      3     3.18152639620209327E+02,     3.22663499126726177E+02,
00108      4     3.27185287703775217E+02,     3.31717887196928473E+02,
00109      5     3.36261181979198477E+02,     3.40815058870799018E+02,
00110      6     3.45379407062266854E+02,     3.49954118040770237E+02,
00111      7     3.54539085519440809E+02,     3.59134205369575399E+02/
00112 C             COEFFICIENTS OF ASYMPTOTIC EXPANSION
00113       DATA CF(1), CF(2), CF(3), CF(4), CF(5), CF(6), CF(7), CF(8),
00114      1     CF(9), CF(10), CF(11), CF(12), CF(13), CF(14), CF(15),
00115      2     CF(16), CF(17), CF(18), CF(19), CF(20), CF(21), CF(22)/
00116      3     8.33333333333333333E-02,    -2.77777777777777778E-03,
00117      4     7.93650793650793651E-04,    -5.95238095238095238E-04,
00118      5     8.41750841750841751E-04,    -1.91752691752691753E-03,
00119      6     6.41025641025641026E-03,    -2.95506535947712418E-02,
00120      7     1.79644372368830573E-01,    -1.39243221690590112E+00,
00121      8     1.34028640441683920E+01,    -1.56848284626002017E+02,
00122      9     2.19310333333333333E+03,    -3.61087712537249894E+04,
00123      A     6.91472268851313067E+05,    -1.52382215394074162E+07,
00124      B     3.82900751391414141E+08,    -1.08822660357843911E+10,
00125      C     3.47320283765002252E+11,    -1.23696021422692745E+13,
00126      D     4.88788064793079335E+14,    -2.13203339609193739E+16/
00127 C
00128 C             LN(2*PI)
00129       DATA CON                    /     1.83787706640934548E+00/
00130 C
00131 C***FIRST EXECUTABLE STATEMENT  GAMLN
00132       IERR=0
00133       IF (Z.LE.0.0E0) GO TO 70
00134       IF (Z.GT.101.0E0) GO TO 10
00135       NZ = INT(Z)
00136       FZ = Z - FLOAT(NZ)
00137       IF (FZ.GT.0.0E0) GO TO 10
00138       IF (NZ.GT.100) GO TO 10
00139       GAMLN = GLN(NZ)
00140       RETURN
00141    10 CONTINUE
00142       WDTOL = R1MACH(4)
00143       WDTOL = AMAX1(WDTOL,0.5E-18)
00144       I1M = I1MACH(11)
00145       RLN = R1MACH(5)*FLOAT(I1M)
00146       FLN = AMIN1(RLN,20.0E0)
00147       FLN = AMAX1(FLN,3.0E0)
00148       FLN = FLN - 3.0E0
00149       ZM = 1.8000E0 + 0.3875E0*FLN
00150       MZ = INT(ZM) + 1
00151       ZMIN = FLOAT(MZ)
00152       ZDMY = Z
00153       ZINC = 0.0E0
00154       IF (Z.GE.ZMIN) GO TO 20
00155       ZINC = ZMIN - FLOAT(NZ)
00156       ZDMY = Z + ZINC
00157    20 CONTINUE
00158       ZP = 1.0E0/ZDMY
00159       T1 = CF(1)*ZP
00160       S = T1
00161       IF (ZP.LT.WDTOL) GO TO 40
00162       ZSQ = ZP*ZP
00163       TST = T1*WDTOL
00164       DO 30 K=2,22
00165         ZP = ZP*ZSQ
00166         TRM = CF(K)*ZP
00167         IF (ABS(TRM).LT.TST) GO TO 40
00168         S = S + TRM
00169    30 CONTINUE
00170    40 CONTINUE
00171       IF (ZINC.NE.0.0E0) GO TO 50
00172       TLG = ALOG(Z)
00173       GAMLN = Z*(TLG-1.0E0) + 0.5E0*(CON-TLG) + S
00174       RETURN
00175    50 CONTINUE
00176       ZP = 1.0E0
00177       NZ = INT(ZINC)
00178       DO 60 I=1,NZ
00179         ZP = ZP*(Z+FLOAT(I-1))
00180    60 CONTINUE
00181       TLG = ALOG(ZDMY)
00182       GAMLN = ZDMY*(TLG-1.0E0) - ALOG(ZP) + 0.5E0*(CON-TLG) + S
00183       RETURN
00184 C
00185 C
00186    70 CONTINUE
00187       IERR=1
00188       RETURN
00189       END