C ACM ALGORITHM 577
C
C ALGORITHMS FOR INCOMPLETE ELLIPTIC INTEGRALS
C
C BY B.C. CARLSON AND E.M. NOTIS
C
C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, SEPTEMBER, 1981.
C
C
C     THIS FILE CONTAINS FOUR SUBROUTINES FOR COMPUTING INCOMPLETE
C     ELLIPTIC INTEGRALS, FOLLOWED BY SIX DRIVERS FOR TESTING THE
C     SUBROUTINES.  EACH SUBROUTINE AND EACH DRIVER IS PRECEDED BY
C     A COMMENT CARD WITH A LINE OF DOLLAR SIGNS, AND EACH DRIVER IS
C     FOLLOWED BY ITS INPUT DATA IF ANY.  THE FOUR SUBROUTINES HAVE
C     THE NAMES RC, RF, RD, RJ IN THAT ORDER.  THE DRIVERS HAVE NO
C     NAMES BUT BEGIN WITH DESCRIPTIVE COMMENTS.  THE FIRST FOUR
C     DRIVERS TEST RC, RF, RD, RJ IN THAT ORDER.  THE FIFTH DRIVER
C     TESTS RC AGAINST LIBRARY ROUTINES.  THE SIXTH DRIVER TESTS RF
C     AGAINST THE FUNPACK SUBROUTINE DELIKM.
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C          ********************************************************
C
      DOUBLE PRECISION FUNCTION RC(X,Y,ERRTOL,IERR)
C
C          THIS FUNCTION SUBROUTINE COMPUTES THE ELEMENTARY INTEGRAL
C          RC(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 RC.  REFERENCE:
C          B. C. CARLSON, COMPUTING ELLIPTIC INTEGRALS BY DUPLICATION,
C          NUMER. MATH. 33 (1979), 1-16.  CODED BY B. C. CARLSON AND
C          ELAINE M. NOTIS, AMES LABORATORY-DOE, IOWA STATE UNIVERSITY,
C          AMES, IOWA 50011.  MARCH 1, 1980.
C
C          CHECK BY ADDITION THEOREM: RC(X,X+Z) + RC(Y,Y+Z) = RC(0,Z),
C          WHERE X, Y, AND Z ARE POSITIVE AND  X * Y = Z * Z.
C
      INTEGER IERR,PRINTR
      DOUBLE PRECISION C1,C2,ERRTOL,LAMDA,LOLIM
      DOUBLE PRECISION MU,S,SN,UPLIM,X,XN,Y,YN
C          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DSQRT
C
C          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
      DATA PRINTR/6/
C
C          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
C          LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5.
C          UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5.
C
C          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
C          IBM 360/370 SERIES   :   3.D-78     1.D+75
C          CDC 6000/7000 SERIES :   1.D-292    1.D+321
C          UNIVAC 1100 SERIES   :   1.D-307    1.D+307
C
C          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
C          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
C          LOLIM = 1.E-37 AND UPLIM = 1.E+37 BECAUSE THE MACHINE
C          EXTREMA CHANGE WITH THE PRECISION.
C
      DATA LOLIM/3.D-78/, UPLIM/1.D+75/
C
C          ON INPUT:
C
C          X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
C
C          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
C          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
C          16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
C
C          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
C                                    ERROR LESS THAN
C                           1.D-3    2.D-17
C                           3.D-3    2.D-14
C                           1.D-2    2.D-11
C                           3.D-2    2.D-8
C                           1.D-1    2.D-5
C
C          ON OUTPUT:
C
C          X, Y, AND ERRTOL ARE UNALTERED.
C
C          IERR IS THE RETURN ERROR CODE:
C               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
C               IERR = 1 FOR ABNORMAL TERMINATION.
C
C          ********************************************************
C          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
C          EXPENSE OF ROBUSTNESS.
C
      IF (X .LT. 0.D0  .OR.  Y .LE. 0.D0) GO TO 100
      IF ((X + Y) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RC)
      WRITE(PRINTR,108) X,Y
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16)
      IERR = 1
      GO TO 124
C
  112 IERR = 0
      XN = X
      YN = Y
C
  116 MU = (XN + YN + YN) / 3.D0
      SN = (YN + MU) / MU - 2.D0
      IF (DABS(SN) .LT. ERRTOL) GO TO 120
      LAMDA = 2.D0 * DSQRT(XN) * DSQRT(YN) + YN
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      GO TO 116
C
  120 C1 = 1.D0 / 7.D0
      C2 = 9.D0 / 22.D0
      S = SN * SN * (0.3D0 + SN * (C1 + SN * (0.375D0 + SN * C2)))
      RC = (1.D0 + S) / DSQRT(MU)
C
  124 RETURN
      END
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C          ********************************************************
C
      DOUBLE PRECISION FUNCTION RF(X,Y,Z,ERRTOL,IERR)
C
C          THIS FUNCTION SUBROUTINE COMPUTES THE INCOMPLETE ELLIPTIC
C          INTEGRAL OF THE FIRST KIND,
C          RF(X,Y,Z) = INTEGRAL FROM ZERO TO INFINITY OF
C
C                                -1/2     -1/2     -1/2
C                      (1/2)(T+X)    (T+Y)    (T+Z)    DT,
C
C          WHERE X, Y, AND Z ARE NONNEGATIVE AND AT MOST ONE OF THEM
C          IS ZERO.  IF ONE OF THEM IS ZERO, THE INTEGRAL IS COMPLETE.
C          THE DUPLICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE
C          NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR
C          SERIES TO FIFTH ORDER.  REFERENCE: B. C. CARLSON, COMPUTING
C          ELLIPTIC INTEGRALS BY DUPLICATION, NUMER. MATH. 33 (1979),
C          1-16.  CODED BY B. C. CARLSON AND ELAINE M. NOTIS, AMES
C          LABORATORY-DOE, IOWA STATE UNIVERSITY, AMES, IOWA 50011.
C          MARCH 1, 1980.
C
C          CHECK BY ADDITION THEOREM: RF(X,X+Z,X+W) + RF(Y,Y+Z,Y+W)
C          = RF(0,Z,W), WHERE X,Y,Z,W ARE POSITIVE AND X * Y = Z * W.
C
      INTEGER IERR,PRINTR
      DOUBLE PRECISION C1,C2,C3,E2,E3,EPSLON,ERRTOL,LAMDA
      DOUBLE PRECISION LOLIM,MU,S,UPLIM,X,XN,XNDEV,XNROOT
      DOUBLE PRECISION Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
C          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DMIN1,DSQRT
C
C          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
      DATA PRINTR/6/
C
C          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
C          LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5.
C          UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5.
C
C          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
C          IBM 360/370 SERIES   :   3.D-78     1.D+75
C          CDC 6000/7000 SERIES :   1.D-292    1.D+321
C          UNIVAC 1100 SERIES   :   1.D-307    1.D+307
C
C          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
C          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
C          LOLIM = 1.E-37 AND UPLIM = 1.E+37 BECAUSE THE MACHINE
C          EXTREMA CHANGE WITH THE PRECISION.
C
      DATA LOLIM/3.D-78/, UPLIM/1.D+75/
C
C          ON INPUT:
C
C          X, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RF(X,Y,Z).
C
C          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
C          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
C          ERRTOL ** 6 / (4 * (1 - ERRTOL)).
C
C          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
C                                    ERROR LESS THAN
C                           1.D-3    3.D-19
C                           3.D-3    2.D-16
C                           1.D-2    3.D-13
C                           3.D-2    2.D-10
C                           1.D-1    3.D-7
C
C          ON OUTPUT:
C
C          X, Y, Z, AND ERRTOL ARE UNALTERED.
C
C          IERR IS THE RETURN ERROR CODE:
C               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
C               IERR = 1 FOR ABNORMAL TERMINATION.
C
C          ********************************************************
C          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
C          EXPENSE OF ROBUSTNESS.
C
      IF (DMIN1(X,Y,Z) .LT. 0.D0) GO TO 100
      IF (DMIN1(X+Y,X+Z,Y+Z) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y,Z) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RF)
      WRITE(PRINTR,108) X,Y,Z
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16,4X,4HZ = ,D23.16)
      IERR = 1
      GO TO 124
C
  112 IERR = 0
      XN = X
      YN = Y
      ZN = Z
C
  116 MU = (XN + YN + ZN) / 3.D0
      XNDEV = 2.D0 - (MU + XN) / MU
      YNDEV = 2.D0 - (MU + YN) / MU
      ZNDEV = 2.D0 - (MU + ZN) / MU
      EPSLON = DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = DSQRT(XN)
      YNROOT = DSQRT(YN)
      ZNROOT = DSQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      ZN = (ZN + LAMDA) * 0.25D0
      GO TO 116
C
  120 C1 = 1.D0 / 24.D0
      C2 = 3.D0 / 44.D0
      C3 = 1.D0 / 14.D0
      E2 = XNDEV * YNDEV - ZNDEV * ZNDEV
      E3 = XNDEV * YNDEV * ZNDEV
      S = 1.D0 + (C1 * E2 - 0.1D0 - C2 * E3) * E2 + C3 * E3
      RF = S / DSQRT(MU)
C
  124 RETURN
      END
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C          ********************************************************
C
      DOUBLE PRECISION FUNCTION RD(X,Y,Z,ERRTOL,IERR)
C
C          THIS FUNCTION SUBROUTINE COMPUTES AN INCOMPLETE ELLIPTIC
C          INTEGRAL OF THE SECOND KIND,
C          RD(X,Y,Z) = INTEGRAL FROM ZERO TO INFINITY OF
C
C                                -1/2     -1/2     -3/2
C                      (3/2)(T+X)    (T+Y)    (T+Z)    DT,
C
C          WHERE X AND Y ARE NONNEGATIVE, X + Y IS POSITIVE, AND Z IS
C          POSITIVE.  IF X OR Y IS ZERO, THE INTEGRAL IS COMPLETE.
C          THE DUPLICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE
C          NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR
C          SERIES TO FIFTH ORDER.  REFERENCE: B. C. CARLSON, COMPUTING
C          ELLIPTIC INTEGRALS BY DUPLICATION, NUMER. MATH. 33 (1979),
C          1-16.  CODED BY B. C. CARLSON AND ELAINE M. NOTIS, AMES
C          LABORATORY-DOE, IOWA STATE UNIVERSITY, AMES, IOWA 50011.
C          MARCH 1, 1980..
C
C          CHECK: RD(X,Y,Z) + RD(Y,Z,X) + RD(Z,X,Y)
C          = 3 / DSQRT(X * Y * Z), WHERE X, Y, AND Z ARE POSITIVE.
C
      INTEGER IERR,PRINTR
      DOUBLE PRECISION C1,C2,C3,C4,EA,EB,EC,ED,EF,EPSLON,ERRTOL,LAMDA
      DOUBLE PRECISION LOLIM,MU,POWER4,SIGMA,S1,S2,UPLIM,X,XN,XNDEV
      DOUBLE PRECISION XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
C          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DMIN1,DSQRT
C
C          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
      DATA PRINTR/6/
C
C          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
C          LOLIM IS NOT LESS THAN 2 / (MACHINE MAXIMUM) ** (2/3).
C          UPLIM IS NOT GREATER THAN (0.1 * ERRTOL / MACHINE
C          MINIMUM) ** (2/3), WHERE ERRTOL IS DESCRIBED BELOW.
C          IN THE FOLLOWING TABLE IT IS ASSUMED THAT ERRTOL WILL
C          NEVER BE CHOSEN SMALLER THAN 1.D-5.
C
C          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
C          IBM 360/370 SERIES   :   6.D-51     1.D+48
C          CDC 6000/7000 SERIES :   5.D-215    2.D+191
C          UNIVAC 1100 SERIES   :   1.D-205    2.D+201
C
C          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
C          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
C          LOLIM = 1.E-25 AND UPLIM = 2.E+21 BECAUSE THE MACHINE
C          EXTREMA CHANGE WITH THE PRECISION.
C
      DATA LOLIM/6.D-51/, UPLIM/1.D+48/
C
C          ON INPUT:
C
C          X, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RD(X,Y,Z).
C
C          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
C          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
C          3 * ERRTOL ** 6 / (1-ERRTOL) ** 3/2.
C
C          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
C                                    ERROR LESS THAN
C                           1.D-3    4.D-18
C                           3.D-3    3.D-15
C                           1.D-2    4.D-12
C                           3.D-2    3.D-9
C                           1.D-1    4.D-6
C
C          ON OUTPUT:
C
C          X, Y, Z, AND ERRTOL ARE UNALTERED.
C
C          IERR IS THE RETURN ERROR CODE:
C               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
C               IERR = 1 FOR ABNORMAL TERMINATION.
C
C          ********************************************************
C          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
C          EXPENSE OF ROBUSTNESS.
C
      IF (DMIN1(X,Y) .LT. 0.D0) GO TO 100
      IF (DMIN1(X+Y,Z) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y,Z) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RD)
      WRITE(PRINTR,108) X,Y,Z
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16,4X,4HZ = ,D23.16)
      IERR = 1
      GO TO 124
C
  112 IERR = 0
      XN = X
      YN = Y
      ZN = Z
      SIGMA = 0.D0
      POWER4 = 1.D0
C
  116 MU = (XN + YN + 3.D0 * ZN) * 0.2D0
      XNDEV = (MU - XN) / MU
      YNDEV = (MU - YN) / MU
      ZNDEV = (MU - ZN) / MU
      EPSLON = DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = DSQRT(XN)
      YNROOT = DSQRT(YN)
      ZNROOT = DSQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      SIGMA = SIGMA + POWER4 / (ZNROOT * (ZN + LAMDA))
      POWER4 = POWER4 * 0.25D0
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      ZN = (ZN + LAMDA) * 0.25D0
      GO TO 116
C
  120 C1 = 3.D0 / 14.D0
      C2 = 1.D0 / 6.D0
      C3 = 9.D0 / 22.D0
      C4 = 3.D0 / 26.D0
      EA = XNDEV * YNDEV
      EB = ZNDEV * ZNDEV
      EC = EA - EB
      ED = EA - 6.D0 * EB
      EF = ED + EC + EC
      S1 = ED * (- C1 + 0.25D0 * C3 * ED - 1.5D0 * C4 * ZNDEV * EF)
      S2 = ZNDEV * (C2 * EF + ZNDEV * (- C3 * EC + ZNDEV * C4 * EA))
      RD = 3.D0 * SIGMA + POWER4 * (1.D0 + S1 + S2) / (MU * DSQRT(MU))
C
  124 RETURN
      END
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C          ********************************************************
C
      DOUBLE PRECISION FUNCTION RJ(X,Y,Z,P,ERRTOL,IERR)
C
C          THIS FUNCTION SUBROUTINE COMPUTES AN INCOMPLETE ELLIPTIC
C          INTEGRAL OF THE THIRD KIND,
C          RJ(X,Y,Z,P) = INTEGRAL FROM ZERO TO INFINITY OF
C
C                                  -1/2     -1/2     -1/2     -1
C                        (3/2)(T+X)    (T+Y)    (T+Z)    (T+P)  DT,
C
C          WHERE X, Y, AND Z ARE NONNEGATIVE, AT MOST ONE OF THEM IS
C          ZERO, AND P IS POSITIVE.  IF X OR Y OR Z IS ZERO, THE
C          INTEGRAL IS COMPLETE.  THE DUPLICATION THEOREM IS ITERATED
C          UNTIL THE VARIABLES ARE NEARLY EQUAL, AND THE FUNCTION IS
C          THEN EXPANDED IN TAYLOR SERIES TO FIFTH ORDER.  REFERENCE:
C          B. C. CARLSON, COMPUTING ELLIPTIC INTEGRALS BY DUPLICATION,
C          NUMER. MATH. 33 (1979), 1-16.  CODED BY B. C. CARLSON AND
C          ELAINE M. NOTIS, AMES LABORATORY-DOE, IOWA STATE UNIVERSITY,
C          AMES, IOWA 50011.  MARCH 1, 1980.
C
C          CHECK BY ADDITION THEOREM: RJ(X,X+Z,X+W,X+P)
C          + RJ(Y,Y+Z,Y+W,Y+P) + (A-B) * RJ(A,B,B,A) + 3 / DSQRT(A)
C          = RJ(0,Z,W,P), WHERE X,Y,Z,W,P ARE POSITIVE AND X * Y
C          = Z * W,  A = P * P * (X+Y+Z+W),  B = P * (P+X) * (P+Y),
C          AND B - A = P * (P-Z) * (P-W).  THE SUM OF THE THIRD AND
C          FOURTH TERMS ON THE LEFT SIDE IS 3 * RC(A,B).
C
      INTEGER IERR,PRINTR
      DOUBLE PRECISION ALFA,BETA,C1,C2,C3,C4,EA,EB,EC,E2,E3
      DOUBLE PRECISION EPSLON,ERRTOL,ETOLRC,LAMDA,LOLIM,MU,P,PN,PNDEV
      DOUBLE PRECISION POWER4,RC,SIGMA,S1,S2,S3,UPLIM,X,XN,XNDEV
      DOUBLE PRECISION XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
C          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DMIN1,DSQRT
C
C          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
      DATA PRINTR/6/
C
C          RC IS A FUNCTION COMPUTED BY AN EXTERNAL SUBROUTINE.
C
C          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
C          LOLIM IS NOT LESS THAN THE CUBE ROOT OF THE VALUE
C          OF LOLIM USED IN THE SUBROUTINE FOR RC.
C          UPLIM IS NOT GREATER THAN 0.3 TIMES THE CUBE ROOT OF
C          THE VALUE OF UPLIM USED IN THE SUBROUTINE FOR RC.
C
C          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
C          IBM 360/370 SERIES   :   2.D-26     3.D+24
C          CDC 6000/7000 SERIES :   5.D-98     3.D+106
C          UNIVAC 1100 SERIES   :   5.D-103    6.D+101
C
C          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
C          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
C          LOLIM = 5.E-13 AND UPLIM = 6.E+11 BECAUSE THE MACHINE
C          EXTREMA CHANGE WITH THE PRECISION.
C
      DATA LOLIM/2.D-26/, UPLIM/3.D+24/
C
C          ON INPUT:
C
C          X, Y, Z, AND P ARE THE VARIABLES IN THE INTEGRAL RJ(X,Y,Z,P).
C
C          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
C          RELATIVE ERROR DUE TO TRUNCATION OF THE SERIES FOR RJ
C          IS LESS THAN 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
C          AN ERROR TOLERANCE (ETOLRC) WILL BE PASSED TO THE SUBROUTINE
C          FOR RC TO MAKE THE TRUNCATION ERROR FOR RC LESS THAN FOR RJ.
C
C          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
C                                    ERROR LESS THAN
C                           1.D-3    4.D-18
C                           3.D-3    3.D-15
C                           1.D-2    4.D-12
C                           3.D-2    3.D-9
C                           1.D-1    4.D-6
C
C          ON OUTPUT:
C
C          X, Y, Z, P, AND ERRTOL ARE UNALTERED.
C
C          IERR IS THE RETURN ERROR CODE:
C               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
C               IERR = 1 FOR ABNORMAL TERMINATION.
C
C          ********************************************************
C          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
C          EXPENSE OF ROBUSTNESS.
C
      IF (DMIN1(X,Y,Z) .LT. 0.D0) GO TO 100
      IF (DMIN1(X+Y,X+Z,Y+Z,P) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y,Z,P) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RJ)
      WRITE(PRINTR,108) X,Y,Z,P
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16,4X,4HZ = ,D23.16,
     1       4X,4HP = ,D23.16)
      IERR = 1
      GO TO 124
C
  112 IERR = 0
      XN = X
      YN = Y
      ZN = Z
      PN = P
      SIGMA = 0.D0
      POWER4 = 1.D0
      ETOLRC = 0.5D0 * ERRTOL
C
  116 MU = (XN + YN + ZN + PN + PN) * 0.2D0
      XNDEV = (MU - XN) / MU
      YNDEV = (MU - YN) / MU
      ZNDEV = (MU - ZN) / MU
      PNDEV = (MU - PN) / MU
      EPSLON = DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV),DABS(PNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = DSQRT(XN)
      YNROOT = DSQRT(YN)
      ZNROOT = DSQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      ALFA = PN * (XNROOT + YNROOT + ZNROOT) + XNROOT * YNROOT * ZNROOT
      ALFA = ALFA * ALFA
      BETA = PN * (PN + LAMDA) * (PN + LAMDA)
      SIGMA = SIGMA + POWER4 * RC(ALFA,BETA,ETOLRC,IERR)
      IF (IERR .NE. 0) GO TO 124
      POWER4 = POWER4 * 0.25D0
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      ZN = (ZN + LAMDA) * 0.25D0
      PN = (PN + LAMDA) * 0.25D0
      GO TO 116
C
  120 C1 = 3.D0 / 14.D0
      C2 = 1.D0 / 3.D0
      C3 = 3.D0 / 22.D0
      C4 = 3.D0 / 26.D0
      EA = XNDEV * (YNDEV + ZNDEV) + YNDEV * ZNDEV
      EB = XNDEV * YNDEV * ZNDEV
      EC = PNDEV * PNDEV
      E2 = EA - 3.D0 * EC
      E3 = EB + 2.D0 * PNDEV * (EA - EC)
      S1 = 1.D0 + E2 * (- C1 + 0.75D0 * C3 * E2 - 1.5D0 * C4 * E3)
      S2 = EB * (0.5D0 * C2 + PNDEV * (- C3 - C3 + PNDEV * C4))
      S3 = PNDEV * EA * (C2 - PNDEV * C3) - C2 * PNDEV * EC
      RJ = 3.D0 * SIGMA + POWER4 * (S1 + S2 + S3) / (MU * DSQRT(MU))
C
  124 RETURN
      END
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C THIS DRIVER TESTS THE DOUBLE PRECISION FUNCTION SUBROUTINE FOR THE
C INTEGRAL RC(X,Y).  THE FIRST SIX SETS OF VALUES OF X AND Y ARE
C EXTREME POINTS OF THE REGION OF VALID ARGUMENTS DEFINED BY THE
C MACHINE-DEPENDENT CONSTANTS LOLIM AND UPLIM.  THE VALUES OF LOLIM,
C UPLIM, X, Y, AND ERRTOL (SEE COMMENTS IN SUBROUTINE) MAY BE USED ON
C MOST MACHINES BUT PROVIDE A SEVERE TEST OF ROBUSTNESS ONLY ON THE
C IBM 360/370 SERIES.  THE SEVENTH SET TESTS THE FAILURE EXIT.  THE
C NEXT THREE SETS ARE CHECK VALUES: RC(0,0.25) = RC(0.0625,0.125) = PI
C AND RC(2.25,2) = LN(2).  THE REMAINING SETS SHOW THE DEPENDENCE ON X
C WHEN Y = 1.  FIXING Y ENTAILS NO LOSS HERE BECAUSE RC IS HOMOGENEOUS.
C
      INTEGER IERR,PRINTR,READR
      DOUBLE PRECISION X,Y,ELIPTC,ERRTOL
      DOUBLE PRECISION RC
C
      DATA PRINTR/6/, READR/5/
C
      WRITE(PRINTR,4)
    4 FORMAT(1H1,24HTEST OF ROBUSTNESS OF RC)
      WRITE(PRINTR,8)
    8 FORMAT(1H0,14X,1HX,26X,1HY,25X,2HRC)
C
      ERRTOL = 1.D-3
C
      DO 20 I = 1,43
        READ(READR,12) X,Y
   12   FORMAT(2D11.2)
        ELIPTC = RC(X,Y,ERRTOL,IERR)
        IF (IERR .NE. 0) GO TO 18
        WRITE(PRINTR,16) X,Y,ELIPTC
   16   FORMAT(1H0,3D27.16)
        GO TO 20
   18   WRITE(PRINTR,19) X,Y
   19   FORMAT(1H0,2D27.16)
   20   CONTINUE
C
      STOP
      END
C
   1.51D-78   1.51D-78
   3.01D-78   0.55D-78
   0.00D+00   3.01D-78
   0.99D+75   0.55D-78
   0.00D+00   0.99D+75
   0.99D+75   0.99D+75
   0.00D+00   2.00D-78
   0.00D+00   2.50D-01
   6.25D-02   1.25D-01
   2.25D+00   2.00D+00
   0.01D+00   1.00D+00
   0.02D+00   1.00D+00
   0.05D+00   1.00D+00
   0.10D+00   1.00D+00
   0.20D+00   1.00D+00
   0.40D+00   1.00D+00
   0.60D+00   1.00D+00
   0.80D+00   1.00D+00
   1.00D+00   1.00D+00
   1.20D+00   1.00D+00
   1.50D+00   1.00D+00
   2.00D+00   1.00D+00
   3.00D+00   1.00D+00
   4.00D+00   1.00D+00
   5.00D+00   1.00D+00
   1.00D+01   1.00D+00
   2.00D+01   1.00D+00
   5.00D+01   1.00D+00
   1.00D+02   1.00D+00
   1.00D+03   1.00D+00
   1.00D+04   1.00D+00
   1.00D+05   1.00D+00
   1.00D+06   1.00D+00
   1.00D+07   1.00D+00
   1.00D+08   1.00D+00
   1.00D+09   1.00D+00
   1.00D+10   1.00D+00
   1.00D+12   1.00D+00
   1.00D+15   1.00D+00
   1.00D+20   1.00D+00
   1.00D+30   1.00D+00
   1.00D+40   1.00D+00
   1.00D+50   1.00D+00
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C THIS DRIVER TESTS THE DOUBLE PRECISION FUNCTION SUBROUTINE FOR THE
C INTEGRAL RF(X,Y,Z), WHICH IS SYMMETRIC IN X,Y,Z.  THE FIRST NINE
C SETS OF VALUES OF X,Y,Z ARE EXTREME POINTS OF THE REGION OF VALID
C ARGUMENTS DEFINED BY THE MACHINE-DEPENDENT CONSTANTS LOLIM AND
C UPLIM.  THE VALUES OF LOLIM, UPLIM, X, Y, Z, AND ERRTOL (SEE
C COMMENTS IN SUBROUTINE) MAY BE USED ON MOST MACHINES BUT PROVIDE A
C SEVERE TEST OF ROBUSTNESS ONLY ON THE IBM 360/370 SERIES.  THE
C TENTH SET TESTS THE FAILURE EXIT.  THE ELEVENTH SET IS A CHECK
C VALUE: RF(0,1,2) = A, WHERE A IS THE FIRST LEMNISCATE CONSTANT.
C THE REMAINING SETS SHOW THE DEPENDENCE ON Z WHEN Y = 1 (NO LOSS OF
C GENERALITY BECAUSE OF HOMOGENEITY) AND X = 0.5 (MIDWAY BETWEEN THE
C COMPLETE CASE X = 0 AND THE DEGENERATE CASE X = Y).
C
      INTEGER IERR,PRINTR,READR
      DOUBLE PRECISION X,Y,Z,ELIPTC,ERRTOL
      DOUBLE PRECISION RF
C
      DATA PRINTR/6/, READR/5/
C
      WRITE(PRINTR,4)
    4 FORMAT(1H1,24HTEST OF ROBUSTNESS OF RF)
      WRITE(PRINTR,8)
    8 FORMAT(1H0,14X,1HX,26X,1HY,26X,1HZ,25X,2HRF)
C
      ERRTOL = 1.D-3
C
      DO 20 I = 1,55
        READ(READR,12) X,Y,Z
   12   FORMAT(3D11.2)
        ELIPTC = RF(X,Y,Z,ERRTOL,IERR)
        IF (IERR .NE. 0) GO TO 18
        WRITE(PRINTR,16) X,Y,Z,ELIPTC
   16   FORMAT(1H0,4D27.16)
        GO TO 20
   18   WRITE(PRINTR,19) X,Y,Z
   19   FORMAT(1H0,3D27.16)
   20   CONTINUE
C
      STOP
      END
C
   1.51D-78   1.51D-78   1.51D-78
   1.51D-78   1.51D-78   0.99D+75
   0.00D+00   3.01D-78   3.01D-78
   0.00D+00   3.01D-78   0.99D+75
   0.00D+00   0.99D+75   0.99D+75
   0.99D+75   0.99D+75   0.99D+75
   0.55D-78   3.01D-78   3.01D-78
   0.55D-78   3.01D-78   0.99D+75
   0.55D-78   0.99D+75   0.99D+75
   0.00D+00   2.00D-78   1.00D+00
   0.00D+00   1.00D+00   2.00D+00
   0.50D+00   1.00D+00   1.00D+00
   0.50D+00   1.00D+00   1.10D+00
   0.50D+00   1.00D+00   1.20D+00
   0.50D+00   1.00D+00   1.30D+00
   0.50D+00   1.00D+00   1.40D+00
   0.50D+00   1.00D+00   1.50D+00
   0.50D+00   1.00D+00   1.60D+00
   0.50D+00   1.00D+00   1.70D+00
   0.50D+00   1.00D+00   1.80D+00
   0.50D+00   1.00D+00   1.90D+00
   0.50D+00   1.00D+00   2.00D+00
   0.50D+00   1.00D+00   2.20D+00
   0.50D+00   1.00D+00   2.40D+00
   0.50D+00   1.00D+00   2.60D+00
   0.50D+00   1.00D+00   2.80D+00
   0.50D+00   1.00D+00   3.00D+00
   0.50D+00   1.00D+00   3.50D+00
   0.50D+00   1.00D+00   4.00D+00
   0.50D+00   1.00D+00   4.50D+00
   0.50D+00   1.00D+00   5.00D+00
   0.50D+00   1.00D+00   6.00D+00
   0.50D+00   1.00D+00   7.00D+00
   0.50D+00   1.00D+00   8.00D+00
   0.50D+00   1.00D+00   9.00D+00
   0.50D+00   1.00D+00   1.00D+01
   0.50D+00   1.00D+00   2.00D+01
   0.50D+00   1.00D+00   3.00D+01
   0.50D+00   1.00D+00   4.00D+01
   0.50D+00   1.00D+00   5.00D+01
   0.50D+00   1.00D+00   1.00D+02
   0.50D+00   1.00D+00   2.00D+02
   0.50D+00   1.00D+00   5.00D+02
   0.50D+00   1.00D+00   1.00D+03
   0.50D+00   1.00D+00   1.00D+04
   0.50D+00   1.00D+00   1.00D+05
   0.50D+00   1.00D+00   1.00D+06
   0.50D+00   1.00D+00   1.00D+08
   0.50D+00   1.00D+00   1.00D+10
   0.50D+00   1.00D+00   1.00D+12
   0.50D+00   1.00D+00   1.00D+15
   0.50D+00   1.00D+00   1.00D+20
   0.50D+00   1.00D+00   1.00D+30
   0.50D+00   1.00D+00   1.00D+40
   0.50D+00   1.00D+00   1.00D+50
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C THIS DRIVER TESTS THE DOUBLE PRECISION FUNCTION SUBROUTINE FOR THE
C INTEGRAL RD(X,Y,Z), WHICH IS SYMMETRIC IN X AND Y.  THE FIRST
C TWELVE SETS OF VALUES OF X,Y,Z ARE EXTREME POINTS OF THE REGION OF
C VALID ARGUMENTS DEFINED BY THE MACHINE-DEPENDENT CONSTANTS LOLIM
C AND UPLIM.  THE VALUES OF LOLIM, UPLIM, X, Y, Z, AND ERRTOL (SEE
C COMMENTS IN SUBROUTINE) MAY BE USED ON MOST MACHINES BUT PROVIDE A
C SEVERE TEST OF ROBUSTNESS ONLY ON THE IBM 360/370 SERIES.  THE
C THIRTEENTH SET TESTS THE FAILURE EXIT.  THE FOURTEENTH SET IS A
C CHECK VALUE: RD(0,2,1) = 3B = 3(PI)/4A, WHERE A AND B ARE THE
C LEMNISCATE CONSTANTS.  THE REMAINING SETS SHOW THE DEPENDENCE
C ON Z WHEN Y = 1 (NO LOSS OF GENERALITY BECAUSE OF HOMOGENEITY)
C AND X = 0.5 (MIDWAY BETWEEN THE COMPLETE CASE X = 0 AND THE
C DEGENERATE CASE X = Y).
C
      INTEGER IERR,PRINTR,READR
      DOUBLE PRECISION X,Y,Z,ELIPTC,ERRTOL
      DOUBLE PRECISION RD
C
      DATA PRINTR/6/, READR/5/
C
      WRITE(PRINTR,4)
    4 FORMAT(1H0,24HTEST OF ROBUSTNESS OF RD)
      WRITE(PRINTR,8)
    8 FORMAT(1H0,15X,1HX,26X,1HY,26X,1HZ,25X,2HRD)
C
      ERRTOL = 1.D-3
C
      DO 20 I = 1,27
        READ(READR,12) X,Y,Z
   12   FORMAT(3D11.2)
        ELIPTC = RD(X,Y,Z,ERRTOL,IERR)
        IF (IERR .NE. 0) GO TO 18
        WRITE(PRINTR,16) X,Y,Z,ELIPTC
   16   FORMAT(1H0,4D27.16)
        GO TO 20
   18   WRITE(PRINTR,19) X,Y,Z
   19   FORMAT(1H0,3D27.16)
   20   CONTINUE
C
      STOP
      END
C
   0.00D+00   6.01D-51   6.01D-51
   0.55D-78   6.01D-51   6.01D-51
   0.00D+00   6.01D-51   0.99D+48
   0.55D-78   6.01D-51   0.99D+48
   0.00D+00   0.99D+48   6.01D-51
   0.55D-78   0.99D+48   6.01D-51
   0.00D+00   0.99D+48   0.99D+48
   0.55D-78   0.99D+48   0.99D+48
   3.01D-51   3.01D-51   6.01D-51
   3.01D-51   3.01D-51   0.99D+48
   0.99D+48   0.99D+48   6.01D-51
   0.99D+48   0.99D+48   0.99D+48
   0.00D+00   3.01D-51   1.00D+00
   0.00D+00   2.00D+00   1.00D+00
   0.50D+00   1.00D+00   1.00D-10
   0.50D+00   1.00D+00   1.00D-05
   0.5D+00    1.00D+00   1.00D-02
   0.50D+00   1.00D+00   1.00D-01
   0.50D+00   1.00D+00   2.00D-01
   0.50D+00   1.00D+00   5.00D-01
   0.50D+00   1.00D+00   1.00D+00
   0.50D+00   1.00D+00   2.00D+00
   0.50D+00   1.00D+00   5.00D+00
   0.50D+00   1.00D+00   1.00D+01
   0.50D+00   1.00D+00   1.00D+02
   0.50D+00   1.00D+00   1.00D+05
   0.50D+00   1.00D+00   1.00D+10
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C THIS DRIVER TESTS THE DOUBLE PRECISION FUNCTION SUBROUTINE FOR THE
C INTEGRAL RJ(X,Y,Z,P), WHICH IS SYMMETRIC IN X,Y,Z.  THE FIRST
C TWENTY SETS OF VALUES OF X,Y,Z,P ARE EXTREME POINTS OF THE REGION
C OF VALID ARGUMENTS DEFINED BY THE MACHINE-DEPENDENT CONSTANTS
C LOLIM AND UPLIM.  THE VALUES OF LOLIM, UPLIM, X, Y, Z, P, AND
C ERRTOL (SEE COMMENTS IN SUBROUTINE) MAY BE USED ON MOST MACHINES
C BUT PROVIDE A SEVERE TEST OF ROBUSTNESS ONLY ON THE IBM 360/370
C SERIES.  THE TWENTY-FIRST SET TESTS THE FAILURE EXIT.  THE TWENTY-
C SECOND SET IS A CHECK VALUE: RJ(2,3,4,5) = 0.14297 57966 71567
C 53833 23308.  THE REMAINING SETS SHOW THE DEPENDENCE ON Z AND P
C WHEN Y = 1 (NO LOSS OF GENERALITY BECAUSE OF HOMOGENEITY) AND
C X = 0.5 (MIDWAY BETWEEN THE COMPLETE CASE X = 0 AND THE DEGENERATE
C CASE X = Y).
C
      INTEGER IERR,PRINTR,READR
      DOUBLE PRECISION X,Y,Z,P,ELIPTC,ERRTOL
      DOUBLE PRECISION RJ
C
      DATA PRINTR/6/, READR/5/
C
      WRITE(PRINTR,4)
    4 FORMAT(1H0,24HTEST OF ROBUSTNESS OF RJ)
      WRITE(PRINTR,8)
    8 FORMAT(1H0,12X,1HX,25X,1HY,25X,1HZ,25X,1HP,25X,2HRJ)
C
      ERRTOL = 1.D-3
C
      DO 20 I = 1,42
        READ(READR,12) X,Y,Z,P
   12   FORMAT(4D11.2)
        ELIPTC = RJ(X,Y,Z,P,ERRTOL,IERR)
        IF (IERR .NE. 0) GO TO 18
        WRITE(PRINTR,16) X,Y,Z,P,ELIPTC
   16   FORMAT(1H0,5D26.16)
        GO TO 20
   18   WRITE(PRINTR,19) X,Y,Z,P
   19   FORMAT(1H0,4D26.16)
   20   CONTINUE
C
      STOP
      END
C
   1.01D-26   1.01D-26   1.01D-26   2.01D-26
   1.01D-26   1.01D-26   2.99D+24   2.01D-26
   0.00D+00   2.01D-26   2.01D-26   2.01D-26
   0.00D+00   2.01D-26   2.99D+24   2.01D-26
   0.00D+00   2.99D+24   2.99D+24   2.01D-26
   2.99D+24   2.99D+24   2.99D+24   2.01D-26
   0.55D-78   2.01D-26   2.01D-26   2.01D-26
   0.55D-78   2.01D-26   2.99D+24   2.01D-26
   0.55D-78   2.99D+24   2.99D+24   2.01D-26
   2.01D-26   2.01D-26   2.01D-26   2.01D-26
   1.01D-26   1.01D-26   1.01D-26   2.99D+24
   1.01D-26   1.01D-26   2.99D+24   2.99D+24
   0.00D+00   2.01D-26   2.01D-26   2.99D+24
   0.00D+00   2.01D-26   2.99D+24   2.99D+24
   0.00D+00   2.99D+24   2.99D+24   2.99D+24
   2.99D+24   2.99D+24   2.99D+24   2.99D+24
   0.55D-78   2.01D-26   2.01D-26   2.99D+24
   0.55D-78   2.01D-26   2.99D+24   2.99D+24
   0.55D-78   2.99D+24   2.99D+24   2.99D+24
   2.01D-26   2.01D-26   2.01D-26   2.99D+24
   0.00D+00   1.90D-26   1.90D-26   1.00D+00
   2.00D+00   3.00D+00   4.00D+00   5.00D+00
   0.50D+00   1.00D+00   1.00D+00   0.25D+00
   0.50D+00   1.00D+00   1.00D+00   0.75D+00
   0.50D+00   1.00D+00   1.00D+00   1.00D+00
   0.50D+00   1.00D+00   1.00D+00   2.00D+00
   0.50D+00   1.00D+00   2.00D+00   0.25D+00
   0.50D+00   1.00D+00   2.00D+00   0.75D+00
   0.50D+00   1.00D+00   2.00D+00   1.50D+00
   0.50D+00   1.00D+00   2.00D+00   4.00D+00
   0.50D+00   1.00D+00   5.00D+00   0.25D+00
   0.50D+00   1.00D+00   5.00D+00   0.75D+00
   0.50D+00   1.00D+00   5.00D+00   3.00D+00
   0.50D+00   1.00D+00   5.00D+00   1.00D+01
   0.50D+00   1.00D+00   1.00D+01   0.25D+00
   0.50D+00   1.00D+00   1.00D+01   0.75D+00
   0.50D+00   1.00D+00   1.00D+01   5.00D+00
   0.50D+00   1.00D+00   1.00D+01   2.00D+01
   0.50D+00   1.00D+00   1.00D+02   0.25D+00
   0.50D+00   1.00D+00   1.00D+02   0.75D+00
   0.50D+00   1.00D+00   1.00D+02   5.00D+01
   0.50D+00   1.00D+00   1.00D+02   2.00D+02
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C THIS DRIVER COMPARES VALUES OF (LOG X)/(X-1) AND ARCTAN(X)
C CALCULATED ON ONE HAND FROM THE SUBROUTINE RC AND ON THE OTHER
C FROM LIBRARY LOG AND ARCTAN ROUTINES.  TO AVOID OVER/UNDERFLOWS
C FOR EXTREME VALUES OF X, WE WRITE (LOG X)/(X-1) = RC(Y,X/Y)/SQRT(Y),
C WHERE Y = (1+X)/2, AND ARCTAN(X) = SQRT(X)*RC(Z,Z+X), WHERE Z = 1/X.
      INTEGER IERR,PRINTR,READR,IPOWER
      DOUBLE PRECISION ERRTOL,IBMARC,IBMLOG,MYARC,MYLOG,RC,V,W,X,Y,Z
C
      DATA PRINTR/6/, READR/5/
C
      WRITE(PRINTR,4)
    4 FORMAT(1H1,5X,1HX,15X,8HFROM LOG,19X,7HFROM RC)
      ERRTOL = 1.D-3
      DO 20 J=1,10
        X = 0.2D0 * J
        Y = (1.D0 + X) / 2.D0
        V = X / Y
        MYLOG = RC(Y,V,ERRTOL,IERR) / DSQRT(Y)
        IF (J .NE. 5) GO TO 12
        WRITE(PRINTR,8) X,MYLOG
    8   FORMAT(1H0,1PD9.1,5X,22H**** ZERO DIVIDE *****,0PD27.16)
        GO TO 20
   12   IBMLOG = DLOG(X) / (X - 1.D0)
        WRITE(PRINTR,16) X,IBMLOG,MYLOG
   16   FORMAT(1H0,1PD9.1,0P2D27.16)
   20   CONTINUE
C
C*************************************************************
C
      WRITE(PRINTR,24)
   24 FORMAT(1H-,30X,19HEXTREME VALUES OF X)
      WRITE(PRINTR,28)
   28 FORMAT(1H0,5X,1HX,15X,8HFROM LOG,19X,7HFROM RC)
      DO 40 I=1,16
        IPOWER = -75 + 10 * (I-1)
        X = 10.D0 ** IPOWER
        Y = (1.D0 + X) / 2.D0
        V = X / Y
        MYLOG = RC(Y,V,ERRTOL,IERR) / DSQRT(Y)
        IBMLOG = DLOG(X) / (X - 1.D0)
        WRITE(PRINTR,30) X,IBMLOG,MYLOG
   30   FORMAT(1H0,1PD9.1,0P2D27.16)
   40   CONTINUE
C
C*************************************************************
C
      WRITE(PRINTR,104)
  104 FORMAT(1H1,5X,1HX,14X,11HFROM ARCTAN,17X,7HFROM RC)
      DO 116 M=1,13
        READ(READR,108) X
  108   FORMAT(1D9.0)
        Z = 1.D0 / X
        W = Z + X
        MYARC = DSQRT(X) * RC(Z,W,ERRTOL,IERR)
        IBMARC = DATAN(X)
        WRITE(PRINTR,112) X,IBMARC,MYARC
  112   FORMAT(1H0,1PD9.1,0P2D27.16)
  116   CONTINUE
      STOP
      END
C
   1.D-75
   1.D-15
   1.D-03
   1.D-01
   2.D-01
   5.D-01
   1.D+00
   2.D+00
   5.D+00
   1.D+01
   1.D+03
   1.D+15
   1.D+75
C
C
C     $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C
C THIS DRIVER COMPARES VALUES OF LEGENDRE'S COMPLETE ELLIPTIC
C INTEGRAL K COMPUTED FROM THE SUBROUTINE RF AND VALUES COMPUTED FROM
C THE FUNPACK SUBROUTINE DELIKM.  THE FIRST TWENTY VALUES OF THE
C SQUARED COMPLEMENTARY MODULUS ETA = KC * KC ARE IN THE UNIT INTERVAL,
C AND THE VALUES OF K MAY BE COMPARED WITH 16S VALUES IN M. ABRAMOWITZ
C AND I. A. STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS (AMS 55).
C THE NEXT TWENTY-NINE VALUES OF KC RANGE FROM 1.D-70 TO 1.D+70.
C DELIKM CAN ACCEPT ONLY EIGHT OF THESE ON THE IBM 360/370 SERIES,
C PARTLY BECAUSE IT DOES NOT ADMIT VALUES OF ETA GREATER THAN UNITY
C AND PARTLY BECAUSE KC MUST BE SQUARED, BUT RF CAN ACCEPT ALL OF
C THEM IF HOMOGENEITY IS USED TO AVOID THE SQUARE.
C
      INTEGER I,IERR,IPOWER,PRINTR
      DOUBLE PRECISION DUPLIC,ERRTOL,FUNPK,RKC,X,Y,Z
      DOUBLE PRECISION DELIKM,RF
C
      DATA PRINTR/6/
C
      WRITE(PRINTR,8)
    8 FORMAT(1H1,47HCOMPARISON OF FUNPACK AND DUPLICATION ALGORITHM)
      WRITE(PRINTR,10)
   10 FORMAT(1H0,22HK(1-ETA) = RF(0,ETA,1))
      WRITE(PRINTR,12)
   12 FORMAT(1H0,12X,3HETA,17X,21HDUPLICATION ALGORITHM,13X,7HFUNPACK)
C
      X = 0.D0
      Z = 1.D0
      ERRTOL = 1.D-3
      DO 20 I = 1,20
        Y = I * .05D0
        FUNPK = DELIKM(Y)
        DUPLIC = RF(X,Y,Z,ERRTOL,IERR)
        IF (IERR .EQ. 0) GO TO 15
        WRITE(PRINTR,14) Y,FUNPK
   14   FORMAT(1H0,D27.16,27X,D27.16)
        GO TO 20
   15   WRITE(PRINTR,16) Y,DUPLIC,FUNPK
   16   FORMAT(1H0,3D27.16)
   20   CONTINUE
C
C*****************************************************************
C
      WRITE(PRINTR,104)
  104 FORMAT(1H1,48HCOMPARISON OF FUNPACK AND DUPLICATION ALGORITHMS)
      WRITE(PRINTR,108)
  108 FORMAT(1H0,46HK((1-KC**2)**0.5) = 1/SQRT(KC) * RF(0,KC,1/KC))
      WRITE(PRINTR,109)
  109 FORMAT(1H0,13X,2HKC,17X,21HDUPLICATION ALGORITHM,13X,7HFUNPACK)
C
      X = 0.D0
      ERRTOL = 1.D-3
      DO 130 I = 1,29
        IPOWER = -(70 - (5 * (I - 1)))
        RKC = 10.D0 ** IPOWER
        Y = RKC
        Z = 1.D0 / RKC
        DUPLIC = 1.D0 / DSQRT(RKC) * RF(X,Y,Z,ERRTOL,IERR)
        IF (IPOWER .GE. -35 .AND. IPOWER .LE. 1) GO TO 114
        IF (IERR .NE. 0) GO TO 112
        WRITE(PRINTR,110) RKC,DUPLIC
  110   FORMAT(1H0,2D27.16)
        GO TO 130
  112   WRITE(PRINTR,113) RKC
  113   FORMAT(1H0,D27.16)
        GO TO 130
  114   Y = RKC * RKC
        FUNPK = DELIKM(Y)
        IF (IERR .NE. 0) GO TO 118
        WRITE(PRINTR,116) RKC,DUPLIC,FUNPK
  116   FORMAT(1H0,3D27.16)
        GO TO 130
  118   WRITE(PRINTR,120) RKC,FUNPK
  120   FORMAT(1H0,D27.16,27X,D27.16)
  130   CONTINUE
C
      STOP
      END