home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C!T ROM 2
/
ctrom_ii_b.zip
/
ctrom_ii_b
/
PROGRAM
/
PASCAL
/
TPL60N19
/
TESTPRGS
/
MACHARIT.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1993-02-15
|
14KB
|
368 lines
{$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
UNIT MachArit; { converted from Fortran original 05-01-92 Norbert Juffa }
INTERFACE
PROCEDURE MACHAR (VAR IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
MAXEXP: LONGINT; VAR EPS,EPSNEG,XMIN,XMAX: REAL);
PROCEDURE PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
MAXEXP: LONGINT; EPS,EPSNEG,XMIN,XMAX: REAL);
IMPLEMENTATION
PROCEDURE MACHAR(VAR IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
MAXEXP: LONGINT; VAR EPS,EPSNEG,XMIN,XMAX: REAL);
{-----------------------------------------------------------------------
C This Fortran 77 subroutine is intended to determine the parameters
C of the floating-point arithmetic system specified below. The
C determination of the first three uses an extension of an algorithm
C due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
C but not all, of the improvements suggested by M. Gentleman and S.
C Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this
C program was published in the book Software Manual for the
C Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
C Englewood Cliffs, NJ, 1980. The present version is documented in
C W. J. Cody, "MACHAR: A subroutine to dynamically determine machine
C parameters," TOMS 14, December, 1988.
C
C The program as given here must be modified before compiling. If
C a single (double) precision version is desired, change all
C occurrences of CS (CD) in columns 1 and 2 to blanks.
C
C Parameter values reported are as follows:
C
C IBETA - the radix for the floating-point representation
C IT - the number of base IBETA digits in the floating-point
C significand
C IRND - 0 if floating-point addition chops
C 1 if floating-point addition rounds, but not in the
C IEEE style
C 2 if floating-point addition rounds in the IEEE style
C 3 if floating-point addition chops, and there is
C partial underflow
C 4 if floating-point addition rounds, but not in the
C IEEE style, and there is partial underflow
C 5 if floating-point addition rounds in the IEEE style,
C and there is partial underflow
C NGRD - the number of guard digits for multiplication with
C truncating arithmetic. It is
C 0 if floating-point arithmetic rounds, or if it
C truncates and only IT base IBETA digits
C participate in the post-normalization shift of the
C floating-point significand in multiplication;
C 1 if floating-point arithmetic truncates and more
C than IT base IBETA digits participate in the
C post-normalization shift of the floating-point
C significand in multiplication.
C MACHEP - the largest negative integer such that
C 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
C MACHEP is bounded below by -(IT+3)
C NEGEPS - the largest negative integer such that
C 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
C NEGEPS is bounded below by -(IT+3)
C IEXP - the number of bits (decimal places if IBETA = 10)
C reserved for the representation of the exponent
C (including the bias or sign) of a floating-point
C number
C MINEXP - the largest in magnitude negative integer such that
C FLOAT(IBETA)**MINEXP is positive and normalized
C MAXEXP - the smallest positive power of BETA that overflows
C EPS - the smallest positive floating-point number such
C that 1.0+EPS .NE. 1.0. In particular, if either
C IBETA = 2 or IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
C Otherwise, EPS = (FLOAT(IBETA)**MACHEP)/2
C EPSNEG - A small positive floating-point number such that
C 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2
C or IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
C Otherwise, EPSNEG = (IBETA**NEGEPS)/2. Because
C NEGEPS is bounded below by -(IT+3), EPSNEG may not
C be the smallest number that can alter 1.0 by
C subtraction.
C XMIN - the smallest non-vanishing normalized floating-point
C power of the radix, i.e., XMIN = FLOAT(IBETA)**MINEXP
C XMAX - the largest finite floating-point number. In
C particular XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
C Note - on some machines XMAX will be only the
C second, or perhaps third, largest number, being
C too small by 1 or 2 units in the last digit of
C the significand.
C
C Latest revision - December 4, 1987
C
C Author - W. J. Cody
C Argonne National Laboratory
C
C-----------------------------------------------------------------------}
VAR I, L, ITEMP, IZ,
J, K, MX, NXRES: LONGINT;
A, B, BETA, BETAIN,
BETAH, ONE, T, TEMP,
TEMPA, TEMP1, TWO,
Y, Z, ZERO: REAL;
CONV: ARRAY [0..10] OF REAL;
LABEL 1, 2, 10, 21, 22, 30, 32, 40,
41, 42, 43, 44, 45, 46, 50, 52;
BEGIN
{-----------------------------------------------------------------------}
FOR L := 1 TO 10 DO BEGIN
CONV [L] := L;
END;
ONE := CONV [1];
TWO := ONE + ONE;
ZERO := ONE - ONE;
{-----------------------------------------------------------------------}
{ Determine IBETA, BETA ala Malcolm. }
{-----------------------------------------------------------------------}
A := ONE;
1: A := A + A;
TEMP := A+ONE;
TEMP1 := TEMP-A;
IF TEMP1-ONE = ZERO THEN
GOTO 1;
B := ONE;
2: B := B + B;
TEMP := A+B;
ITEMP := TRUNC (TEMP-A);
IF ITEMP = 0 THEN
GOTO 2;
IBETA := ITEMP;
BETA := CONV [IBETA];
{-----------------------------------------------------------------------}
{ Determine IT, IRND. }
{-----------------------------------------------------------------------}
IT := 0;
B := ONE;
10:IT := IT + 1;
B := B * BETA;
TEMP := B+ONE;
TEMP1 := TEMP-B;
IF TEMP1-ONE = ZERO THEN
GOTO 10;
IRND := 0;
BETAH := BETA / TWO;
TEMP := A+BETAH;
IF TEMP-A <> ZERO THEN
IRND := 1;
TEMPA := A + BETA;
TEMP := TEMPA+BETAH;
IF (IRND = 0) AND (TEMP-TEMPA <> ZERO) THEN
IRND := 2;
{-----------------------------------------------------------------------}
{ Determine NEGEP, EPSNEG. }
{-----------------------------------------------------------------------}
NEGEP := IT + 3;
BETAIN := ONE / BETA;
A := ONE;
FOR I := 1 TO NEGEP DO BEGIN
A := A * BETAIN;
END;
B := A;
21:TEMP := ONE-A;
IF TEMP-ONE <> ZERO THEN
GOTO 22;
A := A * BETA;
NEGEP := NEGEP - 1;
GOTO 21;
22:NEGEP := -NEGEP;
EPSNEG := A;
{-----------------------------------------------------------------------}
{ Determine MACHEP, EPS. }
{-----------------------------------------------------------------------}
MACHEP := -IT - 3;
A := B;
30:TEMP := ONE+A;
IF TEMP-ONE <> ZERO THEN
GOTO 32;
A := A * BETA;
MACHEP := MACHEP + 1;
GOTO 30;
32:EPS := A;
{-----------------------------------------------------------------------}
{ Determine NGRD. }
{-----------------------------------------------------------------------}
NGRD := 0;
TEMP := ONE+EPS;
IF (IRND = 0) AND (TEMP*ONE-ONE <> ZERO) THEN
NGRD := 1;
{-----------------------------------------------------------------------
C Determine IEXP, MINEXP, XMIN.
C
C Loop to determine largest I and K = 2**I such that
C (1/BETA) ** (2**(I))
C does not underflow.
C Exit from loop is signaled by an underflow.
C-----------------------------------------------------------------------}
I := 0;
K := 1;
Z := BETAIN;
T := ONE + EPS;
NXRES := 0;
40:Y := Z;
Z := Y * Y;
{-----------------------------------------------------------------------}
{ Check for underflow here. }
{-----------------------------------------------------------------------}
A := Z * ONE;
TEMP := Z * T;
IF ((A+A = ZERO) OR (ABS(Z) >= Y)) THEN
GOTO 41;
TEMP1 := TEMP * BETAIN;
IF TEMP1*BETA = Z THEN
GOTO 41;
I := I + 1;
K := K + K;
GOTO 40;
41:IF IBETA = 10 THEN
GOTO 42;
IEXP := I + 1;
MX := K + K;
GOTO 45;
{-----------------------------------------------------------------------}
{ This segment is for decimal machines only. }
{-----------------------------------------------------------------------}
42:IEXP := 2;
IZ := IBETA;
43:IF K < IZ THEN
GOTO 44;
IZ := IZ * IBETA;
IEXP := IEXP + 1;
GOTO 43;
44:MX := IZ + IZ - 1;
{-----------------------------------------------------------------------}
{ Loop to determine MINEXP, XMIN. }
{ Exit from loop is signaled by an underflow. }
{-----------------------------------------------------------------------}
45:XMIN := Y;
Y := Y * BETAIN;
{-----------------------------------------------------------------------}
{ Check for underflow here. }
{-----------------------------------------------------------------------}
A := Y * ONE;
TEMP := Y * T;
IF (A+A = ZERO) OR (ABS(Y) >= XMIN) THEN
GOTO 46;
K := K + 1;
TEMP1 := TEMP * BETAIN;
IF (TEMP1*BETA <> Y) OR (TEMP = Y) THEN
GOTO 45
ELSE BEGIN
NXRES := 3;
XMIN := Y;
END;
46:MINEXP := -K;
{-----------------------------------------------------------------------}
{ Determine MAXEXP, XMAX. }
{-----------------------------------------------------------------------}
IF (MX > K+K-3) OR (IBETA = 10) THEN
GOTO 50;
MX := MX + MX;
IEXP := IEXP + 1;
50:MAXEXP := MX + MINEXP;
{-----------------------------------------------------------------}
{ Adjust IRND to reflect partial underflow. }
{-----------------------------------------------------------------}
IRND := IRND + NXRES;
{-----------------------------------------------------------------}
{ Adjust for IEEE-style machines. }
{-----------------------------------------------------------------}
IF IRND >= 2 THEN
MAXEXP := MAXEXP - 2;
{-----------------------------------------------------------------}
{ Adjust for machines with implicit leading bit in binary }
{ significand, and machines with radix point at extreme }
{ right of significand. }
{-----------------------------------------------------------------}
I := MAXEXP + MINEXP;
IF (IBETA = 2) AND (I = 0) THEN
MAXEXP := MAXEXP - 1;
IF I > 20 THEN
MAXEXP := MAXEXP - 1;
IF A <> Y THEN
MAXEXP := MAXEXP - 2;
XMAX := ONE - EPSNEG;
IF XMAX*ONE <> XMAX THEN
XMAX := ONE - BETA * EPSNEG;
XMAX := XMAX / (BETA * BETA * BETA * XMIN);
I := MAXEXP + MINEXP + 3;
IF I <= 0 THEN
GOTO 52;
FOR L := 1 TO I DO BEGIN
IF IBETA = 2 THEN
XMAX := XMAX + XMAX;
IF IBETA <> 2 THEN
XMAX := XMAX * BETA;
END;
52:
END; { MachAr }
PROCEDURE PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
MAXEXP: LONGINT; EPS,EPSNEG,XMIN,XMAX: REAL);
BEGIN
WRITELN ('MACHAR DETERMINED THE FOLLOWING PARAMETERS OF THE FLOATING-POINT ARITHMETIC');
WRITELN;
WRITELN ('RADIX FOR FLOATING-POINT REPRESENTATION: ', IBETA);
WRITELN ('NUMBER OF BASE ', IBETA:2, ' DIGITS IN SIGNIFICAND: ', IT);
IF IBETA <> 10 THEN
WRITELN ('NUMBER OF BITS USED TO REPRESENT THE EXPONENT: ', IEXP)
ELSE
WRITELN ('NUMBER OF DECIMAL PLACES USED FOR THE EXPONENT: ', IEXP);
WRITELN ('SMALLEST POSITIVE EPS SUCH THAT 1+EPS <> 1 : ', EPS:14, ' = ', IBETA, ' ** -', ABS(MACHEP));
WRITELN ('SMALLEST POSITIVE EPS SUCH THAT 1-EPS <> 1 : ', EPSNEG:14, ' = ', IBETA, ' ** -', ABS(NEGEP));
WRITELN ('SMALLEST POS. NORMALIZED FLOATING-POINT NUMBER:', XMIN:14, ' = ', IBETA, ' ** -', ABS(MINEXP)) ;
WRITELN ('LARGEST FLOATING-POINT NUMBER: ', XMAX:14, ' = ', IBETA, ' ** +', ABS(MAXEXP), ' - 1');
WRITELN;
WRITE ('FLOATING-POINT ARITHMETIC ');
CASE IRND MOD 3 OF
0: WRITELN ('CHOPS');
1: WRITELN ('ROUNDS, BUT NOT IEEE STYLE,');
2: WRITELN ('ROUNDS IEEE STYLE');
END;
IF IRND DIV 3 = 1 THEN
WRITELN ('AND SUPPORTS GRADUAL UNDERFLOW')
ELSE
WRITELN ('AND DOES NOT SUPPORT GRADUAL UNDERFLOW');
END; { PrintParam }
END. { MachArit }