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