1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
5 -- S Y S T E M . F A T _ G E N --
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 -- As a special exception, if other files instantiate generics from this --
24 -- unit, or you link this unit with other files to produce an executable, --
25 -- this unit does not by itself cause the resulting executable to be --
26 -- covered by the GNU General Public License. This exception does not --
27 -- however invalidate any other reasons why the executable file might be --
28 -- covered by the GNU Public License. --
30 -- GNAT was originally developed by the GNAT team at New York University. --
31 -- Extensive contributions were provided by Ada Core Technologies Inc. --
33 ------------------------------------------------------------------------------
35 -- The implementation here is portable to any IEEE implementation. It does
36 -- not handle non-binary radix, and also assumes that model numbers and
37 -- machine numbers are basically identical, which is not true of all possible
38 -- floating-point implementations. On a non-IEEE machine, this body must be
39 -- specialized appropriately, or better still, its generic instantiations
40 -- should be replaced by efficient machine-specific code.
42 with Ada
.Unchecked_Conversion
;
44 package body System
.Fat_Gen
is
46 Float_Radix
: constant T
:= T
(T
'Machine_Radix);
47 Float_Radix_Inv
: constant T
:= 1.0 / Float_Radix
;
48 Radix_To_M_Minus_1
: constant T
:= Float_Radix
** (T
'Machine_Mantissa - 1);
50 pragma Assert
(T
'Machine_Radix = 2);
51 -- This version does not handle radix 16
53 -- Constants for Decompose and Scaling
55 Rad
: constant T
:= T
(T
'Machine_Radix);
56 Invrad
: constant T
:= 1.0 / Rad
;
58 subtype Expbits
is Integer range 0 .. 6;
59 -- 2 ** (2 ** 7) might overflow. how big can radix-16 exponents get?
61 Log_Power
: constant array (Expbits
) of Integer := (1, 2, 4, 8, 16, 32, 64);
63 R_Power
: constant array (Expbits
) of T
:=
72 R_Neg_Power
: constant array (Expbits
) of T
:=
81 -----------------------
82 -- Local Subprograms --
83 -----------------------
85 procedure Decompose
(XX
: T
; Frac
: out T
; Expo
: out UI
);
86 -- Decomposes a floating-point number into fraction and exponent parts
88 function Gradual_Scaling
(Adjustment
: UI
) return T
;
89 -- Like Scaling with a first argument of 1.0, but returns the smallest
90 -- denormal rather than zero when the adjustment is smaller than
91 -- Machine_Emin. Used for Succ and Pred.
97 function Adjacent
(X
, Towards
: T
) return T
is
102 elsif Towards
> X
then
114 function Ceiling
(X
: T
) return T
is
115 XT
: constant T
:= Truncation
(X
);
133 function Compose
(Fraction
: T
; Exponent
: UI
) return T
is
138 Decompose
(Fraction
, Arg_Frac
, Arg_Exp
);
139 return Scaling
(Arg_Frac
, Exponent
);
146 function Copy_Sign
(Value
, Sign
: T
) return T
is
149 function Is_Negative
(V
: T
) return Boolean;
150 pragma Import
(Intrinsic
, Is_Negative
);
155 if Is_Negative
(Sign
) then
166 procedure Decompose
(XX
: T
; Frac
: out T
; Expo
: out UI
) is
167 X
: T
:= T
'Machine (XX
);
174 -- More useful would be defining Expo to be T'Machine_Emin - 1 or
175 -- T'Machine_Emin - T'Machine_Mantissa, which would preserve
176 -- monotonicity of the exponent function ???
178 -- Check for infinities, transfinites, whatnot.
180 elsif X
> T
'Safe_Last then
182 Expo
:= T
'Machine_Emax + 1;
184 elsif X
< T
'Safe_First then
186 Expo
:= T
'Machine_Emax + 2; -- how many extra negative values?
189 -- Case of nonzero finite x. Essentially, we just multiply
190 -- by Rad ** (+-2**N) to reduce the range.
196 -- Ax * Rad ** Ex is invariant.
200 while Ax
>= R_Power
(Expbits
'Last) loop
201 Ax
:= Ax
* R_Neg_Power
(Expbits
'Last);
202 Ex
:= Ex
+ Log_Power
(Expbits
'Last);
207 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
208 if Ax
>= R_Power
(N
) then
209 Ax
:= Ax
* R_Neg_Power
(N
);
210 Ex
:= Ex
+ Log_Power
(N
);
224 while Ax
< R_Neg_Power
(Expbits
'Last) loop
225 Ax
:= Ax
* R_Power
(Expbits
'Last);
226 Ex
:= Ex
- Log_Power
(Expbits
'Last);
229 -- Rad ** -64 <= Ax < 1
231 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
232 if Ax
< R_Neg_Power
(N
) then
233 Ax
:= Ax
* R_Power
(N
);
234 Ex
:= Ex
- Log_Power
(N
);
237 -- R_Neg_Power (N) <= Ax < 1
256 function Exponent
(X
: T
) return UI
is
261 Decompose
(X
, X_Frac
, X_Exp
);
269 function Floor
(X
: T
) return T
is
270 XT
: constant T
:= Truncation
(X
);
288 function Fraction
(X
: T
) return T
is
293 Decompose
(X
, X_Frac
, X_Exp
);
297 ---------------------
298 -- Gradual_Scaling --
299 ---------------------
301 function Gradual_Scaling
(Adjustment
: UI
) return T
is
304 Ex
: UI
:= Adjustment
;
307 if Adjustment
< T
'Machine_Emin then
308 Y
:= 2.0 ** T
'Machine_Emin;
310 Ex
:= Ex
- T
'Machine_Emin;
313 Y
:= T
'Machine (Y
/ 2.0);
326 return Scaling
(1.0, Adjustment
);
334 function Leading_Part
(X
: T
; Radix_Digits
: UI
) return T
is
339 if Radix_Digits
>= T
'Machine_Mantissa then
343 L
:= Exponent
(X
) - Radix_Digits
;
344 Y
:= Truncation
(Scaling
(X
, -L
));
355 -- The trick with Machine is to force the compiler to store the result
356 -- in memory so that we do not have extra precision used. The compiler
357 -- is clever, so we have to outwit its possible optimizations! We do
358 -- this by using an intermediate pragma Volatile location.
360 function Machine
(X
: T
) return T
is
362 pragma Volatile
(Temp
);
373 -- We treat Model as identical to Machine. This is true of IEEE and other
374 -- nice floating-point systems, but not necessarily true of all systems.
376 function Model
(X
: T
) return T
is
385 -- Subtract from the given number a number equivalent to the value of its
386 -- least significant bit. Given that the most significant bit represents
387 -- a value of 1.0 * radix ** (exp - 1), the value we want is obtained by
388 -- shifting this by (mantissa-1) bits to the right, i.e. decreasing the
389 -- exponent by that amount.
391 -- Zero has to be treated specially, since its exponent is zero
393 function Pred
(X
: T
) return T
is
402 Decompose
(X
, X_Frac
, X_Exp
);
404 -- A special case, if the number we had was a positive power of
405 -- two, then we want to subtract half of what we would otherwise
406 -- subtract, since the exponent is going to be reduced.
408 if X_Frac
= 0.5 and then X
> 0.0 then
409 return X
- Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa - 1);
411 -- Otherwise the exponent stays the same
414 return X
- Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa);
423 function Remainder
(X
, Y
: T
) return T
is
451 P_Exp
:= Exponent
(P
);
454 Decompose
(Arg
, Arg_Frac
, Arg_Exp
);
455 Decompose
(P
, P_Frac
, P_Exp
);
457 P
:= Compose
(P_Frac
, Arg_Exp
);
458 K
:= Arg_Exp
- P_Exp
;
462 for Cnt
in reverse 0 .. K
loop
463 if IEEE_Rem
>= P
then
465 IEEE_Rem
:= IEEE_Rem
- P
;
474 -- That completes the calculation of modulus remainder. The final
475 -- step is get the IEEE remainder. Here we need to compare Rem with
476 -- (abs Y) / 2. We must be careful of unrepresentable Y/2 value
477 -- caused by subnormal numbers
488 if A
> B
or else (A
= B
and then not P_Even
) then
489 IEEE_Rem
:= IEEE_Rem
- abs Y
;
492 return Sign_X
* IEEE_Rem
;
500 function Rounding
(X
: T
) return T
is
505 Result
:= Truncation
(abs X
);
506 Tail
:= abs X
- Result
;
509 Result
:= Result
+ 1.0;
518 -- For zero case, make sure sign of zero is preserved
530 -- Return x * rad ** adjustment quickly,
531 -- or quietly underflow to zero, or overflow naturally.
533 function Scaling
(X
: T
; Adjustment
: UI
) return T
is
535 if X
= 0.0 or else Adjustment
= 0 then
539 -- Nonzero x. essentially, just multiply repeatedly by Rad ** (+-2**n).
543 Ex
: UI
:= Adjustment
;
545 -- Y * Rad ** Ex is invariant
549 while Ex
<= -Log_Power
(Expbits
'Last) loop
550 Y
:= Y
* R_Neg_Power
(Expbits
'Last);
551 Ex
:= Ex
+ Log_Power
(Expbits
'Last);
556 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
557 if Ex
<= -Log_Power
(N
) then
558 Y
:= Y
* R_Neg_Power
(N
);
559 Ex
:= Ex
+ Log_Power
(N
);
562 -- -Log_Power (N) < Ex <= 0
570 while Ex
>= Log_Power
(Expbits
'Last) loop
571 Y
:= Y
* R_Power
(Expbits
'Last);
572 Ex
:= Ex
- Log_Power
(Expbits
'Last);
577 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
578 if Ex
>= Log_Power
(N
) then
579 Y
:= Y
* R_Power
(N
);
580 Ex
:= Ex
- Log_Power
(N
);
583 -- 0 <= Ex < Log_Power (N)
596 -- Similar computation to that of Pred: find value of least significant
597 -- bit of given number, and add. Zero has to be treated specially since
598 -- the exponent can be zero, and also we want the smallest denormal if
599 -- denormals are supported.
601 function Succ
(X
: T
) return T
is
608 X1
:= 2.0 ** T
'Machine_Emin;
610 -- Following loop generates smallest denormal
613 X2
:= T
'Machine (X1
/ 2.0);
621 Decompose
(X
, X_Frac
, X_Exp
);
623 -- A special case, if the number we had was a negative power of
624 -- two, then we want to add half of what we would otherwise add,
625 -- since the exponent is going to be reduced.
627 if X_Frac
= 0.5 and then X
< 0.0 then
628 return X
+ Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa - 1);
630 -- Otherwise the exponent stays the same
633 return X
+ Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa);
642 -- The basic approach is to compute
644 -- T'Machine (RM1 + N) - RM1.
646 -- where N >= 0.0 and RM1 = radix ** (mantissa - 1)
648 -- This works provided that the intermediate result (RM1 + N) does not
649 -- have extra precision (which is why we call Machine). When we compute
650 -- RM1 + N, the exponent of N will be normalized and the mantissa shifted
651 -- shifted appropriately so the lower order bits, which cannot contribute
652 -- to the integer part of N, fall off on the right. When we subtract RM1
653 -- again, the significant bits of N are shifted to the left, and what we
654 -- have is an integer, because only the first e bits are different from
655 -- zero (assuming binary radix here).
657 function Truncation
(X
: T
) return T
is
663 if Result
>= Radix_To_M_Minus_1
then
667 Result
:= Machine
(Radix_To_M_Minus_1
+ Result
) - Radix_To_M_Minus_1
;
669 if Result
> abs X
then
670 Result
:= Result
- 1.0;
679 -- For zero case, make sure sign of zero is preserved
688 -----------------------
689 -- Unbiased_Rounding --
690 -----------------------
692 function Unbiased_Rounding
(X
: T
) return T
is
693 Abs_X
: constant T
:= abs X
;
698 Result
:= Truncation
(Abs_X
);
699 Tail
:= Abs_X
- Result
;
702 Result
:= Result
+ 1.0;
704 elsif Tail
= 0.5 then
705 Result
:= 2.0 * Truncation
((Result
/ 2.0) + 0.5);
714 -- For zero case, make sure sign of zero is preserved
720 end Unbiased_Rounding
;
726 function Valid
(X
: access T
) return Boolean is
728 IEEE_Emin
: constant Integer := T
'Machine_Emin - 1;
729 IEEE_Emax
: constant Integer := T
'Machine_Emax - 1;
731 IEEE_Bias
: constant Integer := -(IEEE_Emin
- 1);
733 subtype IEEE_Exponent_Range
is
734 Integer range IEEE_Emin
- 1 .. IEEE_Emax
+ 1;
736 -- The implementation of this floating point attribute uses
737 -- a representation type Float_Rep that allows direct access to
738 -- the exponent and mantissa parts of a floating point number.
740 -- The Float_Rep type is an array of Float_Word elements. This
741 -- representation is chosen to make it possible to size the
742 -- type based on a generic parameter.
744 -- The following conditions must be met for all possible
745 -- instantiations of the attributes package:
747 -- - T'Size is an integral multiple of Float_Word'Size
749 -- - The exponent and sign are completely contained in a single
750 -- component of Float_Rep, named Most_Significant_Word (MSW).
752 -- - The sign occupies the most significant bit of the MSW
753 -- and the exponent is in the following bits.
754 -- Unused bits (if any) are in the least significant part.
756 type Float_Word
is mod 2**32;
757 type Rep_Index
is range 0 .. 7;
759 Rep_Last
: constant Rep_Index
:= (T
'Size - 1) / Float_Word
'Size;
761 type Float_Rep
is array (Rep_Index
range 0 .. Rep_Last
) of Float_Word
;
763 Most_Significant_Word
: constant Rep_Index
:=
764 Rep_Last
* Standard
'Default_Bit_Order;
765 -- Finding the location of the Exponent_Word is a bit tricky.
766 -- In general we assume Word_Order = Bit_Order.
767 -- This expression needs to be refined for VMS.
769 Exponent_Factor
: constant Float_Word
:=
770 2**(Float_Word
'Size - 1) /
771 Float_Word
(IEEE_Emax
- IEEE_Emin
+ 3) *
772 Boolean'Pos (T
'Size /= 96) +
773 Boolean'Pos (T
'Size = 96);
774 -- Factor that the extracted exponent needs to be divided by
775 -- to be in range 0 .. IEEE_Emax - IEEE_Emin + 2.
776 -- Special kludge: Exponent_Factor is 0 for x86 double extended
777 -- as GCC adds 16 unused bits to the type.
779 Exponent_Mask
: constant Float_Word
:=
780 Float_Word
(IEEE_Emax
- IEEE_Emin
+ 2) *
782 -- Value needed to mask out the exponent field.
783 -- This assumes that the range IEEE_Emin - 1 .. IEEE_Emax + 1
784 -- contains 2**N values, for some N in Natural.
786 function To_Float
is new Ada
.Unchecked_Conversion
(Float_Rep
, T
);
788 type Float_Access
is access all T
;
789 function To_Address
is
790 new Ada
.Unchecked_Conversion
(Float_Access
, System
.Address
);
792 XA
: constant System
.Address
:= To_Address
(Float_Access
(X
));
795 pragma Import
(Ada
, R
);
796 for R
'Address use XA
;
797 -- R is a view of the input floating-point parameter. Note that we
798 -- must avoid copying the actual bits of this parameter in float
799 -- form (since it may be a signalling NaN.
801 E
: constant IEEE_Exponent_Range
:=
802 Integer ((R
(Most_Significant_Word
) and Exponent_Mask
) /
805 -- Mask/Shift T to only get bits from the exponent
806 -- Then convert biased value to integer value.
809 -- Float_Rep representation of significant of X.all
814 -- All denormalized numbers are valid, so only invalid numbers
815 -- are overflows and NaN's, both with exponent = Emax + 1.
817 return E
/= IEEE_Emax
+ 1;
821 -- All denormalized numbers except 0.0 are invalid
823 -- Set exponent of X to zero, so we end up with the significand, which
824 -- definitely is a valid number and can be converted back to a float.
827 SR
(Most_Significant_Word
) :=
828 (SR
(Most_Significant_Word
)
829 and not Exponent_Mask
) + Float_Word
(IEEE_Bias
) * Exponent_Factor
;
831 return (E
in IEEE_Emin
.. IEEE_Emax
) or else
832 ((E
= IEEE_Emin
- 1) and then abs To_Float
(SR
) = 1.0);