/*
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----------------------------------------------------------------------
*/

#include "celefunt.h"

int IFLAG, RFLAG, PFLAG;
int KX1, KX2, KX3, KY1, KY2, KY3;
#ifdef TDOUBLE
double AX, AY, BX, BY, RX6, RX7, RY6, RY7;
double complex CF1, CF2, CF3, CV, CX, CX1, CX2, CX3, CY1, CY2, CY3, CZ,
  CZZ;
#endif
#ifdef TSINGLE
float AX, AY, BX, BY, RX6, RX7, RY6, RY7;
float complex CF1, CF2, CF3, CV, CX, CX1, CX2, CX3, CY1, CY2, CY3, CZ,
  CZZ;
#endif
/*
C----------------------------------------------------------------------
C  Declare local variables
C----------------------------------------------------------------------
*/
#ifdef TDOUBLE
double SCALE;
#endif
#ifdef TSINGLE
float SCALE;
#endif

/*
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
*/
/*
C----------------------------------------------------------------------
C  Statement functions for conversion between integer or complex
C  and float
C----------------------------------------------------------------------
*/
/*
C----------------------------------------------------------------------
*/
int
main ()
{
  /*
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
  */
  IBETA = 2;
  IT = IT_VAL;
  XMIN = XMIN_VAL;
  XMAX = XMAX_VAL;
  BETA = CONV (IBETA);
  ALBETA = LOG (BETA);
  AIT = CONV (IT);
  I1 = IT / 2 + 4;
  /* SCALE = BETA ** I1; */
#ifdef TDOUBLE
  SCALE = POW (BETA, (double) I1);
#endif
#ifdef TSINGLE
  SCALE = POW (BETA, (float) I1);
#endif
  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*/
  for (J = 1; J <= 2; J++)
    {
      RESET ();
      DELX = (BX - AX) / XN;
      DELY = BY - AY;
      XL = AX;
      /*         DO 200 I = 1, N */
      for (I = 1; I <= N; I++)
 {
   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;
   CZZ = CONVC (X, Y);
   Z = X * X - Y * Y;
   W = X * Y;
   /* For TABLAT, the function argument must be named CX.  */
   CX = CONVC (Z, W + W);
   if (J == 1)
     CZ = FSQRT (CX);
   else
     {
       CZ = -FSQRT (CX);
#if 0
       /* This did not seem right in the Fortran program.
   The question is, Is CZ a square root of CX?
 Fix so TABLAT compares the square of CZ with CX.  */
       CZZ = CX;
       CX = CZ;
       CZ = CZ * CZ;
#endif
     }
/*
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
*/
   TABLAT ();
   XL = XL + DELX;
   /*  200    CONTINUE */
 }
/*
C-----------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------
*/
      if (J == 1)
 printf (" Test of SQRT(Z*Z) vs Z\n");
      else
 printf (" Test of -SQRT(Z*Z) vs Z\n");

      REPORT (IOUT);
      BX = -BX * BX;
      BY = -BX;
      /*  300 CONTINUE */
    }
/*
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
*/
  printf (" Test of special arguments\n");
  CX = CONVC (ONE, ZERO);
  CZ = FSQRT (CX);
  printf (" SQRT(1.0,0.0) = (%15.7E,%15.7E)\n", creal (CZ), cimag (CZ));
  printf (" Tests near extreme arguments.  These should not\n");
  printf (" but may\n trigger error messages\n");
  CX = CONVC (XMIN, XMIN);
  printf (" SQRT will be called with the argument (%15.4E,%15.4E)\n",
   creal (CX), cimag (CX));
  CZ = FSQRT (CX);
  printf (" SQRT(XMIN,XMIN) = SQRT(%15.7E,%15.7E) = (%15.7E,%15.7E)\n",
   creal (CX), cimag (CX), creal (CZ), cimag (CZ));
  CX = CONVC (XMAX, XMAX);
  printf (" SQRT will be called with the argument (%15.4E,%15.4E)\n",
   creal (CX), cimag (CX));
  CZ = FSQRT (CX);
  printf (" SQRT(XMAX,XMAX) = SQRT(%15.7E,%15.7E) = (%15.7E,%15.7E)\n",
   creal (CX), cimag (CX), creal (CZ), cimag (CZ));
  printf (" This concludes the tests\n");
  exit (0);
}

/*
C---------------------------------------------------------------------
*/
/*
C--------- Last card of CSQRT Test Program ----------
*/

#include "tablat.c"