C     PROGRAM TO TEST ATAN, ATAN2
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 SIX
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                 IRND   - 0 IF FLOATING-POINT ADDITION CHOPS,
C                          1 IF FLOATING-POINT ADDITION ROUNDS
C                 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE
C                          INTEGER SUCH THAT FLOAT(IBETA)**MINEXP
C                          IS A POSITIVE 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     STANDARD FORTRAN SUBPROGRAMS REQUIRED
C
C         ABS, ALOG, AMAX1, ATAN, ATAN2, FLOAT, SQRT
C
C
C     LATEST REVISION - DECEMBER 6, 1979
C
C     AUTHOR - W. J. CODY
C              ARGONNE NATIONAL LABORATORY
C
C
      INTEGER I,IBETA,IEXP,IOUT,IRND,II,IT,I1,J,K1,K2,K3,MACHEP,
     1        MAXEXP,MINEXP,N,NEGEP,NGRD
      REAL A,AIT,ALBETA,B,BETA,BETAP,DEL,EM,EPS,EPSNEG,EXPON,HALF,OB32,
     1   ONE,RAN,R6,R7,SUM,TWO,W,X,XL,XMAX,XMIN,XN,XSQ,X1,Y,Z,ZERO,ZZ
C
      IOUT = 6
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
      BETA = FLOAT(IBETA)
      ALBETA = ALOG(BETA)
      AIT = FLOAT(IT)
      ONE = 1.0E0
      HALF = 0.5E0
      TWO = 2.0E0
      ZERO = 0.0E0
      A = -0.0625E0
      B = -A
      OB32 = B * HALF
      N = 2000
      XN = FLOAT(N)
      I1 = 0
C-----------------------------------------------------------------
C     RANDOM ARGUMENT ACCURACY TESTS
C-----------------------------------------------------------------
      DO 300 J = 1, 4
         K1 = 0
         K3 = 0
         X1 = ZERO
         R6 = ZERO
         R7 = ZERO
         DEL = (B - A) / XN
         XL = A
C
         DO 200 I = 1, N
            X = DEL * REN(I1) + XL
            IF (J .EQ. 2) X = ((1.0E0+X*A)-ONE)*16.0E0
            Z = ATAN(X)
            IF (J .NE. 1) GO TO 100
            XSQ = X * X
            EM = 17.0E0
            SUM = XSQ / EM
C
            DO 80 II = 1, 7
               EM = EM - TWO
               SUM = (ONE/EM - SUM) * XSQ
   80       CONTINUE
C
            SUM = -X * SUM
            ZZ = X + SUM
            SUM = (X - ZZ) + SUM
            IF (IRND .EQ. 0) ZZ = ZZ + (SUM + SUM)
            GO TO 110
  100       IF (J .NE. 2) GO TO 105
            Y = X - .0625E0
            Y = Y / (ONE + X*A)
            ZZ = (ATAN(Y) - 8.1190004042651526021E-5) + OB32
            ZZ = ZZ + OB32
            GO TO 110
  105       Z = Z + Z
            Y = X / ((HALF + X * HALF)*((HALF - X) + HALF))
            ZZ = ATAN(Y)
  110       W = ONE
            IF (Z .NE. ZERO) W = (Z - ZZ) / Z
            IF (W .GT. ZERO) K1 = K1 + 1
            IF (W .LT. ZERO) K3 = K3 + 1
            W = ABS(W)
            IF (W .LE. R6) GO TO 120
            R6 = W
            X1 = X
  120       R7 = R7 + W * W
            XL = XL + DEL
  200    CONTINUE
C
         K2 = N - K3 - K1
         R7 = SQRT(R7/XN)
         IF (J .EQ. 1) WRITE (IOUT,1000)
         IF (J .EQ. 2) WRITE (IOUT,1001)
         IF (J .GT. 2) WRITE (IOUT,1002)
         WRITE (IOUT,1010) N,A,B
         WRITE (IOUT,1011) K1,K2,K3
         WRITE (IOUT,1020) IT,IBETA
         W = -999.0E0
         IF (R6 .NE. ZERO) W = ALOG(ABS(R6))/ALBETA
         WRITE (IOUT,1021) R6,IBETA,W,X1
         W = AMAX1(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         W = -999.0E0
         IF (R7 .NE. ZERO) W = ALOG(ABS(R7))/ALBETA
         WRITE (IOUT,1023) R7,IBETA,W
         W = AMAX1(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         A = B
         IF (J .EQ. 1) B = TWO - SQRT(3.0E0)
         IF (J .EQ. 2) B = SQRT(TWO) - ONE
         IF (J .EQ. 3) B = ONE
  300 CONTINUE
C-----------------------------------------------------------------
C     SPECIAL TESTS
C-----------------------------------------------------------------
      WRITE (IOUT,1025)
      WRITE (IOUT,1030)
      A = 5.0E0
C
      DO 320 I = 1, 5
         X = REN(I1) * A
         Z = ATAN(X) + ATAN(-X)
         WRITE (IOUT,1060) X, Z
  320 CONTINUE
C
      WRITE (IOUT,1031)
      BETAP = BETA ** IT
      X = REN(I1) / BETAP
C
      DO 330 I = 1, 5
         Z = X - ATAN(X)
         WRITE (IOUT,1060) X, Z
         X = X / BETA
  330 CONTINUE
C
      WRITE (IOUT,1032)
      A = -TWO
      B = 4.0E0
C
      DO 340 I = 1, 5
         X = REN(I1) * B + A
         Y = REN(I1)
         W = -Y
         Z = ATAN(X/Y) - ATAN2(X,Y)
         ZZ = ATAN(X/W) - ATAN2(X,W)
         WRITE (IOUT,1059) X, Y, Z, ZZ
  340 CONTINUE
C
      WRITE (IOUT,1035)
      EXPON = FLOAT(MINEXP) * 0.75E0
      X = BETA ** EXPON
      Y = ATAN(X)
      WRITE (IOUT,1061) X, Y
C-----------------------------------------------------------------
C     TEST OF ERROR RETURNS
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      WRITE (IOUT,1051) XMAX
      Z = ATAN(XMAX)
      WRITE (IOUT,1061) XMAX, Z
      X = ONE
      Y = ZERO
      WRITE (IOUT,1053) X, Y
      Z = ATAN2(X,Y)
      WRITE (IOUT,1062) X, Y, Z
      WRITE (IOUT,1053) XMIN, XMAX
      Z = ATAN2(XMIN,XMAX)
      WRITE (IOUT,1062) XMIN, XMAX, Z
      WRITE (IOUT,1053) XMAX, XMIN
      Z = ATAN2(XMAX,XMIN)
      WRITE (IOUT,1062) XMAX, XMIN, Z
      X = ZERO
      WRITE (IOUT,1054) X, Y
      Z = ATAN2(X,Y)
      WRITE (IOUT,1062) X, Y, Z
      WRITE (IOUT,1100)
      STOP
 1000 FORMAT(43H1TEST OF ATAN(X) VS TRUNCATED TAYLOR SERIES  //)
 1001 FORMAT(20H1TEST OF ATAN(X) VS ,
     1  40HATAN(1/16) + ATAN((X-1/16)/(1+X/16))    //)
 1002 FORMAT(40H1TEST OF 2*ATAN(X) VS ATAN(2X/(1-X*X))    //)
 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL /
     1 6X,1H(,E15.4,1H,,E15.4,1H)//)
 1011 FORMAT(19H ATAN(X) WAS LARGER,I6,7H TIMES,/
     1    12X,7H AGREED,I6,11H TIMES, AND /
     1   8X,11HWAS SMALLER,I6,7H TIMES.//)
 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4,
     1    46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER  //)
 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **,
     1  F7.2/4X,16HOCCURRED FOR X =,E17.6)
 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4,
     1  22H SIGNIFICANT DIGITS IS,F7.2//)
 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,E15.4,
     1    3H = ,I4,3H **,F7.2)
 1025 FORMAT(14H1SPECIAL TESTS//)
 1030 FORMAT(53H THE IDENTITY   ATAN(-X) = -ATAN(X)   WILL BE TESTED.//
     1      8X,1HX,9X,12HF(X) + F(-X)/)
 1031 FORMAT(52H THE IDENTITY ATAN(X) = X , X SMALL, WILL BE TESTED.//
     1       8X,1HX,9X,8HX - F(X)/)
 1032 FORMAT(51H THE IDENTITY ATAN(X/Y) = ATAN2(X,Y) WILL BE TESTED /
     1 57H THE FIRST COLUMN OF RESULTS SHOULD BE 0, THE SECOND +-PI//
     2 8X,1HX,13X,1HY,5X,15HF1(X/Y)-F2(X,Y),18HF1(X/Y)-F2(X/(-Y)) /)
 1035 FORMAT(43H TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT. /)
 1050 FORMAT(22H1TEST OF ERROR RETURNS //)
 1051 FORMAT(38H ATAN WILL BE CALLED WITH THE ARGUMENT,E15.7/
     1       41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE/)
 1053 FORMAT(40H ATAN2 WILL BE CALLED WITH THE ARGUMENTS/2E15.7/
     1       41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE/)
 1054 FORMAT(40H ATAN2 WILL BE CALLED WITH THE ARGUMENTS/2E15.7/
     1       37H THIS SHOULD TRIGGER AN ERROR MESSAGE//)
 1059 FORMAT(4E15.7/)
 1060 FORMAT(2E15.7/)
 1061 FORMAT(6X,6H ATAN(,E13.6,3H) =,E13.6/)
 1062 FORMAT(6X,7H ATAN2(,2E13.6,3H) =,E13.6/)
 1100 FORMAT(25H THIS CONCLUDES THE TESTS )
C     ---------- LAST CARD OF ATAN/ATAN2 TEST PROGRAM ----------
      END
      REAL FUNCTION REN(K)
C
C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
C      VOL. 8, NO. 10, OCTOBER 1965.
C
C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
C      29 BITS.
C
      INTEGER IY,J,K
      DATA IY/100001/
C
      J = K
      IY = IY * 125
      IY = IY - (IY/2796203) * 2796203
      REN = FLOAT(IY) / 2796203.0E0
      RETURN
C     ---------- LAST CARD OF RAN ----------
      END