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