Geek Gadgets 1

home *** CD-ROM | disk | FTP | other *** search

/ Geek Gadgets 1 / ADE-1.bin / ade-bin / lib / gcc-lib / m68k-amigaos / 2.7.2.1 / adainclude / s-fatgen.adb < prev next >

Wrap

Text File | 1995-10-17 | 17.0 KB | 659 lines

------------------------------------------------------------------------------ -- -- -- GNAT COMPILER COMPONENTS -- -- -- -- S Y S T E M . F A T _ G E N -- -- -- -- B o d y -- -- -- -- $Revision: 1.5 $ -- -- -- -- The GNAT library is free software; you can redistribute it and/or modify -- -- it under terms of the GNU Library General Public License as published by -- -- the Free Software Foundation; either version 2, or (at your option) any -- -- later version. The GNAT library is distributed in the hope that it will -- -- be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -- -- of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -- -- Library General Public License for more details. You should have -- -- received a copy of the GNU Library General Public License along with -- -- the GNAT library; see the file COPYING.LIB. If not, write to the Free -- -- Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -- -- -- ------------------------------------------------------------------------------ -- The implementation here is portable to any IEEE implementation. It does -- not handle non-binary radix, and also assumes that model numbers and -- machine numbers are basically identical, which is not true of all possible -- floating-point implementations. On a non-IEEE machine, this body must be -- specialized appropriately, or better still, its generic instantiations -- should be replaced by efficient machine-specific code. package body System.Fat_Gen is Float_Radix : constant T := T (T'Machine_Radix); Float_Radix_Inv : constant T := 1.0 / Float_Radix; Radix_To_M_Minus_1 : constant T := Float_Radix ** (T'Machine_Mantissa - 1); pragma Assert (T'Machine_Radix = 2); -- This version does not handle radix 16 -- Constants for Decompose and Scaling Rad : constant T := T (T'Machine_Radix); Invrad : constant T := 1.0 / Rad; subtype Expbits is Integer range 0 .. 6; -- 2 ** (2 ** 7) might overflow. how big can radix-16 exponents get? Log_Power : constant array (Expbits) of Integer := (1, 2, 4, 8, 16, 32, 64); R_Power : constant array (Expbits) of T := (Rad ** 1, Rad ** 2, Rad ** 4, Rad ** 8, Rad ** 16, Rad ** 32, Rad ** 64); R_Neg_Power : constant array (Expbits) of T := (Invrad ** 1, Invrad ** 2, Invrad ** 4, Invrad ** 8, Invrad ** 16, Invrad ** 32, Invrad ** 64); ----------------------- -- Local Subprograms -- ----------------------- procedure Decompose (XX : T; Frac : out T; Expo : out UI); -- Decomposes a floating-point number into fraction and exponent parts -------------- -- Adjacent -- -------------- -- Adjacent( X, 0.0) will be incorrect for abs(X) any power of two -- because of a bug in Succ and Pred. An extra test must be made -- for this case ???. function Adjacent (X, Towards : T) return T is begin if Towards = X then return X; elsif Towards > X then return Succ (X); else return Pred (X); end if; end Adjacent; ------------- -- Ceiling -- ------------- function Ceiling (X : T) return T is XT : constant T := Truncation (X); begin if X <= 0.0 then return XT; elsif X = XT then return X; else return XT + 1.0; end if; end Ceiling; ------------- -- Compose -- ------------- function Compose (Fraction : T; Exponent : UI) return T is Arg_Frac : T; Arg_Exp : UI; begin Decompose (Fraction, Arg_Frac, Arg_Exp); return Scaling (Arg_Frac, Exponent); end Compose; --------------- -- Copy_Sign -- --------------- -- Configuration requirement: If this unit is used to implement the -- floating-point attributes, then if Signed_Zeros is true, we must -- have Machine_Overflows False, and we must generate infinities for -- overflow (otherwise we can't implement signed zeroes properly). function Copy_Sign (Value, Sign : T) return T is Result : T; begin Result := abs Value; -- If we have signed zeroes true, then we know that machine overflows -- is false, and that infinities are generated for overflow, so we can -- test for minus zero by looking at the sign of the corresponding -- infinity value (neat trick eh?) if Sign = 0.0 and then T'Signed_Zeros then if 1.0 / Sign > 0.0 then return Result; else return -Result; end if; end if; -- Handle non-zero cases, and also the zero case if signed zeroes -- is false. In the latter case we always treat zero as positive. if Sign >= 0.0 then return Result; else return -Result; end if; end Copy_Sign; --------------- -- Decompose -- --------------- procedure Decompose (XX : T; Frac : out T; Expo : out UI) is X : T := T'Machine (XX); begin if X = 0.0 then Frac := X; Expo := 0; -- More useful would be defining Expo to be T'Machine_Emin - 1 or -- T'Machine_Emin - T'Machine_Mantissa, which would preserve -- monotonicity of the exponent fuction.. -- Check for infinities, transfinites, whatnot. elsif X > T'Safe_Last then Frac := Invrad; Expo := T'Machine_Emax + 1; elsif X < T'Safe_First then Frac := -Invrad; Expo := T'Machine_Emax + 2; -- how many extra negative values? else -- Case of nonzero finite x. Essentially, we just multiply -- by Rad ** (+-2**N) to reduce the range. declare Ax : T := abs X; Ex : UI := 0; -- Ax * Rad ** Ex is invariant. begin if Ax >= 1.0 then while Ax >= R_Power (Expbits'Last) loop Ax := Ax * R_Neg_Power (Expbits'Last); Ex := Ex + Log_Power (Expbits'Last); end loop; -- Ax < Rad ** 64 for N in reverse Expbits'First .. Expbits'Last - 1 loop if Ax >= R_Power (N) then Ax := Ax * R_Neg_Power (N); Ex := Ex + Log_Power (N); end if; -- Ax < R_Power (N) end loop; -- 1 <= Ax < Rad Ax := Ax * Invrad; Ex := Ex + 1; else -- 0 < ax < 1 while Ax < R_Neg_Power (Expbits'Last) loop Ax := Ax * R_Power (Expbits'Last); Ex := Ex - Log_Power (Expbits'Last); end loop; -- Rad ** -64 <= Ax < 1 for N in reverse Expbits'First .. Expbits'Last - 1 loop if Ax < R_Neg_Power (N) then Ax := Ax * R_Power (N); Ex := Ex - Log_Power (N); end if; -- R_Neg_Power (N) <= Ax < 1 end loop; end if; if X > 0.0 then Frac := Ax; else Frac := -Ax; end if; Expo := Ex; end; end if; end Decompose; -------------- -- Exponent -- -------------- function Exponent (X : T) return UI is X_Frac : T; X_Exp : UI; begin Decompose (X, X_Frac, X_Exp); return X_Exp; end Exponent; ----------- -- Floor -- ----------- function Floor (X : T) return T is XT : constant T := Truncation (X); begin if X >= 0.0 then return XT; elsif XT = X then return X; else return XT - 1.0; end if; end Floor; -------------- -- Fraction -- -------------- function Fraction (X : T) return T is X_Frac : T; X_Exp : UI; begin Decompose (X, X_Frac, X_Exp); return X_Frac; end Fraction; ------------------ -- Leading_Part -- ------------------ function Leading_Part (X : T; Radix_Digits : UI) return T is L : UI; Y, Z : T; begin if Radix_Digits >= T'Machine_Mantissa then return X; else L := Exponent (X) - Radix_Digits; Y := Truncation (Scaling (X, -L)); Z := Scaling (Y, L); return Z; end if; end Leading_Part; ------------- -- Machine -- ------------- -- The trick with Machine is to force the compiler to store the result -- in memory so that we do not have extra precision used. The compiler -- is clever, so we have to outwit its possible optimizations! -- This is achieved by using the following array. In fact only the first -- element is used, and Machine_Index is always 1, but we make sure the -- compiler can't figure this out. Machine_Array : array (1 .. 2) of T; Machine_Index : Integer; function Machine (X : T) return T is begin Machine_Array (1) := X; return Machine_Array (Machine_Index); end Machine; ----------- -- Model -- ----------- -- We treat Model as identical to Machine. This is true of IEEE and other -- nice floating-point systems, but not necessarily true of all systems. function Model (X : T) return T is begin return Machine (X); end Model; ---------- -- Pred -- ---------- -- Subtract from the given number a number equivalent to the value of its -- least significant bit. Given that the most significant bit represents a -- value of 1.0 * radix ** (exp - 1), the value we want is obtained by -- shifting this by (mantissa-1) bits to the right, i.e. decreasing the -- exponent by that amount. function Pred (X : T) return T is X_Frac : T; X_Exp : UI; begin Decompose (X, X_Frac, X_Exp); return X - Scaling (1.0, X_Exp - T'Machine_Mantissa); end Pred; --------------- -- Remainder -- --------------- function Remainder (X, Y : T) return T is A : T; B : T; Arg : T; P : T; Arg_Frac : T; P_Frac : T; Sign_X : T; IEEE_Rem : T; Arg_Exp : UI; P_Exp : UI; K : UI; P_Even : Boolean; begin if X > 0.0 then Sign_X := 1.0; else Sign_X := -1.0; end if; Arg := abs X; P := abs Y; if Arg < P then P_Even := True; IEEE_Rem := Arg; P_Exp := Exponent (P); else Decompose (Arg, Arg_Frac, Arg_Exp); Decompose (P, P_Frac, P_Exp); P := Compose (P_Frac, Arg_Exp); K := UI (Arg_Exp) - UI (P_Exp); P_Even := True; IEEE_Rem := Arg; for Cnt in reverse 0 .. K loop if IEEE_Rem >= P then P_Even := False; IEEE_Rem := IEEE_Rem - P; else P_Even := True; end if; P := P * 0.5; end loop; end if; -- That completes the calculation of modulus remainder. The final -- step is get the IEEE remainder. Here we need to compare Rem with -- (abs Y) / 2. We must be careful of unrepresentable Y/2 value -- caused by subnormal numbers if P_Exp >= 0 then A := IEEE_Rem; B := abs Y * 0.5; else A := IEEE_Rem * 2.0; B := abs Y; end if; if A > B or else (A = B and then not P_Even) then IEEE_Rem := IEEE_Rem - abs Y; end if; return Sign_X * IEEE_Rem; end Remainder; -------------- -- Rounding -- -------------- function Rounding (X : T) return T is Result : T; Tail : T; begin Result := Truncation (abs X); Tail := abs X - Result; if Tail >= 0.5 then Result := Result + 1.0; end if; if X > 0.0 then return Result; elsif X < 0.0 then return -Result; -- For zero case, make sure sign of zero is preserved else return X; end if; end Rounding; ------------- -- Scaling -- ------------- -- Return x * rad ** adjustment quickly, -- or quietly underflow to zero, or overflow naturally. function Scaling (X : T; Adjustment : UI) return T is begin if X = 0.0 or else Adjustment = 0 then return X; end if; -- Nonzero x. essentially, just multiply repeatedly by Rad ** (+-2**n). declare Y : T := X; Ex : UI := Adjustment; -- Y * Rad ** Ex is invariant begin if Ex < 0 then while Ex <= -Log_Power (Expbits'Last) loop Y := Y * R_Neg_Power (Expbits'Last); Ex := Ex + Log_Power (Expbits'Last); end loop; -- -64 < Ex <= 0 for N in reverse Expbits'First .. Expbits'Last - 1 loop if Ex <= -Log_Power (N) then Y := Y * R_Neg_Power (N); Ex := Ex + Log_Power (N); end if; -- -Log_Power (N) < Ex <= 0 end loop; -- Ex = 0 else -- Ex >= 0 while Ex >= Log_Power (Expbits'Last) loop Y := Y * R_Power (Expbits'Last); Ex := Ex - Log_Power (Expbits'Last); end loop; -- 0 <= Ex < 64 for N in reverse Expbits'First .. Expbits'Last - 1 loop if Ex >= Log_Power (N) then Y := Y * R_Power (N); Ex := Ex - Log_Power (N); end if; -- 0 <= Ex < Log_Power (N) end loop; -- Ex = 0 end if; return Y; end; end Scaling; ---------- -- Succ -- ---------- -- Similar computation to that of Pred: find value of least significant -- bit of given number, and add. function Succ (X : T) return T is X_Frac : T; X_Exp : UI; begin Decompose (X, X_Frac, X_Exp); return X + Scaling (1.0, X_Exp - T'Machine_Mantissa); end Succ; ---------------- -- Truncation -- ---------------- -- The basic approach is to compute -- T'Machine (RM1 + N) - RM1. -- where N >= 0.0 and RM1 = radix ** (mantissa - 1) -- This works provided that the intermediate result (RM1 + N) does not -- have extra precision (which is why we call Machine). When we compute -- RM1 + N, the epxonent of N will be normalized and the mantissa shifted -- shifted appropriately so the lower order bits, which cannot contribute -- to the integer part of N, fall off on the right. When we subtract RM1 -- again, the significant bits of N are shifted to the left, and what we -- have is an integer, because only the first e bits are different from -- zero (assuming binary radix here). function Truncation (X : T) return T is Result : T; begin Result := abs X; if Result >= Radix_To_M_Minus_1 then return Machine (X); else Result := Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1; if Result > abs X then Result := Result - 1.0; end if; if X > 0.0 then return Result; elsif X < 0.0 then return -Result; -- For zero case, make sure sign of zero is preserved else return X; end if; end if; end Truncation; ----------------------- -- Unbiased_Rounding -- ----------------------- function Unbiased_Rounding (X : T) return T is Result : T; Tail : T; begin Result := Truncation (abs X); Tail := abs X - Result; if Tail > 0.5 then Result := Result + 1.0; elsif Tail = 0.5 then Result := 2.0 * Truncation ((Result / 2.0) + 0.5); end if; if X > 0.0 then return Result; elsif X < 0.0 then return -Result; -- For zero case, make sure sign of zero is preserved else return X; end if; end Unbiased_Rounding; ---------------------------- -- Package Initialization -- ---------------------------- begin -- Set Machine_Index to 1, but not in an obvious manner (see function -- Machine to understand why we are behaving in this secretive manner) Machine_Index := 2 ** 3; for J in 1 .. 3 loop Machine_Index := Machine_Index / 2; end loop; end System.Fat_Gen;