home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / TPL60N11.ZIP / TESTPRGS.ZIP / DATAN.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1992-03-28  |  22.4 KB  |  671 lines
  1. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  2.  
  3.  
  4. (*
  5. C     PROGRAM TO TEST DATAN, DATAN2
  6. C
  7. C     DATA REQUIRED
  8. C
  9. C        NONE
  10. C
  11. C     SUBPROGRAMS REQUIRED FROM THIS PACKAGE
  12. C
  13. C        MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
  14. C                 INFORMATION ON THE FLOATING-POINT ARITHMETIC
  15. C                 SYSTEM.  NOTE THAT THE CALL TO MACHAR CAN
  16. C                 BE DELETED PROVIDED THE FOLLOWING SIX
  17. C                 PARAMETERS ARE ASSIGNED THE VALUES INDICATED
  18. C
  19. C                 IBETA  - THE RADIX OF THE FLOATING-POINT SYSTEM
  20. C                 IT     - THE NUMBER OF BASE-IBETA DIGITS IN THE
  21. C                          SIGNIFICAND OF A FLOATING-POINT NUMBER
  22. C                 IRND   - 0 IF FLOATING-POINT ADDITION CHOPS,
  23. C                          1 IF FLOATING-POINT ADDITION ROUNDS
  24. C                 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE
  25. C                          INTEGER SUCH THAT DFLOAT(IBETA)**MINEXP
  26. C                          IS A POSITIVE FLOATING-POINT NUMBER
  27. C                 XMIN   - THE SMALLEST NON-VANISHING FLOATING-POINT
  28. C                          POWER OF THE RADIX
  29. C                 XMAX   - THE LARGEST FINITE FLOATING-POINT NO.
  30. C
  31. C        REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
  32. C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
  33. C
  34. C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
  35. C
  36. C         DABS, DLOG, DMAX1, DATAN, DATAN2, DFLOAT, DSQRT
  37. C
  38. C
  39. C     LATEST REVISION - DECEMBER 6, 1979
  40. C
  41. C     AUTHOR - W. J. CODY
  42. C              ARGONNE NATIONAL LABORATORY
  43. C
  44. C
  45. *)
  46.  
  47.  
  48.  
  49. FUNCTION REN (K: LONGINT): REAL;
  50.  
  51. (*
  52.       DOUBLE PRECISION FUNCTION REN(K)
  53. C
  54. C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
  55. C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
  56. C      VOL. 8, NO. 10, OCTOBER 1965.
  57. C
  58. C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
  59. C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
  60. C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
  61. C      29 BITS.
  62. C
  63. *)
  64.  
  65.  
  66.  
  67. VAR   J: LONGINT;
  68. CONST IY: LONGINT = 100001;
  69.  
  70. BEGIN
  71.       J  := K;
  72.       IY := IY * 125;
  73.       IY := IY - (IY DIV 2796203) * 2796203;
  74.       REN:= 1.0 * (IY) / 2796203.0e0 * (1.0e0 + 1.0e-6 + 1.0e-12);
  75. END;
  76.  
  77.  
  78.  
  79.  
  80.    function Power (X, Y: real): real;
  81.  
  82.   var z:real;
  83.     iy:integer;
  84.     reciproc:boolean;
  85.       (* Crude power function:  see Cody & Waite for a better one. *)
  86.  
  87.       begin (* Power *)
  88.       if Y = 0.0 then
  89.          Power := 1.0
  90.       else if (X = 0.0) and (Y > 0.0) then
  91.          Power := 0.0
  92.       else if (Y = Int (Y)) AND (Abs(Y) <= MAXINT) THEN BEGIN
  93.          IY := TRUNC (Y);
  94.          reciproc := iy < 0;
  95.          iy := abs (iy);
  96.          z := 1.0;
  97.          while iy > 0 do begin
  98.            if odd(iy) then begin
  99.              z := z*x;
  100.  
  101.            end else begin
  102.            end; {endif}
  103.            iy := iy div 2;
  104.            if iy >0 then begin
  105.               x := sqr (x);
  106.            end else begin
  107.            end; {endif}
  108.  
  109.          end; {endwhile}
  110.          if reciproc then begin
  111.            power := 1 / z;
  112.  
  113.          end else begin
  114.            power := z;
  115.          end {endif}
  116.          end
  117.       else if (X < 0.0) and (Y = Int(Y)) then
  118.          if odd(trunc(Y)) then Power := -exp (Y * ln (-X))
  119.          else Power := exp (Y * ln (-X))
  120.       else
  121.          Power := exp (Y * ln (X));
  122.       end (* Power *);
  123.  
  124.  
  125. FUNCTION MAX1 (A, B:REAL): REAL;
  126. BEGIN
  127.    IF A > B THEN
  128.       MAX1 := A
  129.    ELSE
  130.       MAX1 := B;
  131. END;
  132.  
  133.  
  134.  
  135.       PROCEDURE MACHAR(VAR IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  136.                       MAXEXP: LONGINT; VAR EPS,EPSNEG,XMIN,XMAX: REAL);
  137. {
  138. C-----------------------------------------------------------------------
  139. C  This Fortran 77 subroutine is intended to determine the parameters
  140. C   of the floating-point arithmetic system specified below.  The
  141. C   determination of the first three uses an extension of an algorithm
  142. C   due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
  143. C   but not all, of the improvements suggested by M. Gentleman and S.
  144. C   Marovich, CACM 17 (1974), pp. 276-277.  An earlier version of this
  145. C   program was published in the book Software Manual for the
  146. C   Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
  147. C   Englewood Cliffs, NJ, 1980.  The present version is documented in
  148. C   W. J. Cody, "MACHAR: A subroutine to dynamically determine machine
  149. C   parameters," TOMS 14, December, 1988.
  150. C
  151. C  The program as given here must be modified before compiling.  If
  152. C   a single (double) precision version is desired, change all
  153. C   occurrences of CS (CD) in columns 1 and 2 to blanks.
  154. C
  155. C  Parameter values reported are as follows:
  156. C
  157. C       IBETA   - the radix for the floating-point representation
  158. C       IT      - the number of base IBETA digits in the floating-point
  159. C                 significand
  160. C       IRND    - 0 if floating-point addition chops
  161. C                 1 if floating-point addition rounds, but not in the
  162. C                   IEEE style
  163. C                 2 if floating-point addition rounds in the IEEE style
  164. C                 3 if floating-point addition chops, and there is
  165. C                   partial underflow
  166. C                 4 if floating-point addition rounds, but not in the
  167. C                   IEEE style, and there is partial underflow
  168. C                 5 if floating-point addition rounds in the IEEE style,
  169. C                   and there is partial underflow
  170. C       NGRD    - the number of guard digits for multiplication with
  171. C                 truncating arithmetic.  It is
  172. C                 0 if floating-point arithmetic rounds, or if it
  173. C                   truncates and only  IT  base  IBETA digits
  174. C                   participate in the post-normalization shift of the
  175. C                   floating-point significand in multiplication;
  176. C                 1 if floating-point arithmetic truncates and more
  177. C                   than  IT  base  IBETA  digits participate in the
  178. C                   post-normalization shift of the floating-point
  179. C                   significand in multiplication.
  180. C       MACHEP  - the largest negative integer such that
  181. C                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
  182. C                 MACHEP is bounded below by  -(IT+3)
  183. C       NEGEPS  - the largest negative integer such that
  184. C                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
  185. C                 NEGEPS is bounded below by  -(IT+3)
  186. C       IEXP    - the number of bits (decimal places if IBETA = 10)
  187. C                 reserved for the representation of the exponent
  188. C                 (including the bias or sign) of a floating-point
  189. C                 number
  190. C       MINEXP  - the largest in magnitude negative integer such that
  191. C                 FLOAT(IBETA)**MINEXP is positive and normalized
  192. C       MAXEXP  - the smallest positive power of  BETA  that overflows
  193. C       EPS     - the smallest positive floating-point number such
  194. C                 that  1.0+EPS .NE. 1.0. In particular, if either
  195. C                 IBETA = 2  or  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
  196. C                 Otherwise,  EPS = (FLOAT(IBETA)**MACHEP)/2
  197. C       EPSNEG  - A small positive floating-point number such that
  198. C                 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2
  199. C                 or  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
  200. C                 Otherwise,  EPSNEG = (IBETA**NEGEPS)/2.  Because
  201. C                 NEGEPS is bounded below by -(IT+3), EPSNEG may not
  202. C                 be the smallest number that can alter 1.0 by
  203. C                 subtraction.
  204. C       XMIN    - the smallest non-vanishing normalized floating-point
  205. C                 power of the radix, i.e.,  XMIN = FLOAT(IBETA)**MINEXP
  206. C       XMAX    - the largest finite floating-point number.  In
  207. C                 particular  XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
  208. C                 Note - on some machines  XMAX  will be only the
  209. C                 second, or perhaps third, largest number, being
  210. C                 too small by 1 or 2 units in the last digit of
  211. C                 the significand.
  212. C
  213. C     Latest revision - December 4, 1987
  214. C
  215. C     Author - W. J. Cody
  216. C              Argonne National Laboratory
  217. C
  218. C-----------------------------------------------------------------------
  219. }
  220. VAR
  221.        I,L, ITEMP,IZ,J,K,
  222.        MX,NXRES: LONGINT;
  223.        A,B,BETA,BETAIN,BETAH,ONE,T,TEMP,TEMPA,
  224.        TEMP1,TWO,Y,Z,ZERO: REAL;
  225.        CONV: ARRAY [0..10] OF REAL;
  226.  
  227. LABEL 10, 20,100, 210, 220, 300, 320, 400, 410, 420, 430, 440, 450, 460, 500, 520;
  228.  
  229. BEGIN
  230. {-----------------------------------------------------------------------}
  231. for l:= 1 to 10 do
  232.     CONV [l] := l;
  233.  
  234.       ONE  := CONV [1];
  235.       TWO  := ONE + ONE;
  236.       ZERO := ONE - ONE;
  237. {-----------------------------------------------------------------------
  238. C  Determine IBETA, BETA ala Malcolm.
  239. C-----------------------------------------------------------------------}
  240.       A := ONE;
  241.    10: A := A + A;
  242.          TEMP  := A+ONE;
  243.          TEMP1 := TEMP-A;
  244.          IF (TEMP1-ONE = ZERO) THEN
  245.             GOTO 10;
  246.       B := ONE;
  247.    20: B := B + B;
  248.          TEMP := A+B;
  249.          ITEMP := TRUNC (TEMP-A);
  250.          IF (ITEMP = 0) THEN
  251.             GOTO 20;
  252.       IBETA := ITEMP;
  253.       BETA  := CONV[IBETA];
  254. {-----------------------------------------------------------------------
  255. C  Determine IT, IRND.
  256. C-----------------------------------------------------------------------}
  257.       IT := 0;
  258.       B  := ONE;
  259.   100: IT := IT + 1;
  260.          B := B * BETA;
  261.          TEMP := B+ONE;
  262.          TEMP1 := TEMP-B;
  263.          IF (TEMP1-ONE = ZERO) THEN
  264.             GOTO 100;
  265.       IRND := 0;
  266.       BETAH := BETA / TWO;
  267.       TEMP := A+BETAH;
  268.       IF (TEMP-A <> ZERO) THEN
  269.          IRND := 1;
  270.       TEMPA := A + BETA;
  271.       TEMP := TEMPA+BETAH;
  272.       IF ((IRND = 0) AND (TEMP-TEMPA <> ZERO))
  273.          THEN IRND := 2;
  274. {-----------------------------------------------------------------------
  275. C  Determine NEGEP, EPSNEG.
  276. C-----------------------------------------------------------------------}
  277.       NEGEP := IT + 3;
  278.       BETAIN := ONE / BETA;
  279.       A := ONE;
  280.       FOR I := 1 TO NEGEP DO BEGIN
  281.          A := A * BETAIN;
  282.       END;
  283.       B := A;
  284.   210:TEMP := ONE-A;
  285.          IF (TEMP-ONE <> ZERO) THEN
  286.            GOTO 220;
  287.          A := A * BETA;
  288.          NEGEP := NEGEP - 1;
  289.       GOTO 210;
  290.   220:NEGEP := -NEGEP;
  291.       EPSNEG:= A;
  292. {-----------------------------------------------------------------------
  293. C  Determine MACHEP, EPS.
  294. C-----------------------------------------------------------------------}
  295.  
  296.        MACHEP := -IT - 3;
  297.        A := B;
  298.   300: TEMP := ONE+A;
  299.          IF (TEMP-ONE <> ZERO) THEN
  300.             GOTO 320;
  301.          A := A * BETA;
  302.          MACHEP := MACHEP + 1;
  303.       GOTO 300;
  304.   320: EPS := A;
  305. {-----------------------------------------------------------------------
  306. C  Determine NGRD.
  307. C-----------------------------------------------------------------------}
  308.       NGRD := 0;
  309.       TEMP := ONE+EPS;
  310.       IF ((IRND = 0) AND (TEMP*ONE-ONE <> ZERO)) THEN
  311.          NGRD := 1;
  312. {-----------------------------------------------------------------------
  313. C  Determine IEXP, MINEXP, XMIN.
  314. C
  315. C  Loop to determine largest I and K = 2**I such that
  316. C         (1/BETA) ** (2**(I))
  317. C  does not underflow.
  318. C  Exit from loop is signaled by an underflow.
  319. C-----------------------------------------------------------------------}
  320.       I := 0;
  321.       K := 1;
  322.       Z := BETAIN;
  323.       T := ONE + EPS;
  324.       NXRES := 0;
  325.   400:   Y := Z;
  326.          Z := Y * Y;
  327. {-----------------------------------------------------------------------
  328. C  Check for underflow here.
  329. C-----------------------------------------------------------------------}
  330.  
  331.          A := Z * ONE;
  332.          TEMP := Z * T;
  333.          IF ((A+A = ZERO) OR (ABS(Z) >= Y)) THEN
  334.             GOTO 410;
  335.          TEMP1 := TEMP * BETAIN;
  336.          IF (TEMP1*BETA = Z) THEN
  337.             GOTO 410;
  338.          I := I + 1;
  339.          K := K + K;
  340.       GOTO 400;
  341.   410: IF (IBETA = 10) THEN
  342.          GOTO 420;
  343.       IEXP := I + 1;
  344.       MX := K + K;
  345.       GOTO 450;
  346. {-----------------------------------------------------------------------
  347. C  This segment is for decimal machines only.
  348. C-----------------------------------------------------------------------}
  349.   420: IEXP := 2;
  350.       IZ := IBETA;
  351.   430: IF (K < IZ) THEN
  352.          GOTO 440;
  353.          IZ := IZ * IBETA;
  354.          IEXP := IEXP + 1;
  355.       GOTO 430;
  356.   440: MX := IZ + IZ - 1;
  357.  
  358. {-----------------------------------------------------------------------
  359. C  Loop to determine MINEXP, XMIN.
  360. C  Exit from loop is signaled by an underflow.
  361. C-----------------------------------------------------------------------}
  362.  
  363.   450:   XMIN := Y;
  364.          Y := Y * BETAIN;
  365. {-----------------------------------------------------------------------
  366. C  Check for underflow here.
  367. C-----------------------------------------------------------------------}
  368.          A    := Y * ONE;
  369.          TEMP := Y * T;
  370.          IF (((A+A) = ZERO) OR (ABS(Y) >= XMIN)) THEN
  371.              GOTO 460;
  372.          K := K + 1;
  373.          TEMP1 := TEMP * BETAIN;
  374.          IF ((TEMP1*BETA <> Y) OR (TEMP = Y)) THEN
  375.                GOTO 450
  376.             ELSE BEGIN
  377.                NXRES := 3;
  378.                XMIN  := Y;
  379.             END;
  380.   460: MINEXP := -K;
  381.  
  382. {-----------------------------------------------------------------------
  383. C  Determine MAXEXP, XMAX.
  384. C-----------------------------------------------------------------------}
  385.       IF ((MX > K+K-3) OR (IBETA = 10)) THEN
  386.          GOTO 500;
  387.       MX := MX + MX;
  388.       IEXP := IEXP + 1;
  389. 500:  MAXEXP := MX + MINEXP;
  390.  
  391. {-----------------------------------------------------------------
  392. C  Adjust IRND to reflect partial underflow.
  393. C-----------------------------------------------------------------}
  394.  
  395.       IRND := IRND + NXRES;
  396.  
  397. {-----------------------------------------------------------------
  398. C  Adjust for IEEE-style machines.
  399. C-----------------------------------------------------------------}
  400.  
  401.       IF (IRND >= 2) THEN
  402.          MAXEXP := MAXEXP - 2;
  403.  
  404. {-----------------------------------------------------------------
  405. C  Adjust for machines with implicit leading bit in binary
  406. C  significand, and machines with radix point at extreme
  407. C  right of significand.
  408. C-----------------------------------------------------------------}
  409.       I := MAXEXP + MINEXP;
  410.       IF ((IBETA = 2) AND (I = 0)) THEN
  411.          MAXEXP := MAXEXP - 1;
  412.       IF (I > 20) THEN
  413.          MAXEXP := MAXEXP - 1;
  414.       IF (A <> Y) THEN
  415.          MAXEXP := MAXEXP - 2;
  416.       XMAX := ONE - EPSNEG;
  417.       IF (XMAX*ONE <> XMAX) THEN
  418.          XMAX := ONE - BETA * EPSNEG;
  419.       XMAX := XMAX / (BETA * BETA * BETA * XMIN);
  420.       I := MAXEXP + MINEXP + 3;
  421.       IF (I <= 0) THEN
  422.          GOTO 520;
  423.       FOR L := 1 TO I DO BEGIN
  424.           IF (IBETA = 2) THEN XMAX := XMAX + XMAX;
  425.           IF (IBETA <>2) THEN XMAX := XMAX * BETA;
  426.       END;
  427.   520:
  428.   END;
  429.  
  430.  
  431. VAR
  432.      I,IBETA,IEXP,IOUT,IRND,II,IT,I1,J,K1,K2,K3,MACHEP,
  433.      MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
  434.      A,AIT,ALBETA,B,BETA,BETAP,DEL,EM,EPS,EPSNEG,EXPON,HALF,
  435.      OB32,ONE,R6,R7,SUM,TWO,W,X,XL,XMAX,XMIN,XN,XSQ,X1,Y,T,
  436.      Z,ZERO,ZZ: REAL;
  437.  
  438. LABEL 100,105,110,120,200;
  439.  
  440. BEGIN
  441.      MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  442.             MAXEXP,EPS,EPSNEG,XMIN,XMAX);
  443.  
  444.    WriteLn ('MACHAR DETERMINED THE FOLLOWING PARAMETERS OF THE FLOATING-POINT ARITHEMTIC');
  445.    WriteLn;
  446.    WriteLn ('Radix for floating-point representation:        ', IBETA);
  447.    WriteLn ('Number of base ', IBETA:2, ' digits in significand:        ', IT);
  448.    IF IBETA <> 10 THEN
  449.       WriteLn ('Number of bits used to represent the exponent:  ', IEXP)
  450.    ELSE
  451.       WriteLn ('Number of decimal places used for the exponent:  ', IEXP);
  452.    WriteLn ('Smallest positive eps such that 1+eps <> 1 :   ', EPS:14, ' = ', IBETA, ' ** -', ABS(MACHEP));
  453.    WriteLn ('Smallest positive eps such that 1-eps <> 1 :   ', EPSNEG:14, ' = ', IBETA, ' ** -', ABS(NEGEP));
  454.    WriteLn ('Smallest pos. normalized floating-point number:', XMIN:14, ' = ', IBETA, ' ** -', ABS(MINEXP)) ;
  455.    WriteLn ('Largest floating-point number:                 ', XMAX:14, ' = ', IBETA, ' ** +', ABS(MAXEXP), ' - 1');
  456.    WriteLn;
  457.    Write   ('Floating-point arithmetic ');
  458.    CASE IRND MOD 3 OF
  459.        0: WriteLn ('chops');
  460.        1: WriteLn ('rounds, but not IEEE style,');
  461.        2: WriteLn ('rounds IEEE style');
  462.    END;
  463.    IF IRND DIV 3 = 1 THEN
  464.       WriteLn ('and supports gradual underflow')
  465.    ELSE
  466.       WriteLn ('and does not support gradual underflow');
  467.  
  468.       BETA  := IBETA;
  469.       ALBETA:= LN(BETA);
  470.       AIT := IT;
  471.       ONE := 1.0;
  472.       HALF:= 0.5;
  473.       TWO := 2.0;
  474.       ZERO:= 0.0;
  475.       A   := -0.0625;
  476.       B   := -A;
  477.       OB32:= B * HALF;
  478.       N := 10000;
  479.       XN:= N;
  480.       I1:= 0;
  481. {*********************************************************************}
  482. {*{----------------------------------------------------------------- *}
  483. {*C     RANDOM ARGUMENT ACCURACY TESTS                               *}
  484. {*C----------------------------------------------------------------- *}
  485. {*********************************************************************}
  486.  
  487.       FOR J := 1 TO 4 DO BEGIN
  488.          K1 := 0;
  489.          K3 := 0;
  490.          X1 := ZERO;
  491.          R6 := ZERO;
  492.          R7 := ZERO;
  493.          DEL:= (B - A) / XN;
  494.          XL := A;
  495.  
  496.          FOR I := 1 TO N DO BEGIN
  497.             X := DEL * REN(I1) + XL;
  498.             IF (J = 2) THEN
  499.                X := ((1.0+X*A)-ONE)*16.0;
  500.             Z := ARCTAN(X);
  501.             IF (J <> 1) THEN
  502.                GOTO 100;
  503.             XSQ := X * X;
  504.             EM  := 17.0;
  505.             SUM := XSQ / EM;
  506.  
  507.             FOR II := 1 TO 7 DO BEGIN
  508.                EM := EM - TWO;
  509.                SUM:= (ONE/EM - SUM) * XSQ;
  510.             END;
  511.  
  512.             SUM := -X * SUM;
  513.             ZZ  :=  X + SUM;
  514.             SUM := (X - ZZ) + SUM;
  515.             IF (IRND = 0) THEN
  516.                ZZ := ZZ + (SUM + SUM);
  517.             GOTO 110;
  518.   100:      IF (J <> 2) THEN
  519.                GOTO 105;
  520.             Y := (X - 0.0625) / (ONE + X * A);
  521.             ZZ:= (ARCTAN(Y));
  522.             ZZ:= ZZ - 8.1190004042651526021e-5 + OB32;
  523.             ZZ:= ZZ + OB32;
  524.             GOTO 110;
  525.   105:      Z := Z + Z;
  526.             Y := X / ((HALF + X * HALF)*((HALF - X) + HALF));
  527.             ZZ:= ARCTAN(Y);
  528.   110:      IF (Z <> ZERO) THEN
  529.                W := (Z - ZZ) / Z
  530.             ELSE IF ZZ <> ZERO THEN
  531.                W := ONE;
  532.             IF (W > ZERO) THEN
  533.                K1 := K1 + 1;
  534.             IF (W < ZERO) THEN
  535.                K3 := K3 + 1;
  536.             W := ABS(W);
  537.             IF (W <= R6) THEN
  538.                GOTO 120;
  539.             R6 := W;
  540.             X1 := X;
  541.   120:      R7 := R7 + W * W;
  542.             XL := XL + DEL;
  543.   200:
  544.          END;
  545.  
  546.          K2 := N - K3 - K1;
  547.          R7 := SQRT(R7/XN);
  548.  
  549.          IF (J = 1) THEN BEGIN
  550.             WRITELN;
  551.             WRITELN ;
  552.             WRITELN ('TEST OF ARCTAN(X) VS TRUNCATED TAYLOR SERIES');
  553.             WRITELN;
  554.             END;
  555.          IF (J = 2) THEN BEGIN
  556.             WRITELN;
  557.             WRITELN;
  558.             WRITELN ('TEST OF ARCTAN(X) VS ARCTAN(1/16) + ARCTAN((X-1/16)/(1+X/16))');
  559.             WRITELN;
  560.             END;
  561.          IF (J > 2) THEN BEGIN
  562.             WRITELN;
  563.             WRITELN;
  564.             WRITELN ('TEST OF 2*ARCTAN(X) VS ARCTAN(2X((1-X*X))');
  565.             WRITELN;
  566.             END;
  567.  
  568.          WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  569.          WRITELN ('(', A, ',', B, ')');
  570.          WRITELN;
  571.  
  572.          WRITELN ('ARCTAN (X) WAS LARGER', K1:6, ' TIMES');
  573.          WRITELN ('               AGREED', K2:6, ' TIMES');
  574.          WRITELN ('      AND WAS SMALLER', K3:6, ' TIMES');
  575.  
  576.          WRITELN;
  577.          WRITELN ('THERE ARE ', IT, ' BASE ', IBETA,
  578.                   ' SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER');
  579.          WRITELN;
  580.  
  581.          W := -999;
  582.          IF (R6 <> ZERO) THEN
  583.             W := LN (ABS(R6))/ALBETA;
  584.          WRITELN ('THE MAXIMUM RELATIVE ERROR OF          ', R6:12,
  585.                   ' = ', IBETA, ' **', W:7:2);
  586.          WRITELN ('OCCURED FOR X = ', X1);
  587.          W := MAX1 (AIT+W,ZERO);
  588.          WRITELN;
  589.          WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  590.                   ' SIGNIFICANT DIGITS IS        ', W:7:2);
  591.          W := -999.0;
  592.          IF (R7 <> ZERO) THEN
  593.             W := LN (ABS(R7))/ALBETA;
  594.          WRITELN;
  595.          WRITELN ('THE ROOT MEAN SQUARE RELATIVE ERROR WAS', R7:12,
  596.                   ' = ', IBETA, ' **' , W:7:2);
  597.          W := MAX1 (AIT+W,ZERO);
  598.          WRITELN;
  599.          WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  600.                   ' SIGNIFICANT DIGITS IS        ', W:7:2);
  601.          A := B;
  602.          IF (J = 1) THEN
  603.             B := TWO - SQRT(3.0);
  604.          IF (J = 2) THEN
  605.             B := SQRT(TWO) - ONE;
  606.          IF (J = 3) THEN
  607.             B := ONE;
  608.      END;
  609.  
  610. {-----------------------------------------------------------------
  611. C     SPECIAL TESTS
  612. C-----------------------------------------------------------------}
  613.  
  614.       WRITELN;
  615.       WRITELN;
  616.       WRITELN ('SPECIAL TESTS');
  617.       WRITELN;
  618.       WRITELN ('THE IDENTITY  ARCTAN(-X) = -ARCTAN(X)  WILL BE TESTED');
  619.       WRITELN;
  620.       WRITELN ('          X           F(X) + F(-X)');
  621.       WRITELN;
  622.  
  623.       A := 5.0;
  624.  
  625.       FOR I := 1 TO 5 DO BEGIN
  626.          X := REN(I1) * A;
  627.          Z := ARCTAN(X) + ARCTAN(-X);
  628.          WRITELN (X:18, Z:18);
  629.       END;
  630.  
  631.       WRITELN;
  632.       WRITELN;
  633.       WRITELN ('THE IDENTITY ARCTAN(X) = X , X SMALL, WILL BE TESTED.');
  634.       WRITELN;
  635.       WRITELN ('          X                 X-F(X)');
  636.       WRITELN;
  637.       BETAP := POWER (BETA,IT);
  638.       X := REN(I1) / BETAP;
  639.  
  640.       FOR I := 1 TO 5 DO BEGIN
  641.          Z := X - ARCTAN (X);
  642.          WRITELN (X:18, Z:18);
  643.          X := X / BETA;
  644.       END;
  645.  
  646.  
  647.       WRITELN;
  648.       WRITELN;
  649.       WRITELN ('TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT');
  650.       WRITELN;
  651.       EXPON := MINEXP * 0.75;
  652.       X := POWER( BETA,EXPON);
  653.       Y := ARCTAN(X);
  654.       WRITELN ('ARCTAN (', X:18, ') = ', Y:18);
  655.  
  656. {C-----------------------------------------------------------------
  657. C     TEST OF ERROR RETURNS
  658. C-----------------------------------------------------------------}
  659.  
  660.       WRITELN;
  661.       WRITELN;
  662.       WRITELN ('TEST OF ERROR RETURNS');
  663.       WRITELN;
  664.       WRITELN ('ARCTAN WILL BE CALLED WITH THE ARGUMENT ', XMAX:18);
  665.       WRITELN ('THIS SHOULD NOT TRIGGER AN ERROR MESSAGE');
  666.       Z := ARCTAN(XMAX);
  667.       WRITELN ('ARCTAN (', XMAX:18, ') = ', Z:18);
  668.       WRITELN;
  669.       WRITELN ('THIS CONCLUDES THE TESTS');
  670. END.
  671.