← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/drj.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/drj.f
 
1
*DECK DRJ
 
2
      DOUBLE PRECISION FUNCTION DRJ (X, Y, Z, P, IER)
 
3
C***BEGIN PROLOGUE  DRJ
 
4
C***PURPOSE  Compute the incomplete or complete (X or Y or Z is zero)
 
5
C            elliptic integral of the 3rd kind.  For X, Y, and Z non-
 
6
C            negative, at most one of them zero, and P positive,
 
7
C             RJ(X,Y,Z,P) = Integral from zero to infinity of
 
8
C                              -1/2     -1/2     -1/2     -1
 
9
C                    (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt.
 
10
C***LIBRARY   SLATEC
 
11
C***CATEGORY  C14
 
12
C***TYPE      DOUBLE PRECISION (RJ-S, DRJ-D)
 
13
C***KEYWORDS  COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM,
 
14
C             INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE THIRD 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.     DRJ
 
30
C          Standard FORTRAN function routine
 
31
C          Double precision version
 
32
C          The routine calculates an approximation result to
 
33
C          DRJ(X,Y,Z,P) = Integral from zero to infinity of
 
34
C
 
35
C                                -1/2     -1/2     -1/2     -1
 
36
C                      (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt,
 
37
C
 
38
C          where X, Y, and Z are nonnegative, at most one of them is
 
39
C          zero, and P is positive.  If X or Y or Z is zero, the
 
40
C          integral is COMPLETE.  The duplication theorem is iterated
 
41
C          until the variables are nearly equal, and the function is
 
42
C          then expanded in Taylor series to fifth order.
 
43
C
 
44
C
 
45
C   2.     Calling Sequence
 
46
C          DRJ( X, Y, Z, P, IER )
 
47
C
 
48
C          Parameters on Entry
 
49
C          Values assigned by the calling routine
 
50
C
 
51
C          X      - Double precision, nonnegative variable
 
52
C
 
53
C          Y      - Double precision, nonnegative variable
 
54
C
 
55
C          Z      - Double precision, nonnegative variable
 
56
C
 
57
C          P      - Double precision, positive variable
 
58
C
 
59
C
 
60
C          On  Return    (values assigned by the DRJ routine)
 
61
C
 
62
C          DRJ     - Double 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
 
73
C          X, Y, Z, P are unaltered.
 
74
C
 
75
C
 
76
C   3.    Error Messages
 
77
C
 
78
C         Value of IER assigned by the DRJ routine
 
79
C
 
80
C              Value assigned         Error Message printed
 
81
C              IER = 1                MIN(X,Y,Z) .LT. 0.0D0
 
82
C                  = 2                MIN(X+Y,X+Z,Y+Z,P) .LT. LOLIM
 
83
C                  = 3                MAX(X,Y,Z,P) .GT. UPLIM
 
84
C
 
85
C
 
86
C
 
87
C   4.     Control Parameters
 
88
C
 
89
C                  Values of LOLIM, UPLIM, and ERRTOL are set by the
 
90
C                  routine.
 
91
C
 
92
C
 
93
C          LOLIM and UPLIM determine the valid range of X, Y, Z, and P
 
94
C
 
95
C          LOLIM is not less than the cube root of the value
 
96
C          of LOLIM used in the routine for DRC.
 
97
C
 
98
C          UPLIM is not greater than 0.3 times the cube root of
 
99
C          the value of UPLIM used in the routine for DRC.
 
100
C
 
101
C
 
102
C                     Acceptable values for:   LOLIM      UPLIM
 
103
C                     IBM 360/370 SERIES   :   2.0D-26     3.0D+24
 
104
C                     CDC 6000/7000 SERIES :   5.0D-98     3.0D+106
 
105
C                     UNIVAC 1100 SERIES   :   5.0D-103    6.0D+101
 
106
C                     CRAY                 :   1.32D-822   1.4D+821
 
107
C                     VAX 11 SERIES        :   2.5D-13     9.0D+11
 
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
 
120
C          Relative error due to truncation of the series for DRJ
 
121
C          is less than 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
 
122
C
 
123
C
 
124
C
 
125
C        The accuracy of the computed approximation to the integral
 
126
C        can be controlled by choosing the value of ERRTOL.
 
127
C        Truncation of a Taylor series after terms of fifth order
 
128
C        introduces an error less than the amount shown in the
 
129
C        second column of the following table for each value of
 
130
C        ERRTOL in the first column.  In addition to the truncation
 
131
C        error there will be round-off error, but in practice the
 
132
C        total error from both sources is usually less than the
 
133
C        amount given in the table.
 
134
C
 
135
C
 
136
C
 
137
C          Sample choices:  ERRTOL   Relative truncation
 
138
C                                    error less than
 
139
C                           1.0D-3    4.0D-18
 
140
C                           3.0D-3    3.0D-15
 
141
C                           1.0D-2    4.0D-12
 
142
C                           3.0D-2    3.0D-9
 
143
C                           1.0D-1    4.0D-6
 
144
C
 
145
C                    Decreasing ERRTOL by a factor of 10 yields six more
 
146
C                    decimal digits of accuracy at the expense of one or
 
147
C                    two more iterations of the duplication theorem.
 
148
C
 
149
C *Long Description:
 
150
C
 
151
C   DRJ Special Comments
 
152
C
 
153
C
 
154
C     Check by addition theorem: DRJ(X,X+Z,X+W,X+P)
 
155
C     + DRJ(Y,Y+Z,Y+W,Y+P) + (A-B) * DRJ(A,B,B,A) + 3.0D0 / SQRT(A)
 
156
C     = DRJ(0,Z,W,P), where X,Y,Z,W,P are positive and X * Y
 
157
C     = Z * W,  A = P * P * (X+Y+Z+W),  B = P * (P+X) * (P+Y),
 
158
C     and B - A = P * (P-Z) * (P-W).  The sum of the third and
 
159
C     fourth terms on the left side is 3.0D0 * DRC(A,B).
 
160
C
 
161
C
 
162
C          On Input:
 
163
C
 
164
C     X, Y, Z, and P are the variables in the integral DRJ(X,Y,Z,P).
 
165
C
 
166
C
 
167
C          On Output:
 
168
C
 
169
C
 
170
C          X, Y, Z, P are unaltered.
 
171
C
 
172
C          ********************************************************
 
173
C
 
174
C          WARNING: Changes in the program may improve speed at the
 
175
C                   expense of robustness.
 
176
C
 
177
C    -------------------------------------------------------------------
 
178
C
 
179
C
 
180
C   Special double precision functions via DRJ and DRF
 
181
C
 
182
C
 
183
C                  Legendre form of ELLIPTIC INTEGRAL of 3rd kind
 
184
C                  -----------------------------------------
 
185
C
 
186
C
 
187
C                          PHI         2         -1
 
188
C             P(PHI,K,N) = INT (1+N SIN (THETA) )   *
 
189
C                           0
 
190
C
 
191
C
 
192
C                                  2    2         -1/2
 
193
C                             *(1-K  SIN (THETA) )     D THETA
 
194
C
 
195
C
 
196
C                                           2          2   2
 
197
C                        = SIN (PHI) DRF(COS (PHI), 1-K SIN (PHI),1)
 
198
C
 
199
C                                   3             2         2   2
 
200
C                         -(N/3) SIN (PHI) DRJ(COS (PHI),1-K SIN (PHI),
 
201
C
 
202
C                                  2
 
203
C                         1,1+N SIN (PHI))
 
204
C
 
205
C
 
206
C
 
207
C                  Bulirsch form of ELLIPTIC INTEGRAL of 3rd kind
 
208
C                  -----------------------------------------
 
209
C
 
210
C
 
211
C                                            2 2    2
 
212
C                  EL3(X,KC,P) = X DRF(1,1+KC X ,1+X ) +
 
213
C
 
214
C                                            3           2 2    2     2
 
215
C                               +(1/3)(1-P) X  DRJ(1,1+KC X ,1+X ,1+PX )
 
216
C
 
217
C
 
218
C                                           2
 
219
C                  CEL(KC,P,A,B) = A RF(0,KC ,1) +
 
220
C
 
221
C
 
222
C                                                      2
 
223
C                                 +(1/3)(B-PA) DRJ(0,KC ,1,P)
 
224
C
 
225
C
 
226
C                  Heuman's LAMBDA function
 
227
C                  -----------------------------------------
 
228
C
 
229
C
 
230
C                                2                      2      2    1/2
 
231
C                  L(A,B,P) =(COS (A)SIN(B)COS(B)/(1-COS (A)SIN (B))   )
 
232
C
 
233
C                                            2         2       2
 
234
C                            *(SIN(P) DRF(COS (P),1-SIN (A) SIN (P),1)
 
235
C
 
236
C                                 2       3            2       2
 
237
C                            +(SIN (A) SIN (P)/(3(1-COS (A) SIN (B))))
 
238
C
 
239
C                                    2         2       2
 
240
C                            *DRJ(COS (P),1-SIN (A) SIN (P),1,1-
 
241
C
 
242
C                                2       2          2       2
 
243
C                            -SIN (A) SIN (P)/(1-COS (A) SIN (B))))
 
244
C
 
245
C
 
246
C
 
247
C                  (PI/2) LAMBDA0(A,B) =L(A,B,PI/2) =
 
248
C
 
249
C                   2                         2       2    -1/2
 
250
C              = COS (A)  SIN(B) COS(B) (1-COS (A) SIN (B))
 
251
C
 
252
C                           2                  2       2
 
253
C                 *DRF(0,COS (A),1) + (1/3) SIN (A) COS (A)
 
254
C
 
255
C                                      2       2    -3/2
 
256
C                 *SIN(B) COS(B) (1-COS (A) SIN (B))
 
257
C
 
258
C                           2         2       2          2       2
 
259
C                 *DRJ(0,COS (A),1,COS (A) COS (B)/(1-COS (A) SIN (B)))
 
260
C
 
261
C
 
262
C                  Jacobi ZETA function
 
263
C                  -----------------------------------------
 
264
C
 
265
C                        2                     2   2    1/2
 
266
C             Z(B,K) = (K/3) SIN(B) COS(B) (1-K SIN (B))
 
267
C
 
268
C
 
269
C                                  2      2   2                 2
 
270
C                        *DRJ(0,1-K ,1,1-K SIN (B)) / DRF (0,1-K ,1)
 
271
C
 
272
C
 
273
C  ---------------------------------------------------------------------
 
274
C
 
275
C***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
 
276
C                 elliptic integrals, ACM Transactions on Mathematical
 
277
C                 Software 7, 3 (September 1981), pp. 398-403.
 
278
C               B. C. Carlson, Computing elliptic integrals by
 
279
C                 duplication, Numerische Mathematik 33, (1979),
 
280
C                 pp. 1-16.
 
281
C               B. C. Carlson, Elliptic integrals of the first kind,
 
282
C                 SIAM Journal of Mathematical Analysis 8, (1977),
 
283
C                 pp. 231-242.
 
284
C***ROUTINES CALLED  D1MACH, DRC, XERMSG
 
285
C***REVISION HISTORY  (YYMMDD)
 
286
C   790801  DATE WRITTEN
 
287
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
288
C   891009  Removed unreferenced statement labels.  (WRB)
 
289
C   891009  REVISION DATE from Version 3.2
 
290
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
291
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
292
C   900326  Removed duplicate information from DESCRIPTION section.
 
293
C           (WRB)
 
294
C   900510  Changed calls to XERMSG to standard form, and some
 
295
C           editorial changes.  (RWC)).
 
296
C   920501  Reformatted the REFERENCES section.  (WRB)
 
297
C***END PROLOGUE  DRJ
 
298
      INTEGER IER
 
299
      CHARACTER*16 XERN3, XERN4, XERN5, XERN6, XERN7
 
300
      DOUBLE PRECISION ALFA, BETA, C1, C2, C3, C4, EA, EB, EC, E2, E3
 
301
      DOUBLE PRECISION LOLIM, UPLIM, EPSLON, ERRTOL, D1MACH
 
302
      DOUBLE PRECISION LAMDA, MU, P, PN, PNDEV
 
303
      DOUBLE PRECISION POWER4, DRC, SIGMA, S1, S2, S3, X, XN, XNDEV
 
304
      DOUBLE PRECISION XNROOT, Y, YN, YNDEV, YNROOT, Z, ZN, ZNDEV,
 
305
     * ZNROOT
 
306
      LOGICAL FIRST
 
307
      SAVE ERRTOL,LOLIM,UPLIM,C1,C2,C3,C4,FIRST
 
308
      DATA FIRST /.TRUE./
 
309
C
 
310
C***FIRST EXECUTABLE STATEMENT  DRJ
 
311
      IF (FIRST) THEN
 
312
         ERRTOL = (D1MACH(3)/3.0D0)**(1.0D0/6.0D0)
 
313
         LOLIM  = (5.0D0 * D1MACH(1))**(1.0D0/3.0D0)
 
314
         UPLIM  = 0.30D0*( D1MACH(2) / 5.0D0)**(1.0D0/3.0D0)
 
315
C
 
316
         C1 = 3.0D0/14.0D0
 
317
         C2 = 1.0D0/3.0D0
 
318
         C3 = 3.0D0/22.0D0
 
319
         C4 = 3.0D0/26.0D0
 
320
      ENDIF
 
321
      FIRST = .FALSE.
 
322
C
 
323
C         CALL ERROR HANDLER IF NECESSARY.
 
324
C
 
325
      DRJ = 0.0D0
 
326
      IF (MIN(X,Y,Z).LT.0.0D0) THEN
 
327
         IER = 1
 
328
         WRITE (XERN3, '(1PE15.6)') X
 
329
         WRITE (XERN4, '(1PE15.6)') Y
 
330
         WRITE (XERN5, '(1PE15.6)') Z
 
331
         CALL XERMSG ('SLATEC', 'DRJ',
 
332
     *      'MIN(X,Y,Z).LT.0 WHERE X = ' // XERN3 // ' Y = ' // XERN4 //
 
333
     *      ' AND Z = ' // XERN5, 1, 1)
 
334
         RETURN
 
335
      ENDIF
 
336
C
 
337
      IF (MAX(X,Y,Z,P).GT.UPLIM) THEN
 
338
         IER = 3
 
339
         WRITE (XERN3, '(1PE15.6)') X
 
340
         WRITE (XERN4, '(1PE15.6)') Y
 
341
         WRITE (XERN5, '(1PE15.6)') Z
 
342
         WRITE (XERN6, '(1PE15.6)') P
 
343
         WRITE (XERN7, '(1PE15.6)') UPLIM
 
344
         CALL XERMSG ('SLATEC', 'DRJ',
 
345
     *      'MAX(X,Y,Z,P).GT.UPLIM WHERE X = ' // XERN3 // ' Y = ' //
 
346
     *      XERN4 // ' Z = ' // XERN5 // ' P = ' // XERN6 //
 
347
     *      ' AND UPLIM = ' // XERN7, 3, 1)
 
348
         RETURN
 
349
      ENDIF
 
350
C
 
351
      IF (MIN(X+Y,X+Z,Y+Z,P).LT.LOLIM) THEN
 
352
         IER = 2
 
353
         WRITE (XERN3, '(1PE15.6)') X
 
354
         WRITE (XERN4, '(1PE15.6)') Y
 
355
         WRITE (XERN5, '(1PE15.6)') Z
 
356
         WRITE (XERN6, '(1PE15.6)') P
 
357
         WRITE (XERN7, '(1PE15.6)') LOLIM
 
358
         CALL XERMSG ('SLATEC', 'RJ',
 
359
     *      'MIN(X+Y,X+Z,Y+Z,P).LT.LOLIM WHERE X = ' // XERN3 //
 
360
     *      ' Y = ' // XERN4 // ' Z = '  // XERN5 // ' P = ' // XERN6 //
 
361
     *      ' AND LOLIM = ', 2, 1)
 
362
         RETURN
 
363
      ENDIF
 
364
C
 
365
      IER = 0
 
366
      XN = X
 
367
      YN = Y
 
368
      ZN = Z
 
369
      PN = P
 
370
      SIGMA = 0.0D0
 
371
      POWER4 = 1.0D0
 
372
C
 
373
   30 MU = (XN+YN+ZN+PN+PN)*0.20D0
 
374
      XNDEV = (MU-XN)/MU
 
375
      YNDEV = (MU-YN)/MU
 
376
      ZNDEV = (MU-ZN)/MU
 
377
      PNDEV = (MU-PN)/MU
 
378
      EPSLON = MAX(ABS(XNDEV), ABS(YNDEV), ABS(ZNDEV), ABS(PNDEV))
 
379
      IF (EPSLON.LT.ERRTOL) GO TO 40
 
380
      XNROOT =  SQRT(XN)
 
381
      YNROOT =  SQRT(YN)
 
382
      ZNROOT =  SQRT(ZN)
 
383
      LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
 
384
      ALFA = PN*(XNROOT+YNROOT+ZNROOT) + XNROOT*YNROOT*ZNROOT
 
385
      ALFA = ALFA*ALFA
 
386
      BETA = PN*(PN+LAMDA)*(PN+LAMDA)
 
387
      SIGMA = SIGMA + POWER4*DRC(ALFA,BETA,IER)
 
388
      POWER4 = POWER4*0.250D0
 
389
      XN = (XN+LAMDA)*0.250D0
 
390
      YN = (YN+LAMDA)*0.250D0
 
391
      ZN = (ZN+LAMDA)*0.250D0
 
392
      PN = (PN+LAMDA)*0.250D0
 
393
      GO TO 30
 
394
C
 
395
   40 EA = XNDEV*(YNDEV+ZNDEV) + YNDEV*ZNDEV
 
396
      EB = XNDEV*YNDEV*ZNDEV
 
397
      EC = PNDEV*PNDEV
 
398
      E2 = EA - 3.0D0*EC
 
399
      E3 = EB + 2.0D0*PNDEV*(EA-EC)
 
400
      S1 = 1.0D0 + E2*(-C1+0.750D0*C3*E2-1.50D0*C4*E3)
 
401
      S2 = EB*(0.50D0*C2+PNDEV*(-C3-C3+PNDEV*C4))
 
402
      S3 = PNDEV*EA*(C2-PNDEV*C3) - C2*PNDEV*EC
 
403
      DRJ = 3.0D0*SIGMA + POWER4*(S1+S2+S3)/(MU* SQRT(MU))
 
404
      RETURN
 
405
      END