1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
9 -- Copyright (C) 1992-2003 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
86 function Machine_Mantissa
(RT
: R
) return Nat
;
87 -- Return value of the Machine_Mantissa attribute
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.
402 -- This rounding mode should not be used for static
403 -- expressions, but only for compile-time evaluation
404 -- of non-static expressions.
406 if (Even
and then N
* 2 > D
)
408 (not Even
and then N
* 2 >= D
)
410 Fraction
:= Fraction
+ 1;
415 -- Do not round to even as is done with IEEE arithmetic,
416 -- but instead round away from zero when the result is
417 -- exactly between two machine numbers. See RM 4.9(38).
420 Fraction
:= Fraction
+ 1;
425 Fraction
:= Fraction
+ 1;
431 -- The result must be normalized to [1.0/Radix, 1.0),
432 -- so adjust if the result is 1.0 because of rounding.
434 if Fraction
= Most_Significant_Digit
* Radix
then
435 Fraction
:= Most_Significant_Digit
;
436 Exponent
:= Exponent
+ 1;
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
471 if UR_Is_Zero
(X
) then
474 Decompose_Int
(RT
, X
, X_Frac
, X_Exp
, Round_Even
);
483 function Floor
(RT
: R
; X
: T
) return T
is
484 XT
: constant T
:= Truncation
(RT
, X
);
487 if UR_Is_Positive
(X
) then
502 function Fraction
(RT
: R
; X
: T
) return T
is
507 if UR_Is_Zero
(X
) then
510 Decompose
(RT
, X
, X_Frac
, X_Exp
);
519 function Leading_Part
(RT
: R
; X
: T
; Radix_Digits
: UI
) return T
is
524 if Radix_Digits
>= Machine_Mantissa
(RT
) then
528 L
:= Exponent
(RT
, X
) - Radix_Digits
;
529 Y
:= Truncation
(RT
, Scaling
(RT
, X
, -L
));
530 Z
:= Scaling
(RT
, Y
, L
);
542 Mode
: Rounding_Mode
;
546 pragma Warnings
(Off
, Enode
); -- not yet referenced
550 Emin
: constant UI
:= UI_From_Int
(Machine_Emin
(RT
));
553 if UR_Is_Zero
(X
) then
557 Decompose
(RT
, X
, X_Frac
, X_Exp
, Mode
);
559 -- Case of denormalized number or (gradual) underflow
561 -- A denormalized number is one with the minimum exponent Emin, but
562 -- that breaks the assumption that the first digit of the mantissa
563 -- is a one. This allows the first non-zero digit to be in any
564 -- of the remaining Mant - 1 spots. The gap between subsequent
565 -- denormalized numbers is the same as for the smallest normalized
566 -- numbers. However, the number of significant digits left decreases
567 -- as a result of the mantissa now having leading seros.
571 Emin_Den
: constant UI
:=
573 (Machine_Emin
(RT
) - Machine_Mantissa
(RT
) + 1);
575 if X_Exp
< Emin_Den
or not Denorm_On_Target
then
576 if UR_Is_Negative
(X
) then
578 ("floating-point value underflows to -0.0?", Enode
);
583 ("floating-point value underflows to 0.0?", Enode
);
587 elsif Denorm_On_Target
then
589 -- Emin - Mant <= X_Exp < Emin, so result is denormal.
590 -- Handle gradual underflow by first computing the
591 -- number of significant bits still available for the
592 -- mantissa and then truncating the fraction to this
595 -- If this value is different from the original
596 -- fraction, precision is lost due to gradual underflow.
598 -- We probably should round here and prevent double
599 -- rounding as a result of first rounding to a model
600 -- number and then to a machine number. However, this
601 -- is an extremely rare case that is not worth the extra
602 -- complexity. In any case, a warning is issued in cases
603 -- where gradual underflow occurs.
606 Denorm_Sig_Bits
: constant UI
:= X_Exp
- Emin_Den
+ 1;
608 X_Frac_Denorm
: constant T
:= UR_From_Components
609 (UR_Trunc
(Scaling
(RT
, abs X_Frac
, Denorm_Sig_Bits
)),
615 if X_Frac_Denorm
/= X_Frac
then
617 ("gradual underflow causes loss of precision?",
619 X_Frac
:= X_Frac_Denorm
;
626 return Scaling
(RT
, X_Frac
, X_Exp
);
634 function Machine_Emin
(RT
: R
) return Int
is
635 Digs
: constant UI
:= Digits_Value
(RT
);
639 if Vax_Float
(RT
) then
640 if Digs
= VAXFF_Digits
then
641 Emin
:= VAXFF_Machine_Emin
;
643 elsif Digs
= VAXDF_Digits
then
644 Emin
:= VAXDF_Machine_Emin
;
647 pragma Assert
(Digs
= VAXGF_Digits
);
648 Emin
:= VAXGF_Machine_Emin
;
651 elsif Is_AAMP_Float
(RT
) then
652 if Digs
= AAMPS_Digits
then
653 Emin
:= AAMPS_Machine_Emin
;
656 pragma Assert
(Digs
= AAMPL_Digits
);
657 Emin
:= AAMPL_Machine_Emin
;
661 if Digs
= IEEES_Digits
then
662 Emin
:= IEEES_Machine_Emin
;
664 elsif Digs
= IEEEL_Digits
then
665 Emin
:= IEEEL_Machine_Emin
;
668 pragma Assert
(Digs
= IEEEX_Digits
);
669 Emin
:= IEEEX_Machine_Emin
;
676 ----------------------
677 -- Machine_Mantissa --
678 ----------------------
680 function Machine_Mantissa
(RT
: R
) return Nat
is
681 Digs
: constant UI
:= Digits_Value
(RT
);
685 if Vax_Float
(RT
) then
686 if Digs
= VAXFF_Digits
then
687 Mant
:= VAXFF_Machine_Mantissa
;
689 elsif Digs
= VAXDF_Digits
then
690 Mant
:= VAXDF_Machine_Mantissa
;
693 pragma Assert
(Digs
= VAXGF_Digits
);
694 Mant
:= VAXGF_Machine_Mantissa
;
697 elsif Is_AAMP_Float
(RT
) then
698 if Digs
= AAMPS_Digits
then
699 Mant
:= AAMPS_Machine_Mantissa
;
702 pragma Assert
(Digs
= AAMPL_Digits
);
703 Mant
:= AAMPL_Machine_Mantissa
;
707 if Digs
= IEEES_Digits
then
708 Mant
:= IEEES_Machine_Mantissa
;
710 elsif Digs
= IEEEL_Digits
then
711 Mant
:= IEEEL_Machine_Mantissa
;
714 pragma Assert
(Digs
= IEEEX_Digits
);
715 Mant
:= IEEEX_Machine_Mantissa
;
720 end Machine_Mantissa
;
726 function Model
(RT
: R
; X
: T
) return T
is
731 Decompose
(RT
, X
, X_Frac
, X_Exp
);
732 return Compose
(RT
, X_Frac
, X_Exp
);
739 function Pred
(RT
: R
; X
: T
) return T
is
744 if abs X
< Eps_Model
(RT
) then
745 if Denorm_On_Target
then
746 return X
- Eps_Denorm
(RT
);
748 elsif X
> Ureal_0
then
750 -- Target does not support denorms, so predecessor is 0.0
755 -- Target does not support denorms, and X is 0.0
756 -- or at least bigger than -Eps_Model (RT)
758 return -Eps_Model
(RT
);
762 Decompose_Int
(RT
, X
, Result_F
, Result_X
, Ceiling
);
763 return UR_From_Components
764 (Num
=> Result_F
- 1,
765 Den
=> Machine_Mantissa
(RT
) - Result_X
,
768 -- Result_F may be false, but this is OK as UR_From_Components
769 -- handles that situation.
777 function Remainder
(RT
: R
; X
, Y
: T
) return T
is
792 if UR_Is_Positive
(X
) then
804 P_Exp
:= Exponent
(RT
, P
);
807 -- ??? what about zero cases?
808 Decompose
(RT
, Arg
, Arg_Frac
, Arg_Exp
);
809 Decompose
(RT
, P
, P_Frac
, P_Exp
);
811 P
:= Compose
(RT
, P_Frac
, Arg_Exp
);
812 K
:= Arg_Exp
- P_Exp
;
816 for Cnt
in reverse 0 .. UI_To_Int
(K
) loop
817 if IEEE_Rem
>= P
then
819 IEEE_Rem
:= IEEE_Rem
- P
;
828 -- That completes the calculation of modulus remainder. The final step
829 -- is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
833 B
:= abs Y
* Ureal_Half
;
836 A
:= IEEE_Rem
* Ureal_2
;
840 if A
> B
or else (A
= B
and then not P_Even
) then
841 IEEE_Rem
:= IEEE_Rem
- abs Y
;
844 return Sign_X
* IEEE_Rem
;
851 function Rounding
(RT
: R
; X
: T
) return T
is
856 Result
:= Truncation
(RT
, abs X
);
857 Tail
:= abs X
- Result
;
859 if Tail
>= Ureal_Half
then
860 Result
:= Result
+ Ureal_1
;
863 if UR_Is_Negative
(X
) then
874 function Scaling
(RT
: R
; X
: T
; Adjustment
: UI
) return T
is
875 pragma Warnings
(Off
, RT
);
878 if Rbase
(X
) = Radix
then
879 return UR_From_Components
880 (Num
=> Numerator
(X
),
881 Den
=> Denominator
(X
) - Adjustment
,
883 Negative
=> UR_Is_Negative
(X
));
885 elsif Adjustment
>= 0 then
886 return X
* Radix
** Adjustment
;
888 return X
/ Radix
** (-Adjustment
);
896 function Succ
(RT
: R
; X
: T
) return T
is
901 if abs X
< Eps_Model
(RT
) then
902 if Denorm_On_Target
then
903 return X
+ Eps_Denorm
(RT
);
905 elsif X
< Ureal_0
then
906 -- Target does not support denorms, so successor is 0.0
910 -- Target does not support denorms, and X is 0.0
911 -- or at least smaller than Eps_Model (RT)
913 return Eps_Model
(RT
);
917 Decompose_Int
(RT
, X
, Result_F
, Result_X
, Floor
);
918 return UR_From_Components
919 (Num
=> Result_F
+ 1,
920 Den
=> Machine_Mantissa
(RT
) - Result_X
,
923 -- Result_F may be false, but this is OK as UR_From_Components
924 -- handles that situation.
932 function Truncation
(RT
: R
; X
: T
) return T
is
933 pragma Warnings
(Off
, RT
);
936 return UR_From_Uint
(UR_Trunc
(X
));
939 -----------------------
940 -- Unbiased_Rounding --
941 -----------------------
943 function Unbiased_Rounding
(RT
: R
; X
: T
) return T
is
944 Abs_X
: constant T
:= abs X
;
949 Result
:= Truncation
(RT
, Abs_X
);
950 Tail
:= Abs_X
- Result
;
952 if Tail
> Ureal_Half
then
953 Result
:= Result
+ Ureal_1
;
955 elsif Tail
= Ureal_Half
then
957 Truncation
(RT
, (Result
/ Ureal_2
) + Ureal_Half
);
960 if UR_Is_Negative
(X
) then
962 elsif UR_Is_Positive
(X
) then
965 -- For zero case, make sure sign of zero is preserved
970 end Unbiased_Rounding
;