------------------------------------------------------------------------------ | |

-- -- | |

-- GNAT COMPILER COMPONENTS -- | |

-- -- | |

-- S Y S T E M . F A T _ G E N -- | |

-- -- | |

-- B o d y -- | |

-- -- | |

-- Copyright (C) 1992-2021, Free Software Foundation, Inc. -- | |

-- -- | |

-- GNAT is free software; you can redistribute it and/or modify it under -- | |

-- terms of the GNU General Public License as published by the Free Soft- -- | |

-- ware Foundation; either version 3, or (at your option) any later ver- -- | |

-- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- | |

-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- | |

-- or FITNESS FOR A PARTICULAR PURPOSE. -- | |

-- -- | |

-- As a special exception under Section 7 of GPL version 3, you are granted -- | |

-- additional permissions described in the GCC Runtime Library Exception, -- | |

-- version 3.1, as published by the Free Software Foundation. -- | |

-- -- | |

-- You should have received a copy of the GNU General Public License and -- | |

-- a copy of the GCC Runtime Library Exception along with this program; -- | |

-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- | |

-- <http://www.gnu.org/licenses/>. -- | |

-- -- | |

-- GNAT was originally developed by the GNAT team at New York University. -- | |

-- Extensive contributions were provided by Ada Core Technologies Inc. -- | |

-- -- | |

------------------------------------------------------------------------------ | |

-- This implementation is portable to any IEEE implementation. It does not | |

-- handle nonbinary radix, and also assumes that model numbers and machine | |

-- numbers are basically identical, which is not true of all possible | |

-- floating-point implementations. | |

with Ada.Unchecked_Conversion; | |

with System.Unsigned_Types; | |

pragma Warnings (Off, "non-static constant in preelaborated unit"); | |

-- Every constant is static given our instantiation model | |

package body System.Fat_Gen is | |

pragma Assert (T'Machine_Radix = 2); | |

-- This version does not handle radix 16 | |

Rad : constant T := T (T'Machine_Radix); | |

-- Renaming for the machine radix | |

Mantissa : constant Integer := T'Machine_Mantissa; | |

-- Renaming for the machine mantissa | |

Invrad : constant T := 1.0 / Rad; | |

-- Smallest positive mantissa in the canonical form (RM A.5.3(4)) | |

-- Small : constant T := Rad ** (T'Machine_Emin - 1); | |

-- Smallest positive normalized number | |

-- Tiny : constant T := Rad ** (T'Machine_Emin - Mantissa); | |

-- Smallest positive denormalized number | |

RM1 : constant T := Rad ** (Mantissa - 1); | |

-- Smallest positive member of the large consecutive integers. It is equal | |

-- to the ratio Small / Tiny, which means that multiplying by it normalizes | |

-- any nonzero denormalized number. | |

IEEE_Emin : constant Integer := T'Machine_Emin - 1; | |

IEEE_Emax : constant Integer := T'Machine_Emax - 1; | |

-- The mantissa is a fraction with first digit set in Ada whereas it is | |

-- shifted by 1 digit to the left in the IEEE floating-point format. | |

subtype IEEE_Erange is Integer range IEEE_Emin - 1 .. IEEE_Emax + 1; | |

-- The IEEE floating-point format extends the machine range by 1 to the | |

-- left for denormalized numbers and 1 to the right for infinities/NaNs. | |

IEEE_Ebias : constant Integer := -(IEEE_Emin - 1); | |

-- The exponent is biased such that denormalized numbers have it zero | |

-- The implementation uses a representation type Float_Rep that allows | |

-- direct access to exponent and mantissa of the floating point number. | |

-- The Float_Rep type is a simple array of Float_Word elements. This | |

-- representation is chosen to make it possible to size the type based | |

-- on a generic parameter. Since the array size is known at compile | |

-- time, efficient code can still be generated. The size of Float_Word | |

-- elements should be large enough to allow accessing the exponent in | |

-- one read, but small enough so that all floating-point object sizes | |

-- are a multiple of Float_Word'Size. | |

-- The following conditions must be met for all possible instantiations | |

-- of the attribute package: | |

-- - T'Size is an integral multiple of Float_Word'Size | |

-- - The exponent and sign are completely contained in a single | |

-- component of Float_Rep, named Most Significant Word (MSW). | |

-- - The sign occupies the most significant bit of the MSW and the | |

-- exponent is in the following bits. | |

-- The low-level primitives Copy_Sign, Decompose, Finite_Succ, Scaling and | |

-- Valid are implemented by accessing the bit pattern of the floating-point | |

-- number. Only the normalization of denormalized numbers, if any, and the | |

-- gradual underflow are left to the hardware, mainly because there is some | |

-- leeway for the hardware implementation in this area: for example the MSB | |

-- of the mantissa, that is 1 for normalized numbers and 0 for denormalized | |

-- numbers, may or may not be stored by the hardware. | |

Siz : constant := 16; | |

type Float_Word is mod 2**Siz; | |

-- We use the GCD of the size of all the supported floating-point formats | |

N : constant Natural := (T'Size + Siz - 1) / Siz; | |

NR : constant Natural := (Mantissa + 16 + Siz - 1) / Siz; | |

Rep_Last : constant Natural := Natural'Min (N, NR) - 1; | |

-- Determine the number of Float_Words needed for representing the | |

-- entire floating-point value. Do not take into account excessive | |

-- padding, as occurs on IA-64 where 80 bits floats get padded to 128 | |

-- bits. In general, the exponent field cannot be larger than 15 bits, | |

-- even for 128-bit floating-point types, so the final format size | |

-- won't be larger than Mantissa + 16. | |

type Float_Rep is array (Natural range 0 .. N - 1) of Float_Word; | |

pragma Suppress_Initialization (Float_Rep); | |

-- This pragma suppresses the generation of an initialization procedure | |

-- for type Float_Rep when operating in Initialize/Normalize_Scalars mode. | |

MSW : constant Natural := Rep_Last * Standard'Default_Bit_Order; | |

-- Finding the location of the Exponent_Word is a bit tricky. In general | |

-- we assume Word_Order = Bit_Order. | |

Exp_Factor : constant Float_Word := | |

2**(Siz - 1) / Float_Word (IEEE_Emax - IEEE_Emin + 3); | |

-- Factor that the extracted exponent needs to be divided by to be in | |

-- range 0 .. IEEE_Emax - IEEE_Emin + 2 | |

Exp_Mask : constant Float_Word := | |

Float_Word (IEEE_Emax - IEEE_Emin + 2) * Exp_Factor; | |

-- Value needed to mask out the exponent field. This assumes that the | |

-- range 0 .. IEEE_Emax - IEEE_Emin + 2 contains 2**N values, for some | |

-- N in Natural. | |

Sign_Mask : constant Float_Word := 2**(Siz - 1); | |

-- Value needed to mask out the sign field | |

----------------------- | |

-- Local Subprograms -- | |

----------------------- | |

procedure Decompose (XX : T; Frac : out T; Expo : out UI); | |

-- Decomposes a floating-point number into fraction and exponent parts. | |

-- Both results are signed, with Frac having the sign of XX, and UI has | |

-- the sign of the exponent. The absolute value of Frac is in the range | |

-- 0.0 <= Frac < 1.0. If Frac = 0.0 or -0.0, then Expo is always zero. | |

function Finite_Succ (X : T) return T; | |

-- Return the successor of X, a finite number not equal to T'Last | |

-------------- | |

-- Adjacent -- | |

-------------- | |

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; | |

pragma Unreferenced (Arg_Exp); | |

begin | |

Decompose (Fraction, Arg_Frac, Arg_Exp); | |

return Scaling (Arg_Frac, Exponent); | |

end Compose; | |

--------------- | |

-- Copy_Sign -- | |

--------------- | |

function Copy_Sign (Value, Sign : T) return T is | |

S : constant T := T'Machine (Sign); | |

Rep_S : Float_Rep; | |

for Rep_S'Address use S'Address; | |

-- Rep_S is a view of the Sign parameter | |

V : T := T'Machine (Value); | |

Rep_V : Float_Rep; | |

for Rep_V'Address use V'Address; | |

-- Rep_V is a view of the Value parameter | |

begin | |

Rep_V (MSW) := | |

(Rep_V (MSW) and not Sign_Mask) or (Rep_S (MSW) and Sign_Mask); | |

return V; | |

end Copy_Sign; | |

--------------- | |

-- Decompose -- | |

--------------- | |

procedure Decompose (XX : T; Frac : out T; Expo : out UI) is | |

X : T := T'Machine (XX); | |

Rep : Float_Rep; | |

for Rep'Address use X'Address; | |

-- Rep is a view of the input floating-point parameter | |

Exp : constant IEEE_Erange := | |

Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias; | |

-- Mask/Shift X to only get bits from the exponent. Then convert biased | |

-- value to final value. | |

Minus : constant Boolean := (Rep (MSW) and Sign_Mask) /= 0; | |

-- Mask/Shift X to only get bit from the sign | |

begin | |

-- The normalized exponent of zero is zero, see RM A.5.3(15) | |

if X = 0.0 then | |

Expo := 0; | |

Frac := X; | |

-- Check for infinities and NaNs | |

elsif Exp = IEEE_Emax + 1 then | |

Expo := T'Machine_Emax + 1; | |

Frac := (if Minus then -Invrad else Invrad); | |

-- Check for nonzero denormalized numbers | |

elsif Exp = IEEE_Emin - 1 then | |

-- Normalize by multiplying by Radix ** (Mantissa - 1) | |

Decompose (X * RM1, Frac, Expo); | |

Expo := Expo - (Mantissa - 1); | |

-- Case of normalized numbers | |

else | |

-- The Ada exponent is the IEEE exponent plus 1, see above | |

Expo := Exp + 1; | |

-- Set Ada exponent of X to zero, so we end up with the fraction | |

Rep (MSW) := (Rep (MSW) and not Exp_Mask) + | |

Float_Word (IEEE_Ebias - 1) * Exp_Factor; | |

Frac := X; | |

end if; | |

end Decompose; | |

-------------- | |

-- Exponent -- | |

-------------- | |

function Exponent (X : T) return UI is | |

X_Frac : T; | |

X_Exp : UI; | |

pragma Unreferenced (X_Frac); | |

begin | |

Decompose (X, X_Frac, X_Exp); | |

return X_Exp; | |

end Exponent; | |

----------------- | |

-- Finite_Succ -- | |

----------------- | |

function Finite_Succ (X : T) return T is | |

XX : T := T'Machine (X); | |

Rep : Float_Rep; | |

for Rep'Address use XX'Address; | |

-- Rep is a view of the input floating-point parameter | |

begin | |

-- If the floating-point type does not support denormalized numbers, | |

-- there is a couple of problematic values, namely -Small and Zero, | |

-- because the increment is equal to Small in these cases. | |

if not T'Denorm then | |

declare | |

Small : constant T := Rad ** (T'Machine_Emin - 1); | |

-- Smallest positive normalized number declared here and not at | |

-- library level for the sake of the CCG compiler, which cannot | |

-- currently compile the constant because the target is C90. | |

begin | |

if X = -Small then | |

XX := 0.0; | |

return -XX; | |

elsif X = 0.0 then | |

return Small; | |

end if; | |

end; | |

end if; | |

-- In all the other cases, the increment is equal to 1 in the binary | |

-- integer representation of the number if X is nonnegative and equal | |

-- to -1 if X is negative. | |

if XX >= 0.0 then | |

-- First clear the sign of negative Zero | |

Rep (MSW) := Rep (MSW) and not Sign_Mask; | |

-- Deal with big endian | |

if MSW = 0 then | |

for J in reverse 0 .. Rep_Last loop | |

Rep (J) := Rep (J) + 1; | |

-- For 80-bit IEEE Extended, the MSB of the mantissa is stored | |

-- so, when it has been flipped, its status must be reanalyzed. | |

if Mantissa = 64 and then J = 1 then | |

-- If the MSB changed from denormalized to normalized, then | |

-- keep it normalized since the exponent will be bumped. | |

if Rep (J) = 2**(Siz - 1) then | |

null; | |

-- If the MSB changed from normalized, restore it since we | |

-- cannot denormalize in this context. | |

elsif Rep (J) = 0 then | |

Rep (J) := 2**(Siz - 1); | |

else | |

exit; | |

end if; | |

-- In other cases, stop if there is no carry | |

else | |

exit when Rep (J) > 0; | |

end if; | |

end loop; | |

-- Deal with little endian | |

else | |

for J in 0 .. Rep_Last loop | |

Rep (J) := Rep (J) + 1; | |

-- For 80-bit IEEE Extended, the MSB of the mantissa is stored | |

-- so, when it has been flipped, its status must be reanalyzed. | |

if Mantissa = 64 and then J = Rep_Last - 1 then | |

-- If the MSB changed from denormalized to normalized, then | |

-- keep it normalized since the exponent will be bumped. | |

if Rep (J) = 2**(Siz - 1) then | |

null; | |

-- If the MSB changed from normalized, restore it since we | |

-- cannot denormalize in this context. | |

elsif Rep (J) = 0 then | |

Rep (J) := 2**(Siz - 1); | |

else | |

exit; | |

end if; | |

-- In other cases, stop if there is no carry | |

else | |

exit when Rep (J) > 0; | |

end if; | |

end loop; | |

end if; | |

else | |

if MSW = 0 then | |

for J in reverse 0 .. Rep_Last loop | |

Rep (J) := Rep (J) - 1; | |

-- For 80-bit IEEE Extended, the MSB of the mantissa is stored | |

-- so, when it has been flipped, its status must be reanalyzed. | |

if Mantissa = 64 and then J = 1 then | |

-- If the MSB changed from normalized to denormalized, then | |

-- keep it normalized if the exponent is not 1. | |

if Rep (J) = 2**(Siz - 1) - 1 then | |

if Rep (0) /= 2**(Siz - 1) + 1 then | |

Rep (J) := 2**Siz - 1; | |

end if; | |

else | |

exit; | |

end if; | |

-- In other cases, stop if there is no borrow | |

else | |

exit when Rep (J) < 2**Siz - 1; | |

end if; | |

end loop; | |

else | |

for J in 0 .. Rep_Last loop | |

Rep (J) := Rep (J) - 1; | |

-- For 80-bit IEEE Extended, the MSB of the mantissa is stored | |

-- so, when it has been flipped, its status must be reanalyzed. | |

if Mantissa = 64 and then J = Rep_Last - 1 then | |

-- If the MSB changed from normalized to denormalized, then | |

-- keep it normalized if the exponent is not 1. | |

if Rep (J) = 2**(Siz - 1) - 1 then | |

if Rep (Rep_Last) /= 2**(Siz - 1) + 1 then | |

Rep (J) := 2**Siz - 1; | |

end if; | |

else | |

exit; | |

end if; | |

-- In other cases, stop if there is no borrow | |

else | |

exit when Rep (J) < 2**Siz - 1; | |

end if; | |

end loop; | |

end if; | |

end if; | |

return XX; | |

end Finite_Succ; | |

----------- | |

-- 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; | |

pragma Unreferenced (X_Exp); | |

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 >= Mantissa then | |

return X; | |

elsif Radix_Digits <= 0 then | |

raise Constraint_Error; | |

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. We do | |

-- this by using an intermediate pragma Volatile location. | |

function Machine (X : T) return T is | |

Temp : T; | |

pragma Volatile (Temp); | |

begin | |

Temp := X; | |

return Temp; | |

end Machine; | |

---------------------- | |

-- Machine_Rounding -- | |

---------------------- | |

-- For now, the implementation is identical to that of Rounding, which is | |

-- a permissible behavior, but is not the most efficient possible approach. | |

function Machine_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 Machine_Rounding; | |

----------- | |

-- 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 T'Machine (X); | |

end Model; | |

---------- | |

-- Pred -- | |

---------- | |

function Pred (X : T) return T is | |

begin | |

-- Special treatment for largest negative number: raise Constraint_Error | |

if X = T'First then | |

raise Constraint_Error with "Pred of largest negative number"; | |

-- For finite numbers, use the symmetry around zero of floating point | |

elsif X > T'First and then X <= T'Last then | |

pragma Annotate (CodePeer, Intentional, "test always true", | |

"Check for invalid float"); | |

pragma Annotate (CodePeer, Intentional, "condition predetermined", | |

"Check for invalid float"); | |

return -Finite_Succ (-X); | |

-- For infinities and NaNs, return unchanged | |

else | |

return X; | |

pragma Annotate (CodePeer, Intentional, "dead code", | |

"Check float range."); | |

end if; | |

end Pred; | |

--------------- | |

-- Remainder -- | |

--------------- | |

function Remainder (X, Y : T) return T is | |

A : T; | |

B : T; | |

Arg : T; | |

P : T; | |

P_Frac : T; | |

Sign_X : T; | |

IEEE_Rem : T; | |

Arg_Exp : UI; | |

P_Exp : UI; | |

K : UI; | |

P_Even : Boolean; | |

Arg_Frac : T; | |

pragma Unreferenced (Arg_Frac); | |

begin | |

if Y = 0.0 then | |

raise Constraint_Error; | |

end if; | |

if X > 0.0 then | |

Sign_X := 1.0; | |

Arg := X; | |

else | |

Sign_X := -1.0; | |

Arg := -X; | |

end if; | |

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 := Arg_Exp - 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 -- | |

------------- | |

function Scaling (X : T; Adjustment : UI) return T is | |

pragma Assert (Mantissa <= 64); | |

-- This implementation handles only 80-bit IEEE Extended or smaller | |

package UST renames System.Unsigned_Types; | |

use type UST.Long_Long_Unsigned; | |

XX : T := T'Machine (X); | |

Rep : Float_Rep; | |

for Rep'Address use XX'Address; | |

-- Rep is a view of the input floating-point parameter | |

Exp : constant IEEE_Erange := | |

Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias; | |

-- Mask/Shift X to only get bits from the exponent. Then convert biased | |

-- value to final value. | |

Minus : constant Boolean := (Rep (MSW) and Sign_Mask) /= 0; | |

-- Mask/Shift X to only get bit from the sign | |

Expi, Expf : IEEE_Erange; | |

begin | |

-- Check for zero, infinities, NaNs as well as no adjustment | |

if X = 0.0 or else Exp = IEEE_Emax + 1 or else Adjustment = 0 then | |

return X; | |

-- Check for nonzero denormalized numbers | |

elsif Exp = IEEE_Emin - 1 then | |

-- Check for zero result to protect the subtraction below | |

if Adjustment < -(Mantissa - 1) then | |

XX := 0.0; | |

return (if Minus then -XX else XX); | |

-- Normalize by multiplying by Radix ** (Mantissa - 1) | |

else | |

return Scaling (XX * RM1, Adjustment - (Mantissa - 1)); | |

end if; | |

-- Case of normalized numbers | |

else | |

-- Check for overflow | |

if Adjustment > IEEE_Emax - Exp then | |

-- Optionally raise Constraint_Error as per RM A.5.3(29) | |

if T'Machine_Overflows then | |

raise Constraint_Error with "Too large exponent"; | |

else | |

XX := 0.0; | |

return (if Minus then -1.0 / XX else 1.0 / XX); | |

pragma Annotate (CodePeer, Intentional, "overflow check", | |

"Infinity produced"); | |

pragma Annotate (CodePeer, Intentional, "divide by zero", | |

"Infinity produced"); | |

end if; | |

-- Check for underflow | |

elsif Adjustment < IEEE_Emin - Exp then | |

-- Check for possibly gradual underflow (up to the hardware) | |

if Adjustment >= IEEE_Emin - Mantissa - Exp then | |

Expf := IEEE_Emin; | |

Expi := Exp + Adjustment - Expf; | |

-- Case of zero result | |

else | |

XX := 0.0; | |

return (if Minus then -XX else XX); | |

end if; | |

-- Case of normalized results | |

else | |

Expf := Exp + Adjustment; | |

Expi := 0; | |

end if; | |

Rep (MSW) := (Rep (MSW) and not Exp_Mask) + | |

Float_Word (IEEE_Ebias + Expf) * Exp_Factor; | |

if Expi < 0 then | |

-- Given that Expi >= -Mantissa, only -64 is problematic | |

if Expi = -64 then | |

pragma Annotate | |

(CodePeer, Intentional, "test always false", | |

"test always false in some instantiations"); | |

XX := XX / 2.0; | |

Expi := -63; | |

end if; | |

XX := XX / T (UST.Long_Long_Unsigned (2) ** (-Expi)); | |

end if; | |

return XX; | |

end if; | |

end Scaling; | |

---------- | |

-- Succ -- | |

---------- | |

function Succ (X : T) return T is | |

begin | |

-- Special treatment for largest positive number: raise Constraint_Error | |

if X = T'Last then | |

raise Constraint_Error with "Succ of largest positive number"; | |

-- For finite numbers, call the specific routine | |

elsif X >= T'First and then X < T'Last then | |

pragma Annotate (CodePeer, Intentional, "test always true", | |

"Check for invalid float"); | |

pragma Annotate (CodePeer, Intentional, "condition predetermined", | |

"Check for invalid float"); | |

return Finite_Succ (X); | |

-- For infinities and NaNs, return unchanged | |

else | |

return X; | |

pragma Annotate (CodePeer, Intentional, "dead code", | |

"Check float range."); | |

end if; | |

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 exponent of N will be normalized and the mantissa 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 >= RM1 then | |

return T'Machine (X); | |

else | |

Result := T'Machine (RM1 + Result) - RM1; | |

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 | |

Abs_X : constant T := abs X; | |

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; | |

----------- | |

-- Valid -- | |

----------- | |

function Valid (X : not null access T) return Boolean is | |

type Access_T is access all T; | |

function To_Address is | |

new Ada.Unchecked_Conversion (Access_T, System.Address); | |

Rep : Float_Rep; | |

for Rep'Address use To_Address (Access_T (X)); | |

-- Rep is a view of the input floating-point parameter. Note that we | |

-- must avoid reading the actual bits of this parameter in float form | |

-- since it may be a signalling NaN. | |

Exp : constant IEEE_Erange := | |

Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias; | |

-- Mask/Shift X to only get bits from the exponent. Then convert biased | |

-- value to final value. | |

begin | |

if Exp = IEEE_Emax + 1 then | |

-- This is an infinity or a NaN, i.e. always invalid | |

return False; | |

elsif Exp in IEEE_Emin .. IEEE_Emax then | |

-- This is a normalized number, i.e. always valid | |

return True; | |

else pragma Assert (Exp = IEEE_Emin - 1); | |

-- This is a denormalized number, valid if T'Denorm is True or 0.0 | |

if T'Denorm then | |

return True; | |

-- Note that we cannot do a direct comparison with 0.0 because the | |

-- hardware may evaluate it to True for all denormalized numbers. | |

else | |

-- First clear the sign bit (the exponent is already zero) | |

Rep (MSW) := Rep (MSW) and not Sign_Mask; | |

return (for all J in 0 .. Rep_Last => Rep (J) = 0); | |

end if; | |

end if; | |

end Valid; | |

end System.Fat_Gen; |