1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
9 -- Copyright (C) 1992-2023, 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 Einfo
.Utils
; use Einfo
.Utils
;
28 with Errout
; use Errout
;
30 with Sem_Util
; use Sem_Util
;
32 package body Eval_Fat
is
34 Radix
: constant Int
:= 2;
35 -- This code is currently only correct for the radix 2 case. We use the
36 -- symbolic value Radix where possible to help in the unlikely case of
37 -- anyone ever having to adjust this code for another value, and for
38 -- documentation purposes.
40 -- Another assumption is that the range of the floating-point type is
41 -- symmetric around zero.
43 type Radix_Power_Table
is array (Int
range 1 .. 4) of Int
;
45 Radix_Powers
: constant Radix_Power_Table
:=
46 (Radix
** 1, Radix
** 2, Radix
** 3, Radix
** 4);
48 -----------------------
49 -- Local Subprograms --
50 -----------------------
57 Mode
: Rounding_Mode
:= Round
);
58 -- Decomposes a non-zero floating-point number into fraction and exponent
59 -- parts. The fraction is in the interval 1.0 / Radix .. T'Pred (1.0) and
60 -- uses Rbase = Radix. The result is rounded to a nearest machine number.
66 function Adjacent
(RT
: R
; X
, Towards
: T
) return T
is
70 elsif Towards
> X
then
81 function Ceiling
(RT
: R
; X
: T
) return T
is
82 XT
: constant T
:= Truncation
(RT
, X
);
84 if UR_Is_Negative
(X
) then
97 function Compose
(RT
: R
; Fraction
: T
; Exponent
: UI
) return T
is
100 pragma Warnings
(Off
, Arg_Exp
);
102 Decompose
(RT
, Fraction
, Arg_Frac
, Arg_Exp
);
103 return Scaling
(RT
, Arg_Frac
, Exponent
);
110 function Copy_Sign
(RT
: R
; Value
, Sign
: T
) return T
is
111 pragma Warnings
(Off
, RT
);
117 if UR_Is_Negative
(Sign
) then
133 Mode
: Rounding_Mode
:= Round
)
138 Decompose_Int
(RT
, abs X
, Int_F
, Exponent
, Mode
);
140 Fraction
:= UR_From_Components
142 Den
=> Machine_Mantissa_Value
(RT
),
146 if UR_Is_Negative
(X
) then
147 Fraction
:= -Fraction
;
157 -- This procedure should be modified with care, as there are many non-
158 -- obvious details that may cause problems that are hard to detect. For
159 -- zero arguments, Fraction and Exponent are set to zero. Note that sign
160 -- of zero cannot be preserved.
162 procedure Decompose_Int
167 Mode
: Rounding_Mode
)
169 Base
: Int
:= Rbase
(X
);
170 N
: UI
:= abs Numerator
(X
);
171 D
: UI
:= Denominator
(X
);
176 -- True iff Fraction is even
178 Most_Significant_Digit
: constant UI
:=
179 Radix
** (Machine_Mantissa_Value
(RT
) - 1);
181 Uintp_Mark
: Uintp
.Save_Mark
;
182 -- The code is divided into blocks that systematically release
183 -- intermediate values (this routine generates lots of junk).
192 Calculate_D_And_Exponent_1
: begin
196 -- In cases where Base > 1, the actual denominator is Base**D. For
197 -- cases where Base is a power of Radix, use the value 1 for the
198 -- Denominator and adjust the exponent.
200 -- Note: Exponent has different sign from D, because D is a divisor
202 for Power
in 1 .. Radix_Powers
'Last loop
203 if Base
= Radix_Powers
(Power
) then
204 Exponent
:= -D
* Power
;
211 Release_And_Save
(Uintp_Mark
, D
, Exponent
);
212 end Calculate_D_And_Exponent_1
;
215 Calculate_Exponent
: begin
218 -- For bases that are a multiple of the Radix, divide the base by
219 -- Radix and adjust the Exponent. This will help because D will be
220 -- much smaller and faster to process.
222 -- This occurs for decimal bases on machines with binary floating-
223 -- point for example. When calculating 1E40, with Radix = 2, N
224 -- will be 93 bits instead of 133.
232 -- = -------------------------- * Radix
234 -- (Base/Radix) * Radix
237 -- = --------------- * Radix
241 -- This code is commented out, because it causes numerous
242 -- failures in the regression suite. To be studied ???
244 while False and then Base
> 0 and then Base
mod Radix
= 0 loop
245 Base
:= Base
/ Radix
;
246 Exponent
:= Exponent
+ D
;
249 Release_And_Save
(Uintp_Mark
, Exponent
);
250 end Calculate_Exponent
;
252 -- For remaining bases we must actually compute the exponentiation
254 -- Because the exponentiation can be negative, and D must be integer,
255 -- the numerator is corrected instead.
257 Calculate_N_And_D
: begin
261 N
:= N
* Base
** (-D
);
267 Release_And_Save
(Uintp_Mark
, N
, D
);
268 end Calculate_N_And_D
;
273 -- Now scale N and D so that N / D is a value in the interval [1.0 /
274 -- Radix, 1.0) and adjust Exponent accordingly, so the value N / D *
275 -- Radix ** Exponent remains unchanged.
277 -- Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0
279 -- N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.
280 -- As this scaling is not possible for N is Uint_0, zero is handled
281 -- explicitly at the start of this subprogram.
283 Calculate_N_And_Exponent
: begin
286 N_Times_Radix
:= N
* Radix
;
287 while not (N_Times_Radix
>= D
) loop
289 Exponent
:= Exponent
- 1;
290 N_Times_Radix
:= N
* Radix
;
293 Release_And_Save
(Uintp_Mark
, N
, Exponent
);
294 end Calculate_N_And_Exponent
;
296 -- Step 2 - Adjust D so N / D < 1
298 -- Scale up D so N / D < 1, so N < D
300 Calculate_D_And_Exponent_2
: begin
303 while not (N
< D
) loop
305 -- As N / D >= 1, N / (D * Radix) will be at least 1 / Radix, so
306 -- the result of Step 1 stays valid
309 Exponent
:= Exponent
+ 1;
312 Release_And_Save
(Uintp_Mark
, D
, Exponent
);
313 end Calculate_D_And_Exponent_2
;
315 -- Here the value N / D is in the range [1.0 / Radix .. 1.0)
317 -- Now find the fraction by doing a very simple-minded division until
318 -- enough digits have been computed.
320 -- This division works for all radices, but is only efficient for a
321 -- binary radix. It is just like a manual division algorithm, but
322 -- instead of moving the denominator one digit right, we move the
323 -- numerator one digit left so the numerator and denominator remain
329 Calculate_Fraction_And_N
: begin
335 Fraction
:= Fraction
+ 1;
339 -- Stop when the result is in [1.0 / Radix, 1.0)
341 exit when Fraction
>= Most_Significant_Digit
;
344 Fraction
:= Fraction
* Radix
;
348 Release_And_Save
(Uintp_Mark
, Fraction
, N
);
349 end Calculate_Fraction_And_N
;
351 Calculate_Fraction_And_Exponent
: begin
354 -- Determine correct rounding based on the remainder which is in
355 -- N and the divisor D. The rounding is performed on the absolute
356 -- value of X, so Ceiling and Floor need to check for the sign of
362 -- This rounding mode corresponds to the unbiased rounding
363 -- method that is used at run time. When the real value is
364 -- exactly between two machine numbers, choose the machine
365 -- number with its least significant bit equal to zero.
367 -- The recommendation advice in RM 4.9(38) is that static
368 -- expressions are rounded to machine numbers in the same
369 -- way as the target machine does.
371 if (Even
and then N
* 2 > D
)
373 (not Even
and then N
* 2 >= D
)
375 Fraction
:= Fraction
+ 1;
380 -- Do not round to even as is done with IEEE arithmetic, but
381 -- instead round away from zero when the result is exactly
382 -- between two machine numbers. This biased rounding method
383 -- should not be used to convert static expressions to
384 -- machine numbers, see AI95-268.
387 Fraction
:= Fraction
+ 1;
391 if N
> Uint_0
and then not UR_Is_Negative
(X
) then
392 Fraction
:= Fraction
+ 1;
396 if N
> Uint_0
and then UR_Is_Negative
(X
) then
397 Fraction
:= Fraction
+ 1;
401 -- The result must be normalized to [1.0/Radix, 1.0), so adjust if
402 -- the result is 1.0 because of rounding.
404 if Fraction
= Most_Significant_Digit
* Radix
then
405 Fraction
:= Most_Significant_Digit
;
406 Exponent
:= Exponent
+ 1;
409 -- Put back sign after applying the rounding
411 if UR_Is_Negative
(X
) then
412 Fraction
:= -Fraction
;
415 Release_And_Save
(Uintp_Mark
, Fraction
, Exponent
);
416 end Calculate_Fraction_And_Exponent
;
423 function Exponent
(RT
: R
; X
: T
) return UI
is
426 pragma Warnings
(Off
, X_Frac
);
428 Decompose_Int
(RT
, X
, X_Frac
, X_Exp
, Round_Even
);
436 function Floor
(RT
: R
; X
: T
) return T
is
437 XT
: constant T
:= Truncation
(RT
, X
);
440 if UR_Is_Positive
(X
) then
455 function Fraction
(RT
: R
; X
: T
) return T
is
458 pragma Warnings
(Off
, X_Exp
);
460 Decompose
(RT
, X
, X_Frac
, X_Exp
);
468 function Leading_Part
(RT
: R
; X
: T
; Radix_Digits
: UI
) return T
is
469 RD
: constant UI
:= UI_Min
(Radix_Digits
, Machine_Mantissa_Value
(RT
));
473 L
:= Exponent
(RT
, X
) - RD
;
474 Y
:= UR_From_Uint
(UR_Trunc
(Scaling
(RT
, X
, -L
)));
475 return Scaling
(RT
, Y
, L
);
485 Mode
: Rounding_Mode
;
486 Enode
: Node_Id
) return T
490 Emin
: constant UI
:= Machine_Emin_Value
(RT
);
493 Decompose
(RT
, X
, X_Frac
, X_Exp
, Mode
);
495 -- Case of denormalized number or (gradual) underflow
497 -- A denormalized number is one with the minimum exponent Emin, but that
498 -- breaks the assumption that the first digit of the mantissa is a one.
499 -- This allows the first non-zero digit to be in any of the remaining
500 -- Mant - 1 spots. The gap between subsequent denormalized numbers is
501 -- the same as for the smallest normalized numbers. However, the number
502 -- of significant digits left decreases as a result of the mantissa now
503 -- having leading seros.
507 Emin_Den
: constant UI
:= Machine_Emin_Value
(RT
) -
508 Machine_Mantissa_Value
(RT
) + Uint_1
;
511 -- Do not issue warnings about underflows in GNATprove mode,
512 -- as calling Machine as part of interval checking may lead
513 -- to spurious warnings.
515 if X_Exp
< Emin_Den
or not Has_Denormals
(RT
) then
516 if Has_Signed_Zeros
(RT
) and then UR_Is_Negative
(X
) then
517 if not GNATprove_Mode
then
519 ("floating-point value underflows to -0.0??", Enode
);
525 if not GNATprove_Mode
then
527 ("floating-point value underflows to 0.0??", Enode
);
533 elsif Has_Denormals
(RT
) then
535 -- Emin - Mant <= X_Exp < Emin, so result is denormal. Handle
536 -- gradual underflow by first computing the number of
537 -- significant bits still available for the mantissa and
538 -- then truncating the fraction to this number of bits.
540 -- If this value is different from the original fraction,
541 -- precision is lost due to gradual underflow.
543 -- We probably should round here and prevent double rounding as
544 -- a result of first rounding to a model number and then to a
545 -- machine number. However, this is an extremely rare case that
546 -- is not worth the extra complexity. In any case, a warning is
547 -- issued in cases where gradual underflow occurs.
550 Denorm_Sig_Bits
: constant UI
:= X_Exp
- Emin_Den
+ 1;
552 X_Frac_Denorm
: constant T
:= UR_From_Components
553 (UR_Trunc
(Scaling
(RT
, abs X_Frac
, Denorm_Sig_Bits
)),
559 -- Do not issue warnings about loss of precision in
560 -- GNATprove mode, as calling Machine as part of interval
561 -- checking may lead to spurious warnings.
563 if X_Frac_Denorm
/= X_Frac
then
564 if not GNATprove_Mode
then
566 ("gradual underflow causes loss of precision??",
569 X_Frac
:= X_Frac_Denorm
;
576 return Scaling
(RT
, X_Frac
, X_Exp
);
583 function Model
(RT
: R
; X
: T
) return T
is
587 Decompose
(RT
, X
, X_Frac
, X_Exp
);
588 return Compose
(RT
, X_Frac
, X_Exp
);
595 function Pred
(RT
: R
; X
: T
) return T
is
597 return -Succ
(RT
, -X
);
604 function Remainder
(RT
: R
; X
, Y
: T
) return T
is
618 pragma Warnings
(Off
, Arg_Frac
);
621 if UR_Is_Positive
(X
) then
633 P_Exp
:= Exponent
(RT
, P
);
636 -- ??? what about zero cases?
637 Decompose
(RT
, Arg
, Arg_Frac
, Arg_Exp
);
638 Decompose
(RT
, P
, P_Frac
, P_Exp
);
640 P
:= Compose
(RT
, P_Frac
, Arg_Exp
);
641 K
:= Arg_Exp
- P_Exp
;
645 for Cnt
in reverse 0 .. UI_To_Int
(K
) loop
646 if IEEE_Rem
>= P
then
648 IEEE_Rem
:= IEEE_Rem
- P
;
657 -- That completes the calculation of modulus remainder. The final step
658 -- is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
662 B
:= abs Y
* Ureal_Half
;
665 A
:= IEEE_Rem
* Ureal_2
;
669 if A
> B
or else (A
= B
and then not P_Even
) then
670 IEEE_Rem
:= IEEE_Rem
- abs Y
;
673 return Sign_X
* IEEE_Rem
;
680 function Rounding
(RT
: R
; X
: T
) return T
is
685 Result
:= Truncation
(RT
, abs X
);
686 Tail
:= abs X
- Result
;
688 if Tail
>= Ureal_Half
then
689 Result
:= Result
+ Ureal_1
;
692 if UR_Is_Negative
(X
) then
703 function Scaling
(RT
: R
; X
: T
; Adjustment
: UI
) return T
is
704 pragma Warnings
(Off
, RT
);
707 if Rbase
(X
) = Radix
then
708 return UR_From_Components
709 (Num
=> Numerator
(X
),
710 Den
=> Denominator
(X
) - Adjustment
,
712 Negative
=> UR_Is_Negative
(X
));
714 elsif Adjustment
>= 0 then
715 return X
* Radix
** Adjustment
;
717 return X
/ Radix
** (-Adjustment
);
725 function Succ
(RT
: R
; X
: T
) return T
is
726 Emin
: constant UI
:= Machine_Emin_Value
(RT
);
727 Mantissa
: constant UI
:= Machine_Mantissa_Value
(RT
);
728 Exp
: UI
:= UI_Max
(Emin
, Exponent
(RT
, X
));
733 -- Treat zero as a regular denormalized number if they are supported,
734 -- otherwise return the smallest normalized number.
736 if UR_Is_Zero
(X
) then
737 if Has_Denormals
(RT
) then
740 return Scaling
(RT
, Ureal_Half
, Emin
);
744 -- Multiply the number by 2.0**(Mantissa-Exp) so that the radix point
745 -- will be directly following the mantissa after scaling.
747 Exp
:= Exp
- Mantissa
;
748 Frac
:= Scaling
(RT
, X
, -Exp
);
750 -- Round to the neareast integer towards +Inf
752 New_Frac
:= Ceiling
(RT
, Frac
);
754 -- If the rounding was a NOP, add one, except for -2.0**(Mantissa-1)
755 -- because the exponent is going to be reduced.
757 if New_Frac
= Frac
then
758 if New_Frac
= Scaling
(RT
, -Ureal_1
, Mantissa
- 1) then
759 New_Frac
:= New_Frac
+ Ureal_Half
;
761 New_Frac
:= New_Frac
+ Ureal_1
;
765 -- Divide back by 2.0**(Mantissa-Exp) to get the final result
767 return Scaling
(RT
, New_Frac
, Exp
);
774 function Truncation
(RT
: R
; X
: T
) return T
is
775 pragma Warnings
(Off
, RT
);
777 return UR_From_Uint
(UR_Trunc
(X
));
780 -----------------------
781 -- Unbiased_Rounding --
782 -----------------------
784 function Unbiased_Rounding
(RT
: R
; X
: T
) return T
is
785 Abs_X
: constant T
:= abs X
;
790 Result
:= Truncation
(RT
, Abs_X
);
791 Tail
:= Abs_X
- Result
;
793 if Tail
> Ureal_Half
then
794 Result
:= Result
+ Ureal_1
;
796 elsif Tail
= Ureal_Half
then
798 Truncation
(RT
, (Result
/ Ureal_2
) + Ureal_Half
);
801 if UR_Is_Negative
(X
) then
803 elsif UR_Is_Positive
(X
) then
806 -- For zero case, make sure sign of zero is preserved
811 end Unbiased_Rounding
;