home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / TPL60N11.ZIP / TESTPRGS.ZIP / DSIN.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1992-03-28  |  22.9 KB  |  675 lines
  1. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  2.  
  3. (*
  4. C     PROGRAM TO TEST DSIN/DCOS
  5. C
  6. C     DATA REQUIRED
  7. C
  8. C        NONE
  9. C
  10. C     SUBPROGRAMS REQUIRED FROM THIS PACKAGE
  11. C
  12. C        MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
  13. C                 INFORMATION ON THE FLOATING-POINT ARITHMETIC
  14. C                 SYSTEM.  NOTE THAT THE CALL TO MACHAR CAN
  15. C                 BE DELETED PROVIDED THE FOLLOWING FIVE
  16. C                 PARAMETERS ARE ASSIGNED THE VALUES INDICATED
  17. C
  18. C                 IBETA  - THE RADIX OF THE FLOATING-POINT SYSTEM
  19. C                 IT     - THE NUMBER OF BASE-IBETA DIGITS IN THE
  20. C                          SIGNIFICAND OF A FLOATING-POINT NUMBER
  21. C                 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE
  22. C                          INTEGER SUCH THAT  DFLOAT(IBETA)**MINEXP
  23. C                          IS A POSITIVE FLOATING-POINT NUMBER
  24. C                 EPS    - THE SMALLEST POSITIVE FLOATING-POINT
  25. C                          NUMBER SUCH THAT 1.0+EPS .NE. 1.0
  26. C                 EPSNEG - THE SMALLEST POSITIVE FLOATING-POINT
  27. C                          NUMBER SUCH THAT 1.0-EPSNEG .NE. 1.0
  28. C
  29. C        REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
  30. C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
  31. C
  32. C
  33. C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
  34. C
  35. C         DABS, DLOG, DMAX1, DCOS, DLOAT, DSIN, DSQRT
  36. C
  37. C
  38. C     LATEST REVISION - DECEMBER 6, 1979
  39. C
  40. C     AUTHOR - W. J. CODY
  41. C              ARGONNE NATIONAL LABORATORY
  42. C
  43. C
  44. *)
  45.  
  46.  
  47.    function Power (X, Y: real): real;
  48.  
  49.   var z:real;
  50.     iy:integer;
  51.     reciproc:boolean;
  52.       (* Crude power function:  see Cody & Waite for a better one. *)
  53.  
  54.       begin (* Power *)
  55.       if Y = 0.0 then
  56.          Power := 1.0
  57.       else if (X = 0.0) and (Y > 0.0) then
  58.          Power := 0.0
  59.       else if (Y = Int (Y)) AND (Abs(Y) <= MAXINT) THEN BEGIN
  60.          IY := TRUNC (Y);
  61.          reciproc := iy < 0;
  62.          iy := abs (iy);
  63.          z := 1.0;
  64.          while iy > 0 do begin
  65.            if odd(iy) then begin
  66.              z := z*x;
  67.  
  68.            end else begin
  69.            end; {endif}
  70.            iy := iy div 2;
  71.            if iy >0 then begin
  72.               x := sqr (x);
  73.            end else begin
  74.            end; {endif}
  75.  
  76.          end; {endwhile}
  77.          if reciproc then begin
  78.            power := 1 / z;
  79.  
  80.          end else begin
  81.            power := z;
  82.          end {endif}
  83.          end
  84.       else if (X < 0.0) and (Y = Int(Y)) then
  85.          if odd(trunc(Y)) then Power := -exp (Y * ln (-X))
  86.          else Power := exp (Y * ln (-X))
  87.       else
  88.          Power := exp (Y * ln (X));
  89.       end (* Power *);
  90.  
  91.  
  92. FUNCTION REN (VAR K: LONGINT): REAL;
  93.  
  94. (*
  95.       DOUBLE PRECISION FUNCTION REN(K)
  96. C
  97. C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
  98. C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
  99. C      VOL. 8, NO. 10, OCTOBER 1965.
  100. C
  101. C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
  102. C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
  103. C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
  104. C      29 BITS.
  105. C
  106. *)
  107.  
  108. VAR   J: LONGINT;
  109. CONST IY: LONGINT = 100001;
  110.  
  111. BEGIN
  112.       J  := K;
  113.       IY := IY * 125;
  114.       IY := IY - (IY DIV 2796203) * 2796203;
  115.       REN:= 1.0 * (IY) / 2796203.0e0 * (1.0e0 + 1.0e-6 + 1.0e-12);
  116. END;
  117.  
  118.  
  119. FUNCTION LOG (X: REAL): REAL;
  120. BEGIN
  121.    LOG  := LN (X) * 0.43429448190325182765;
  122. END;
  123.  
  124.  
  125.  
  126. FUNCTION MAX1 (A, B:REAL): REAL;
  127. BEGIN
  128.    IF A > B THEN
  129.       MAX1 := A
  130.    ELSE
  131.       MAX1 := B;
  132. END;
  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.  
  432. VAR
  433.            I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,
  434.            MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
  435.  
  436.            A,AIT,ALBETA,B,BETA,C,DEL,EPS,EPSNEG,HALF,ONE,T,THREE,BETAP,
  437.            TEMP,EXPON,R6,R7,TENTH,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ: REAL;
  438.  
  439. LABEL      100, 110, 120, 150, 160, 210, 220, 230, 240, 300;
  440.  
  441. BEGIN
  442.       MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  443.              MAXEXP,EPS,EPSNEG,XMIN,XMAX);
  444.  
  445.    WriteLn ('MACHAR DETERMINED THE FOLLOWING PARAMETERS OF THE FLOATING-POINT ARITHEMTIC');
  446.    WriteLn;
  447.    WriteLn ('Radix for floating-point representation:        ', IBETA);
  448.    WriteLn ('Number of base ', IBETA:2, ' digits in significand:        ', IT);
  449.    IF IBETA <> 10 THEN
  450.       WriteLn ('Number of bits used to represent the exponent:  ', IEXP)
  451.    ELSE
  452.       WriteLn ('Number of decimal places used for the exponent:  ', IEXP);
  453.    WriteLn ('Smallest positive eps such that 1+eps <> 1 :   ', EPS:14, ' = ', IBETA, ' ** -', ABS(MACHEP));
  454.    WriteLn ('Smallest positive eps such that 1-eps <> 1 :   ', EPSNEG:14, ' = ', IBETA, ' ** -', ABS(NEGEP));
  455.    WriteLn ('Smallest pos. normalized floating-point number:', XMIN:14, ' = ', IBETA, ' ** -', ABS(MINEXP)) ;
  456.    WriteLn ('Largest floating-point number:                 ', XMAX:14, ' = ', IBETA, ' ** +', ABS(MAXEXP), ' - 1');
  457.    WriteLn;
  458.    Write   ('Floating-point arithmetic ');
  459.    CASE IRND MOD 3 OF
  460.        0: WriteLn ('chops');
  461.        1: WriteLn ('rounds, but not IEEE style,');
  462.        2: WriteLn ('rounds IEEE style');
  463.    END;
  464.    IF IRND DIV 3 = 1 THEN
  465.       WriteLn ('and supports gradual underflow')
  466.    ELSE
  467.       WriteLn ('and does not support gradual underflow');
  468.  
  469.       BETA  := IBETA;
  470.       ALBETA:= LN(BETA);
  471.       AIT   := IT;
  472.       ONE   := 1.0;
  473.       ZERO  := 0.0;
  474.       THREE := 3.0;
  475.       A := ZERO;
  476.       B := 1.570796327;
  477.       C := B;
  478.       N := 10000;
  479.       XN:= N;
  480.       I1 := 0;
  481. {-----------------------------------------------------------------
  482. C     RANDOM ARGUMENT ACCURACY TESTS
  483. C-----------------------------------------------------------------}
  484.       FOR J := 1 TO 3 DO BEGIN
  485.          K1 := 0;
  486.          K3 := 0;
  487.          X1 := ZERO;
  488.          R6 := ZERO;
  489.          R7 := ZERO;
  490.          DEL:= (B - A) / XN;
  491.          XL := A;
  492.  
  493.          FOR I := 1 TO N DO BEGIN
  494.             X := DEL * REN(I1) + XL;
  495.             Y := X / THREE;
  496.             Y := (X + Y);
  497.             Y := Y - X;
  498.             X := THREE * Y;
  499.             IF J = 3 THEN GOTO 100;
  500.             Z := SIN(X);
  501.             ZZ:= SIN(Y);
  502.             W := ONE;
  503.             IF Z <> ZERO THEN BEGIN
  504.                TEMP := 4.0*ZZ*ZZ;
  505.                TEMP := THREE-TEMP;
  506.                TEMP := ZZ*TEMP;
  507.                TEMP := Z - TEMP;
  508.                W := TEMP / Z;
  509.                END;
  510.             GOTO 110;
  511.   100:      Z  := COS(X);
  512.             ZZ := COS(Y);
  513.             W  := ONE;
  514.             IF Z <> ZERO THEN BEGIN
  515.                TEMP := 4.0*ZZ*ZZ;
  516.                TEMP := THREE-TEMP;
  517.                TEMP := ZZ*TEMP;
  518.                TEMP := Z + TEMP;
  519.                W := TEMP / Z;
  520.             END;
  521.   110:      IF W > ZERO THEN K1 := K1 + 1;
  522.             IF W < ZERO THEN K3 := K3 + 1;
  523.             W := ABS(W);
  524.             IF W <= R6 THEN GOTO 120;
  525.             R6 := W;
  526.             X1 := X;
  527.   120:      R7 := R7 + W * W;
  528.             XL := XL + DEL;
  529.          END;
  530.  
  531.          K2 := N - K3 - K1;
  532.          R7 := SQRT (R7/XN);
  533.          IF (J = 3) THEN GOTO 210;
  534.          WriteLn;
  535.          WriteLn;
  536.          WriteLn ('TEST OF SIN(X) VS 3*SIN(X/3)-4*SIN(X/3)**3');
  537.          WRITELN;
  538.          WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  539.          WRITELN ('(', A, ',', B, ')');
  540.          WRITELN;
  541.          WRITELN ('SIN (X) WAS LARGER', K1:6, ' TIMES');
  542.          WRITELN ('            AGREED', K2:6, ' TIMES');
  543.          WRITELN ('   AND WAS SMALLER', K3:6, ' TIMES');
  544.          GOTO 220;
  545.   210:   WRITELN;
  546.          WRITELN;
  547.          WRITELN ('TEST OF COS(X) VS 4*COS(X/3)**3-3*COS(X/3)');
  548.          WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  549.          WRITELN ('(', A, ',', B, ')');
  550.          WRITELN;
  551.          WRITELN ('COS (X) WAS LARGER', K1:6, ' TIMES');
  552.          WRITELN ('            AGREED', K2:6, ' TIMES');
  553.          WRITELN ('   AND WAS SMALLER', K3:6, ' TIMES');
  554.   220:   WRITELN;
  555.          WRITELN ('THERE ARE ', IT, ' BASE ', IBETA,
  556.                   ' SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER');
  557.          WRITELN;
  558.          W := -999.0;
  559.          IF (R6 <> ZERO) THEN
  560.             W := LN (ABS(R6))/ALBETA;
  561.          WRITELN ('THE MAXIMUM RELATIVE ERROR OF          ', R6:12,
  562.                   ' = ', IBETA, ' **', W:7:2);
  563.          WRITELN ('OCCURED FOR X = ', X1);
  564.          W := MAX1 (AIT+W,ZERO);
  565.          WRITELN;
  566.          WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  567.                   ' SIGNIFICANT DIGITS IS        ', W:7:2);
  568.          W := -999.0;
  569.          IF (R7 <> ZERO) THEN
  570.             W := LN (ABS(R7))/ALBETA;
  571.          WRITELN;
  572.          WRITELN ('THE ROOT MEAN SQUARE RELATIVE ERROR WAS', R7:12,
  573.                   ' = ', IBETA, ' **' , W:7:2);
  574.          W := MAX1 (AIT+W,ZERO);
  575.          WRITELN;
  576.          WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  577.                   ' SIGNIFICANT DIGITS IS        ', W:7:2);
  578.          A := 18.84955592;
  579.          IF (J = 2) THEN A := B + C;
  580.          B := A + C;
  581.       END;
  582.  
  583. {-----------------------------------------------------------------
  584. C     SPECIAL TESTS
  585. C-----------------------------------------------------------------}
  586.  
  587.       WRITELN;
  588.       WRITELN;
  589.       WRITELN ('SPECIAL TESTS');
  590.       WRITELN;
  591.       C := ONE / POWER (BETA, (IT/2));
  592.       Z := SIN (A+C);
  593.       Z := Z - SIN (A-C);
  594.       TEMP := C+C;
  595.       Z := Z / TEMP;
  596.       WRITELN ('IF ',Z:18,' IS NOT ALMOST 1.0, SIN HAS THE WRONG PERIOD.');
  597.  
  598.       WRITELN;
  599.       WRITELN;
  600.       WRITELN ('THE IDENTITY   SIN(-X) = -SIN(X)   WILL BE TESTED.');
  601.       WRITELN;
  602.       WRITELN ('        F(X)           F(X)+F(-X)');
  603.       WRITELN;
  604.  
  605.       FOR I := 1 TO 5 DO BEGIN
  606.          X := REN(I1) * A;
  607.          Z := SIN(X) + SIN(-X);
  608.          WRITELN (X:18, Z:18);
  609.       END;
  610.  
  611.       WRITELN;
  612.       WRITELN;
  613.       WRITELN ('THE IDENTITY SIN(X) = X , X SMALL, WILL BE TESTED.');
  614.       WRITELN;
  615.       WRITELN ('         X                 X-F(X)');
  616.       WRITELN;
  617.       BETAP := POWER (BETA,IT);
  618.       X := REN(I1) / BETAP;
  619.  
  620.       FOR I := 1 TO 5 DO BEGIN
  621.          Z := X - SIN(X);
  622.          WRITELN (X:18, Z:18);
  623.          X := X / BETA;
  624.       END;
  625.  
  626.       WRITELN;
  627.       WRITELN;
  628.       WRITELN ('THE IDENTITY   COS(-X) = COS(X)   WILL BE TESTED.');
  629.       WRITELN;
  630.       WRITELN ('         X              F(X)-F(-X)');
  631.       WRITELN;
  632.       FOR I := 1 TO 5 DO BEGIN
  633.          X := REN (I1) * A;
  634.          Z := COS (X) - COS (-X);
  635.          WRITELN (X:18, Z:18);
  636.       END;
  637.  
  638.       WRITELN;
  639.       WRITELN;
  640.       WRITELN ('TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT.');
  641.       WRITELN;
  642.       EXPON := MINEXP * 0.75;
  643.       X := POWER (BETA, EXPON);
  644.       Y := SIN(X);
  645.       WRITELN ('SIN (', X:15, ') =                   ', Y:15);
  646.       WRITELN;
  647.       WRITELN ('THE FOLLOWING THREE LINES ILLUSTRATE THE LOSS IN SIGNIFICANCE');
  648.       WRITELN ('FOR LARGE ARGUMENTS.  THE ARGUMENTS ARE CONSECUTIVE.');
  649.       WRITELN;
  650.       Z := SQRT (BETAP);
  651.       X := Z * (ONE - EPSNEG);
  652.       Y := SIN(X);
  653.       WRITELN ('SIN (', X:15, ') =                   ', Y:15);
  654.       Y := SIN(Z);
  655.       WRITELN ('SIN (', X:15, ') =                   ', Y:15);
  656.       X := Z * (ONE + EPS);
  657.       Y := SIN(X);
  658.       WRITELN ('SIN (', X:15, ') =                   ', Y:15);
  659.  
  660. {-----------------------------------------------------------------
  661. C     TEST OF ERROR RETURNS
  662. C-----------------------------------------------------------------}
  663.  
  664.       WRITELN;
  665.       WRITELN;
  666.       WRITELN ('TEST OF ERROR RETURNS');
  667.       WRITELN;
  668.       X := BETAP;
  669.       WRITELN ('SIN WILL BE CALLED WITH THE ARGUMENT ',  X:15);
  670.       WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE');
  671.       Y := SIN(X);
  672.       WRITELN ('SIN RETURNED THE VALUE ', Y:15);
  673.       WRITELN ('THIS CONCLUDES THE TESTS');
  674. END.
  675.