C----------------------------------------------------------------------
C  Program to test CSQRT
C
C  Accuracy tests are based on the identity
C
C            z = sqrt(z*z)
C
C  Data required:
C
C     None
C
C  Subprograms required from this package:
C
C     MACHAR - An environmental inquiry program providing
C              information on the floating point arithmetic
C              system.  Note that the call to MACHAR can
C              be deleted provided the following four
C              parameters are assigned the values indicated:
C
C              IBETA  - the radix of the floating point system;
C              IT     - the number of base-IBETA digits in the
C                       significand of a floating-point number;
C              XMIN   - the smallest non-vanishing floating-point
C                       power of the radix;
C              XMAX   - the largest finite floating-point no.
C
C     REN(K) - a function subprogram returning random real
C              numbers uniformly distributed over (0,1)
C
C     RESET, TABLAT, REPORT - programs to report results.
C
C
C  Generic Fortran subprograms required:
C
C         ABS, CMPLX (or DCMPLX), LOG, REAL (or DBLE), SQRT
C
C  The program as given here must be modified before compiling.
C   If a single (double) precision version is desired, change all
C   occurrences of CS (CD) in columns 1 and 2 to blanks.  For a
C   calibration version, change all occurrences of CS and of CC
C   in columns 1 and 2 to blanks.
C
C  Latest revision - November 15, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
C  Declare named common
C----------------------------------------------------------------------
      COMMON /ELEFUN/ CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,
     1       CZ,CZZ,AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN,
     2       IBETA,IT,KX1,KX2,KX3,KY1,KY2,KY3,K1,K2,K3,N,
     3       IFLAG,RFLAG,PFLAG
      LOGICAL IFLAG,RFLAG,PFLAG
      INTEGER IBETA,IT,KX1,KX2,KX3,KY1,KY2,KY3,K1,K2,K3,N
CS    REAL
      DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
      COMPLEX*16
     1       CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,CZ,CZZ
C----------------------------------------------------------------------
C  Declare local variables
C----------------------------------------------------------------------
      INTEGER I,IEXP,IOUT,IRND,I1,J,MACHEP,MAXEXP,MINEXP,NDUM,
     1        NEGEP,NGRD
CS    REAL
      DOUBLE PRECISION
     1   BETA,CONV,DELX,DELY,EPS,EPSNEG,ONE,REN,SCALE,TEN,
     2   W,X,XL,XMAX,XMIN,Y,Z,ZERO
CS    COMPLEX
      COMPLEX*16
     1   CONVC,CY
CC   2   ,FSQRT
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ZERO,ONE,TEN/0.0E0,1.0E0,10.0E0/
      DATA ZERO,ONE,TEN/0.0D0,1.0D0,10.0D0/
C----------------------------------------------------------------------
C  Statement functions for conversion between integer or complex
C  and float
C----------------------------------------------------------------------
CS    CONV(NDUM) = REAL(NDUM)
CS    CONVC(X,Y) = CMPLX(X,Y)
      CONV(NDUM) = DBLE(NDUM)
      CONVC(X,Y) = DCMPLX(X,Y)
C----------------------------------------------------------------------
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
      BETA = CONV(IBETA)
      ALBETA = LOG(BETA)
      AIT = CONV(IT)
      I1 = IT/2 + 4
      SCALE = BETA ** I1
      AX = ZERO
      BX = TEN
      AY = AX
      BY = BX
      N = 2000
      XN = CONV(N)
      I1 = 0
      RFLAG = .TRUE.
      IFLAG = .TRUE.
      PFLAG = .FALSE.
C---------------------------------------------------------------------
C  Random argument accuracy tests
C---------------------------------------------------------------------
      DO 300 J = 1, 2
         CALL RESET
         DELX = (BX - AX) / XN
         DELY = BY - AY
         XL = AX
         DO 200 I = 1, N
            X = DELX * REN(I1) + XL
            Y = DELY * REN(I1)
C-----------------------------------------------------------------
C  Carefully purify arguments and evaluate identities
C-----------------------------------------------------------------
            Z = X * SCALE
            W = Z + X
            X = W - Z
            Z = Y * SCALE
            W = Z + Y
            Y = W - Z
            CX = CONVC(X,Y)
            Z = X*X - Y*Y
            W = X*Y
            CY = CONVC(Z,W+W)
            CZZ = CX
            IF (J .EQ. 1) THEN
                  CZ = SQRT(CY)
CC                CZ = FSQRT(CY)
               ELSE
                  CZ = -SQRT(CY)
CC                CZ = -FSQRT(CY)
            END IF
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
            CALL TABLAT
            XL = XL + DELX
  200    CONTINUE
C-----------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------
         IF (J .EQ. 1) THEN
               WRITE (IOUT,1000)
            ELSE
               WRITE (IOUT,1001)
         END IF
         CALL REPORT(IOUT)
         BX = -BX * BX
         BY = -BX
  300 CONTINUE
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
      WRITE (IOUT,1040)
      CX = CONVC(ONE,ZERO)
      CZ = SQRT(CX)
      WRITE (IOUT,1041) CZ
      WRITE (IOUT,1050)
      CX = CONVC(XMIN,XMIN)
      WRITE (IOUT,1061) CX
      CZ = SQRT(CX)
      WRITE (IOUT,1051) CX,CZ
      CX = CONVC(XMAX,XMAX)
      WRITE (IOUT,1061) CX
      CZ = SQRT(CX)
      WRITE (IOUT,1052) CX,CZ
      WRITE (IOUT,1100)
      STOP
C---------------------------------------------------------------------
 1000 FORMAT(' Test of SQRT(Z*Z) vs Z'//)
 1001 FORMAT(' Test of -SQRT(Z*Z) vs Z'//)
 1040 FORMAT(//' Test of special arguments'//)
 1041 FORMAT(' SQRT(1.0,0.0) = (',(E15.7,E15.7),')'//)
 1050 FORMAT(//' Tests near extreme arguments.  These should not,',
     1         ' but may'/' trigger error messages'//)
 1051 FORMAT(' SQRT(XMIN,XMIN) = SQRT(',(E15.7,E15.7),') ='/40X,'(',
     1         (E15.7,E15.7),')'/)
 1052 FORMAT(' SQRT(XMAX,XMAX) = SQRT(',(E15.7,E15.7),') ='/40X,'(',
     1         (E15.7,E15.7),')'/)
 1061 FORMAT(/' SQRT will be called with the argument (',
     1          (E15.4,E15.4),')'//)
 1100 FORMAT(' This concludes the tests')
C--------- Last card of CSQRT Test Program ----------
      END
