1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
5 -- S Y S T E M . F A T _ G E N --
9 -- Copyright (C) 1992-2007, 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, 51 Franklin Street, Fifth Floor, --
20 -- Boston, MA 02110-1301, USA. --
22 -- As a special exception, if other files instantiate generics from this --
23 -- unit, or you link this unit with other files to produce an executable, --
24 -- this unit does not by itself cause the resulting executable to be --
25 -- covered by the GNU General Public License. This exception does not --
26 -- however invalidate any other reasons why the executable file might be --
27 -- covered by the GNU Public License. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
34 -- The implementation here is portable to any IEEE implementation. It does
35 -- not handle non-binary radix, and also assumes that model numbers and
36 -- machine numbers are basically identical, which is not true of all possible
37 -- floating-point implementations. On a non-IEEE machine, this body must be
38 -- specialized appropriately, or better still, its generic instantiations
39 -- should be replaced by efficient machine-specific code.
41 with Ada
.Unchecked_Conversion
;
43 package body System
.Fat_Gen
is
45 Float_Radix
: constant T
:= T
(T
'Machine_Radix);
46 Radix_To_M_Minus_1
: constant T
:= Float_Radix
** (T
'Machine_Mantissa - 1);
48 pragma Assert
(T
'Machine_Radix = 2);
49 -- This version does not handle radix 16
51 -- Constants for Decompose and Scaling
53 Rad
: constant T
:= T
(T
'Machine_Radix);
54 Invrad
: constant T
:= 1.0 / Rad
;
56 subtype Expbits
is Integer range 0 .. 6;
57 -- 2 ** (2 ** 7) might overflow. How big can radix-16 exponents get?
59 Log_Power
: constant array (Expbits
) of Integer := (1, 2, 4, 8, 16, 32, 64);
61 R_Power
: constant array (Expbits
) of T
:=
70 R_Neg_Power
: constant array (Expbits
) of T
:=
79 -----------------------
80 -- Local Subprograms --
81 -----------------------
83 procedure Decompose
(XX
: T
; Frac
: out T
; Expo
: out UI
);
84 -- Decomposes a floating-point number into fraction and exponent parts.
85 -- Both results are signed, with Frac having the sign of XX, and UI has
86 -- the sign of the exponent. The absolute value of Frac is in the range
87 -- 0.0 <= Frac < 1.0. If Frac = 0.0 or -0.0, then Expo is always zero.
89 function Gradual_Scaling
(Adjustment
: UI
) return T
;
90 -- Like Scaling with a first argument of 1.0, but returns the smallest
91 -- denormal rather than zero when the adjustment is smaller than
92 -- Machine_Emin. Used for Succ and Pred.
98 function Adjacent
(X
, Towards
: T
) return T
is
102 elsif Towards
> X
then
113 function Ceiling
(X
: T
) return T
is
114 XT
: constant T
:= Truncation
(X
);
129 function Compose
(Fraction
: T
; Exponent
: UI
) return T
is
132 pragma Unreferenced
(Arg_Exp
);
134 Decompose
(Fraction
, Arg_Frac
, Arg_Exp
);
135 return Scaling
(Arg_Frac
, Exponent
);
142 function Copy_Sign
(Value
, Sign
: T
) return T
is
145 function Is_Negative
(V
: T
) return Boolean;
146 pragma Import
(Intrinsic
, Is_Negative
);
151 if Is_Negative
(Sign
) then
162 procedure Decompose
(XX
: T
; Frac
: out T
; Expo
: out UI
) is
163 X
: constant T
:= T
'Machine (XX
);
170 -- More useful would be defining Expo to be T'Machine_Emin - 1 or
171 -- T'Machine_Emin - T'Machine_Mantissa, which would preserve
172 -- monotonicity of the exponent function ???
174 -- Check for infinities, transfinites, whatnot
176 elsif X
> T
'Safe_Last then
178 Expo
:= T
'Machine_Emax + 1;
180 elsif X
< T
'Safe_First then
182 Expo
:= T
'Machine_Emax + 2; -- how many extra negative values?
185 -- Case of nonzero finite x. Essentially, we just multiply
186 -- by Rad ** (+-2**N) to reduce the range.
192 -- Ax * Rad ** Ex is invariant
196 while Ax
>= R_Power
(Expbits
'Last) loop
197 Ax
:= Ax
* R_Neg_Power
(Expbits
'Last);
198 Ex
:= Ex
+ Log_Power
(Expbits
'Last);
203 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
204 if Ax
>= R_Power
(N
) then
205 Ax
:= Ax
* R_Neg_Power
(N
);
206 Ex
:= Ex
+ Log_Power
(N
);
220 while Ax
< R_Neg_Power
(Expbits
'Last) loop
221 Ax
:= Ax
* R_Power
(Expbits
'Last);
222 Ex
:= Ex
- Log_Power
(Expbits
'Last);
225 -- Rad ** -64 <= Ax < 1
227 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
228 if Ax
< R_Neg_Power
(N
) then
229 Ax
:= Ax
* R_Power
(N
);
230 Ex
:= Ex
- Log_Power
(N
);
233 -- R_Neg_Power (N) <= Ax < 1
252 function Exponent
(X
: T
) return UI
is
255 pragma Unreferenced
(X_Frac
);
257 Decompose
(X
, X_Frac
, X_Exp
);
265 function Floor
(X
: T
) return T
is
266 XT
: constant T
:= Truncation
(X
);
281 function Fraction
(X
: T
) return T
is
284 pragma Unreferenced
(X_Exp
);
286 Decompose
(X
, X_Frac
, X_Exp
);
290 ---------------------
291 -- Gradual_Scaling --
292 ---------------------
294 function Gradual_Scaling
(Adjustment
: UI
) return T
is
297 Ex
: UI
:= Adjustment
;
300 if Adjustment
< T
'Machine_Emin - 1 then
301 Y
:= 2.0 ** T
'Machine_Emin;
303 Ex
:= Ex
- T
'Machine_Emin;
305 Y
:= T
'Machine (Y
/ 2.0);
318 return Scaling
(1.0, Adjustment
);
326 function Leading_Part
(X
: T
; Radix_Digits
: UI
) return T
is
331 if Radix_Digits
>= T
'Machine_Mantissa then
334 elsif Radix_Digits
<= 0 then
335 raise Constraint_Error
;
338 L
:= Exponent
(X
) - Radix_Digits
;
339 Y
:= Truncation
(Scaling
(X
, -L
));
349 -- The trick with Machine is to force the compiler to store the result
350 -- in memory so that we do not have extra precision used. The compiler
351 -- is clever, so we have to outwit its possible optimizations! We do
352 -- this by using an intermediate pragma Volatile location.
354 function Machine
(X
: T
) return T
is
356 pragma Volatile
(Temp
);
362 ----------------------
363 -- Machine_Rounding --
364 ----------------------
366 -- For now, the implementation is identical to that of Rounding, which is
367 -- a permissible behavior, but is not the most efficient possible approach.
369 function Machine_Rounding
(X
: T
) return T
is
374 Result
:= Truncation
(abs X
);
375 Tail
:= abs X
- Result
;
378 Result
:= Result
+ 1.0;
387 -- For zero case, make sure sign of zero is preserved
392 end Machine_Rounding
;
398 -- We treat Model as identical to Machine. This is true of IEEE and other
399 -- nice floating-point systems, but not necessarily true of all systems.
401 function Model
(X
: T
) return T
is
410 -- Subtract from the given number a number equivalent to the value of its
411 -- least significant bit. Given that the most significant bit represents
412 -- a value of 1.0 * radix ** (exp - 1), the value we want is obtained by
413 -- shifting this by (mantissa-1) bits to the right, i.e. decreasing the
414 -- exponent by that amount.
416 -- Zero has to be treated specially, since its exponent is zero
418 function Pred
(X
: T
) return T
is
427 Decompose
(X
, X_Frac
, X_Exp
);
429 -- A special case, if the number we had was a positive power of
430 -- two, then we want to subtract half of what we would otherwise
431 -- subtract, since the exponent is going to be reduced.
433 -- Note that X_Frac has the same sign as X, so if X_Frac is 0.5,
434 -- then we know that we have a positive number (and hence a
435 -- positive power of 2).
438 return X
- Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa - 1);
440 -- Otherwise the exponent is unchanged
443 return X
- Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa);
452 function Remainder
(X
, Y
: T
) return T
is
466 pragma Unreferenced
(Arg_Frac
);
470 raise Constraint_Error
;
486 P_Exp
:= Exponent
(P
);
489 Decompose
(Arg
, Arg_Frac
, Arg_Exp
);
490 Decompose
(P
, P_Frac
, P_Exp
);
492 P
:= Compose
(P_Frac
, Arg_Exp
);
493 K
:= Arg_Exp
- P_Exp
;
497 for Cnt
in reverse 0 .. K
loop
498 if IEEE_Rem
>= P
then
500 IEEE_Rem
:= IEEE_Rem
- P
;
509 -- That completes the calculation of modulus remainder. The final
510 -- step is get the IEEE remainder. Here we need to compare Rem with
511 -- (abs Y) / 2. We must be careful of unrepresentable Y/2 value
512 -- caused by subnormal numbers
523 if A
> B
or else (A
= B
and then not P_Even
) then
524 IEEE_Rem
:= IEEE_Rem
- abs Y
;
527 return Sign_X
* IEEE_Rem
;
534 function Rounding
(X
: T
) return T
is
539 Result
:= Truncation
(abs X
);
540 Tail
:= abs X
- Result
;
543 Result
:= Result
+ 1.0;
552 -- For zero case, make sure sign of zero is preserved
563 -- Return x * rad ** adjustment quickly,
564 -- or quietly underflow to zero, or overflow naturally.
566 function Scaling
(X
: T
; Adjustment
: UI
) return T
is
568 if X
= 0.0 or else Adjustment
= 0 then
572 -- Nonzero x essentially, just multiply repeatedly by Rad ** (+-2**n)
576 Ex
: UI
:= Adjustment
;
578 -- Y * Rad ** Ex is invariant
582 while Ex
<= -Log_Power
(Expbits
'Last) loop
583 Y
:= Y
* R_Neg_Power
(Expbits
'Last);
584 Ex
:= Ex
+ Log_Power
(Expbits
'Last);
589 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
590 if Ex
<= -Log_Power
(N
) then
591 Y
:= Y
* R_Neg_Power
(N
);
592 Ex
:= Ex
+ Log_Power
(N
);
595 -- -Log_Power (N) < Ex <= 0
603 while Ex
>= Log_Power
(Expbits
'Last) loop
604 Y
:= Y
* R_Power
(Expbits
'Last);
605 Ex
:= Ex
- Log_Power
(Expbits
'Last);
610 for N
in reverse Expbits
'First .. Expbits
'Last - 1 loop
611 if Ex
>= Log_Power
(N
) then
612 Y
:= Y
* R_Power
(N
);
613 Ex
:= Ex
- Log_Power
(N
);
616 -- 0 <= Ex < Log_Power (N)
631 -- Similar computation to that of Pred: find value of least significant
632 -- bit of given number, and add. Zero has to be treated specially since
633 -- the exponent can be zero, and also we want the smallest denormal if
634 -- denormals are supported.
636 function Succ
(X
: T
) return T
is
643 X1
:= 2.0 ** T
'Machine_Emin;
645 -- Following loop generates smallest denormal
648 X2
:= T
'Machine (X1
/ 2.0);
656 Decompose
(X
, X_Frac
, X_Exp
);
658 -- A special case, if the number we had was a negative power of
659 -- two, then we want to add half of what we would otherwise add,
660 -- since the exponent is going to be reduced.
662 -- Note that X_Frac has the same sign as X, so if X_Frac is -0.5,
663 -- then we know that we have a negative number (and hence a
664 -- negative power of 2).
666 if X_Frac
= -0.5 then
667 return X
+ Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa - 1);
669 -- Otherwise the exponent is unchanged
672 return X
+ Gradual_Scaling
(X_Exp
- T
'Machine_Mantissa);
681 -- The basic approach is to compute
683 -- T'Machine (RM1 + N) - RM1
685 -- where N >= 0.0 and RM1 = radix ** (mantissa - 1)
687 -- This works provided that the intermediate result (RM1 + N) does not
688 -- have extra precision (which is why we call Machine). When we compute
689 -- RM1 + N, the exponent of N will be normalized and the mantissa shifted
690 -- shifted appropriately so the lower order bits, which cannot contribute
691 -- to the integer part of N, fall off on the right. When we subtract RM1
692 -- again, the significant bits of N are shifted to the left, and what we
693 -- have is an integer, because only the first e bits are different from
694 -- zero (assuming binary radix here).
696 function Truncation
(X
: T
) return T
is
702 if Result
>= Radix_To_M_Minus_1
then
706 Result
:= Machine
(Radix_To_M_Minus_1
+ Result
) - Radix_To_M_Minus_1
;
708 if Result
> abs X
then
709 Result
:= Result
- 1.0;
718 -- For zero case, make sure sign of zero is preserved
726 -----------------------
727 -- Unbiased_Rounding --
728 -----------------------
730 function Unbiased_Rounding
(X
: T
) return T
is
731 Abs_X
: constant T
:= abs X
;
736 Result
:= Truncation
(Abs_X
);
737 Tail
:= Abs_X
- Result
;
740 Result
:= Result
+ 1.0;
742 elsif Tail
= 0.5 then
743 Result
:= 2.0 * Truncation
((Result
/ 2.0) + 0.5);
752 -- For zero case, make sure sign of zero is preserved
757 end Unbiased_Rounding
;
763 -- Note: this routine does not work for VAX float. We compensate for this
764 -- in Exp_Attr by using the Valid functions in Vax_Float_Operations rather
765 -- than the corresponding instantiation of this function.
767 function Valid
(X
: not null access T
) return Boolean is
769 IEEE_Emin
: constant Integer := T
'Machine_Emin - 1;
770 IEEE_Emax
: constant Integer := T
'Machine_Emax - 1;
772 IEEE_Bias
: constant Integer := -(IEEE_Emin
- 1);
774 subtype IEEE_Exponent_Range
is
775 Integer range IEEE_Emin
- 1 .. IEEE_Emax
+ 1;
777 -- The implementation of this floating point attribute uses a
778 -- representation type Float_Rep that allows direct access to the
779 -- exponent and mantissa parts of a floating point number.
781 -- The Float_Rep type is an array of Float_Word elements. This
782 -- representation is chosen to make it possible to size the type based
783 -- on a generic parameter. Since the array size is known at compile
784 -- time, efficient code can still be generated. The size of Float_Word
785 -- elements should be large enough to allow accessing the exponent in
786 -- one read, but small enough so that all floating point object sizes
787 -- are a multiple of the Float_Word'Size.
789 -- The following conditions must be met for all possible
790 -- instantiations of the attributes package:
792 -- - T'Size is an integral multiple of Float_Word'Size
794 -- - The exponent and sign are completely contained in a single
795 -- component of Float_Rep, named Most_Significant_Word (MSW).
797 -- - The sign occupies the most significant bit of the MSW and the
798 -- exponent is in the following bits. Unused bits (if any) are in
799 -- the least significant part.
801 type Float_Word
is mod 2**Positive'Min (System
.Word_Size
, 32);
802 type Rep_Index
is range 0 .. 7;
804 Rep_Words
: constant Positive :=
805 (T
'Size + Float_Word
'Size - 1) / Float_Word
'Size;
806 Rep_Last
: constant Rep_Index
:= Rep_Index
'Min
807 (Rep_Index
(Rep_Words
- 1), (T
'Mantissa + 16) / Float_Word
'Size);
808 -- Determine the number of Float_Words needed for representing the
809 -- entire floating-point value. Do not take into account excessive
810 -- padding, as occurs on IA-64 where 80 bits floats get padded to 128
811 -- bits. In general, the exponent field cannot be larger than 15 bits,
812 -- even for 128-bit floating-point types, so the final format size
813 -- won't be larger than T'Mantissa + 16.
816 array (Rep_Index
range 0 .. Rep_Index
(Rep_Words
- 1)) of Float_Word
;
818 pragma Suppress_Initialization
(Float_Rep
);
819 -- This pragma suppresses the generation of an initialization procedure
820 -- for type Float_Rep when operating in Initialize/Normalize_Scalars
821 -- mode. This is not just a matter of efficiency, but of functionality,
822 -- since Valid has a pragma Inline_Always, which is not permitted if
823 -- there are nested subprograms present.
825 Most_Significant_Word
: constant Rep_Index
:=
826 Rep_Last
* Standard
'Default_Bit_Order;
827 -- Finding the location of the Exponent_Word is a bit tricky. In general
828 -- we assume Word_Order = Bit_Order. This expression needs to be refined
831 Exponent_Factor
: constant Float_Word
:=
832 2**(Float_Word
'Size - 1) /
833 Float_Word
(IEEE_Emax
- IEEE_Emin
+ 3) *
834 Boolean'Pos (Most_Significant_Word
/= 2) +
835 Boolean'Pos (Most_Significant_Word
= 2);
836 -- Factor that the extracted exponent needs to be divided by to be in
837 -- range 0 .. IEEE_Emax - IEEE_Emin + 2. Special kludge: Exponent_Factor
838 -- is 1 for x86/IA64 double extended as GCC adds unused bits to the
841 Exponent_Mask
: constant Float_Word
:=
842 Float_Word
(IEEE_Emax
- IEEE_Emin
+ 2) *
844 -- Value needed to mask out the exponent field. This assumes that the
845 -- range IEEE_Emin - 1 .. IEEE_Emax + contains 2**N values, for some N
848 function To_Float
is new Ada
.Unchecked_Conversion
(Float_Rep
, T
);
850 type Float_Access
is access all T
;
851 function To_Address
is
852 new Ada
.Unchecked_Conversion
(Float_Access
, System
.Address
);
854 XA
: constant System
.Address
:= To_Address
(Float_Access
(X
));
857 pragma Import
(Ada
, R
);
858 for R
'Address use XA
;
859 -- R is a view of the input floating-point parameter. Note that we
860 -- must avoid copying the actual bits of this parameter in float
861 -- form (since it may be a signalling NaN.
863 E
: constant IEEE_Exponent_Range
:=
864 Integer ((R
(Most_Significant_Word
) and Exponent_Mask
) /
867 -- Mask/Shift T to only get bits from the exponent. Then convert biased
868 -- value to integer value.
871 -- Float_Rep representation of significant of X.all
876 -- All denormalized numbers are valid, so the only invalid numbers
877 -- are overflows and NaNs, both with exponent = Emax + 1.
879 return E
/= IEEE_Emax
+ 1;
883 -- All denormalized numbers except 0.0 are invalid
885 -- Set exponent of X to zero, so we end up with the significand, which
886 -- definitely is a valid number and can be converted back to a float.
889 SR
(Most_Significant_Word
) :=
890 (SR
(Most_Significant_Word
)
891 and not Exponent_Mask
) + Float_Word
(IEEE_Bias
) * Exponent_Factor
;
893 return (E
in IEEE_Emin
.. IEEE_Emax
) or else
894 ((E
= IEEE_Emin
- 1) and then abs To_Float
(SR
) = 1.0);
897 ---------------------
898 -- Unaligned_Valid --
899 ---------------------
901 function Unaligned_Valid
(A
: System
.Address
) return Boolean is
902 subtype FS
is String (1 .. T
'Size / Character'Size);
903 type FSP
is access FS
;
905 function To_FSP
is new Ada
.Unchecked_Conversion
(Address
, FSP
);
910 -- Note that we have to be sure that we do not load the value into a
911 -- floating-point register, since a signalling NaN may cause a trap.
912 -- The following assignment is what does the actual alignment, since
913 -- we know that the target Local_T is aligned.
915 To_FSP
(Local_T
'Address).all := To_FSP
(A
).all;
917 -- Now that we have an aligned value, we can use the normal aligned
918 -- version of Valid to obtain the required result.
920 return Valid
(Local_T
'Access);