← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/rf.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/rf.f
 
1
*DECK RF
 
2
      REAL FUNCTION RF (X, Y, Z, IER)
 
3
C***BEGIN PROLOGUE  RF
 
4
C***PURPOSE  Compute the incomplete or complete elliptic integral of the
 
5
C            1st kind.  For X, Y, and Z non-negative and at most one of
 
6
C            them zero, RF(X,Y,Z) = Integral from zero to infinity of
 
7
C                                -1/2     -1/2     -1/2
 
8
C                      (1/2)(t+X)    (t+Y)    (t+Z)    dt.
 
9
C            If X, Y or Z is zero, the integral is complete.
 
10
C***LIBRARY   SLATEC
 
11
C***CATEGORY  C14
 
12
C***TYPE      SINGLE PRECISION (RF-S, DRF-D)
 
13
C***KEYWORDS  COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM,
 
14
C             INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE FIRST KIND,
 
15
C             TAYLOR SERIES
 
16
C***AUTHOR  Carlson, B. C.
 
17
C             Ames Laboratory-DOE
 
18
C             Iowa State University
 
19
C             Ames, IA  50011
 
20
C           Notis, E. M.
 
21
C             Ames Laboratory-DOE
 
22
C             Iowa State University
 
23
C             Ames, IA  50011
 
24
C           Pexton, R. L.
 
25
C             Lawrence Livermore National Laboratory
 
26
C             Livermore, CA  94550
 
27
C***DESCRIPTION
 
28
C
 
29
C   1.     RF
 
30
C          Evaluate an INCOMPLETE (or COMPLETE) ELLIPTIC INTEGRAL
 
31
C          of the first kind
 
32
C          Standard FORTRAN function routine
 
33
C          Single precision version
 
34
C          The routine calculates an approximation result to
 
35
C          RF(X,Y,Z) = Integral from zero to infinity of
 
36
C
 
37
C                               -1/2     -1/2     -1/2
 
38
C                     (1/2)(t+X)    (t+Y)    (t+Z)    dt,
 
39
C
 
40
C          where X, Y, and Z are nonnegative and at most one of them
 
41
C          is zero.  If one of them is zero, the integral is COMPLETE.
 
42
C          The duplication theorem is iterated until the variables are
 
43
C          nearly equal, and the function is then expanded in Taylor
 
44
C          series to fifth order.
 
45
C
 
46
C   2.     Calling Sequence
 
47
C          RF( X, Y, Z, IER )
 
48
C
 
49
C          Parameters on Entry
 
50
C          Values assigned by the calling routine
 
51
C
 
52
C          X      - Single precision, nonnegative variable
 
53
C
 
54
C          Y      - Single precision, nonnegative variable
 
55
C
 
56
C          Z      - Single precision, nonnegative variable
 
57
C
 
58
C
 
59
C
 
60
C          On Return     (values assigned by the RF routine)
 
61
C
 
62
C          RF     - Single precision approximation to the integral
 
63
C
 
64
C          IER    - Integer
 
65
C
 
66
C                   IER = 0 Normal and reliable termination of the
 
67
C                           routine.  It is assumed that the requested
 
68
C                           accuracy has been achieved.
 
69
C
 
70
C                   IER >  0 Abnormal termination of the routine
 
71
C
 
72
C          X, Y, Z are unaltered.
 
73
C
 
74
C
 
75
C   3.    Error Messages
 
76
C
 
77
C         Value of IER assigned by the RF routine
 
78
C
 
79
C                  Value assigned         Error Message Printed
 
80
C                  IER = 1                MIN(X,Y,Z) .LT. 0.0E0
 
81
C                      = 2                MIN(X+Y,X+Z,Y+Z) .LT. LOLIM
 
82
C                      = 3                MAX(X,Y,Z) .GT. UPLIM
 
83
C
 
84
C
 
85
C
 
86
C   4.     Control Parameters
 
87
C
 
88
C                  Values of LOLIM, UPLIM, and ERRTOL are set by the
 
89
C                  routine.
 
90
C
 
91
C          LOLIM and UPLIM determine the valid range of X, Y and Z
 
92
C
 
93
C          LOLIM  - Lower limit of valid arguments
 
94
C
 
95
C                   Not less than 5 * (machine minimum).
 
96
C
 
97
C          UPLIM  - Upper limit of valid arguments
 
98
C
 
99
C                   Not greater than (machine maximum) / 5.
 
100
C
 
101
C
 
102
C                     Acceptable Values For:   LOLIM      UPLIM
 
103
C                     IBM 360/370 SERIES   :   3.0E-78     1.0E+75
 
104
C                     CDC 6000/7000 SERIES :   1.0E-292    1.0E+321
 
105
C                     UNIVAC 1100 SERIES   :   1.0E-37     1.0E+37
 
106
C                     CRAY                 :   2.3E-2466   1.09E+2465
 
107
C                     VAX 11 SERIES        :   1.5E-38     3.0E+37
 
108
C
 
109
C
 
110
C
 
111
C          ERRTOL determines the accuracy of the answer
 
112
C
 
113
C                 The value assigned by the routine will result
 
114
C                 in solution precision within 1-2 decimals of
 
115
C                 "machine precision".
 
116
C
 
117
C
 
118
C
 
119
C          ERRTOL - Relative error due to truncation is less than
 
120
C                   ERRTOL ** 6 / (4 * (1-ERRTOL)  .
 
121
C
 
122
C
 
123
C
 
124
C              The accuracy of the computed approximation to the inte-
 
125
C              gral can be controlled by choosing the value of ERRTOL.
 
126
C              Truncation of a Taylor series after terms of fifth order
 
127
C              introduces an error less than the amount shown in the
 
128
C              second column of the following table for each value of
 
129
C              ERRTOL in the first column.  In addition to the trunca-
 
130
C              tion error there will be round-off error, but in prac-
 
131
C              tice the total error from both sources is usually less
 
132
C              than the amount given in the table.
 
133
C
 
134
C
 
135
C
 
136
C
 
137
C
 
138
C          Sample Choices:  ERRTOL   Relative Truncation
 
139
C                                    error less than
 
140
C                           1.0E-3    3.0E-19
 
141
C                           3.0E-3    2.0E-16
 
142
C                           1.0E-2    3.0E-13
 
143
C                           3.0E-2    2.0E-10
 
144
C                           1.0E-1    3.0E-7
 
145
C
 
146
C
 
147
C                    Decreasing ERRTOL by a factor of 10 yields six more
 
148
C                    decimal digits of accuracy at the expense of one or
 
149
C                    two more iterations of the duplication theorem.
 
150
C
 
151
C *Long Description:
 
152
C
 
153
C   RF Special Comments
 
154
C
 
155
C
 
156
C
 
157
C          Check by addition theorem: RF(X,X+Z,X+W) + RF(Y,Y+Z,Y+W)
 
158
C          = RF(0,Z,W), where X,Y,Z,W are positive and X * Y = Z * W.
 
159
C
 
160
C
 
161
C          On Input:
 
162
C
 
163
C          X, Y, and Z are the variables in the integral RF(X,Y,Z).
 
164
C
 
165
C
 
166
C          On Output:
 
167
C
 
168
C
 
169
C          X, Y, and Z are unaltered.
 
170
C
 
171
C
 
172
C
 
173
C          ********************************************************
 
174
C
 
175
C          Warning: Changes in the program may improve speed at the
 
176
C                   expense of robustness.
 
177
C
 
178
C
 
179
C
 
180
C   Special Functions via RF
 
181
C
 
182
C
 
183
C                  Legendre form of ELLIPTIC INTEGRAL of 1st kind
 
184
C                  ----------------------------------------------
 
185
C
 
186
C
 
187
C                                            2         2   2
 
188
C                  F(PHI,K) = SIN(PHI) RF(COS (PHI),1-K SIN (PHI),1)
 
189
C
 
190
C
 
191
C                                 2
 
192
C                  K(K) = RF(0,1-K ,1)
 
193
C
 
194
C                         PI/2     2   2      -1/2
 
195
C                       = INT  (1-K SIN (PHI) )   D PHI
 
196
C                          0
 
197
C
 
198
C
 
199
C
 
200
C
 
201
C
 
202
C                  Bulirsch form of ELLIPTIC INTEGRAL of 1st kind
 
203
C                  ----------------------------------------------
 
204
C
 
205
C
 
206
C                                         2 2    2
 
207
C                  EL1(X,KC) = X RF(1,1+KC X ,1+X )
 
208
C
 
209
C
 
210
C
 
211
C
 
212
C                  Lemniscate constant A
 
213
C                  ---------------------
 
214
C
 
215
C
 
216
C                       1      4 -1/2
 
217
C                  A = INT (1-S )    DS = RF(0,1,2) = RF(0,2,1)
 
218
C                       0
 
219
C
 
220
C
 
221
C    -------------------------------------------------------------------
 
222
C
 
223
C***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
 
224
C                 elliptic integrals, ACM Transactions on Mathematical
 
225
C                 Software 7, 3 (September 1981), pp. 398-403.
 
226
C               B. C. Carlson, Computing elliptic integrals by
 
227
C                 duplication, Numerische Mathematik 33, (1979),
 
228
C                 pp. 1-16.
 
229
C               B. C. Carlson, Elliptic integrals of the first kind,
 
230
C                 SIAM Journal of Mathematical Analysis 8, (1977),
 
231
C                 pp. 231-242.
 
232
C***ROUTINES CALLED  R1MACH, XERMSG
 
233
C***REVISION HISTORY  (YYMMDD)
 
234
C   790801  DATE WRITTEN
 
235
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
236
C   891009  Removed unreferenced statement labels.  (WRB)
 
237
C   891009  REVISION DATE from Version 3.2
 
238
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
239
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
240
C   900326  Removed duplicate information from DESCRIPTION section.
 
241
C           (WRB)
 
242
C   900510  Changed calls to XERMSG to standard form, and some
 
243
C           editorial changes.  (RWC))
 
244
C   920501  Reformatted the REFERENCES section.  (WRB)
 
245
C***END PROLOGUE  RF
 
246
      CHARACTER*16 XERN3, XERN4, XERN5, XERN6
 
247
      INTEGER IER
 
248
      REAL LOLIM, UPLIM, EPSLON, ERRTOL
 
249
      REAL C1, C2, C3, E2, E3, LAMDA
 
250
      REAL MU, S, X, XN, XNDEV
 
251
      REAL XNROOT, Y, YN, YNDEV, YNROOT, Z, ZN, ZNDEV,
 
252
     * ZNROOT
 
253
      LOGICAL FIRST
 
254
      SAVE ERRTOL,LOLIM,UPLIM,C1,C2,C3,FIRST
 
255
      DATA FIRST /.TRUE./
 
256
C
 
257
C***FIRST EXECUTABLE STATEMENT  RF
 
258
C
 
259
      IF (FIRST) THEN
 
260
         ERRTOL = (4.0E0*R1MACH(3))**(1.0E0/6.0E0)
 
261
         LOLIM  = 5.0E0 * R1MACH(1)
 
262
         UPLIM  = R1MACH(2)/5.0E0
 
263
C
 
264
         C1 = 1.0E0/24.0E0
 
265
         C2 = 3.0E0/44.0E0
 
266
         C3 = 1.0E0/14.0E0
 
267
      ENDIF
 
268
      FIRST = .FALSE.
 
269
C
 
270
C         CALL ERROR HANDLER IF NECESSARY.
 
271
C
 
272
      RF = 0.0E0
 
273
      IF (MIN(X,Y,Z).LT.0.0E0) THEN
 
274
         IER = 1
 
275
         WRITE (XERN3, '(1PE15.6)') X
 
276
         WRITE (XERN4, '(1PE15.6)') Y
 
277
         WRITE (XERN5, '(1PE15.6)') Z
 
278
         CALL XERMSG ('SLATEC', 'RF',
 
279
     *      'MIN(X,Y,Z).LT.0 WHERE X = ' // XERN3 // ' Y = ' // XERN4 //
 
280
     *      ' AND Z = ' // XERN5, 1, 1)
 
281
         RETURN
 
282
      ENDIF
 
283
C
 
284
      IF (MAX(X,Y,Z).GT.UPLIM) THEN
 
285
         IER = 3
 
286
         WRITE (XERN3, '(1PE15.6)') X
 
287
         WRITE (XERN4, '(1PE15.6)') Y
 
288
         WRITE (XERN5, '(1PE15.6)') Z
 
289
         WRITE (XERN6, '(1PE15.6)') UPLIM
 
290
         CALL XERMSG ('SLATEC', 'RF',
 
291
     *      'MAX(X,Y,Z).GT.UPLIM WHERE X = '  // XERN3 // ' Y = ' //
 
292
     *      XERN4 // ' Z = ' // XERN5 // ' AND UPLIM = ' // XERN6, 3, 1)
 
293
         RETURN
 
294
      ENDIF
 
295
C
 
296
      IF (MIN(X+Y,X+Z,Y+Z).LT.LOLIM) THEN
 
297
         IER = 2
 
298
         WRITE (XERN3, '(1PE15.6)') X
 
299
         WRITE (XERN4, '(1PE15.6)') Y
 
300
         WRITE (XERN5, '(1PE15.6)') Z
 
301
         WRITE (XERN6, '(1PE15.6)') LOLIM
 
302
         CALL XERMSG ('SLATEC', 'RF',
 
303
     *      'MIN(X+Y,X+Z,Y+Z).LT.LOLIM WHERE X = ' // XERN3 //
 
304
     *      ' Y = ' // XERN4 // ' Z = ' // XERN5 // ' AND LOLIM = ' //
 
305
     *      XERN6, 2, 1)
 
306
         RETURN
 
307
      ENDIF
 
308
C
 
309
      IER = 0
 
310
      XN = X
 
311
      YN = Y
 
312
      ZN = Z
 
313
C
 
314
   30 MU = (XN+YN+ZN)/3.0E0
 
315
      XNDEV = 2.0E0 - (MU+XN)/MU
 
316
      YNDEV = 2.0E0 - (MU+YN)/MU
 
317
      ZNDEV = 2.0E0 - (MU+ZN)/MU
 
318
      EPSLON = MAX(ABS(XNDEV), ABS(YNDEV), ABS(ZNDEV))
 
319
      IF (EPSLON.LT.ERRTOL) GO TO 40
 
320
      XNROOT =  SQRT(XN)
 
321
      YNROOT =  SQRT(YN)
 
322
      ZNROOT =  SQRT(ZN)
 
323
      LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
 
324
      XN = (XN+LAMDA)*0.250E0
 
325
      YN = (YN+LAMDA)*0.250E0
 
326
      ZN = (ZN+LAMDA)*0.250E0
 
327
      GO TO 30
 
328
C
 
329
   40 E2 = XNDEV*YNDEV - ZNDEV*ZNDEV
 
330
      E3 = XNDEV*YNDEV*ZNDEV
 
331
      S  = 1.0E0 + (C1*E2-0.10E0-C2*E3)*E2 + C3*E3
 
332
      RF = S/SQRT(MU)
 
333
C
 
334
      RETURN
 
335
      END