1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
9 -- Copyright (C) 1992-2004 Free Software Foundation, Inc. --
11 -- GNAT is free software; you can redistribute it and/or modify it under --
12 -- terms of the GNU General Public License as published by the Free Soft- --
13 -- ware Foundation; either version 2, or (at your option) any later ver- --
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
16 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
17 -- for more details. You should have received a copy of the GNU General --
18 -- Public License distributed with GNAT; see file COPYING. If not, write --
19 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
20 -- MA 02111-1307, USA. --
22 -- GNAT was originally developed by the GNAT team at New York University. --
23 -- Extensive contributions were provided by Ada Core Technologies Inc. --
25 ------------------------------------------------------------------------------
27 with Einfo; use Einfo;
28 with Errout; use Errout;
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).
77 function Eps_Model (RT : R) return T;
78 -- Return the smallest model number of R.
80 function Eps_Denorm (RT : R) return T;
81 -- Return the smallest denormal of type R.
83 function Machine_Emin (RT : R) return Int;
84 -- Return value of the Machine_Emin attribute
90 function Adjacent (RT : R; X, Towards : T) return T is
95 elsif Towards > X then
107 function Ceiling (RT : R; X : T) return T is
108 XT : constant T := Truncation (RT, X);
111 if UR_Is_Negative (X) then
126 function Compose (RT : R; Fraction : T; Exponent : UI) return T is
130 if UR_Is_Zero (Fraction) then
133 Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
134 return Scaling (RT, Arg_Frac, Exponent);
142 function Copy_Sign (RT : R; Value, Sign : T) return T is
143 pragma Warnings (Off, RT);
149 if UR_Is_Negative (Sign) then
165 Mode : Rounding_Mode := Round)
170 Decompose_Int (RT, abs X, Int_F, Exponent, Mode);
172 Fraction := UR_From_Components
174 Den => UI_From_Int (Machine_Mantissa (RT)),
178 if UR_Is_Negative (X) then
179 Fraction := -Fraction;
189 -- This procedure should be modified with care, as there are many
190 -- non-obvious details that may cause problems that are hard to
191 -- detect. The cases of positive and negative zeroes are also
192 -- special and should be verified separately.
194 procedure Decompose_Int
199 Mode : Rounding_Mode)
201 Base : Int := Rbase (X);
202 N : UI := abs Numerator (X);
203 D : UI := Denominator (X);
208 -- True iff Fraction is even
210 Most_Significant_Digit : constant UI :=
211 Radix ** (Machine_Mantissa (RT) - 1);
213 Uintp_Mark : Uintp.Save_Mark;
214 -- The code is divided into blocks that systematically release
215 -- intermediate values (this routine generates lots of junk!)
218 Calculate_D_And_Exponent_1 : begin
222 -- In cases where Base > 1, the actual denominator is
223 -- Base**D. For cases where Base is a power of Radix, use
224 -- the value 1 for the Denominator and adjust the exponent.
226 -- Note: Exponent has different sign from D, because D is a divisor
228 for Power in 1 .. Radix_Powers'Last loop
229 if Base = Radix_Powers (Power) then
230 Exponent := -D * Power;
237 Release_And_Save (Uintp_Mark, D, Exponent);
238 end Calculate_D_And_Exponent_1;
241 Calculate_Exponent : begin
244 -- For bases that are a multiple of the Radix, divide
245 -- the base by Radix and adjust the Exponent. This will
246 -- help because D will be much smaller and faster to process.
248 -- This occurs for decimal bases on a machine with binary
249 -- floating-point for example. When calculating 1E40,
250 -- with Radix = 2, N will be 93 bits instead of 133.
258 -- = -------------------------- * Radix
260 -- (Base/Radix) * Radix
263 -- = --------------- * Radix
267 -- This code is commented out, because it causes numerous
268 -- failures in the regression suite. To be studied ???
270 while False and then Base > 0 and then Base mod Radix = 0 loop
271 Base := Base / Radix;
272 Exponent := Exponent + D;
275 Release_And_Save (Uintp_Mark, Exponent);
276 end Calculate_Exponent;
278 -- For remaining bases we must actually compute
279 -- the exponentiation.
281 -- Because the exponentiation can be negative, and D must
282 -- be integer, the numerator is corrected instead.
284 Calculate_N_And_D : begin
288 N := N * Base ** (-D);
294 Release_And_Save (Uintp_Mark, N, D);
295 end Calculate_N_And_D;
300 -- Now scale N and D so that N / D is a value in the
301 -- interval [1.0 / Radix, 1.0) and adjust Exponent accordingly,
302 -- so the value N / D * Radix ** Exponent remains unchanged.
304 -- Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0
306 -- N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.
307 -- This scaling is not possible for N is Uint_0 as there
308 -- is no way to scale Uint_0 so the first digit is non-zero.
310 Calculate_N_And_Exponent : begin
313 N_Times_Radix := N * Radix;
316 while not (N_Times_Radix >= D) loop
318 Exponent := Exponent - 1;
320 N_Times_Radix := N * Radix;
324 Release_And_Save (Uintp_Mark, N, Exponent);
325 end Calculate_N_And_Exponent;
327 -- Step 2 - Adjust D so N / D < 1
329 -- Scale up D so N / D < 1, so N < D
331 Calculate_D_And_Exponent_2 : begin
334 while not (N < D) loop
336 -- As N / D >= 1, N / (D * Radix) will be at least 1 / Radix,
337 -- so the result of Step 1 stays valid
340 Exponent := Exponent + 1;
343 Release_And_Save (Uintp_Mark, D, Exponent);
344 end Calculate_D_And_Exponent_2;
346 -- Here the value N / D is in the range [1.0 / Radix .. 1.0)
348 -- Now find the fraction by doing a very simple-minded
349 -- division until enough digits have been computed.
351 -- This division works for all radices, but is only efficient for
352 -- a binary radix. It is just like a manual division algorithm,
353 -- but instead of moving the denominator one digit right, we move
354 -- the numerator one digit left so the numerator and denominator
360 Calculate_Fraction_And_N : begin
366 Fraction := Fraction + 1;
370 -- Stop when the result is in [1.0 / Radix, 1.0)
372 exit when Fraction >= Most_Significant_Digit;
375 Fraction := Fraction * Radix;
379 Release_And_Save (Uintp_Mark, Fraction, N);
380 end Calculate_Fraction_And_N;
382 Calculate_Fraction_And_Exponent : begin
385 -- Determine correct rounding based on the remainder
386 -- which is in N and the divisor D. The rounding is
387 -- performed on the absolute value of X, so Ceiling
388 -- and Floor need to check for the sign of X explicitly.
393 -- This rounding mode should not be used for static
394 -- expressions, but only for compile-time evaluation
395 -- of non-static expressions.
397 if (Even and then N * 2 > D)
399 (not Even and then N * 2 >= D)
401 Fraction := Fraction + 1;
406 -- Do not round to even as is done with IEEE arithmetic,
407 -- but instead round away from zero when the result is
408 -- exactly between two machine numbers. See RM 4.9(38).
411 Fraction := Fraction + 1;
415 if N > Uint_0 and then not UR_Is_Negative (X) then
416 Fraction := Fraction + 1;
420 if N > Uint_0 and then UR_Is_Negative (X) then
421 Fraction := Fraction + 1;
425 -- The result must be normalized to [1.0/Radix, 1.0),
426 -- so adjust if the result is 1.0 because of rounding.
428 if Fraction = Most_Significant_Digit * Radix then
429 Fraction := Most_Significant_Digit;
430 Exponent := Exponent + 1;
433 -- Put back sign after applying the rounding.
435 if UR_Is_Negative (X) then
436 Fraction := -Fraction;
439 Release_And_Save (Uintp_Mark, Fraction, Exponent);
440 end Calculate_Fraction_And_Exponent;
447 function Eps_Denorm (RT : R) return T is
449 return Float_Radix ** UI_From_Int
450 (Machine_Emin (RT) - Machine_Mantissa (RT));
457 function Eps_Model (RT : R) return T is
459 return Float_Radix ** UI_From_Int (Machine_Emin (RT));
466 function Exponent (RT : R; X : T) return UI is
470 if UR_Is_Zero (X) then
473 Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);
482 function Floor (RT : R; X : T) return T is
483 XT : constant T := Truncation (RT, X);
486 if UR_Is_Positive (X) then
501 function Fraction (RT : R; X : T) return T is
505 if UR_Is_Zero (X) then
508 Decompose (RT, X, X_Frac, X_Exp);
517 function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
518 RD : constant UI := UI_Min (Radix_Digits, Machine_Mantissa (RT));
522 L := Exponent (RT, X) - RD;
523 Y := UR_From_Uint (UR_Trunc (Scaling (RT, X, -L)));
524 return Scaling (RT, Y, L);
534 Mode : Rounding_Mode;
535 Enode : Node_Id) return T
539 Emin : constant UI := UI_From_Int (Machine_Emin (RT));
542 if UR_Is_Zero (X) then
546 Decompose (RT, X, X_Frac, X_Exp, Mode);
548 -- Case of denormalized number or (gradual) underflow
550 -- A denormalized number is one with the minimum exponent Emin, but
551 -- that breaks the assumption that the first digit of the mantissa
552 -- is a one. This allows the first non-zero digit to be in any
553 -- of the remaining Mant - 1 spots. The gap between subsequent
554 -- denormalized numbers is the same as for the smallest normalized
555 -- numbers. However, the number of significant digits left decreases
556 -- as a result of the mantissa now having leading seros.
560 Emin_Den : constant UI :=
562 (Machine_Emin (RT) - Machine_Mantissa (RT) + 1);
564 if X_Exp < Emin_Den or not Denorm_On_Target then
565 if UR_Is_Negative (X) then
567 ("floating-point value underflows to -0.0?", Enode);
572 ("floating-point value underflows to 0.0?", Enode);
576 elsif Denorm_On_Target then
578 -- Emin - Mant <= X_Exp < Emin, so result is denormal.
579 -- Handle gradual underflow by first computing the
580 -- number of significant bits still available for the
581 -- mantissa and then truncating the fraction to this
584 -- If this value is different from the original
585 -- fraction, precision is lost due to gradual underflow.
587 -- We probably should round here and prevent double
588 -- rounding as a result of first rounding to a model
589 -- number and then to a machine number. However, this
590 -- is an extremely rare case that is not worth the extra
591 -- complexity. In any case, a warning is issued in cases
592 -- where gradual underflow occurs.
595 Denorm_Sig_Bits : constant UI := X_Exp - Emin_Den + 1;
597 X_Frac_Denorm : constant T := UR_From_Components
598 (UR_Trunc (Scaling (RT, abs X_Frac, Denorm_Sig_Bits)),
604 if X_Frac_Denorm /= X_Frac then
606 ("gradual underflow causes loss of precision?",
608 X_Frac := X_Frac_Denorm;
615 return Scaling (RT, X_Frac, X_Exp);
623 function Machine_Emin (RT : R) return Int is
624 Digs : constant UI := Digits_Value (RT);
628 if Vax_Float (RT) then
629 if Digs = VAXFF_Digits then
630 Emin := VAXFF_Machine_Emin;
632 elsif Digs = VAXDF_Digits then
633 Emin := VAXDF_Machine_Emin;
636 pragma Assert (Digs = VAXGF_Digits);
637 Emin := VAXGF_Machine_Emin;
640 elsif Is_AAMP_Float (RT) then
641 if Digs = AAMPS_Digits then
642 Emin := AAMPS_Machine_Emin;
645 pragma Assert (Digs = AAMPL_Digits);
646 Emin := AAMPL_Machine_Emin;
650 if Digs = IEEES_Digits then
651 Emin := IEEES_Machine_Emin;
653 elsif Digs = IEEEL_Digits then
654 Emin := IEEEL_Machine_Emin;
657 pragma Assert (Digs = IEEEX_Digits);
658 Emin := IEEEX_Machine_Emin;
665 ----------------------
666 -- Machine_Mantissa --
667 ----------------------
669 function Machine_Mantissa (RT : R) return Nat is
670 Digs : constant UI := Digits_Value (RT);
674 if Vax_Float (RT) then
675 if Digs = VAXFF_Digits then
676 Mant := VAXFF_Machine_Mantissa;
678 elsif Digs = VAXDF_Digits then
679 Mant := VAXDF_Machine_Mantissa;
682 pragma Assert (Digs = VAXGF_Digits);
683 Mant := VAXGF_Machine_Mantissa;
686 elsif Is_AAMP_Float (RT) then
687 if Digs = AAMPS_Digits then
688 Mant := AAMPS_Machine_Mantissa;
691 pragma Assert (Digs = AAMPL_Digits);
692 Mant := AAMPL_Machine_Mantissa;
696 if Digs = IEEES_Digits then
697 Mant := IEEES_Machine_Mantissa;
699 elsif Digs = IEEEL_Digits then
700 Mant := IEEEL_Machine_Mantissa;
703 pragma Assert (Digs = IEEEX_Digits);
704 Mant := IEEEX_Machine_Mantissa;
709 end Machine_Mantissa;
715 function Machine_Radix (RT : R) return Nat is
716 pragma Warnings (Off, RT);
725 function Model (RT : R; X : T) return T is
729 Decompose (RT, X, X_Frac, X_Exp);
730 return Compose (RT, X_Frac, X_Exp);
737 function Pred (RT : R; X : T) return T is
742 if abs X < Eps_Model (RT) then
743 if Denorm_On_Target then
744 return X - Eps_Denorm (RT);
746 elsif X > Ureal_0 then
748 -- Target does not support denorms, so predecessor is 0.0
753 -- Target does not support denorms, and X is 0.0
754 -- or at least bigger than -Eps_Model (RT)
756 return -Eps_Model (RT);
760 Decompose_Int (RT, X, Result_F, Result_X, Ceiling);
761 return UR_From_Components
762 (Num => Result_F - 1,
763 Den => Machine_Mantissa (RT) - Result_X,
766 -- Result_F may be false, but this is OK as UR_From_Components
767 -- handles that situation.
775 function Remainder (RT : R; X, Y : T) return T is
790 if UR_Is_Positive (X) then
802 P_Exp := Exponent (RT, P);
805 -- ??? what about zero cases?
806 Decompose (RT, Arg, Arg_Frac, Arg_Exp);
807 Decompose (RT, P, P_Frac, P_Exp);
809 P := Compose (RT, P_Frac, Arg_Exp);
810 K := Arg_Exp - P_Exp;
814 for Cnt in reverse 0 .. UI_To_Int (K) loop
815 if IEEE_Rem >= P then
817 IEEE_Rem := IEEE_Rem - P;
826 -- That completes the calculation of modulus remainder. The final step
827 -- is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
831 B := abs Y * Ureal_Half;
834 A := IEEE_Rem * Ureal_2;
838 if A > B or else (A = B and then not P_Even) then
839 IEEE_Rem := IEEE_Rem - abs Y;
842 return Sign_X * IEEE_Rem;
849 function Rounding (RT : R; X : T) return T is
854 Result := Truncation (RT, abs X);
855 Tail := abs X - Result;
857 if Tail >= Ureal_Half then
858 Result := Result + Ureal_1;
861 if UR_Is_Negative (X) then
872 function Scaling (RT : R; X : T; Adjustment : UI) return T is
873 pragma Warnings (Off, RT);
876 if Rbase (X) = Radix then
877 return UR_From_Components
878 (Num => Numerator (X),
879 Den => Denominator (X) - Adjustment,
881 Negative => UR_Is_Negative (X));
883 elsif Adjustment >= 0 then
884 return X * Radix ** Adjustment;
886 return X / Radix ** (-Adjustment);
894 function Succ (RT : R; X : T) return T is
899 if abs X < Eps_Model (RT) then
900 if Denorm_On_Target then
901 return X + Eps_Denorm (RT);
903 elsif X < Ureal_0 then
904 -- Target does not support denorms, so successor is 0.0
908 -- Target does not support denorms, and X is 0.0
909 -- or at least smaller than Eps_Model (RT)
911 return Eps_Model (RT);
915 Decompose_Int (RT, X, Result_F, Result_X, Floor);
916 return UR_From_Components
917 (Num => Result_F + 1,
918 Den => Machine_Mantissa (RT) - Result_X,
921 -- Result_F may be false, but this is OK as UR_From_Components
922 -- handles that situation.
930 function Truncation (RT : R; X : T) return T is
931 pragma Warnings (Off, RT);
934 return UR_From_Uint (UR_Trunc (X));
937 -----------------------
938 -- Unbiased_Rounding --
939 -----------------------
941 function Unbiased_Rounding (RT : R; X : T) return T is
942 Abs_X : constant T := abs X;
947 Result := Truncation (RT, Abs_X);
948 Tail := Abs_X - Result;
950 if Tail > Ureal_Half then
951 Result := Result + Ureal_1;
953 elsif Tail = Ureal_Half then
955 Truncation (RT, (Result / Ureal_2) + Ureal_Half);
958 if UR_Is_Negative (X) then
960 elsif UR_Is_Positive (X) then
963 -- For zero case, make sure sign of zero is preserved
968 end Unbiased_Rounding;