with Ada.Unchecked_Conversion;
with System;
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);
Rad : constant T := T (T'Machine_Radix);
Invrad : constant T := 1.0 / Rad;
subtype Expbits is Integer range 0 .. 6;
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);
procedure Decompose (XX : T; Frac : out T; Expo : out UI);
function Gradual_Scaling (Adjustment : UI) return T;
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;
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;
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;
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;
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;
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;
else
declare
Ax : T := abs X;
Ex : UI := 0;
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;
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;
end loop;
Ax := Ax * Invrad;
Ex := Ex + 1;
else
while Ax < R_Neg_Power (Expbits'Last) loop
Ax := Ax * R_Power (Expbits'Last);
Ex := Ex - Log_Power (Expbits'Last);
end loop;
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;
end loop;
end if;
if X > 0.0 then
Frac := Ax;
else
Frac := -Ax;
end if;
Expo := Ex;
end;
end if;
end Decompose;
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;
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;
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;
function Gradual_Scaling (Adjustment : UI) return T is
Y : T;
Y1 : T;
Ex : UI := Adjustment;
begin
if Adjustment < T'Machine_Emin 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;
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;
function Machine (X : T) return T is
Temp : T;
pragma Volatile (Temp);
begin
Temp := X;
return Temp;
end Machine;
function Model (X : T) return T is
begin
return Machine (X);
end Model;
function Pred (X : T) return T is
X_Frac : T;
X_Exp : UI;
begin
if X = 0.0 then
return -Succ (X);
else
Decompose (X, X_Frac, X_Exp);
if X_Frac = 0.5 and then X > 0.0 then
return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
else
return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa);
end if;
end if;
end Pred;
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;
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;
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;
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;
else
return X;
end if;
end Rounding;
function Scaling (X : T; Adjustment : UI) return T is
begin
if X = 0.0 or else Adjustment = 0 then
return X;
end if;
declare
Y : T := X;
Ex : UI := Adjustment;
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;
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;
end loop;
else
while Ex >= Log_Power (Expbits'Last) loop
Y := Y * R_Power (Expbits'Last);
Ex := Ex - Log_Power (Expbits'Last);
end loop;
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;
end loop;
end if;
return Y;
end;
end Scaling;
function Succ (X : T) return T is
X_Frac : T;
X_Exp : UI;
X1, X2 : T;
begin
if X = 0.0 then
X1 := 2.0 ** T'Machine_Emin;
loop
X2 := T'Machine (X1 / 2.0);
exit when X2 = 0.0;
X1 := X2;
end loop;
return X1;
else
Decompose (X, X_Frac, X_Exp);
if X_Frac = 0.5 and then X < 0.0 then
return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
else
return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa);
end if;
end if;
end Succ;
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;
else
return X;
end if;
end if;
end Truncation;
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;
else
return X;
end if;
end Unbiased_Rounding;
function Valid (X : 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;
type Float_Word is mod 2**32;
type Rep_Index is range 0 .. 7;
Rep_Last : constant Rep_Index := (T'Size - 1) / Float_Word'Size;
type Float_Rep is array (Rep_Index range 0 .. Rep_Last) of Float_Word;
Most_Significant_Word : constant Rep_Index :=
Rep_Last * Standard'Default_Bit_Order;
Exponent_Factor : constant Float_Word :=
2**(Float_Word'Size - 1) /
Float_Word (IEEE_Emax - IEEE_Emin + 3) *
Boolean'Pos (T'Size /= 96) +
Boolean'Pos (T'Size = 96);
Exponent_Mask : constant Float_Word :=
Float_Word (IEEE_Emax - IEEE_Emin + 2) *
Exponent_Factor;
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;
E : constant IEEE_Exponent_Range :=
Integer ((R (Most_Significant_Word) and Exponent_Mask) /
Exponent_Factor)
- IEEE_Bias;
SR : Float_Rep;
begin
if T'Denorm then
return E /= IEEE_Emax + 1;
end if;
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;