1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
9 -- Copyright (C) 1992-2018, 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 3, 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 COPYING3. If not, go to --
19 -- http://www.gnu.org/licenses for a complete copy of the license. --
21 -- GNAT was originally developed by the GNAT team at New York University. --
22 -- Extensive contributions were provided by Ada Core Technologies Inc. --
24 ------------------------------------------------------------------------------
26 with Einfo; use Einfo;
27 with Errout; use Errout;
29 with Sem_Util; use Sem_Util;
31 package body Eval_Fat is
33 Radix : constant Int := 2;
34 -- This code is currently only correct for the radix 2 case. We use the
35 -- symbolic value Radix where possible to help in the unlikely case of
36 -- anyone ever having to adjust this code for another value, and for
37 -- documentation purposes.
39 -- Another assumption is that the range of the floating-point type is
40 -- symmetric around zero.
42 type Radix_Power_Table is array (Int range 1 .. 4) of Int;
44 Radix_Powers : constant Radix_Power_Table :=
45 (Radix ** 1, Radix ** 2, Radix ** 3, Radix ** 4);
47 -----------------------
48 -- Local Subprograms --
49 -----------------------
56 Mode : Rounding_Mode := Round);
57 -- Decomposes a non-zero floating-point number into fraction and exponent
58 -- parts. The fraction is in the interval 1.0 / Radix .. T'Pred (1.0) and
59 -- uses Rbase = Radix. The result is rounded to a nearest machine number.
65 function Adjacent (RT : R; X, Towards : T) return T is
69 elsif Towards > X then
80 function Ceiling (RT : R; X : T) return T is
81 XT : constant T := Truncation (RT, X);
83 if UR_Is_Negative (X) then
96 function Compose (RT : R; Fraction : T; Exponent : UI) return T is
99 pragma Warnings (Off, Arg_Exp);
101 Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
102 return Scaling (RT, Arg_Frac, Exponent);
109 function Copy_Sign (RT : R; Value, Sign : T) return T is
110 pragma Warnings (Off, RT);
116 if UR_Is_Negative (Sign) then
132 Mode : Rounding_Mode := Round)
137 Decompose_Int (RT, abs X, Int_F, Exponent, Mode);
139 Fraction := UR_From_Components
141 Den => Machine_Mantissa_Value (RT),
145 if UR_Is_Negative (X) then
146 Fraction := -Fraction;
156 -- This procedure should be modified with care, as there are many non-
157 -- obvious details that may cause problems that are hard to detect. For
158 -- zero arguments, Fraction and Exponent are set to zero. Note that sign
159 -- of zero cannot be preserved.
161 procedure Decompose_Int
166 Mode : Rounding_Mode)
168 Base : Int := Rbase (X);
169 N : UI := abs Numerator (X);
170 D : UI := Denominator (X);
175 -- True iff Fraction is even
177 Most_Significant_Digit : constant UI :=
178 Radix ** (Machine_Mantissa_Value (RT) - 1);
180 Uintp_Mark : Uintp.Save_Mark;
181 -- The code is divided into blocks that systematically release
182 -- intermediate values (this routine generates lots of junk).
191 Calculate_D_And_Exponent_1 : begin
195 -- In cases where Base > 1, the actual denominator is Base**D. For
196 -- cases where Base is a power of Radix, use the value 1 for the
197 -- Denominator and adjust the exponent.
199 -- Note: Exponent has different sign from D, because D is a divisor
201 for Power in 1 .. Radix_Powers'Last loop
202 if Base = Radix_Powers (Power) then
203 Exponent := -D * Power;
210 Release_And_Save (Uintp_Mark, D, Exponent);
211 end Calculate_D_And_Exponent_1;
214 Calculate_Exponent : begin
217 -- For bases that are a multiple of the Radix, divide the base by
218 -- Radix and adjust the Exponent. This will help because D will be
219 -- much smaller and faster to process.
221 -- This occurs for decimal bases on machines with binary floating-
222 -- point for example. When calculating 1E40, with Radix = 2, N
223 -- will be 93 bits instead of 133.
231 -- = -------------------------- * Radix
233 -- (Base/Radix) * Radix
236 -- = --------------- * Radix
240 -- This code is commented out, because it causes numerous
241 -- failures in the regression suite. To be studied ???
243 while False and then Base > 0 and then Base mod Radix = 0 loop
244 Base := Base / Radix;
245 Exponent := Exponent + D;
248 Release_And_Save (Uintp_Mark, Exponent);
249 end Calculate_Exponent;
251 -- For remaining bases we must actually compute the exponentiation
253 -- Because the exponentiation can be negative, and D must be integer,
254 -- the numerator is corrected instead.
256 Calculate_N_And_D : begin
260 N := N * Base ** (-D);
266 Release_And_Save (Uintp_Mark, N, D);
267 end Calculate_N_And_D;
272 -- Now scale N and D so that N / D is a value in the interval [1.0 /
273 -- Radix, 1.0) and adjust Exponent accordingly, so the value N / D *
274 -- Radix ** Exponent remains unchanged.
276 -- Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0
278 -- N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.
279 -- As this scaling is not possible for N is Uint_0, zero is handled
280 -- explicitly at the start of this subprogram.
282 Calculate_N_And_Exponent : begin
285 N_Times_Radix := N * Radix;
286 while not (N_Times_Radix >= D) loop
288 Exponent := Exponent - 1;
289 N_Times_Radix := N * Radix;
292 Release_And_Save (Uintp_Mark, N, Exponent);
293 end Calculate_N_And_Exponent;
295 -- Step 2 - Adjust D so N / D < 1
297 -- Scale up D so N / D < 1, so N < D
299 Calculate_D_And_Exponent_2 : begin
302 while not (N < D) loop
304 -- As N / D >= 1, N / (D * Radix) will be at least 1 / Radix, so
305 -- the result of Step 1 stays valid
308 Exponent := Exponent + 1;
311 Release_And_Save (Uintp_Mark, D, Exponent);
312 end Calculate_D_And_Exponent_2;
314 -- Here the value N / D is in the range [1.0 / Radix .. 1.0)
316 -- Now find the fraction by doing a very simple-minded division until
317 -- enough digits have been computed.
319 -- This division works for all radices, but is only efficient for a
320 -- binary radix. It is just like a manual division algorithm, but
321 -- instead of moving the denominator one digit right, we move the
322 -- numerator one digit left so the numerator and denominator remain
328 Calculate_Fraction_And_N : begin
334 Fraction := Fraction + 1;
338 -- Stop when the result is in [1.0 / Radix, 1.0)
340 exit when Fraction >= Most_Significant_Digit;
343 Fraction := Fraction * Radix;
347 Release_And_Save (Uintp_Mark, Fraction, N);
348 end Calculate_Fraction_And_N;
350 Calculate_Fraction_And_Exponent : begin
353 -- Determine correct rounding based on the remainder which is in
354 -- N and the divisor D. The rounding is performed on the absolute
355 -- value of X, so Ceiling and Floor need to check for the sign of
361 -- This rounding mode corresponds to the unbiased rounding
362 -- method that is used at run time. When the real value is
363 -- exactly between two machine numbers, choose the machine
364 -- number with its least significant bit equal to zero.
366 -- The recommendation advice in RM 4.9(38) is that static
367 -- expressions are rounded to machine numbers in the same
368 -- way as the target machine does.
370 if (Even and then N * 2 > D)
372 (not Even and then N * 2 >= D)
374 Fraction := Fraction + 1;
379 -- Do not round to even as is done with IEEE arithmetic, but
380 -- instead round away from zero when the result is exactly
381 -- between two machine numbers. This biased rounding method
382 -- should not be used to convert static expressions to
383 -- machine numbers, see AI95-268.
386 Fraction := Fraction + 1;
390 if N > Uint_0 and then not UR_Is_Negative (X) then
391 Fraction := Fraction + 1;
395 if N > Uint_0 and then UR_Is_Negative (X) then
396 Fraction := Fraction + 1;
400 -- The result must be normalized to [1.0/Radix, 1.0), so adjust if
401 -- the result is 1.0 because of rounding.
403 if Fraction = Most_Significant_Digit * Radix then
404 Fraction := Most_Significant_Digit;
405 Exponent := Exponent + 1;
408 -- Put back sign after applying the rounding
410 if UR_Is_Negative (X) then
411 Fraction := -Fraction;
414 Release_And_Save (Uintp_Mark, Fraction, Exponent);
415 end Calculate_Fraction_And_Exponent;
422 function Exponent (RT : R; X : T) return UI is
425 pragma Warnings (Off, X_Frac);
427 Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);
435 function Floor (RT : R; X : T) return T is
436 XT : constant T := Truncation (RT, X);
439 if UR_Is_Positive (X) then
454 function Fraction (RT : R; X : T) return T is
457 pragma Warnings (Off, X_Exp);
459 Decompose (RT, X, X_Frac, X_Exp);
467 function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
468 RD : constant UI := UI_Min (Radix_Digits, Machine_Mantissa_Value (RT));
472 L := Exponent (RT, X) - RD;
473 Y := UR_From_Uint (UR_Trunc (Scaling (RT, X, -L)));
474 return Scaling (RT, Y, L);
484 Mode : Rounding_Mode;
485 Enode : Node_Id) return T
489 Emin : constant UI := Machine_Emin_Value (RT);
492 Decompose (RT, X, X_Frac, X_Exp, Mode);
494 -- Case of denormalized number or (gradual) underflow
496 -- A denormalized number is one with the minimum exponent Emin, but that
497 -- breaks the assumption that the first digit of the mantissa is a one.
498 -- This allows the first non-zero digit to be in any of the remaining
499 -- Mant - 1 spots. The gap between subsequent denormalized numbers is
500 -- the same as for the smallest normalized numbers. However, the number
501 -- of significant digits left decreases as a result of the mantissa now
502 -- having leading seros.
506 Emin_Den : constant UI := Machine_Emin_Value (RT) -
507 Machine_Mantissa_Value (RT) + Uint_1;
510 -- Do not issue warnings about underflows in GNATprove mode,
511 -- as calling Machine as part of interval checking may lead
512 -- to spurious warnings.
514 if X_Exp < Emin_Den or not Has_Denormals (RT) then
515 if Has_Signed_Zeros (RT) and then UR_Is_Negative (X) then
516 if not GNATprove_Mode then
518 ("floating-point value underflows to -0.0??", Enode);
524 if not GNATprove_Mode then
526 ("floating-point value underflows to 0.0??", Enode);
532 elsif Has_Denormals (RT) then
534 -- Emin - Mant <= X_Exp < Emin, so result is denormal. Handle
535 -- gradual underflow by first computing the number of
536 -- significant bits still available for the mantissa and
537 -- then truncating the fraction to this number of bits.
539 -- If this value is different from the original fraction,
540 -- precision is lost due to gradual underflow.
542 -- We probably should round here and prevent double rounding as
543 -- a result of first rounding to a model number and then to a
544 -- machine number. However, this is an extremely rare case that
545 -- is not worth the extra complexity. In any case, a warning is
546 -- issued in cases where gradual underflow occurs.
549 Denorm_Sig_Bits : constant UI := X_Exp - Emin_Den + 1;
551 X_Frac_Denorm : constant T := UR_From_Components
552 (UR_Trunc (Scaling (RT, abs X_Frac, Denorm_Sig_Bits)),
558 -- Do not issue warnings about loss of precision in
559 -- GNATprove mode, as calling Machine as part of interval
560 -- checking may lead to spurious warnings.
562 if X_Frac_Denorm /= X_Frac then
563 if not GNATprove_Mode then
565 ("gradual underflow causes loss of precision??",
568 X_Frac := X_Frac_Denorm;
575 return Scaling (RT, X_Frac, X_Exp);
582 function Model (RT : R; X : T) return T is
586 Decompose (RT, X, X_Frac, X_Exp);
587 return Compose (RT, X_Frac, X_Exp);
594 function Pred (RT : R; X : T) return T is
596 return -Succ (RT, -X);
603 function Remainder (RT : R; X, Y : T) return T is
617 pragma Warnings (Off, Arg_Frac);
620 if UR_Is_Positive (X) then
632 P_Exp := Exponent (RT, P);
635 -- ??? what about zero cases?
636 Decompose (RT, Arg, Arg_Frac, Arg_Exp);
637 Decompose (RT, P, P_Frac, P_Exp);
639 P := Compose (RT, P_Frac, Arg_Exp);
640 K := Arg_Exp - P_Exp;
644 for Cnt in reverse 0 .. UI_To_Int (K) loop
645 if IEEE_Rem >= P then
647 IEEE_Rem := IEEE_Rem - P;
656 -- That completes the calculation of modulus remainder. The final step
657 -- is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
661 B := abs Y * Ureal_Half;
664 A := IEEE_Rem * Ureal_2;
668 if A > B or else (A = B and then not P_Even) then
669 IEEE_Rem := IEEE_Rem - abs Y;
672 return Sign_X * IEEE_Rem;
679 function Rounding (RT : R; X : T) return T is
684 Result := Truncation (RT, abs X);
685 Tail := abs X - Result;
687 if Tail >= Ureal_Half then
688 Result := Result + Ureal_1;
691 if UR_Is_Negative (X) then
702 function Scaling (RT : R; X : T; Adjustment : UI) return T is
703 pragma Warnings (Off, RT);
706 if Rbase (X) = Radix then
707 return UR_From_Components
708 (Num => Numerator (X),
709 Den => Denominator (X) - Adjustment,
711 Negative => UR_Is_Negative (X));
713 elsif Adjustment >= 0 then
714 return X * Radix ** Adjustment;
716 return X / Radix ** (-Adjustment);
724 function Succ (RT : R; X : T) return T is
725 Emin : constant UI := Machine_Emin_Value (RT);
726 Mantissa : constant UI := Machine_Mantissa_Value (RT);
727 Exp : UI := UI_Max (Emin, Exponent (RT, X));
732 if UR_Is_Zero (X) then
736 -- Set exponent such that the radix point will be directly following the
737 -- mantissa after scaling.
739 if Has_Denormals (RT) or Exp /= Emin then
740 Exp := Exp - Mantissa;
745 Frac := Scaling (RT, X, -Exp);
746 New_Frac := Ceiling (RT, Frac);
748 if New_Frac = Frac then
749 if New_Frac = Scaling (RT, -Ureal_1, Mantissa - 1) then
750 New_Frac := New_Frac + Scaling (RT, Ureal_1, Uint_Minus_1);
752 New_Frac := New_Frac + Ureal_1;
756 return Scaling (RT, New_Frac, Exp);
763 function Truncation (RT : R; X : T) return T is
764 pragma Warnings (Off, RT);
766 return UR_From_Uint (UR_Trunc (X));
769 -----------------------
770 -- Unbiased_Rounding --
771 -----------------------
773 function Unbiased_Rounding (RT : R; X : T) return T is
774 Abs_X : constant T := abs X;
779 Result := Truncation (RT, Abs_X);
780 Tail := Abs_X - Result;
782 if Tail > Ureal_Half then
783 Result := Result + Ureal_1;
785 elsif Tail = Ureal_Half then
787 Truncation (RT, (Result / Ureal_2) + Ureal_Half);
790 if UR_Is_Negative (X) then
792 elsif UR_Is_Positive (X) then
795 -- For zero case, make sure sign of zero is preserved
800 end Unbiased_Rounding;