with Einfo; use Einfo;
with Sem_Util; use Sem_Util;
with Ttypef; use Ttypef;
with Targparm; use Targparm;
package body Eval_Fat is
Radix : constant Int := 2;
type Radix_Power_Table is array (Int range 1 .. 4) of Int;
Radix_Powers : constant Radix_Power_Table
:= (Radix**1, Radix**2, Radix**3, Radix**4);
function Float_Radix return T renames Ureal_2;
procedure Decompose
(RT : R;
X : in T;
Fraction : out T;
Exponent : out UI;
Mode : Rounding_Mode := Round);
procedure Decompose_Int
(RT : R;
X : in T;
Fraction : out UI;
Exponent : out UI;
Mode : Rounding_Mode);
function Eps_Model (RT : R) return T;
function Eps_Denorm (RT : R) return T;
function Machine_Mantissa (RT : R) return Nat;
function Adjacent (RT : R; X, Towards : T) return T is
begin
if Towards = X then
return X;
elsif Towards > X then
return Succ (RT, X);
else
return Pred (RT, X);
end if;
end Adjacent;
function Ceiling (RT : R; X : T) return T is
XT : constant T := Truncation (RT, X);
begin
if UR_Is_Negative (X) then
return XT;
elsif X = XT then
return X;
else
return XT + Ureal_1;
end if;
end Ceiling;
function Compose (RT : R; Fraction : T; Exponent : UI) return T is
Arg_Frac : T;
Arg_Exp : UI;
begin
if UR_Is_Zero (Fraction) then
return Fraction;
else
Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
return Scaling (RT, Arg_Frac, Exponent);
end if;
end Compose;
function Copy_Sign (RT : R; Value, Sign : T) return T is
pragma Warnings (Off, RT);
Result : T;
begin
Result := abs Value;
if UR_Is_Negative (Sign) then
return -Result;
else
return Result;
end if;
end Copy_Sign;
procedure Decompose
(RT : R;
X : in T;
Fraction : out T;
Exponent : out UI;
Mode : Rounding_Mode := Round)
is
Int_F : UI;
begin
Decompose_Int (RT, abs X, Int_F, Exponent, Mode);
Fraction := UR_From_Components
(Num => Int_F,
Den => UI_From_Int (Machine_Mantissa (RT)),
Rbase => Radix,
Negative => False);
if UR_Is_Negative (X) then
Fraction := -Fraction;
end if;
return;
end Decompose;
procedure Decompose_Int
(RT : R;
X : in T;
Fraction : out UI;
Exponent : out UI;
Mode : Rounding_Mode)
is
Base : Int := Rbase (X);
N : UI := abs Numerator (X);
D : UI := Denominator (X);
N_Times_Radix : UI;
Even : Boolean;
Most_Significant_Digit : constant UI :=
Radix ** (Machine_Mantissa (RT) - 1);
Uintp_Mark : Uintp.Save_Mark;
begin
Calculate_D_And_Exponent_1 : begin
Uintp_Mark := Mark;
Exponent := Uint_0;
for Power in 1 .. Radix_Powers'Last loop
if Base = Radix_Powers (Power) then
Exponent := -D * Power;
Base := 0;
D := Uint_1;
exit;
end if;
end loop;
Release_And_Save (Uintp_Mark, D, Exponent);
end Calculate_D_And_Exponent_1;
if Base > 0 then
Calculate_Exponent : begin
Uintp_Mark := Mark;
while False and then Base > 0 and then Base mod Radix = 0 loop
Base := Base / Radix;
Exponent := Exponent + D;
end loop;
Release_And_Save (Uintp_Mark, Exponent);
end Calculate_Exponent;
Calculate_N_And_D : begin
Uintp_Mark := Mark;
if D < 0 then
N := N * Base ** (-D);
D := Uint_1;
else
D := Base ** D;
end if;
Release_And_Save (Uintp_Mark, N, D);
end Calculate_N_And_D;
Base := 0;
end if;
Calculate_N_And_Exponent : begin
Uintp_Mark := Mark;
N_Times_Radix := N * Radix;
if N /= Uint_0 then
while not (N_Times_Radix >= D) loop
N := N_Times_Radix;
Exponent := Exponent - 1;
N_Times_Radix := N * Radix;
end loop;
end if;
Release_And_Save (Uintp_Mark, N, Exponent);
end Calculate_N_And_Exponent;
Calculate_D_And_Exponent_2 : begin
Uintp_Mark := Mark;
while not (N < D) loop
D := D * Radix;
Exponent := Exponent + 1;
end loop;
Release_And_Save (Uintp_Mark, D, Exponent);
end Calculate_D_And_Exponent_2;
Fraction := Uint_0;
Even := True;
Calculate_Fraction_And_N : begin
Uintp_Mark := Mark;
loop
while N >= D loop
N := N - D;
Fraction := Fraction + 1;
Even := not Even;
end loop;
exit when Fraction >= Most_Significant_Digit;
N := N * Radix;
Fraction := Fraction * Radix;
Even := True;
end loop;
Release_And_Save (Uintp_Mark, Fraction, N);
end Calculate_Fraction_And_N;
Calculate_Fraction_And_Exponent : begin
Uintp_Mark := Mark;
if UR_Is_Negative (X) then
Fraction := -Fraction;
end if;
Rounding_Was_Biased := False;
case Mode is
when Round_Even =>
if (Even and then N * 2 > D)
or else
(not Even and then N * 2 >= D)
then
Fraction := Fraction + 1;
end if;
when Round =>
if N * 2 >= D then
Fraction := Fraction + 1;
Rounding_Was_Biased := Even and then N * 2 = D;
end if;
when Ceiling =>
if N > Uint_0 then
Fraction := Fraction + 1;
end if;
when Floor => null;
end case;
if Fraction = Most_Significant_Digit * Radix then
Fraction := Most_Significant_Digit;
Exponent := Exponent + 1;
end if;
Release_And_Save (Uintp_Mark, Fraction, Exponent);
end Calculate_Fraction_And_Exponent;
end Decompose_Int;
function Eps_Denorm (RT : R) return T is
Digs : constant UI := Digits_Value (RT);
Emin : Int;
Mant : Int;
begin
if Vax_Float (RT) then
if Digs = VAXFF_Digits then
Emin := VAXFF_Machine_Emin;
Mant := VAXFF_Machine_Mantissa;
elsif Digs = VAXDF_Digits then
Emin := VAXDF_Machine_Emin;
Mant := VAXDF_Machine_Mantissa;
else
pragma Assert (Digs = VAXGF_Digits);
Emin := VAXGF_Machine_Emin;
Mant := VAXGF_Machine_Mantissa;
end if;
elsif Is_AAMP_Float (RT) then
if Digs = AAMPS_Digits then
Emin := AAMPS_Machine_Emin;
Mant := AAMPS_Machine_Mantissa;
else
pragma Assert (Digs = AAMPL_Digits);
Emin := AAMPL_Machine_Emin;
Mant := AAMPL_Machine_Mantissa;
end if;
else
if Digs = IEEES_Digits then
Emin := IEEES_Machine_Emin;
Mant := IEEES_Machine_Mantissa;
elsif Digs = IEEEL_Digits then
Emin := IEEEL_Machine_Emin;
Mant := IEEEL_Machine_Mantissa;
else
pragma Assert (Digs = IEEEX_Digits);
Emin := IEEEX_Machine_Emin;
Mant := IEEEX_Machine_Mantissa;
end if;
end if;
return Float_Radix ** UI_From_Int (Emin - Mant);
end Eps_Denorm;
function Eps_Model (RT : R) return T is
Digs : constant UI := Digits_Value (RT);
Emin : Int;
begin
if Vax_Float (RT) then
if Digs = VAXFF_Digits then
Emin := VAXFF_Machine_Emin;
elsif Digs = VAXDF_Digits then
Emin := VAXDF_Machine_Emin;
else
pragma Assert (Digs = VAXGF_Digits);
Emin := VAXGF_Machine_Emin;
end if;
elsif Is_AAMP_Float (RT) then
if Digs = AAMPS_Digits then
Emin := AAMPS_Machine_Emin;
else
pragma Assert (Digs = AAMPL_Digits);
Emin := AAMPL_Machine_Emin;
end if;
else
if Digs = IEEES_Digits then
Emin := IEEES_Machine_Emin;
elsif Digs = IEEEL_Digits then
Emin := IEEEL_Machine_Emin;
else
pragma Assert (Digs = IEEEX_Digits);
Emin := IEEEX_Machine_Emin;
end if;
end if;
return Float_Radix ** UI_From_Int (Emin);
end Eps_Model;
function Exponent (RT : R; X : T) return UI is
X_Frac : UI;
X_Exp : UI;
begin
if UR_Is_Zero (X) then
return Uint_0;
else
Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);
return X_Exp;
end if;
end Exponent;
function Floor (RT : R; X : T) return T is
XT : constant T := Truncation (RT, X);
begin
if UR_Is_Positive (X) then
return XT;
elsif XT = X then
return X;
else
return XT - Ureal_1;
end if;
end Floor;
function Fraction (RT : R; X : T) return T is
X_Frac : T;
X_Exp : UI;
begin
if UR_Is_Zero (X) then
return X;
else
Decompose (RT, X, X_Frac, X_Exp);
return X_Frac;
end if;
end Fraction;
function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
L : UI;
Y, Z : T;
begin
if Radix_Digits >= Machine_Mantissa (RT) then
return X;
else
L := Exponent (RT, X) - Radix_Digits;
Y := Truncation (RT, Scaling (RT, X, -L));
Z := Scaling (RT, Y, L);
return Z;
end if;
end Leading_Part;
function Machine (RT : R; X : T; Mode : Rounding_Mode) return T is
X_Frac : T;
X_Exp : UI;
begin
if UR_Is_Zero (X) then
return X;
else
Decompose (RT, X, X_Frac, X_Exp, Mode);
return Scaling (RT, X_Frac, X_Exp);
end if;
end Machine;
function Machine_Mantissa (RT : R) return Nat is
Digs : constant UI := Digits_Value (RT);
Mant : Nat;
begin
if Vax_Float (RT) then
if Digs = VAXFF_Digits then
Mant := VAXFF_Machine_Mantissa;
elsif Digs = VAXDF_Digits then
Mant := VAXDF_Machine_Mantissa;
else
pragma Assert (Digs = VAXGF_Digits);
Mant := VAXGF_Machine_Mantissa;
end if;
elsif Is_AAMP_Float (RT) then
if Digs = AAMPS_Digits then
Mant := AAMPS_Machine_Mantissa;
else
pragma Assert (Digs = AAMPL_Digits);
Mant := AAMPL_Machine_Mantissa;
end if;
else
if Digs = IEEES_Digits then
Mant := IEEES_Machine_Mantissa;
elsif Digs = IEEEL_Digits then
Mant := IEEEL_Machine_Mantissa;
else
pragma Assert (Digs = IEEEX_Digits);
Mant := IEEEX_Machine_Mantissa;
end if;
end if;
return Mant;
end Machine_Mantissa;
function Model (RT : R; X : T) return T is
X_Frac : T;
X_Exp : UI;
begin
Decompose (RT, X, X_Frac, X_Exp);
return Compose (RT, X_Frac, X_Exp);
end Model;
function Pred (RT : R; X : T) return T is
Result_F : UI;
Result_X : UI;
begin
if abs X < Eps_Model (RT) then
if Denorm_On_Target then
return X - Eps_Denorm (RT);
elsif X > Ureal_0 then
return Ureal_0;
else
return -Eps_Model (RT);
end if;
else
Decompose_Int (RT, X, Result_F, Result_X, Ceiling);
return UR_From_Components
(Num => Result_F - 1,
Den => Machine_Mantissa (RT) - Result_X,
Rbase => Radix,
Negative => False);
end if;
end Pred;
function Remainder (RT : R; 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 UR_Is_Positive (X) then
Sign_X := Ureal_1;
else
Sign_X := -Ureal_1;
end if;
Arg := abs X;
P := abs Y;
if Arg < P then
P_Even := True;
IEEE_Rem := Arg;
P_Exp := Exponent (RT, P);
else
Decompose (RT, Arg, Arg_Frac, Arg_Exp);
Decompose (RT, P, P_Frac, P_Exp);
P := Compose (RT, P_Frac, Arg_Exp);
K := Arg_Exp - P_Exp;
P_Even := True;
IEEE_Rem := Arg;
for Cnt in reverse 0 .. UI_To_Int (K) loop
if IEEE_Rem >= P then
P_Even := False;
IEEE_Rem := IEEE_Rem - P;
else
P_Even := True;
end if;
P := P * Ureal_Half;
end loop;
end if;
if P_Exp >= 0 then
A := IEEE_Rem;
B := abs Y * Ureal_Half;
else
A := IEEE_Rem * Ureal_2;
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 (RT : R; X : T) return T is
Result : T;
Tail : T;
begin
Result := Truncation (RT, abs X);
Tail := abs X - Result;
if Tail >= Ureal_Half then
Result := Result + Ureal_1;
end if;
if UR_Is_Negative (X) then
return -Result;
else
return Result;
end if;
end Rounding;
function Scaling (RT : R; X : T; Adjustment : UI) return T is
pragma Warnings (Off, RT);
begin
if Rbase (X) = Radix then
return UR_From_Components
(Num => Numerator (X),
Den => Denominator (X) - Adjustment,
Rbase => Radix,
Negative => UR_Is_Negative (X));
elsif Adjustment >= 0 then
return X * Radix ** Adjustment;
else
return X / Radix ** (-Adjustment);
end if;
end Scaling;
function Succ (RT : R; X : T) return T is
Result_F : UI;
Result_X : UI;
begin
if abs X < Eps_Model (RT) then
if Denorm_On_Target then
return X + Eps_Denorm (RT);
elsif X < Ureal_0 then
return Ureal_0;
else
return Eps_Model (RT);
end if;
else
Decompose_Int (RT, X, Result_F, Result_X, Floor);
return UR_From_Components
(Num => Result_F + 1,
Den => Machine_Mantissa (RT) - Result_X,
Rbase => Radix,
Negative => False);
end if;
end Succ;
function Truncation (RT : R; X : T) return T is
pragma Warnings (Off, RT);
begin
return UR_From_Uint (UR_Trunc (X));
end Truncation;
function Unbiased_Rounding (RT : R; X : T) return T is
Abs_X : constant T := abs X;
Result : T;
Tail : T;
begin
Result := Truncation (RT, Abs_X);
Tail := Abs_X - Result;
if Tail > Ureal_Half then
Result := Result + Ureal_1;
elsif Tail = Ureal_Half then
Result := Ureal_2 *
Truncation (RT, (Result / Ureal_2) + Ureal_Half);
end if;
if UR_Is_Negative (X) then
return -Result;
elsif UR_Is_Positive (X) then
return Result;
else
return X;
end if;
end Unbiased_Rounding;
end Eval_Fat;