← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to scipy/special/amos/dgamln.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/special/amos/dgamln.f
 
1
      DOUBLE PRECISION FUNCTION DGAMLN(Z,IERR)
 
2
C***BEGIN PROLOGUE  DGAMLN
 
3
C***DATE WRITTEN   830501   (YYMMDD)
 
4
C***REVISION DATE  830501   (YYMMDD)
 
5
C***CATEGORY NO.  B5F
 
6
C***KEYWORDS  GAMMA FUNCTION,LOGARITHM OF GAMMA FUNCTION
 
7
C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
 
8
C***PURPOSE  TO COMPUTE THE LOGARITHM OF THE GAMMA FUNCTION
 
9
C***DESCRIPTION
 
10
C
 
11
C               **** A DOUBLE PRECISION ROUTINE ****
 
12
C         DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR
 
13
C         Z.GT.0.  THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES
 
14
C         GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION
 
15
C         G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN.  THE FUNCTION WAS MADE AS
 
16
C         PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE
 
17
C         10 DIGITS IN A WORD, RLN=AMAX1(-ALOG10(R1MACH(4)),0.5E-18)
 
18
C         LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY.
 
19
C
 
20
C         SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100
 
21
C         VALUES IS USED FOR SPEED OF EXECUTION.
 
22
C
 
23
C     DESCRIPTION OF ARGUMENTS
 
24
C
 
25
C         INPUT      Z IS D0UBLE PRECISION
 
26
C           Z      - ARGUMENT, Z.GT.0.0D0
 
27
C
 
28
C         OUTPUT      DGAMLN IS DOUBLE PRECISION
 
29
C           DGAMLN  - NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0
 
30
C           IERR    - ERROR FLAG
 
31
C                     IERR=0, NORMAL RETURN, COMPUTATION COMPLETED
 
32
C                     IERR=1, Z.LE.0.0D0,    NO COMPUTATION
 
33
C
 
34
C
 
35
C***REFERENCES  COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
 
36
C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
 
37
C***ROUTINES CALLED  I1MACH,D1MACH
 
38
C***END PROLOGUE  DGAMLN
 
39
      DOUBLE PRECISION CF, CON, FLN, FZ, GLN, RLN, S, TLG, TRM, TST,
 
40
     * T1, WDTOL, Z, ZDMY, ZINC, ZM, ZMIN, ZP, ZSQ, D1MACH
 
41
      INTEGER I, IERR, I1M, K, MZ, NZ, I1MACH
 
42
      DIMENSION CF(22), GLN(100)
 
43
C           LNGAMMA(N), N=1,100
 
44
      DATA GLN(1), GLN(2), GLN(3), GLN(4), GLN(5), GLN(6), GLN(7),
 
45
     1     GLN(8), GLN(9), GLN(10), GLN(11), GLN(12), GLN(13), GLN(14),
 
46
     2     GLN(15), GLN(16), GLN(17), GLN(18), GLN(19), GLN(20),
 
47
     3     GLN(21), GLN(22)/
 
48
     4     0.00000000000000000D+00,     0.00000000000000000D+00,
 
49
     5     6.93147180559945309D-01,     1.79175946922805500D+00,
 
50
     6     3.17805383034794562D+00,     4.78749174278204599D+00,
 
51
     7     6.57925121201010100D+00,     8.52516136106541430D+00,
 
52
     8     1.06046029027452502D+01,     1.28018274800814696D+01,
 
53
     9     1.51044125730755153D+01,     1.75023078458738858D+01,
 
54
     A     1.99872144956618861D+01,     2.25521638531234229D+01,
 
55
     B     2.51912211827386815D+01,     2.78992713838408916D+01,
 
56
     C     3.06718601060806728D+01,     3.35050734501368889D+01,
 
57
     D     3.63954452080330536D+01,     3.93398841871994940D+01,
 
58
     E     4.23356164607534850D+01,     4.53801388984769080D+01/
 
59
      DATA GLN(23), GLN(24), GLN(25), GLN(26), GLN(27), GLN(28),
 
60
     1     GLN(29), GLN(30), GLN(31), GLN(32), GLN(33), GLN(34),
 
61
     2     GLN(35), GLN(36), GLN(37), GLN(38), GLN(39), GLN(40),
 
62
     3     GLN(41), GLN(42), GLN(43), GLN(44)/
 
63
     4     4.84711813518352239D+01,     5.16066755677643736D+01,
 
64
     5     5.47847293981123192D+01,     5.80036052229805199D+01,
 
65
     6     6.12617017610020020D+01,     6.45575386270063311D+01,
 
66
     7     6.78897431371815350D+01,     7.12570389671680090D+01,
 
67
     8     7.46582363488301644D+01,     7.80922235533153106D+01,
 
68
     9     8.15579594561150372D+01,     8.50544670175815174D+01,
 
69
     A     8.85808275421976788D+01,     9.21361756036870925D+01,
 
70
     B     9.57196945421432025D+01,     9.93306124547874269D+01,
 
71
     C     1.02968198614513813D+02,     1.06631760260643459D+02,
 
72
     D     1.10320639714757395D+02,     1.14034211781461703D+02,
 
73
     E     1.17771881399745072D+02,     1.21533081515438634D+02/
 
74
      DATA GLN(45), GLN(46), GLN(47), GLN(48), GLN(49), GLN(50),
 
75
     1     GLN(51), GLN(52), GLN(53), GLN(54), GLN(55), GLN(56),
 
76
     2     GLN(57), GLN(58), GLN(59), GLN(60), GLN(61), GLN(62),
 
77
     3     GLN(63), GLN(64), GLN(65), GLN(66)/
 
78
     4     1.25317271149356895D+02,     1.29123933639127215D+02,
 
79
     5     1.32952575035616310D+02,     1.36802722637326368D+02,
 
80
     6     1.40673923648234259D+02,     1.44565743946344886D+02,
 
81
     7     1.48477766951773032D+02,     1.52409592584497358D+02,
 
82
     8     1.56360836303078785D+02,     1.60331128216630907D+02,
 
83
     9     1.64320112263195181D+02,     1.68327445448427652D+02,
 
84
     A     1.72352797139162802D+02,     1.76395848406997352D+02,
 
85
     B     1.80456291417543771D+02,     1.84533828861449491D+02,
 
86
     C     1.88628173423671591D+02,     1.92739047287844902D+02,
 
87
     D     1.96866181672889994D+02,     2.01009316399281527D+02,
 
88
     E     2.05168199482641199D+02,     2.09342586752536836D+02/
 
89
      DATA GLN(67), GLN(68), GLN(69), GLN(70), GLN(71), GLN(72),
 
90
     1     GLN(73), GLN(74), GLN(75), GLN(76), GLN(77), GLN(78),
 
91
     2     GLN(79), GLN(80), GLN(81), GLN(82), GLN(83), GLN(84),
 
92
     3     GLN(85), GLN(86), GLN(87), GLN(88)/
 
93
     4     2.13532241494563261D+02,     2.17736934113954227D+02,
 
94
     5     2.21956441819130334D+02,     2.26190548323727593D+02,
 
95
     6     2.30439043565776952D+02,     2.34701723442818268D+02,
 
96
     7     2.38978389561834323D+02,     2.43268849002982714D+02,
 
97
     8     2.47572914096186884D+02,     2.51890402209723194D+02,
 
98
     9     2.56221135550009525D+02,     2.60564940971863209D+02,
 
99
     A     2.64921649798552801D+02,     2.69291097651019823D+02,
 
100
     B     2.73673124285693704D+02,     2.78067573440366143D+02,
 
101
     C     2.82474292687630396D+02,     2.86893133295426994D+02,
 
102
     D     2.91323950094270308D+02,     2.95766601350760624D+02,
 
103
     E     3.00220948647014132D+02,     3.04686856765668715D+02/
 
104
      DATA GLN(89), GLN(90), GLN(91), GLN(92), GLN(93), GLN(94),
 
105
     1     GLN(95), GLN(96), GLN(97), GLN(98), GLN(99), GLN(100)/
 
106
     2     3.09164193580146922D+02,     3.13652829949879062D+02,
 
107
     3     3.18152639620209327D+02,     3.22663499126726177D+02,
 
108
     4     3.27185287703775217D+02,     3.31717887196928473D+02,
 
109
     5     3.36261181979198477D+02,     3.40815058870799018D+02,
 
110
     6     3.45379407062266854D+02,     3.49954118040770237D+02,
 
111
     7     3.54539085519440809D+02,     3.59134205369575399D+02/
 
112
C             COEFFICIENTS OF ASYMPTOTIC EXPANSION
 
113
      DATA CF(1), CF(2), CF(3), CF(4), CF(5), CF(6), CF(7), CF(8),
 
114
     1     CF(9), CF(10), CF(11), CF(12), CF(13), CF(14), CF(15),
 
115
     2     CF(16), CF(17), CF(18), CF(19), CF(20), CF(21), CF(22)/
 
116
     3     8.33333333333333333D-02,    -2.77777777777777778D-03,
 
117
     4     7.93650793650793651D-04,    -5.95238095238095238D-04,
 
118
     5     8.41750841750841751D-04,    -1.91752691752691753D-03,
 
119
     6     6.41025641025641026D-03,    -2.95506535947712418D-02,
 
120
     7     1.79644372368830573D-01,    -1.39243221690590112D+00,
 
121
     8     1.34028640441683920D+01,    -1.56848284626002017D+02,
 
122
     9     2.19310333333333333D+03,    -3.61087712537249894D+04,
 
123
     A     6.91472268851313067D+05,    -1.52382215394074162D+07,
 
124
     B     3.82900751391414141D+08,    -1.08822660357843911D+10,
 
125
     C     3.47320283765002252D+11,    -1.23696021422692745D+13,
 
126
     D     4.88788064793079335D+14,    -2.13203339609193739D+16/
 
127
C
 
128
C             LN(2*PI)
 
129
      DATA CON                    /     1.83787706640934548D+00/
 
130
C
 
131
C***FIRST EXECUTABLE STATEMENT  DGAMLN
 
132
      IERR=0
 
133
      IF (Z.LE.0.0D0) GO TO 70
 
134
      IF (Z.GT.101.0D0) GO TO 10
 
135
      NZ = INT(SNGL(Z))
 
136
      FZ = Z - FLOAT(NZ)
 
137
      IF (FZ.GT.0.0D0) GO TO 10
 
138
      IF (NZ.GT.100) GO TO 10
 
139
      DGAMLN = GLN(NZ)
 
140
      RETURN
 
141
   10 CONTINUE
 
142
      WDTOL = D1MACH(4)
 
143
      WDTOL = DMAX1(WDTOL,0.5D-18)
 
144
      I1M = I1MACH(14)
 
145
      RLN = D1MACH(5)*FLOAT(I1M)
 
146
      FLN = DMIN1(RLN,20.0D0)
 
147
      FLN = DMAX1(FLN,3.0D0)
 
148
      FLN = FLN - 3.0D0
 
149
      ZM = 1.8000D0 + 0.3875D0*FLN
 
150
      MZ = INT(SNGL(ZM)) + 1
 
151
      ZMIN = FLOAT(MZ)
 
152
      ZDMY = Z
 
153
      ZINC = 0.0D0
 
154
      IF (Z.GE.ZMIN) GO TO 20
 
155
      ZINC = ZMIN - FLOAT(NZ)
 
156
      ZDMY = Z + ZINC
 
157
   20 CONTINUE
 
158
      ZP = 1.0D0/ZDMY
 
159
      T1 = CF(1)*ZP
 
160
      S = T1
 
161
      IF (ZP.LT.WDTOL) GO TO 40
 
162
      ZSQ = ZP*ZP
 
163
      TST = T1*WDTOL
 
164
      DO 30 K=2,22
 
165
        ZP = ZP*ZSQ
 
166
        TRM = CF(K)*ZP
 
167
        IF (DABS(TRM).LT.TST) GO TO 40
 
168
        S = S + TRM
 
169
   30 CONTINUE
 
170
   40 CONTINUE
 
171
      IF (ZINC.NE.0.0D0) GO TO 50
 
172
      TLG = DLOG(Z)
 
173
      DGAMLN = Z*(TLG-1.0D0) + 0.5D0*(CON-TLG) + S
 
174
      RETURN
 
175
   50 CONTINUE
 
176
      ZP = 1.0D0
 
177
      NZ = INT(SNGL(ZINC))
 
178
      DO 60 I=1,NZ
 
179
        ZP = ZP*(Z+FLOAT(I-1))
 
180
   60 CONTINUE
 
181
      TLG = DLOG(ZDMY)
 
182
      DGAMLN = ZDMY*(TLG-1.0D0) - DLOG(ZP) + 0.5D0*(CON-TLG) + S
 
183
      RETURN
 
184
C
 
185
C
 
186
   70 CONTINUE
 
187
      IERR=1
 
188
      RETURN
 
189
      END