← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/special/cdflib/alngam.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:
Lib/special/cdflib/alngam.f
1
 
      DOUBLE PRECISION FUNCTION alngam(x)
2
 
C**********************************************************************
3
 
C
4
 
C     DOUBLE PRECISION FUNCTION ALNGAM(X)
5
 
C                 double precision LN of the GAMma function
6
 
C
7
 
C
8
 
C                              Function
9
 
C
10
 
C
11
 
C     Returns the natural logarithm of GAMMA(X).
12
 
C
13
 
C
14
 
C                              Arguments
15
 
C
16
 
C
17
 
C     X --> value at which scaled log gamma is to be returned
18
 
C                    X is DOUBLE PRECISION
19
 
C
20
 
C
21
 
C                              Method
22
 
C
23
 
C
24
 
C     If X .le. 6.0, then use recursion to get X below 3
25
 
C     then apply rational approximation number 5236 of
26
 
C     Hart et al, Computer Approximations, John Wiley and
27
 
C     Sons, NY, 1968.
28
 
C
29
 
C     If X .gt. 6.0, then use recursion to get X to at least 12 and
30
 
C     then use formula 5423 of the same source.
31
 
C
32
 
C**********************************************************************
33
 
C
34
 
C     .. Parameters ..
35
 
      DOUBLE PRECISION hln2pi
36
 
      PARAMETER (hln2pi=0.91893853320467274178D0)
37
 
C     ..
38
 
C     .. Scalar Arguments ..
39
 
      DOUBLE PRECISION x
40
 
C     ..
41
 
C     .. Local Scalars ..
42
 
      DOUBLE PRECISION offset,prod,xx
43
 
      INTEGER i,n
44
 
C     ..
45
 
C     .. Local Arrays ..
46
 
      DOUBLE PRECISION coef(5),scoefd(4),scoefn(9)
47
 
C     ..
48
 
C     .. External Functions ..
49
 
      DOUBLE PRECISION devlpl
50
 
      EXTERNAL devlpl
51
 
C     ..
52
 
C     .. Intrinsic Functions ..
53
 
      INTRINSIC log,dble,int
54
 
C     ..
55
 
C     .. Data statements ..
56
 
      DATA scoefn(1)/0.62003838007127258804D2/,
57
 
     +     scoefn(2)/0.36036772530024836321D2/,
58
 
     +     scoefn(3)/0.20782472531792126786D2/,
59
 
     +     scoefn(4)/0.6338067999387272343D1/,
60
 
     +     scoefn(5)/0.215994312846059073D1/,
61
 
     +     scoefn(6)/0.3980671310203570498D0/,
62
 
     +     scoefn(7)/0.1093115956710439502D0/,
63
 
     +     scoefn(8)/0.92381945590275995D-2/,
64
 
     +     scoefn(9)/0.29737866448101651D-2/
65
 
      DATA scoefd(1)/0.62003838007126989331D2/,
66
 
     +     scoefd(2)/0.9822521104713994894D1/,
67
 
     +     scoefd(3)/-0.8906016659497461257D1/,
68
 
     +     scoefd(4)/0.1000000000000000000D1/
69
 
      DATA coef(1)/0.83333333333333023564D-1/,
70
 
     +     coef(2)/-0.27777777768818808D-2/,
71
 
     +     coef(3)/0.79365006754279D-3/,coef(4)/-0.594997310889D-3/,
72
 
     +     coef(5)/0.8065880899D-3/
73
 
C     ..
74
 
C     .. Executable Statements ..
75
 
      IF (.NOT. (x.LE.6.0D0)) GO TO 70
76
 
      prod = 1.0D0
77
 
      xx = x
78
 
      IF (.NOT. (x.GT.3.0D0)) GO TO 30
79
 
   10 IF (.NOT. (xx.GT.3.0D0)) GO TO 20
80
 
      xx = xx - 1.0D0
81
 
      prod = prod*xx
82
 
      GO TO 10
83
 
 
84
 
   20 CONTINUE
85
 
   30 IF (.NOT. (x.LT.2.0D0)) GO TO 60
86
 
   40 IF (.NOT. (xx.LT.2.0D0)) GO TO 50
87
 
      prod = prod/xx
88
 
      xx = xx + 1.0D0
89
 
      GO TO 40
90
 
 
91
 
   50 CONTINUE
92
 
   60 alngam = devlpl(scoefn,9,xx-2.0D0)/devlpl(scoefd,4,xx-2.0D0)
93
 
C
94
 
C
95
 
C     COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
96
 
C
97
 
C
98
 
      alngam = alngam*prod
99
 
      alngam = log(alngam)
100
 
      GO TO 110
101
 
 
102
 
   70 offset = hln2pi
103
 
C
104
 
C
105
 
C     IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
106
 
C
107
 
C
108
 
      n = int(12.0D0-x)
109
 
      IF (.NOT. (n.GT.0)) GO TO 90
110
 
      prod = 1.0D0
111
 
      DO 80,i = 1,n
112
 
          prod = prod* (x+dble(i-1))
113
 
   80 CONTINUE
114
 
      offset = offset - log(prod)
115
 
      xx = x + dble(n)
116
 
      GO TO 100
117
 
 
118
 
   90 xx = x
119
 
C
120
 
C
121
 
C     COMPUTE POWER SERIES
122
 
C
123
 
C
124
 
  100 alngam = devlpl(coef,5,1.0D0/ (xx**2))/xx
125
 
      alngam = alngam + offset + (xx-0.5D0)*log(xx) - xx
126
 
  110 RETURN
127
 
 
128
 
      END