home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / BPL70N12.ZIP / TESTPRGS.ZIP / DPOWER.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1993-02-14  |  9.7 KB  |  339 lines
  1. PROGRAM DPower;    { converted from Fortran original 05-01-92 Norbert Juffa }
  2.  
  3. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  4.  
  5. USES MachArit, Power;
  6.  
  7. {
  8. C     PROGRAM TO TEST POWER FUNCTION (**)
  9. C
  10. C     DATA REQUIRED
  11. C
  12. C        NONE
  13. C
  14. C     SUBPROGRAMS REQUIRED FROM THIS PACKAGE
  15. C
  16. C        MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
  17. C                 INFORMATION ON THE FLOATING-POINT ARITHMETIC
  18. C                 SYSTEM.  NOTE THAT THE CALL TO MACHAR CAN
  19. C                 BE DELETED PROVIDED THE FOLLOWING SIX
  20. C                 PARAMETERS ARE ASSIGNED THE VALUES INDICATED
  21. C
  22. C                 IBETA  - THE RADIX OF THE FLOATING-POINT SYSTEM
  23. C                 IT     - THE NUMBER OF BASE-IBETA DIGITS IN THE
  24. C                          SIGNIFICAND OF A FLOATING-POINT NUMBER
  25. C                 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE
  26. C                          INTEGER SUCH THAT  DFLOAT(IBETA)**MINEXP
  27. C                          IS A POSITIVE FLOATING-POINT NUMBER
  28. C                 MAXEXP - THE LARGEST POSITIVE INTEGER EXPONENT
  29. C                          FOR A FINITE FLOATING-POINT NUMBER
  30. C                 XMIN   - THE SMALLEST NON-VANISHING FLOATING-POINT
  31. C                          POWER OF THE RADIX
  32. C                 XMAX   - THE LARGEST FINITE FLOATING-POINT
  33. C                          NUMBER
  34. C
  35. C        REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
  36. C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
  37. C
  38. C
  39. C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
  40. C
  41. C         DABS, DLOG, DMAX1, DEXP, DFLOAT, DSQRT
  42. C
  43. C
  44. C     LATEST REVISION - DECEMBER 6, 1979
  45. C
  46. C     AUTHOR - W. J. CODY
  47. C              ARGONNE NATIONAL LABORATORY
  48. C
  49. C
  50. }
  51.  
  52.  
  53. FUNCTION REN (K: LONGINT): REAL;
  54.  
  55. {
  56.       DOUBLE PRECISION FUNCTION REN(K)
  57. C
  58. C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
  59. C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
  60. C      VOL. 8, NO. 10, OCTOBER 1965.
  61. C
  62. C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
  63. C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
  64. C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
  65. C      29 BITS.
  66. C
  67. }
  68.  
  69. VAR   J:  LONGINT;
  70. CONST IY: LONGINT = 100001;
  71.  
  72. BEGIN
  73.    J  := K;
  74.    IY := IY * 125;
  75.    IY := IY - (IY DIV 2796203) * 2796203;
  76.    REN:= 1.0 * (IY) / 2796203.0e0 * (1.0e0 + 1.0e-6 + 1.0e-12);
  77. END;
  78.  
  79.  
  80. FUNCTION MAX1 (A, B:REAL): REAL;
  81. BEGIN
  82.    IF A > B THEN
  83.       MAX1 := A
  84.    ELSE
  85.       MAX1 := B;
  86. END;
  87.  
  88.  
  89.  
  90. VAR   I,IBETA,IEXP,IOUT,IRND,
  91.       IT,I1,J,K1,K2,K3,MACHEP,
  92.       MAXEXP,MINEXP,N,NEGEP,NGRD:  LONGINT;
  93.       A,AIT,ALBETA,ALXMAX,B,BETA,
  94.       C,DEL,DELY,EPS,EPSNEG,ONE,
  95.       ONEP5,R6,R7,SCALE,TWO,W,
  96.       X,XL,XMAX,XMIN,XN,XSQ,X1,
  97.       Y,Y1,Y2,Z,ZERO,TEN,THREE,
  98.       ZZ,ONEHUNDREDTH:             REAL;
  99.  
  100. LABEL 50,70,110,120,210,215,220,300;
  101.  
  102. BEGIN
  103.  
  104.    N := 1000000;  { number of trials }
  105.  
  106.    MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  107.            EPS,EPSNEG,XMIN,XMAX);
  108.    PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  109.                EPS,EPSNEG,XMIN,XMAX);
  110.    BETA        := IBETA;
  111.    ALBETA      := LN (BETA);
  112.    AIT         := IT;
  113.    ALXMAX      := LN (XMAX);
  114.    ZERO        := 0;
  115.    ONE         := 1;
  116.    TWO         := 2;
  117.    THREE       := 3;
  118.    TEN         := 10;
  119.    ONEHUNDREDTH:= 0.01;
  120.    ONEP5       := (TWO + ONE) / TWO;
  121.    SCALE       := ONE;
  122.    J           := (IT+1) DIV 2;
  123.  
  124.    FOR I := 1 TO J DO BEGIN
  125.       SCALE := SCALE * BETA;
  126.    END;
  127.  
  128.    A    := ONE / BETA;
  129.    B    := ONE;
  130.    C    := -MAX1 (ALXMAX, -LN(XMIN))/ LN(100.0);
  131.    DELY := -C - C;
  132.    XN   := N;
  133.    I1   := 0;
  134.    Y1   := ZERO;
  135.  
  136. {-----------------------------------------------------------------
  137. C     RANDOM ARGUMENT ACCURACY TESTS
  138. C-----------------------------------------------------------------}
  139.  
  140.    FOR J := 1 TO 4 DO BEGIN
  141.       K1 := 0;
  142.       K3 := 0;
  143.       X1 := ZERO;
  144.       R6 := ZERO;
  145.       R7 := ZERO;
  146.       DEL:= (B - A) / XN;
  147.       XL := A;
  148.  
  149.       FOR I := 1 TO N DO BEGIN
  150.          X := DEL * REN(I1) + XL;
  151.          IF (J <> 1) THEN
  152.             GOTO 50;
  153.          ZZ := POW (X, ONE);
  154.          Z  := X;
  155.          GOTO 110;
  156. 50:      W  := SCALE * X;
  157.          X  := (X + W);
  158.          X  := X - W;
  159.          XSQ:= X * X;
  160.          IF (J = 4) THEN
  161.             GOTO 70;
  162.          ZZ := POW (XSQ, ONEP5);
  163.          Z  := X * XSQ;
  164.          GOTO 110;
  165. 70:      Y  := DELY * REN(I1) + C;
  166.          Y2 := Y / TWO;
  167.          Y2 := Y2 + Y;
  168.          Y2 := Y2 - Y;
  169.          Y  := Y2 + Y2;
  170.          Z  := Pow (X, Y);
  171.          ZZ := Pow (XSQ, Y2);
  172. 110:     IF Z <> ZERO THEN
  173.             W := (Z - ZZ) / Z
  174.          ELSE IF ZZ <> ZERO THEN
  175.             W := ONE;
  176.          IF W > ZERO THEN
  177.             K1 := K1 + 1;
  178.          IF W < ZERO THEN
  179.             K3 := K3 + 1;
  180.          W := ABS (W);
  181.          IF W <= R6 THEN
  182.             GOTO 120;
  183.          R6 := W;
  184.          X1 := X;
  185.          IF J = 4 THEN
  186.             Y1 := Y;
  187. 120:     R7 := R7 + W * W;
  188.          XL := XL + DEL;
  189.       END;
  190.  
  191.       K2 := N - K3 - K1;
  192.       R7 := SQRT (R7/XN);
  193.  
  194.       IF J > 1 THEN
  195.          GOTO 210;
  196.       WRITELN;
  197.       WRITELN;
  198.       WRITELN ('TEST OF X**1.0 VS X');
  199.       WRITELN;
  200.       WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  201.       WRITELN ('(', A, ',', B, ')');
  202.       WRITELN;
  203.       WRITELN ('X**1.0 WAS LARGER', K1:6, ' TIMES');
  204.       WRITELN ('           AGREED', K2:6, ' TIMES');
  205.       WRITELN ('  AND WAS SMALLER', K3:6, ' TIMES');
  206.       GOTO 220;
  207. 210:  IF J = 4 THEN
  208.          GOTO 215;
  209.       WRITELN;
  210.       WRITELN;
  211.       WRITELN ('TEST OF XSQ**1.5 VS XSQ*X');
  212.       WRITELN;
  213.       WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  214.       WRITELN ('(', A, ',', B, ')');
  215.       WRITELN;
  216.       WRITELN ('X**1.5 WAS LARGER', K1:6, ' TIMES');
  217.       WRITELN ('           AGREED', K2:6, ' TIMES');
  218.       WRITELN ('  AND WAS SMALLER', K3:6, ' TIMES');
  219.       GOTO 220;
  220. 215:  WRITELN;
  221.       WRITELN;
  222.       WRITELN ('TEST OF X**Y VS XSQ**(Y/2)');
  223.       W := C + DELY;
  224.       WRITELN;
  225.       WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE REGION');
  226.       WRITELN ('X IN (', A, ',', B, '),');
  227.       WRITELN ('Y IN (', C, ',', W, ')');
  228.       WRITELN;
  229.       WRITELN (' X**Y  WAS LARGER', K1:6, ' TIMES');
  230.       WRITELN ('           AGREED', K2:6, ' TIMES');
  231.       WRITELN ('  AND WAS SMALLER', K3:6, ' TIMES');
  232. 220:  WRITELN;
  233.       WRITELN ('THERE ARE ', IT, ' BASE ', IBETA,
  234.               ' SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER');
  235.       WRITELN;
  236.       W := -999;
  237.       IF R6 <> ZERO THEN
  238.          W := LN (ABS(R6))/ALBETA;
  239.       WRITELN ('THE MAXIMUM RELATIVE ERROR OF          ', R6:12,
  240.                ' = ', IBETA, ' **', W:7:2);
  241.       WRITELN ('OCCURED FOR X = ', X1);
  242.       W := MAX1 (AIT+W,ZERO);
  243.       WRITELN;
  244.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  245.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  246.       W := -999.0;
  247.       IF (R7 <> ZERO) THEN
  248.          W := LN (ABS(R7))/ALBETA;
  249.       WRITELN;
  250.       WRITELN ('THE ROOT MEAN SQUARE RELATIVE ERROR WAS', R7:12,
  251.                ' = ', IBETA, ' **' , W:7:2);
  252.       W := MAX1 (AIT+W,ZERO);
  253.       WRITELN;
  254.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  255.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  256.       IF J = 1 THEN
  257.          GOTO 300;
  258.       B := TEN;
  259.       A := ONEHUNDREDTH;
  260.       IF J = 3 THEN
  261.          GOTO 300;
  262.       A := ONE;
  263.       B := EXP (ALXMAX/THREE);
  264. 300:
  265.    END;
  266.  
  267. {-----------------------------------------------------------------}
  268. {     SPECIAL TESTS                                               }
  269. {-----------------------------------------------------------------}
  270.  
  271.    WRITELN;
  272.    WRITELN;
  273.    WRITELN ('SPECIAL TESTS');
  274.    WRITELN;
  275.    WRITELN ('THE IDENTITY  X ** Y = (1/X) ** (-Y)  WILL BE TESTED.');
  276.    WRITELN;
  277.    WRITELN ('         X                Y           X**Y-(1/X)**(-Y)/X**Y ');
  278.    B := TEN;
  279.  
  280.    FOR I := 1 TO 5 DO BEGIN
  281.       X := REN(I1) * B + ONE;
  282.       Y := REN(I1) * B + ONE;
  283.       Z := POW (X, Y);
  284.       ZZ:= POW ((ONE/X), (-Y));
  285.       W := (Z - ZZ) / Z;
  286.       WRITELN (X:18, Y:18, W:18);
  287.    END;
  288.  
  289. {-----------------------------------------------------------------}
  290. {     TEST OF ERROR RETURNS                                       }
  291. {-----------------------------------------------------------------}
  292.  
  293.    WRITELN;
  294.    WRITELN;
  295.    WRITELN ('TEST OF ERROR RETURNS');
  296.    X := BETA;
  297.    Y := MINEXP;
  298.    WRITELN;
  299.    WRITELN ('(', X, ') ** (', Y, ') WILL BE COMPUTED.');
  300.    WRITELN ('THIS SHOULD NOT TRIGGER AN ERROR MESSAGE.');
  301.    Z := POW (X, Y);
  302.    WRITELN ('THE VALUE RETURNED IS: ', Z);
  303.    Y := (MAXEXP-1);
  304.    WRITELN;
  305.    WRITELN ('(', X, ') ** (', Y, ') WILL BE COMPUTED.');
  306.    WRITELN ('THIS SHOULD NOT TRIGGER AN ERROR MESSAGE.');
  307.    Z := POW (X, Y);
  308.    WRITELN ('THE VALUE RETURNED IS: ', Z);
  309.    X := ZERO;
  310.    Y := TWO;
  311.    WRITELN;
  312.    WRITELN ('(', X, ') ** (', Y, ') WILL BE COMPUTED.');
  313.    WRITELN ('THIS SHOULD NOT TRIGGER AN ERROR MESSAGE.');
  314.    Z := POW (X, Y);
  315.    WRITELN ('THE VALUE RETURNED IS: ', Z);
  316.    X := -Y;
  317.    Y := ZERO;
  318.    WRITELN;
  319.    WRITELN ('(', X, ') ** (', Y, ') WILL BE COMPUTED.');
  320.    WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE.');
  321.    Z := POW (X, Y);
  322.    WRITELN ('THE VALUE RETURNED IS: ', Z);
  323.    Y := TWO;
  324.    WRITELN;
  325.    WRITELN ('(', X, ') ** (', Y, ') WILL BE COMPUTED.');
  326.    WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE.');
  327.    Z := POW (X, Y);
  328.    WRITELN ('THE VALUE RETURNED IS: ', Z);
  329.    X := ZERO;
  330.    Y := ZERO;
  331.    WRITELN;
  332.    WRITELN ('(', X, ') ** (', Y, ') WILL BE COMPUTED.');
  333.    WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE.');
  334.    Z := POW (X, Y);
  335.    WRITELN ('THE VALUE RETURNED IS: ', Z);
  336.    WRITELN;
  337.    WRITELN ('THIS CONCLUDES THE TESTS');
  338. END. { DPower }
  339.