932 lines
27 KiB
Ada
932 lines
27 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- GNAT COMPILER COMPONENTS --
|
|
-- --
|
|
-- S Y S T E M . F A T _ G E N --
|
|
-- --
|
|
-- B o d y --
|
|
-- --
|
|
-- Copyright (C) 1992-2015, 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. --
|
|
-- --
|
|
------------------------------------------------------------------------------
|
|
|
|
-- The implementation here 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. 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.
|
|
|
|
with Ada.Unchecked_Conversion;
|
|
with System;
|
|
package body System.Fat_Gen is
|
|
|
|
Float_Radix : constant T := T (T'Machine_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.
|
|
-- 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 Gradual_Scaling (Adjustment : UI) return T;
|
|
-- Like Scaling with a first argument of 1.0, but returns the smallest
|
|
-- denormal rather than zero when the adjustment is smaller than
|
|
-- Machine_Emin. Used for Succ and Pred.
|
|
|
|
--------------
|
|
-- 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
|
|
Result : T;
|
|
|
|
function Is_Negative (V : T) return Boolean;
|
|
pragma Import (Intrinsic, Is_Negative);
|
|
|
|
begin
|
|
Result := abs Value;
|
|
|
|
if Is_Negative (Sign) then
|
|
return -Result;
|
|
else
|
|
return Result;
|
|
end if;
|
|
end Copy_Sign;
|
|
|
|
---------------
|
|
-- Decompose --
|
|
---------------
|
|
|
|
procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
|
|
X : constant T := T'Machine (XX);
|
|
|
|
begin
|
|
if X = 0.0 then
|
|
|
|
-- The normalized exponent of zero is zero, see RM A.5.2(15)
|
|
|
|
Frac := X;
|
|
Expo := 0;
|
|
|
|
-- 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;
|
|
|
|
Frac := (if X > 0.0 then Ax else -Ax);
|
|
Expo := Ex;
|
|
end;
|
|
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;
|
|
|
|
-----------
|
|
-- 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;
|
|
|
|
---------------------
|
|
-- Gradual_Scaling --
|
|
---------------------
|
|
|
|
function Gradual_Scaling (Adjustment : UI) return T is
|
|
Y : T;
|
|
Y1 : T;
|
|
Ex : UI := Adjustment;
|
|
|
|
begin
|
|
if Adjustment < T'Machine_Emin - 1 then
|
|
Y := 2.0 ** T'Machine_Emin;
|
|
Y1 := Y;
|
|
Ex := Ex - T'Machine_Emin;
|
|
while Ex < 0 loop
|
|
Y := T'Machine (Y / 2.0);
|
|
|
|
if Y = 0.0 then
|
|
return Y1;
|
|
end if;
|
|
|
|
Ex := Ex + 1;
|
|
Y1 := Y;
|
|
end loop;
|
|
|
|
return Y1;
|
|
|
|
else
|
|
return Scaling (1.0, Adjustment);
|
|
end if;
|
|
end Gradual_Scaling;
|
|
|
|
------------------
|
|
-- 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;
|
|
|
|
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 Machine (X);
|
|
end Model;
|
|
|
|
----------
|
|
-- Pred --
|
|
----------
|
|
|
|
function Pred (X : T) return T is
|
|
X_Frac : T;
|
|
X_Exp : UI;
|
|
|
|
begin
|
|
-- Zero has to be treated specially, since its exponent is zero
|
|
|
|
if X = 0.0 then
|
|
return -Succ (X);
|
|
|
|
-- Special treatment for most negative number
|
|
|
|
elsif X = T'First then
|
|
|
|
-- If not generating infinities, we raise a constraint error
|
|
|
|
if T'Machine_Overflows then
|
|
raise Constraint_Error with "Pred of largest negative number";
|
|
|
|
-- Otherwise generate a negative infinity
|
|
|
|
else
|
|
return X / (X - X);
|
|
end if;
|
|
|
|
-- For infinities, return unchanged
|
|
|
|
elsif X < T'First or else X > T'Last then
|
|
return X;
|
|
|
|
-- 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.
|
|
|
|
else
|
|
Decompose (X, X_Frac, X_Exp);
|
|
|
|
-- A special case, if the number we had was a positive power of
|
|
-- two, then we want to subtract half of what we would otherwise
|
|
-- subtract, since the exponent is going to be reduced.
|
|
|
|
-- Note that X_Frac has the same sign as X, so if X_Frac is 0.5,
|
|
-- then we know that we have a positive number (and hence a
|
|
-- positive power of 2).
|
|
|
|
if X_Frac = 0.5 then
|
|
return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
|
|
|
|
-- Otherwise the exponent is unchanged
|
|
|
|
else
|
|
return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa);
|
|
end if;
|
|
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 --
|
|
-------------
|
|
|
|
-- 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 --
|
|
----------
|
|
|
|
function Succ (X : T) return T is
|
|
X_Frac : T;
|
|
X_Exp : UI;
|
|
X1, X2 : T;
|
|
|
|
begin
|
|
-- Treat zero specially since it has a zero exponent
|
|
|
|
if X = 0.0 then
|
|
X1 := 2.0 ** T'Machine_Emin;
|
|
|
|
-- Following loop generates smallest denormal
|
|
|
|
loop
|
|
X2 := T'Machine (X1 / 2.0);
|
|
exit when X2 = 0.0;
|
|
X1 := X2;
|
|
end loop;
|
|
|
|
return X1;
|
|
|
|
-- Special treatment for largest positive number
|
|
|
|
elsif X = T'Last then
|
|
|
|
-- If not generating infinities, we raise a constraint error
|
|
|
|
if T'Machine_Overflows then
|
|
raise Constraint_Error with "Succ of largest negative number";
|
|
|
|
-- Otherwise generate a positive infinity
|
|
|
|
else
|
|
return X / (X - X);
|
|
end if;
|
|
|
|
-- For infinities, return unchanged
|
|
|
|
elsif X < T'First or else X > T'Last then
|
|
return X;
|
|
|
|
-- Add to 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.
|
|
|
|
else
|
|
Decompose (X, X_Frac, X_Exp);
|
|
|
|
-- A special case, if the number we had was a negative power of two,
|
|
-- then we want to add half of what we would otherwise add, since the
|
|
-- exponent is going to be reduced.
|
|
|
|
-- Note that X_Frac has the same sign as X, so if X_Frac is -0.5,
|
|
-- then we know that we have a negative number (and hence a negative
|
|
-- power of 2).
|
|
|
|
if X_Frac = -0.5 then
|
|
return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
|
|
|
|
-- Otherwise the exponent is unchanged
|
|
|
|
else
|
|
return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa);
|
|
end if;
|
|
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
|
|
-- 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
|
|
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
|
|
IEEE_Emin : constant Integer := T'Machine_Emin - 1;
|
|
IEEE_Emax : constant Integer := T'Machine_Emax - 1;
|
|
|
|
IEEE_Bias : constant Integer := -(IEEE_Emin - 1);
|
|
|
|
subtype IEEE_Exponent_Range is
|
|
Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
|
|
|
|
-- The implementation of this floating point attribute uses a
|
|
-- representation type Float_Rep that allows direct access to the
|
|
-- exponent and mantissa parts of a floating point number.
|
|
|
|
-- The Float_Rep type is an 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 the Float_Word'Size.
|
|
|
|
-- The following conditions must be met for all possible instantiations
|
|
-- of the attributes 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. Unused bits (if any) are in
|
|
-- the least significant part.
|
|
|
|
type Float_Word is mod 2**Positive'Min (System.Word_Size, 32);
|
|
type Rep_Index is range 0 .. 7;
|
|
|
|
Rep_Words : constant Positive :=
|
|
(T'Size + Float_Word'Size - 1) / Float_Word'Size;
|
|
Rep_Last : constant Rep_Index :=
|
|
Rep_Index'Min
|
|
(Rep_Index (Rep_Words - 1),
|
|
(T'Mantissa + 16) / Float_Word'Size);
|
|
-- 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 T'Mantissa + 16.
|
|
|
|
type Float_Rep is
|
|
array (Rep_Index range 0 .. Rep_Index (Rep_Words - 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. This is not just a matter of efficiency, but of functionality,
|
|
-- since Valid has a pragma Inline_Always, which is not permitted if
|
|
-- there are nested subprograms present.
|
|
|
|
Most_Significant_Word : constant Rep_Index :=
|
|
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.
|
|
|
|
Exponent_Factor : constant Float_Word :=
|
|
2**(Float_Word'Size - 1) /
|
|
Float_Word (IEEE_Emax - IEEE_Emin + 3) *
|
|
Boolean'Pos (Most_Significant_Word /= 2) +
|
|
Boolean'Pos (Most_Significant_Word = 2);
|
|
-- Factor that the extracted exponent needs to be divided by to be in
|
|
-- range 0 .. IEEE_Emax - IEEE_Emin + 2. Special case: Exponent_Factor
|
|
-- is 1 for x86/IA64 double extended (GCC adds unused bits to the type).
|
|
|
|
Exponent_Mask : constant Float_Word :=
|
|
Float_Word (IEEE_Emax - IEEE_Emin + 2) *
|
|
Exponent_Factor;
|
|
-- Value needed to mask out the exponent field. This assumes that the
|
|
-- range IEEE_Emin - 1 .. IEEE_Emax + contains 2**N values, for some N
|
|
-- in Natural.
|
|
|
|
function To_Float is new Ada.Unchecked_Conversion (Float_Rep, T);
|
|
|
|
type Float_Access is access all T;
|
|
function To_Address is
|
|
new Ada.Unchecked_Conversion (Float_Access, System.Address);
|
|
|
|
XA : constant System.Address := To_Address (Float_Access (X));
|
|
|
|
R : Float_Rep;
|
|
pragma Import (Ada, R);
|
|
for R'Address use XA;
|
|
-- R is a view of the input floating-point parameter. Note that we
|
|
-- must avoid copying the actual bits of this parameter in float
|
|
-- form (since it may be a signalling NaN).
|
|
|
|
E : constant IEEE_Exponent_Range :=
|
|
Integer ((R (Most_Significant_Word) and Exponent_Mask) /
|
|
Exponent_Factor)
|
|
- IEEE_Bias;
|
|
-- Mask/Shift T to only get bits from the exponent. Then convert biased
|
|
-- value to integer value.
|
|
|
|
SR : Float_Rep;
|
|
-- Float_Rep representation of significant of X.all
|
|
|
|
begin
|
|
if T'Denorm then
|
|
|
|
-- All denormalized numbers are valid, so the only invalid numbers
|
|
-- are overflows and NaNs, both with exponent = Emax + 1.
|
|
|
|
return E /= IEEE_Emax + 1;
|
|
|
|
end if;
|
|
|
|
-- All denormalized numbers except 0.0 are invalid
|
|
|
|
-- Set exponent of X to zero, so we end up with the significand, which
|
|
-- definitely is a valid number and can be converted back to a float.
|
|
|
|
SR := R;
|
|
SR (Most_Significant_Word) :=
|
|
(SR (Most_Significant_Word)
|
|
and not Exponent_Mask) + Float_Word (IEEE_Bias) * Exponent_Factor;
|
|
|
|
return (E in IEEE_Emin .. IEEE_Emax) or else
|
|
((E = IEEE_Emin - 1) and then abs To_Float (SR) = 1.0);
|
|
end Valid;
|
|
|
|
end System.Fat_Gen;
|