← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/special/cdflib/psi_fort.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cdflib/psi_fort.f
 
1
      DOUBLE PRECISION FUNCTION psi(xx)
 
2
C---------------------------------------------------------------------
 
3
C
 
4
C                 EVALUATION OF THE DIGAMMA FUNCTION
 
5
C
 
6
C                           -----------
 
7
C
 
8
C     PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
 
9
C     BE COMPUTED.
 
10
C
 
11
C     THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
 
12
C     APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
 
13
C     CODY, STRECOK AND THACHER.
 
14
C
 
15
C---------------------------------------------------------------------
 
16
C     PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
 
17
C     PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
 
18
C     A.H. MORRIS (NSWC).
 
19
C---------------------------------------------------------------------
 
20
C     .. Scalar Arguments ..
 
21
      DOUBLE PRECISION xx
 
22
C     ..
 
23
C     .. Local Scalars ..
 
24
      DOUBLE PRECISION aug,den,dx0,piov4,sgn,upper,w,x,xmax1,xmx0,
 
25
     +                 xsmall,z
 
26
      INTEGER i,m,n,nq
 
27
C     ..
 
28
C     .. Local Arrays ..
 
29
      DOUBLE PRECISION p1(7),p2(4),q1(6),q2(4)
 
30
C     ..
 
31
C     .. External Functions ..
 
32
      DOUBLE PRECISION spmpar
 
33
      INTEGER ipmpar
 
34
      EXTERNAL spmpar,ipmpar
 
35
C     ..
 
36
C     .. Intrinsic Functions ..
 
37
      INTRINSIC abs,cos,dble,dlog,dmin1,int,sin
 
38
C     ..
 
39
C     .. Data statements ..
 
40
C---------------------------------------------------------------------
 
41
C
 
42
C     PIOV4 = PI/4
 
43
C     DX0 = ZERO OF PSI TO EXTENDED PRECISION
 
44
C
 
45
C---------------------------------------------------------------------
 
46
C---------------------------------------------------------------------
 
47
C
 
48
C     COEFFICIENTS FOR RATIONAL APPROXIMATION OF
 
49
C     PSI(X) / (X - X0),  0.5 .LE. X .LE. 3.0
 
50
C
 
51
C---------------------------------------------------------------------
 
52
C---------------------------------------------------------------------
 
53
C
 
54
C     COEFFICIENTS FOR RATIONAL APPROXIMATION OF
 
55
C     PSI(X) - LN(X) + 1 / (2*X),  X .GT. 3.0
 
56
C
 
57
C---------------------------------------------------------------------
 
58
      DATA piov4/.785398163397448D0/
 
59
      DATA dx0/1.461632144968362341262659542325721325D0/
 
60
      DATA p1(1)/.895385022981970D-02/,p1(2)/.477762828042627D+01/,
 
61
     +     p1(3)/.142441585084029D+03/,p1(4)/.118645200713425D+04/,
 
62
     +     p1(5)/.363351846806499D+04/,p1(6)/.413810161269013D+04/,
 
63
     +     p1(7)/.130560269827897D+04/
 
64
      DATA q1(1)/.448452573429826D+02/,q1(2)/.520752771467162D+03/,
 
65
     +     q1(3)/.221000799247830D+04/,q1(4)/.364127349079381D+04/,
 
66
     +     q1(5)/.190831076596300D+04/,q1(6)/.691091682714533D-05/
 
67
      DATA p2(1)/-.212940445131011D+01/,p2(2)/-.701677227766759D+01/,
 
68
     +     p2(3)/-.448616543918019D+01/,p2(4)/-.648157123766197D+00/
 
69
      DATA q2(1)/.322703493791143D+02/,q2(2)/.892920700481861D+02/,
 
70
     +     q2(3)/.546117738103215D+02/,q2(4)/.777788548522962D+01/
 
71
C     ..
 
72
C     .. Executable Statements ..
 
73
C---------------------------------------------------------------------
 
74
C
 
75
C     MACHINE DEPENDENT CONSTANTS ...
 
76
C
 
77
C        XMAX1  = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
 
78
C                 WITH ENTIRELY INTEGER REPRESENTATION.  ALSO USED
 
79
C                 AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
 
80
C                 ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
 
81
C                 PSI MAY BE REPRESENTED AS ALOG(X).
 
82
C
 
83
C        XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
 
84
C                 MAY BE REPRESENTED BY 1/X.
 
85
C
 
86
C---------------------------------------------------------------------
 
87
      xmax1 = ipmpar(3)
 
88
      xmax1 = dmin1(xmax1,1.0D0/spmpar(1))
 
89
      xsmall = 1.D-9
 
90
C---------------------------------------------------------------------
 
91
      x = xx
 
92
      aug = 0.0D0
 
93
      IF (x.GE.0.5D0) GO TO 50
 
94
C---------------------------------------------------------------------
 
95
C     X .LT. 0.5,  USE REFLECTION FORMULA
 
96
C     PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
 
97
C---------------------------------------------------------------------
 
98
      IF (abs(x).GT.xsmall) GO TO 10
 
99
      IF (x.EQ.0.0D0) GO TO 100
 
100
C---------------------------------------------------------------------
 
101
C     0 .LT. ABS(X) .LE. XSMALL.  USE 1/X AS A SUBSTITUTE
 
102
C     FOR  PI*COTAN(PI*X)
 
103
C---------------------------------------------------------------------
 
104
      aug = -1.0D0/x
 
105
      GO TO 40
 
106
C---------------------------------------------------------------------
 
107
C     REDUCTION OF ARGUMENT FOR COTAN
 
108
C---------------------------------------------------------------------
 
109
   10 w = -x
 
110
      sgn = piov4
 
111
      IF (w.GT.0.0D0) GO TO 20
 
112
      w = -w
 
113
      sgn = -sgn
 
114
C---------------------------------------------------------------------
 
115
C     MAKE AN ERROR EXIT IF X .LE. -XMAX1
 
116
C---------------------------------------------------------------------
 
117
   20 IF (w.GE.xmax1) GO TO 100
 
118
      nq = int(w)
 
119
      w = w - dble(nq)
 
120
      nq = int(w*4.0D0)
 
121
      w = 4.0D0* (w-dble(nq)*.25D0)
 
122
C---------------------------------------------------------------------
 
123
C     W IS NOW RELATED TO THE FRACTIONAL PART OF  4.0 * X.
 
124
C     ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
 
125
C     QUADRANT AND DETERMINE SIGN
 
126
C---------------------------------------------------------------------
 
127
      n = nq/2
 
128
      IF ((n+n).NE.nq) w = 1.0D0 - w
 
129
      z = piov4*w
 
130
      m = n/2
 
131
      IF ((m+m).NE.n) sgn = -sgn
 
132
C---------------------------------------------------------------------
 
133
C     DETERMINE FINAL VALUE FOR  -PI*COTAN(PI*X)
 
134
C---------------------------------------------------------------------
 
135
      n = (nq+1)/2
 
136
      m = n/2
 
137
      m = m + m
 
138
      IF (m.NE.n) GO TO 30
 
139
C---------------------------------------------------------------------
 
140
C     CHECK FOR SINGULARITY
 
141
C---------------------------------------------------------------------
 
142
      IF (z.EQ.0.0D0) GO TO 100
 
143
C---------------------------------------------------------------------
 
144
C     USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
 
145
C     SIN/COS AS A SUBSTITUTE FOR TAN
 
146
C---------------------------------------------------------------------
 
147
      aug = sgn* ((cos(z)/sin(z))*4.0D0)
 
148
      GO TO 40
 
149
 
 
150
   30 aug = sgn* ((sin(z)/cos(z))*4.0D0)
 
151
   40 x = 1.0D0 - x
 
152
   50 IF (x.GT.3.0D0) GO TO 70
 
153
C---------------------------------------------------------------------
 
154
C     0.5 .LE. X .LE. 3.0
 
155
C---------------------------------------------------------------------
 
156
      den = x
 
157
      upper = p1(1)*x
 
158
C
 
159
      DO 60 i = 1,5
 
160
          den = (den+q1(i))*x
 
161
          upper = (upper+p1(i+1))*x
 
162
   60 CONTINUE
 
163
C
 
164
      den = (upper+p1(7))/ (den+q1(6))
 
165
      xmx0 = dble(x) - dx0
 
166
      psi = den*xmx0 + aug
 
167
      RETURN
 
168
C---------------------------------------------------------------------
 
169
C     IF X .GE. XMAX1, PSI = LN(X)
 
170
C---------------------------------------------------------------------
 
171
   70 IF (x.GE.xmax1) GO TO 90
 
172
C---------------------------------------------------------------------
 
173
C     3.0 .LT. X .LT. XMAX1
 
174
C---------------------------------------------------------------------
 
175
      w = 1.0D0/ (x*x)
 
176
      den = w
 
177
      upper = p2(1)*w
 
178
C
 
179
      DO 80 i = 1,3
 
180
          den = (den+q2(i))*w
 
181
          upper = (upper+p2(i+1))*w
 
182
   80 CONTINUE
 
183
C
 
184
      aug = upper/ (den+q2(4)) - 0.5D0/x + aug
 
185
   90 psi = aug + dlog(x)
 
186
      RETURN
 
187
C---------------------------------------------------------------------
 
188
C     ERROR RETURN
 
189
C---------------------------------------------------------------------
 
190
  100 psi = 0.0D0
 
191
      RETURN
 
192
 
 
193
      END