1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
10 -- Copyright (C) 1992-2001 Free Software Foundation, Inc. --
12 -- GNAT is free software; you can redistribute it and/or modify it under --
13 -- terms of the GNU General Public License as published by the Free Soft- --
14 -- ware Foundation; either version 2, or (at your option) any later ver- --
15 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
16 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
17 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
18 -- for more details. You should have received a copy of the GNU General --
19 -- Public License distributed with GNAT; see file COPYING. If not, write --
20 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
21 -- MA 02111-1307, USA. --
23 -- GNAT was originally developed by the GNAT team at New York University. --
24 -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
26 ------------------------------------------------------------------------------
28 with Einfo; use Einfo;
29 with Sem_Util; use Sem_Util;
30 with Ttypef; use Ttypef;
31 with Targparm; use Targparm;
33 package body Eval_Fat is
35 Radix : constant Int := 2;
36 -- This code is currently only correct for the radix 2 case. We use
37 -- the symbolic value Radix where possible to help in the unlikely
38 -- case of anyone ever having to adjust this code for another value,
39 -- and for documentation purposes.
41 type Radix_Power_Table is array (Int range 1 .. 4) of Int;
43 Radix_Powers : constant Radix_Power_Table
44 := (Radix**1, Radix**2, Radix**3, Radix**4);
46 function Float_Radix return T renames Ureal_2;
47 -- Radix expressed in real form
49 -----------------------
50 -- Local Subprograms --
51 -----------------------
58 Mode : Rounding_Mode := Round);
59 -- Decomposes a non-zero floating-point number into fraction and
60 -- exponent parts. The fraction is in the interval 1.0 / Radix ..
61 -- T'Pred (1.0) and uses Rbase = Radix.
62 -- The result is rounded to a nearest machine number.
64 procedure Decompose_Int
69 Mode : Rounding_Mode);
70 -- This is similar to Decompose, except that the Fraction value returned
71 -- is an integer representing the value Fraction * Scale, where Scale is
72 -- the value (Radix ** Machine_Mantissa (RT)). The value is obtained by
73 -- using biased rounding (halfway cases round away from zero), round to
74 -- even, a floor operation or a ceiling operation depending on the setting
75 -- of Mode (see corresponding descriptions in Urealp).
76 -- In case rounding was specified, Rounding_Was_Biased is set True
77 -- if the input was indeed halfway between to machine numbers and
78 -- got rounded away from zero to an odd number.
80 function Eps_Model (RT : R) return T;
81 -- Return the smallest model number of R.
83 function Eps_Denorm (RT : R) return T;
84 -- Return the smallest denormal of type R.
86 function Machine_Mantissa (RT : R) return Nat;
87 -- Get value of machine mantissa
93 function Adjacent (RT : R; X, Towards : T) return T is
98 elsif Towards > X then
110 function Ceiling (RT : R; X : T) return T is
111 XT : constant T := Truncation (RT, X);
114 if UR_Is_Negative (X) then
129 function Compose (RT : R; Fraction : T; Exponent : UI) return T is
134 if UR_Is_Zero (Fraction) then
137 Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
138 return Scaling (RT, Arg_Frac, Exponent);
146 function Copy_Sign (RT : R; Value, Sign : T) return T is
147 pragma Warnings (Off, RT);
153 if UR_Is_Negative (Sign) then
169 Mode : Rounding_Mode := Round)
174 Decompose_Int (RT, abs X, Int_F, Exponent, Mode);
176 Fraction := UR_From_Components
178 Den => UI_From_Int (Machine_Mantissa (RT)),
182 if UR_Is_Negative (X) then
183 Fraction := -Fraction;
193 -- This procedure should be modified with care, as there
194 -- are many non-obvious details that may cause problems
195 -- that are hard to detect. The cases of positive and
196 -- negative zeroes are also special and should be
197 -- verified separately.
199 procedure Decompose_Int
204 Mode : Rounding_Mode)
206 Base : Int := Rbase (X);
207 N : UI := abs Numerator (X);
208 D : UI := Denominator (X);
213 -- True iff Fraction is even
215 Most_Significant_Digit : constant UI :=
216 Radix ** (Machine_Mantissa (RT) - 1);
218 Uintp_Mark : Uintp.Save_Mark;
219 -- The code is divided into blocks that systematically release
220 -- intermediate values (this routine generates lots of junk!)
223 Calculate_D_And_Exponent_1 : begin
227 -- In cases where Base > 1, the actual denominator is
228 -- Base**D. For cases where Base is a power of Radix, use
229 -- the value 1 for the Denominator and adjust the exponent.
231 -- Note: Exponent has different sign from D, because D is a divisor
233 for Power in 1 .. Radix_Powers'Last loop
234 if Base = Radix_Powers (Power) then
235 Exponent := -D * Power;
242 Release_And_Save (Uintp_Mark, D, Exponent);
243 end Calculate_D_And_Exponent_1;
246 Calculate_Exponent : begin
249 -- For bases that are a multiple of the Radix, divide
250 -- the base by Radix and adjust the Exponent. This will
251 -- help because D will be much smaller and faster to process.
253 -- This occurs for decimal bases on a machine with binary
254 -- floating-point for example. When calculating 1E40,
255 -- with Radix = 2, N will be 93 bits instead of 133.
263 -- = -------------------------- * Radix
265 -- (Base/Radix) * Radix
268 -- = --------------- * Radix
272 -- This code is commented out, because it causes numerous
273 -- failures in the regression suite. To be studied ???
275 while False and then Base > 0 and then Base mod Radix = 0 loop
276 Base := Base / Radix;
277 Exponent := Exponent + D;
280 Release_And_Save (Uintp_Mark, Exponent);
281 end Calculate_Exponent;
283 -- For remaining bases we must actually compute
284 -- the exponentiation.
286 -- Because the exponentiation can be negative, and D must
287 -- be integer, the numerator is corrected instead.
289 Calculate_N_And_D : begin
293 N := N * Base ** (-D);
299 Release_And_Save (Uintp_Mark, N, D);
300 end Calculate_N_And_D;
305 -- Now scale N and D so that N / D is a value in the
306 -- interval [1.0 / Radix, 1.0) and adjust Exponent accordingly,
307 -- so the value N / D * Radix ** Exponent remains unchanged.
309 -- Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0
311 -- N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.
312 -- This scaling is not possible for N is Uint_0 as there
313 -- is no way to scale Uint_0 so the first digit is non-zero.
315 Calculate_N_And_Exponent : begin
318 N_Times_Radix := N * Radix;
321 while not (N_Times_Radix >= D) loop
323 Exponent := Exponent - 1;
325 N_Times_Radix := N * Radix;
329 Release_And_Save (Uintp_Mark, N, Exponent);
330 end Calculate_N_And_Exponent;
332 -- Step 2 - Adjust D so N / D < 1
334 -- Scale up D so N / D < 1, so N < D
336 Calculate_D_And_Exponent_2 : begin
339 while not (N < D) loop
341 -- As N / D >= 1, N / (D * Radix) will be at least 1 / Radix,
342 -- so the result of Step 1 stays valid
345 Exponent := Exponent + 1;
348 Release_And_Save (Uintp_Mark, D, Exponent);
349 end Calculate_D_And_Exponent_2;
351 -- Here the value N / D is in the range [1.0 / Radix .. 1.0)
353 -- Now find the fraction by doing a very simple-minded
354 -- division until enough digits have been computed.
356 -- This division works for all radices, but is only efficient for
357 -- a binary radix. It is just like a manual division algorithm,
358 -- but instead of moving the denominator one digit right, we move
359 -- the numerator one digit left so the numerator and denominator
365 Calculate_Fraction_And_N : begin
371 Fraction := Fraction + 1;
375 -- Stop when the result is in [1.0 / Radix, 1.0)
377 exit when Fraction >= Most_Significant_Digit;
380 Fraction := Fraction * Radix;
384 Release_And_Save (Uintp_Mark, Fraction, N);
385 end Calculate_Fraction_And_N;
387 Calculate_Fraction_And_Exponent : begin
390 -- Put back sign before applying the rounding.
392 if UR_Is_Negative (X) then
393 Fraction := -Fraction;
396 -- Determine correct rounding based on the remainder
397 -- which is in N and the divisor D.
399 Rounding_Was_Biased := False; -- Until proven otherwise
404 -- This rounding mode should not be used for static
405 -- expressions, but only for compile-time evaluation
406 -- of non-static expressions.
408 if (Even and then N * 2 > D)
410 (not Even and then N * 2 >= D)
412 Fraction := Fraction + 1;
417 -- Do not round to even as is done with IEEE arithmetic,
418 -- but instead round away from zero when the result is
419 -- exactly between two machine numbers. See RM 4.9(38).
422 Fraction := Fraction + 1;
424 Rounding_Was_Biased := Even and then N * 2 = D;
425 -- Check for the case where the result is actually
426 -- different from Round_Even.
431 Fraction := Fraction + 1;
437 -- The result must be normalized to [1.0/Radix, 1.0),
438 -- so adjust if the result is 1.0 because of rounding.
440 if Fraction = Most_Significant_Digit * Radix then
441 Fraction := Most_Significant_Digit;
442 Exponent := Exponent + 1;
445 Release_And_Save (Uintp_Mark, Fraction, Exponent);
446 end Calculate_Fraction_And_Exponent;
454 function Eps_Denorm (RT : R) return T is
455 Digs : constant UI := Digits_Value (RT);
460 if Vax_Float (RT) then
461 if Digs = VAXFF_Digits then
462 Emin := VAXFF_Machine_Emin;
463 Mant := VAXFF_Machine_Mantissa;
465 elsif Digs = VAXDF_Digits then
466 Emin := VAXDF_Machine_Emin;
467 Mant := VAXDF_Machine_Mantissa;
470 pragma Assert (Digs = VAXGF_Digits);
471 Emin := VAXGF_Machine_Emin;
472 Mant := VAXGF_Machine_Mantissa;
475 elsif Is_AAMP_Float (RT) then
476 if Digs = AAMPS_Digits then
477 Emin := AAMPS_Machine_Emin;
478 Mant := AAMPS_Machine_Mantissa;
481 pragma Assert (Digs = AAMPL_Digits);
482 Emin := AAMPL_Machine_Emin;
483 Mant := AAMPL_Machine_Mantissa;
487 if Digs = IEEES_Digits then
488 Emin := IEEES_Machine_Emin;
489 Mant := IEEES_Machine_Mantissa;
491 elsif Digs = IEEEL_Digits then
492 Emin := IEEEL_Machine_Emin;
493 Mant := IEEEL_Machine_Mantissa;
496 pragma Assert (Digs = IEEEX_Digits);
497 Emin := IEEEX_Machine_Emin;
498 Mant := IEEEX_Machine_Mantissa;
502 return Float_Radix ** UI_From_Int (Emin - Mant);
509 function Eps_Model (RT : R) return T is
510 Digs : constant UI := Digits_Value (RT);
514 if Vax_Float (RT) then
515 if Digs = VAXFF_Digits then
516 Emin := VAXFF_Machine_Emin;
518 elsif Digs = VAXDF_Digits then
519 Emin := VAXDF_Machine_Emin;
522 pragma Assert (Digs = VAXGF_Digits);
523 Emin := VAXGF_Machine_Emin;
526 elsif Is_AAMP_Float (RT) then
527 if Digs = AAMPS_Digits then
528 Emin := AAMPS_Machine_Emin;
531 pragma Assert (Digs = AAMPL_Digits);
532 Emin := AAMPL_Machine_Emin;
536 if Digs = IEEES_Digits then
537 Emin := IEEES_Machine_Emin;
539 elsif Digs = IEEEL_Digits then
540 Emin := IEEEL_Machine_Emin;
543 pragma Assert (Digs = IEEEX_Digits);
544 Emin := IEEEX_Machine_Emin;
548 return Float_Radix ** UI_From_Int (Emin);
555 function Exponent (RT : R; X : T) return UI is
560 if UR_Is_Zero (X) then
563 Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);
572 function Floor (RT : R; X : T) return T is
573 XT : constant T := Truncation (RT, X);
576 if UR_Is_Positive (X) then
591 function Fraction (RT : R; X : T) return T is
596 if UR_Is_Zero (X) then
599 Decompose (RT, X, X_Frac, X_Exp);
608 function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
613 if Radix_Digits >= Machine_Mantissa (RT) then
617 L := Exponent (RT, X) - Radix_Digits;
618 Y := Truncation (RT, Scaling (RT, X, -L));
619 Z := Scaling (RT, Y, L);
629 function Machine (RT : R; X : T; Mode : Rounding_Mode) return T is
634 if UR_Is_Zero (X) then
637 Decompose (RT, X, X_Frac, X_Exp, Mode);
638 return Scaling (RT, X_Frac, X_Exp);
642 ----------------------
643 -- Machine_Mantissa --
644 ----------------------
646 function Machine_Mantissa (RT : R) return Nat is
647 Digs : constant UI := Digits_Value (RT);
651 if Vax_Float (RT) then
652 if Digs = VAXFF_Digits then
653 Mant := VAXFF_Machine_Mantissa;
655 elsif Digs = VAXDF_Digits then
656 Mant := VAXDF_Machine_Mantissa;
659 pragma Assert (Digs = VAXGF_Digits);
660 Mant := VAXGF_Machine_Mantissa;
663 elsif Is_AAMP_Float (RT) then
664 if Digs = AAMPS_Digits then
665 Mant := AAMPS_Machine_Mantissa;
668 pragma Assert (Digs = AAMPL_Digits);
669 Mant := AAMPL_Machine_Mantissa;
673 if Digs = IEEES_Digits then
674 Mant := IEEES_Machine_Mantissa;
676 elsif Digs = IEEEL_Digits then
677 Mant := IEEEL_Machine_Mantissa;
680 pragma Assert (Digs = IEEEX_Digits);
681 Mant := IEEEX_Machine_Mantissa;
686 end Machine_Mantissa;
692 function Model (RT : R; X : T) return T is
697 Decompose (RT, X, X_Frac, X_Exp);
698 return Compose (RT, X_Frac, X_Exp);
705 function Pred (RT : R; X : T) return T is
710 if abs X < Eps_Model (RT) then
711 if Denorm_On_Target then
712 return X - Eps_Denorm (RT);
714 elsif X > Ureal_0 then
715 -- Target does not support denorms, so predecessor is 0.0
719 -- Target does not support denorms, and X is 0.0
720 -- or at least bigger than -Eps_Model (RT)
722 return -Eps_Model (RT);
726 Decompose_Int (RT, X, Result_F, Result_X, Ceiling);
727 return UR_From_Components
728 (Num => Result_F - 1,
729 Den => Machine_Mantissa (RT) - Result_X,
732 -- Result_F may be false, but this is OK as UR_From_Components
733 -- handles that situation.
741 function Remainder (RT : R; X, Y : T) return T is
756 if UR_Is_Positive (X) then
768 P_Exp := Exponent (RT, P);
771 -- ??? what about zero cases?
772 Decompose (RT, Arg, Arg_Frac, Arg_Exp);
773 Decompose (RT, P, P_Frac, P_Exp);
775 P := Compose (RT, P_Frac, Arg_Exp);
776 K := Arg_Exp - P_Exp;
780 for Cnt in reverse 0 .. UI_To_Int (K) loop
781 if IEEE_Rem >= P then
783 IEEE_Rem := IEEE_Rem - P;
792 -- That completes the calculation of modulus remainder. The final step
793 -- is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
797 B := abs Y * Ureal_Half;
800 A := IEEE_Rem * Ureal_2;
804 if A > B or else (A = B and then not P_Even) then
805 IEEE_Rem := IEEE_Rem - abs Y;
808 return Sign_X * IEEE_Rem;
816 function Rounding (RT : R; X : T) return T is
821 Result := Truncation (RT, abs X);
822 Tail := abs X - Result;
824 if Tail >= Ureal_Half then
825 Result := Result + Ureal_1;
828 if UR_Is_Negative (X) then
840 function Scaling (RT : R; X : T; Adjustment : UI) return T is
841 pragma Warnings (Off, RT);
844 if Rbase (X) = Radix then
845 return UR_From_Components
846 (Num => Numerator (X),
847 Den => Denominator (X) - Adjustment,
849 Negative => UR_Is_Negative (X));
851 elsif Adjustment >= 0 then
852 return X * Radix ** Adjustment;
854 return X / Radix ** (-Adjustment);
862 function Succ (RT : R; X : T) return T is
867 if abs X < Eps_Model (RT) then
868 if Denorm_On_Target then
869 return X + Eps_Denorm (RT);
871 elsif X < Ureal_0 then
872 -- Target does not support denorms, so successor is 0.0
876 -- Target does not support denorms, and X is 0.0
877 -- or at least smaller than Eps_Model (RT)
879 return Eps_Model (RT);
883 Decompose_Int (RT, X, Result_F, Result_X, Floor);
884 return UR_From_Components
885 (Num => Result_F + 1,
886 Den => Machine_Mantissa (RT) - Result_X,
889 -- Result_F may be false, but this is OK as UR_From_Components
890 -- handles that situation.
898 function Truncation (RT : R; X : T) return T is
899 pragma Warnings (Off, RT);
902 return UR_From_Uint (UR_Trunc (X));
905 -----------------------
906 -- Unbiased_Rounding --
907 -----------------------
909 function Unbiased_Rounding (RT : R; X : T) return T is
910 Abs_X : constant T := abs X;
915 Result := Truncation (RT, Abs_X);
916 Tail := Abs_X - Result;
918 if Tail > Ureal_Half then
919 Result := Result + Ureal_1;
921 elsif Tail = Ureal_Half then
923 Truncation (RT, (Result / Ureal_2) + Ureal_Half);
926 if UR_Is_Negative (X) then
928 elsif UR_Is_Positive (X) then
931 -- For zero case, make sure sign of zero is preserved
937 end Unbiased_Rounding;