C  Program to test CSIN/CCOS
C
C  Accuracy tests are based on the identities
C
C            sin(z) = sin(z-w)*cos(w) + cos(z-w)*sin(w)
C     and
C            cos(z) = cos(z-w)*cos(w) - sin(z-w)*sin(w)
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 five
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              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, AIMAG (or DIMAG), CMPLX (or DCMPLX), COS,
C
C         LOG, MAX, REAL (or DBLE), SIN
C
C
C  Latest revision - October 31, 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   ALXMAX,BETA,BETAP,C,CONV,CONVR,CONVI,DELX,DELY,ELEVEN,EPS,
     2   EPSNEG,HALF,ONE,REN,SIXTEN,TEN,W,X,XL,XMAX,XMIN,Y,ZERO
CS    COMPLEX
      COMPLEX*16
     1   CC,CONVC,CDEL,SINDEL,COSDEL,CS,CY,SIN11A,SIN11B,SIN22A,
     2   SIN22B,ZDUM
CC   3   ,FSIN,FCOS
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ZERO,HALF,ONE/0.0E0,0.5E0,1.0E0/
CS    DATA TEN,ELEVEN,SIXTEN/10.0E0,11.0E0,16.0E0/
CS    DATA SINDEL/(6.2581348413276935585E-2,6.2418588008436587236E-2)/
CS    DATA COSDEL/(-2.5431314180235545803E-6,-3.9062493377261771826E-3)/
CS    DATA SIN11A/(-134167.0E0,+593.0E0)/
CS    DATA SIN22A/(-1187.0E0,-134163.0E0)/
CS    DATA SIN11B/(-3.2928849557976510994E-1,7.8989452897413630729E-1)/
CS    DATA SIN22B/(-5.6815858848429427210E-1,-3.8738909081835398666E-1)/
CS    DATA CDEL/(0.0625,0.0625)/
      DATA ZERO,HALF,ONE/0.0D0,0.5D0,1.0D0/
      DATA TEN,ELEVEN,SIXTEN/10.0D0,11.0D0,16.0D0/
      DATA SINDEL/(6.2581348413276935585D-2,6.2418588008436587236D-2)/
      DATA COSDEL/(-2.5431314180235545803D-6,-3.9062493377261771826D-3)/
      DATA SIN11A/(-134167.0D0,+593.0D0)/
      DATA SIN22A/(-1187.0D0,-134163.0D0)/
      DATA SIN11B/(-3.2928849557976510994D-1,7.8989452897413630729D-1)/
      DATA SIN22B/(-5.6815858848429427210D-1,-3.8738909081835398666D-1)/
      DATA CDEL/(0.0625,0.0625)/
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)
CS    CONVR(ZDUM) = REAL(ZDUM)
CS    CONVI(ZDUM) = AIMAG(ZDUM)
      CONV(NDUM) = DBLE(NDUM)
      CONVC(X,Y) = DCMPLX(X,Y)
      CONVR(ZDUM) = DBLE(ZDUM)
      CONVI(ZDUM) = DIMAG(ZDUM)
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)
      AX = ONE/SIXTEN
      BX = ONE
      AY = ZERO
      BY = ONE
      N = 2000
      XN = CONV(N)
      I1 = 0
      IFLAG = .FALSE.
      RFLAG = .TRUE.
      PFLAG = .FALSE.
C-----------------------------------------------------------------
C  Random argument accuracy tests
C-----------------------------------------------------------------
      DO 300 J = 1, 3
         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) + AY
            CX = CONVC(X,Y)
C----------------------------------------------------------------
C  Carefully purify arguments and evaluate identities
C----------------------------------------------------------------
            CY = CX - CDEL
            CS = SIN(CY)
            CC = COS(CY)
CC          CS = FSIN(CY)
CC          CC = FCOS(CY)
            IF (J .NE. 3) THEN
                  CZ = SIN(CX)
CC                CZ = FSIN(CX)
                  CZZ = CS + (CS*COSDEL+CC*SINDEL)
               ELSE
                  CZ = COS(CX)
CC                CZ = FCOS(CX)
                  CZZ = CC + (CC*COSDEL-CS*SINDEL)
            END IF
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
            CALL TABLAT
            XL = XL + DELX
  200    CONTINUE
C---------------------------------------------------------------------
C  Process and output statistics
C---------------------------------------------------------------------
         IF (J .NE. 3) THEN
               WRITE (IOUT,1000)
            ELSE
               WRITE (IOUT,1005)
         END IF
         CALL REPORT(IOUT)
         IF (J .EQ. 1) THEN
               AX = SIXTEN
               BX = AX+ONE
         END IF
         AY = AX
         BY = BX
         IFLAG = .TRUE.
  300 CONTINUE
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
      WRITE (IOUT,1025)
      C = ONE / BETA ** (IT/2)
      X = ELEVEN
      Y = X + ONE + HALF
      X = X + X
      CX = CONVC(X+C,Y)
      CY = CONVC(X-C,Y)
      CZ = (SIN(CX) - SIN(CY)) / (C + C)
      WRITE (IOUT,1026) CZ
      WRITE (IOUT,1030)
      DO 320 I = 1, 5
         X = REN(I1) * TEN
         Y = REN(I1) * TEN
         CX = CONVC(X,Y)
         CZ = SIN(CX) + SIN(-CX)
         WRITE (IOUT,1060) CX, CZ
  320 CONTINUE
      WRITE (IOUT,1040)
      X = ELEVEN
      Y = X + ONE + HALF
      CX = CONVC(X,Y)
      CZ = SIN (CX)
      CS = (CZ-SIN11A)-SIN11B
      C = CONVI(CS)/CONVI(CZ)
      IF (C .NE. ZERO) THEN
            W = LOG(ABS(C))/ALBETA+AIT
         ELSE
            W = ZERO
      END IF
      W = MAX(W,ZERO)
      WRITE (IOUT,1042) CX,CZ,IBETA,W
      IF (W .LT. ONE) THEN
            WRITE (IOUT,1046)
         ELSE
            WRITE (IOUT,1048)
      END IF
      X = X + X
      CX = CONVC(X,Y)
      CZ = SIN (CX)
      CS = (CZ-SIN22A)-SIN22B
      C = CONVR(CS)/CONVR(CZ)
      IF (C .NE. ZERO) THEN
            W = LOG(ABS(C))/ALBETA+AIT
         ELSE
            W = ZERO
      END IF
      W = MAX(W,ZERO)
      WRITE (IOUT,1044) CX,CZ,IBETA,W
      IF (W .LT. ONE) THEN
            WRITE (IOUT,1046)
         ELSE
            WRITE (IOUT,1048)
      END IF
C-----------------------------------------------------------------
C  Test of error returns
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      ALXMAX = LOG(XMAX)-HALF
      CX = CONVC(ONE,ALXMAX)
      WRITE (IOUT,1053) CX
      CY = SIN(CX)
      WRITE (IOUT,1055) CY
      BETAP = BETA ** IT
      CX = CONVC(BETAP,ONE)
      WRITE (IOUT,1052) CX
      CY = SIN(CX)
      WRITE (IOUT,1055) CY
      CX = CONVC(ONE,BETAP)
      WRITE (IOUT,1052) CX
      CY = SIN(CX)
      WRITE (IOUT,1055) CY
      WRITE (IOUT,1100)
      STOP
C----------------------------------------------------------------
 1000 FORMAT(' Test of SIN(X) vs SIN(W+D), W = X-D, D = (1/16,1/16)'//)
 1005 FORMAT(' Test of COS(X) vs COS(W+D), W = X-D, D = (1/16,1/16)'//)
 1025 FORMAT(' Special Tests'//)
 1026 FORMAT(' If (',(F9.1,F7.1),') is not almost (-134163,1188),'/
     1   14X,25HSIN has the wrong period.  //)
 1030 FORMAT(' The identity   SIN(-Z) = -SIN(Z)   will be tested.'//
     1      16X,'X',23X,'F(Z) + F(-Z)'/)
 1040 FORMAT(/' Test the significance in the real or imaginary',
     1  ' components near'/' the zeros of the real sin',
     2  ' or cos.'/)
 1042 FORMAT(' SIN(',(E15.7,E15.7),') = (',(E15.7,E15.7),')'/
     1  ' The estimated loss of base ',I4,' digits in the',
     2  ' imaginary component is ',F7.2)
 1044 FORMAT(' SIN(',(E15.7,E15.7),') = (',(E15.7,E15.7),')'/
     1  ' The estimated loss of base ',I4,' digits in the',
     2  ' real component is ',F7.2)
 1046 FORMAT(' This is acceptable'/)
 1048 FORMAT(' This is too large'/)
 1050 FORMAT(' Test of error returns'//)
 1052 FORMAT(' SIN will be called with the argument (',(E15.4,E15.4),
     1       ')'/' This should trigger an error message'//)
 1053 FORMAT(' SIN will be called with the argument (',(E15.4,E15.4),
     1       ')'/' This should not trigger an error message'//)
 1055 FORMAT(' SIN returned the value (',(E15.4,E15.4),')'///)
 1060 FORMAT(2(E15.7,E15.7)/)
 1100 FORMAT(' This concludes the tests')
C---------- Last card of CSIN/CCOS Test Program ----------
      END
 
CC    COMPLEX FUNCTION FSIN(CX)
C----------------------------------------------------------------------
C  Calibration function for csin tests.
C
C  Latest revision - October 31, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
CC    COMPLEX CX
CC    COMPLEX*16 DCX
C----------------------------------------------------------------------
CC    DCX = CX
CC    FSIN = SIN(DCX)
CC    RETURN
C---------- Last line of FSIN ----------
CC    END
 
CC    COMPLEX FUNCTION FCOS(CX)
C----------------------------------------------------------------------
C  Calibration function for ccos tests.
C
C  Latest revision - October 31, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
CC    COMPLEX CX
CC    COMPLEX*16 DCX
C----------------------------------------------------------------------
CC    DCX = CX
CC    FCOS = COS(DCX)
CC    RETURN
C---------- Last line of FCOS ----------
CC    END
