*DECK DRC DOUBLE PRECISION FUNCTION DRC (X, Y, IER) C***BEGIN PROLOGUE DRC C***PURPOSE Calculate a double precision approximation to C DRC(X,Y) = Integral from zero to infinity of C -1/2 -1 C (1/2)(t+X) (t+Y) dt, C where X is nonnegative and Y is positive. C***LIBRARY SLATEC C***CATEGORY C14 C***TYPE DOUBLE PRECISION (RC-S, DRC-D) C***KEYWORDS DUPLICATION THEOREM, ELEMENTARY FUNCTIONS, C ELLIPTIC INTEGRAL, TAYLOR SERIES C***AUTHOR Carlson, B. C. C Ames Laboratory-DOE C Iowa State University C Ames, IA 50011 C Notis, E. M. C Ames Laboratory-DOE C Iowa State University C Ames, IA 50011 C Pexton, R. L. C Lawrence Livermore National Laboratory C Livermore, CA 94550 C***DESCRIPTION C C 1. DRC C Standard FORTRAN function routine C Double precision version C The routine calculates an approximation result to C DRC(X,Y) = integral from zero to infinity of C C -1/2 -1 C (1/2)(t+X) (t+Y) dt, C C where X is nonnegative and Y is positive. The duplication C theorem is iterated until the variables are nearly equal, C and the function is then expanded in Taylor series to fifth C order. Logarithmic, inverse circular, and inverse hyper- C bolic functions can be expressed in terms of DRC. C C 2. Calling Sequence C DRC( X, Y, IER ) C C Parameters On Entry C Values assigned by the calling routine C C X - Double precision, nonnegative variable C C Y - Double precision, positive variable C C C C On Return (values assigned by the DRC routine) C C DRC - Double precision approximation to the integral C C IER - Integer to indicate normal or abnormal termination. C C IER = 0 Normal and reliable termination of the C routine. It is assumed that the requested C accuracy has been achieved. C C IER > 0 Abnormal termination of the routine C C X and Y are unaltered. C C 3. Error messages C C Value of IER assigned by the DRC routine C C Value assigned Error message printed C IER = 1 X.LT.0.0D0.OR.Y.LE.0.0D0 C = 2 X+Y.LT.LOLIM C = 3 MAX(X,Y) .GT. UPLIM C C 4. Control parameters C C Values of LOLIM, UPLIM, and ERRTOL are set by the C routine. C C LOLIM and UPLIM determine the valid range of X and Y C C LOLIM - Lower limit of valid arguments C C Not less than 5 * (machine minimum) . C C UPLIM - Upper limit of valid arguments C C Not greater than (machine maximum) / 5 . C C C Acceptable values for: LOLIM UPLIM C IBM 360/370 SERIES : 3.0D-78 1.0D+75 C CDC 6000/7000 SERIES : 1.0D-292 1.0D+321 C UNIVAC 1100 SERIES : 1.0D-307 1.0D+307 C CRAY : 2.3D-2466 1.0D+2465 C VAX 11 SERIES : 1.5D-38 3.0D+37 C C ERRTOL determines the accuracy of the answer C C The value assigned by the routine will result C in solution precision within 1-2 decimals of C "machine precision". C C C ERRTOL - relative error due to truncation is less than C 16 * ERRTOL ** 6 / (1 - 2 * ERRTOL). C C C The accuracy of the computed approximation to the inte- C gral can be controlled by choosing the value of ERRTOL. C Truncation of a Taylor series after terms of fifth order C introduces an error less than the amount shown in the C second column of the following table for each value of C ERRTOL in the first column. In addition to the trunca- C tion error there will be round-off error, but in prac- C tice the total error from both sources is usually less C than the amount given in the table. C C C C Sample choices: ERRTOL Relative truncation C error less than C 1.0D-3 2.0D-17 C 3.0D-3 2.0D-14 C 1.0D-2 2.0D-11 C 3.0D-2 2.0D-8 C 1.0D-1 2.0D-5 C C C Decreasing ERRTOL by a factor of 10 yields six more C decimal digits of accuracy at the expense of one or C two more iterations of the duplication theorem. C C *Long Description: C C DRC special comments C C C C C Check: DRC(X,X+Z) + DRC(Y,Y+Z) = DRC(0,Z) C C where X, Y, and Z are positive and X * Y = Z * Z C C C On Input: C C X, and Y are the variables in the integral DRC(X,Y). C C On Output: C C X and Y are unaltered. C C C C DRC(0,1/4)=DRC(1/16,1/8)=PI=3.14159... C C DRC(9/4,2)=LN(2) C C C C ******************************************************** C C WARNING: Changes in the program may improve speed at the C expense of robustness. C C C -------------------------------------------------------------------- C C Special functions via DRC C C C C LN X X .GT. 0 C C 2 C LN(X) = (X-1) DRC(((1+X)/2) , X ) C C C -------------------------------------------------------------------- C C ARCSIN X -1 .LE. X .LE. 1 C C 2 C ARCSIN X = X DRC (1-X ,1 ) C C -------------------------------------------------------------------- C C ARCCOS X 0 .LE. X .LE. 1 C C C 2 2 C ARCCOS X = SQRT(1-X ) DRC(X ,1 ) C C -------------------------------------------------------------------- C C ARCTAN X -INF .LT. X .LT. +INF C C 2 C ARCTAN X = X DRC(1,1+X ) C C -------------------------------------------------------------------- C C ARCCOT X 0 .LE. X .LT. INF C C 2 2 C ARCCOT X = DRC(X ,X +1 ) C C -------------------------------------------------------------------- C C ARCSINH X -INF .LT. X .LT. +INF C C 2 C ARCSINH X = X DRC(1+X ,1 ) C C -------------------------------------------------------------------- C C ARCCOSH X X .GE. 1 C C 2 2 C ARCCOSH X = SQRT(X -1) DRC(X ,1 ) C C -------------------------------------------------------------------- C C ARCTANH X -1 .LT. X .LT. 1 C C 2 C ARCTANH X = X DRC(1,1-X ) C C -------------------------------------------------------------------- C C ARCCOTH X X .GT. 1 C C 2 2 C ARCCOTH X = DRC(X ,X -1 ) C C -------------------------------------------------------------------- C C***REFERENCES B. C. Carlson and E. M. Notis, Algorithms for incomplete C elliptic integrals, ACM Transactions on Mathematical C Software 7, 3 (September 1981), pp. 398-403. C B. C. Carlson, Computing elliptic integrals by C duplication, Numerische Mathematik 33, (1979), C pp. 1-16. C B. C. Carlson, Elliptic integrals of the first kind, C SIAM Journal of Mathematical Analysis 8, (1977), C pp. 231-242. C***ROUTINES CALLED D1MACH, XERMSG C***REVISION HISTORY (YYMMDD) C 790801 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 891009 Removed unreferenced statement labels. (WRB) C 891009 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900326 Removed duplicate information from DESCRIPTION section. C (WRB) C 900510 Changed calls to XERMSG to standard form, and some C editorial changes. (RWC)) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DRC CHARACTER*16 XERN3, XERN4, XERN5 INTEGER IER DOUBLE PRECISION C1, C2, ERRTOL, LAMDA, LOLIM, D1MACH DOUBLE PRECISION MU, S, SN, UPLIM, X, XN, Y, YN LOGICAL FIRST SAVE ERRTOL,LOLIM,UPLIM,C1,C2,FIRST DATA FIRST /.TRUE./ C C***FIRST EXECUTABLE STATEMENT DRC IF (FIRST) THEN ERRTOL = (D1MACH(3)/16.0D0)**(1.0D0/6.0D0) LOLIM = 5.0D0 * D1MACH(1) UPLIM = D1MACH(2) / 5.0D0 C C1 = 1.0D0/7.0D0 C2 = 9.0D0/22.0D0 ENDIF FIRST = .FALSE. C C CALL ERROR HANDLER IF NECESSARY. C DRC = 0.0D0 IF (X.LT.0.0D0.OR.Y.LE.0.0D0) THEN IER = 1 WRITE (XERN3, '(1PE15.6)') X WRITE (XERN4, '(1PE15.6)') Y CALL XERMSG ('SLATEC', 'DRC', * 'X.LT.0 .OR. Y.LE.0 WHERE X = ' // XERN3 // ' AND Y = ' // * XERN4, 1, 1) RETURN ENDIF C IF (MAX(X,Y).GT.UPLIM) THEN IER = 3 WRITE (XERN3, '(1PE15.6)') X WRITE (XERN4, '(1PE15.6)') Y WRITE (XERN5, '(1PE15.6)') UPLIM CALL XERMSG ('SLATEC', 'DRC', * 'MAX(X,Y).GT.UPLIM WHERE X = ' // XERN3 // ' Y = ' // * XERN4 // ' AND UPLIM = ' // XERN5, 3, 1) RETURN ENDIF C IF (X+Y.LT.LOLIM) THEN IER = 2 WRITE (XERN3, '(1PE15.6)') X WRITE (XERN4, '(1PE15.6)') Y WRITE (XERN5, '(1PE15.6)') LOLIM CALL XERMSG ('SLATEC', 'DRC', * 'X+Y.LT.LOLIM WHERE X = ' // XERN3 // ' Y = ' // XERN4 // * ' AND LOLIM = ' // XERN5, 2, 1) RETURN ENDIF C IER = 0 XN = X YN = Y C 30 MU = (XN+YN+YN)/3.0D0 SN = (YN+MU)/MU - 2.0D0 IF (ABS(SN).LT.ERRTOL) GO TO 40 LAMDA = 2.0D0*SQRT(XN)*SQRT(YN) + YN XN = (XN+LAMDA)*0.250D0 YN = (YN+LAMDA)*0.250D0 GO TO 30 C 40 S = SN*SN*(0.30D0+SN*(C1+SN*(0.3750D0+SN*C2))) DRC = (1.0D0+S)/SQRT(MU) RETURN END