1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
5 -- S Y S T E M . F A T _ G E N --
9 -- Copyright (C) 1992-2009, 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. --
18 -- As a special exception under Section 7 of GPL version 3, you are granted --
19 -- additional permissions described in the GCC Runtime Library Exception, --
20 -- version 3.1, as published by the Free Software Foundation. --
22 -- You should have received a copy of the GNU General Public License and --
23 -- a copy of the GCC Runtime Library Exception along with this program; --
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
25 -- <http://www.gnu.org/licenses/>. --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 ------------------------------------------------------------------------------
32 -- The implementation here is portable to any IEEE implementation. It does
33 -- not handle non-binary radix, and also assumes that model numbers and
34 -- machine numbers are basically identical, which is not true of all possible
35 -- floating-point implementations. On a non-IEEE machine, this body must be
36 -- specialized appropriately, or better still, its generic instantiations
37 -- should be replaced by efficient machine-specific code.
39 with Ada.Unchecked_Conversion;
41 package body System.Fat_Gen is
43 Float_Radix : constant T := T (T'Machine_Radix);
44 Radix_To_M_Minus_1 : constant T := Float_Radix ** (T'Machine_Mantissa - 1);
46 pragma Assert (T'Machine_Radix = 2);
47 -- This version does not handle radix 16
49 -- Constants for Decompose and Scaling
51 Rad : constant T := T (T'Machine_Radix);
52 Invrad : constant T := 1.0 / Rad;
54 subtype Expbits is Integer range 0 .. 6;
55 -- 2 ** (2 ** 7) might overflow. How big can radix-16 exponents get?
57 Log_Power : constant array (Expbits) of Integer := (1, 2, 4, 8, 16, 32, 64);
59 R_Power : constant array (Expbits) of T :=
68 R_Neg_Power : constant array (Expbits) of T :=
77 -----------------------
78 -- Local Subprograms --
79 -----------------------
81 procedure Decompose (XX : T; Frac : out T; Expo : out UI);
82 -- Decomposes a floating-point number into fraction and exponent parts.
83 -- Both results are signed, with Frac having the sign of XX, and UI has
84 -- the sign of the exponent. The absolute value of Frac is in the range
85 -- 0.0 <= Frac < 1.0. If Frac = 0.0 or -0.0, then Expo is always zero.
87 function Gradual_Scaling (Adjustment : UI) return T;
88 -- Like Scaling with a first argument of 1.0, but returns the smallest
89 -- denormal rather than zero when the adjustment is smaller than
90 -- Machine_Emin. Used for Succ and Pred.
96 function Adjacent (X, Towards : T) return T is
100 elsif Towards > X then
111 function Ceiling (X : T) return T is
112 XT : constant T := Truncation (X);
127 function Compose (Fraction : T; Exponent : UI) return T is
130 pragma Unreferenced (Arg_Exp);
132 Decompose (Fraction, Arg_Frac, Arg_Exp);
133 return Scaling (Arg_Frac, Exponent);
140 function Copy_Sign (Value, Sign : T) return T is
143 function Is_Negative (V : T) return Boolean;
144 pragma Import (Intrinsic, Is_Negative);
149 if Is_Negative (Sign) then
160 procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
161 X : constant T := T'Machine (XX);
168 -- More useful would be defining Expo to be T'Machine_Emin - 1 or
169 -- T'Machine_Emin - T'Machine_Mantissa, which would preserve
170 -- monotonicity of the exponent function ???
172 -- Check for infinities, transfinites, whatnot
174 elsif X > T'Safe_Last then
176 Expo := T'Machine_Emax + 1;
178 elsif X < T'Safe_First then
180 Expo := T'Machine_Emax + 2; -- how many extra negative values?
183 -- Case of nonzero finite x. Essentially, we just multiply
184 -- by Rad ** (+-2**N) to reduce the range.
190 -- Ax * Rad ** Ex is invariant
194 while Ax >= R_Power (Expbits'Last) loop
195 Ax := Ax * R_Neg_Power (Expbits'Last);
196 Ex := Ex + Log_Power (Expbits'Last);
201 for N in reverse Expbits'First .. Expbits'Last - 1 loop
202 if Ax >= R_Power (N) then
203 Ax := Ax * R_Neg_Power (N);
204 Ex := Ex + Log_Power (N);
218 while Ax < R_Neg_Power (Expbits'Last) loop
219 Ax := Ax * R_Power (Expbits'Last);
220 Ex := Ex - Log_Power (Expbits'Last);
223 -- Rad ** -64 <= Ax < 1
225 for N in reverse Expbits'First .. Expbits'Last - 1 loop
226 if Ax < R_Neg_Power (N) then
227 Ax := Ax * R_Power (N);
228 Ex := Ex - Log_Power (N);
231 -- R_Neg_Power (N) <= Ax < 1
235 Frac := (if X > 0.0 then Ax else -Ax);
245 function Exponent (X : T) return UI is
248 pragma Unreferenced (X_Frac);
250 Decompose (X, X_Frac, X_Exp);
258 function Floor (X : T) return T is
259 XT : constant T := Truncation (X);
274 function Fraction (X : T) return T is
277 pragma Unreferenced (X_Exp);
279 Decompose (X, X_Frac, X_Exp);
283 ---------------------
284 -- Gradual_Scaling --
285 ---------------------
287 function Gradual_Scaling (Adjustment : UI) return T is
290 Ex : UI := Adjustment;
293 if Adjustment < T'Machine_Emin - 1 then
294 Y := 2.0 ** T'Machine_Emin;
296 Ex := Ex - T'Machine_Emin;
298 Y := T'Machine (Y / 2.0);
311 return Scaling (1.0, Adjustment);
319 function Leading_Part (X : T; Radix_Digits : UI) return T is
324 if Radix_Digits >= T'Machine_Mantissa then
327 elsif Radix_Digits <= 0 then
328 raise Constraint_Error;
331 L := Exponent (X) - Radix_Digits;
332 Y := Truncation (Scaling (X, -L));
342 -- The trick with Machine is to force the compiler to store the result
343 -- in memory so that we do not have extra precision used. The compiler
344 -- is clever, so we have to outwit its possible optimizations! We do
345 -- this by using an intermediate pragma Volatile location.
347 function Machine (X : T) return T is
349 pragma Volatile (Temp);
355 ----------------------
356 -- Machine_Rounding --
357 ----------------------
359 -- For now, the implementation is identical to that of Rounding, which is
360 -- a permissible behavior, but is not the most efficient possible approach.
362 function Machine_Rounding (X : T) return T is
367 Result := Truncation (abs X);
368 Tail := abs X - Result;
371 Result := Result + 1.0;
380 -- For zero case, make sure sign of zero is preserved
385 end Machine_Rounding;
391 -- We treat Model as identical to Machine. This is true of IEEE and other
392 -- nice floating-point systems, but not necessarily true of all systems.
394 function Model (X : T) return T is
403 -- Subtract from the given number a number equivalent to the value of its
404 -- least significant bit. Given that the most significant bit represents
405 -- a value of 1.0 * radix ** (exp - 1), the value we want is obtained by
406 -- shifting this by (mantissa-1) bits to the right, i.e. decreasing the
407 -- exponent by that amount.
409 -- Zero has to be treated specially, since its exponent is zero
411 function Pred (X : T) return T is
420 Decompose (X, X_Frac, X_Exp);
422 -- A special case, if the number we had was a positive power of
423 -- two, then we want to subtract half of what we would otherwise
424 -- subtract, since the exponent is going to be reduced.
426 -- Note that X_Frac has the same sign as X, so if X_Frac is 0.5,
427 -- then we know that we have a positive number (and hence a
428 -- positive power of 2).
431 return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
433 -- Otherwise the exponent is unchanged
436 return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa);
445 function Remainder (X, Y : T) return T is
459 pragma Unreferenced (Arg_Frac);
463 raise Constraint_Error;
479 P_Exp := Exponent (P);
482 Decompose (Arg, Arg_Frac, Arg_Exp);
483 Decompose (P, P_Frac, P_Exp);
485 P := Compose (P_Frac, Arg_Exp);
486 K := Arg_Exp - P_Exp;
490 for Cnt in reverse 0 .. K loop
491 if IEEE_Rem >= P then
493 IEEE_Rem := IEEE_Rem - P;
502 -- That completes the calculation of modulus remainder. The final
503 -- step is get the IEEE remainder. Here we need to compare Rem with
504 -- (abs Y) / 2. We must be careful of unrepresentable Y/2 value
505 -- caused by subnormal numbers
516 if A > B or else (A = B and then not P_Even) then
517 IEEE_Rem := IEEE_Rem - abs Y;
520 return Sign_X * IEEE_Rem;
527 function Rounding (X : T) return T is
532 Result := Truncation (abs X);
533 Tail := abs X - Result;
536 Result := Result + 1.0;
545 -- For zero case, make sure sign of zero is preserved
556 -- Return x * rad ** adjustment quickly,
557 -- or quietly underflow to zero, or overflow naturally.
559 function Scaling (X : T; Adjustment : UI) return T is
561 if X = 0.0 or else Adjustment = 0 then
565 -- Nonzero x essentially, just multiply repeatedly by Rad ** (+-2**n)
569 Ex : UI := Adjustment;
571 -- Y * Rad ** Ex is invariant
575 while Ex <= -Log_Power (Expbits'Last) loop
576 Y := Y * R_Neg_Power (Expbits'Last);
577 Ex := Ex + Log_Power (Expbits'Last);
582 for N in reverse Expbits'First .. Expbits'Last - 1 loop
583 if Ex <= -Log_Power (N) then
584 Y := Y * R_Neg_Power (N);
585 Ex := Ex + Log_Power (N);
588 -- -Log_Power (N) < Ex <= 0
596 while Ex >= Log_Power (Expbits'Last) loop
597 Y := Y * R_Power (Expbits'Last);
598 Ex := Ex - Log_Power (Expbits'Last);
603 for N in reverse Expbits'First .. Expbits'Last - 1 loop
604 if Ex >= Log_Power (N) then
605 Y := Y * R_Power (N);
606 Ex := Ex - Log_Power (N);
609 -- 0 <= Ex < Log_Power (N)
624 -- Similar computation to that of Pred: find value of least significant
625 -- bit of given number, and add. Zero has to be treated specially since
626 -- the exponent can be zero, and also we want the smallest denormal if
627 -- denormals are supported.
629 function Succ (X : T) return T is
636 X1 := 2.0 ** T'Machine_Emin;
638 -- Following loop generates smallest denormal
641 X2 := T'Machine (X1 / 2.0);
649 Decompose (X, X_Frac, X_Exp);
651 -- A special case, if the number we had was a negative power of
652 -- two, then we want to add half of what we would otherwise add,
653 -- since the exponent is going to be reduced.
655 -- Note that X_Frac has the same sign as X, so if X_Frac is -0.5,
656 -- then we know that we have a negative number (and hence a
657 -- negative power of 2).
659 if X_Frac = -0.5 then
660 return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
662 -- Otherwise the exponent is unchanged
665 return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa);
674 -- The basic approach is to compute
676 -- T'Machine (RM1 + N) - RM1
678 -- where N >= 0.0 and RM1 = radix ** (mantissa - 1)
680 -- This works provided that the intermediate result (RM1 + N) does not
681 -- have extra precision (which is why we call Machine). When we compute
682 -- RM1 + N, the exponent of N will be normalized and the mantissa shifted
683 -- shifted appropriately so the lower order bits, which cannot contribute
684 -- to the integer part of N, fall off on the right. When we subtract RM1
685 -- again, the significant bits of N are shifted to the left, and what we
686 -- have is an integer, because only the first e bits are different from
687 -- zero (assuming binary radix here).
689 function Truncation (X : T) return T is
695 if Result >= Radix_To_M_Minus_1 then
699 Result := Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1;
701 if Result > abs X then
702 Result := Result - 1.0;
711 -- For zero case, make sure sign of zero is preserved
719 -----------------------
720 -- Unbiased_Rounding --
721 -----------------------
723 function Unbiased_Rounding (X : T) return T is
724 Abs_X : constant T := abs X;
729 Result := Truncation (Abs_X);
730 Tail := Abs_X - Result;
733 Result := Result + 1.0;
735 elsif Tail = 0.5 then
736 Result := 2.0 * Truncation ((Result / 2.0) + 0.5);
745 -- For zero case, make sure sign of zero is preserved
750 end Unbiased_Rounding;
756 -- Note: this routine does not work for VAX float. We compensate for this
757 -- in Exp_Attr by using the Valid functions in Vax_Float_Operations rather
758 -- than the corresponding instantiation of this function.
760 function Valid (X : not null access T) return Boolean is
762 IEEE_Emin : constant Integer := T'Machine_Emin - 1;
763 IEEE_Emax : constant Integer := T'Machine_Emax - 1;
765 IEEE_Bias : constant Integer := -(IEEE_Emin - 1);
767 subtype IEEE_Exponent_Range is
768 Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
770 -- The implementation of this floating point attribute uses a
771 -- representation type Float_Rep that allows direct access to the
772 -- exponent and mantissa parts of a floating point number.
774 -- The Float_Rep type is an array of Float_Word elements. This
775 -- representation is chosen to make it possible to size the type based
776 -- on a generic parameter. Since the array size is known at compile
777 -- time, efficient code can still be generated. The size of Float_Word
778 -- elements should be large enough to allow accessing the exponent in
779 -- one read, but small enough so that all floating point object sizes
780 -- are a multiple of the Float_Word'Size.
782 -- The following conditions must be met for all possible
783 -- instantiations of the attributes package:
785 -- - T'Size is an integral multiple of Float_Word'Size
787 -- - The exponent and sign are completely contained in a single
788 -- component of Float_Rep, named Most_Significant_Word (MSW).
790 -- - The sign occupies the most significant bit of the MSW and the
791 -- exponent is in the following bits. Unused bits (if any) are in
792 -- the least significant part.
794 type Float_Word is mod 2**Positive'Min (System.Word_Size, 32);
795 type Rep_Index is range 0 .. 7;
797 Rep_Words : constant Positive :=
798 (T'Size + Float_Word'Size - 1) / Float_Word'Size;
799 Rep_Last : constant Rep_Index := Rep_Index'Min
800 (Rep_Index (Rep_Words - 1), (T'Mantissa + 16) / Float_Word'Size);
801 -- Determine the number of Float_Words needed for representing the
802 -- entire floating-point value. Do not take into account excessive
803 -- padding, as occurs on IA-64 where 80 bits floats get padded to 128
804 -- bits. In general, the exponent field cannot be larger than 15 bits,
805 -- even for 128-bit floating-point types, so the final format size
806 -- won't be larger than T'Mantissa + 16.
809 array (Rep_Index range 0 .. Rep_Index (Rep_Words - 1)) of Float_Word;
811 pragma Suppress_Initialization (Float_Rep);
812 -- This pragma suppresses the generation of an initialization procedure
813 -- for type Float_Rep when operating in Initialize/Normalize_Scalars
814 -- mode. This is not just a matter of efficiency, but of functionality,
815 -- since Valid has a pragma Inline_Always, which is not permitted if
816 -- there are nested subprograms present.
818 Most_Significant_Word : constant Rep_Index :=
819 Rep_Last * Standard'Default_Bit_Order;
820 -- Finding the location of the Exponent_Word is a bit tricky. In general
821 -- we assume Word_Order = Bit_Order. This expression needs to be refined
824 Exponent_Factor : constant Float_Word :=
825 2**(Float_Word'Size - 1) /
826 Float_Word (IEEE_Emax - IEEE_Emin + 3) *
827 Boolean'Pos (Most_Significant_Word /= 2) +
828 Boolean'Pos (Most_Significant_Word = 2);
829 -- Factor that the extracted exponent needs to be divided by to be in
830 -- range 0 .. IEEE_Emax - IEEE_Emin + 2. Special kludge: Exponent_Factor
831 -- is 1 for x86/IA64 double extended as GCC adds unused bits to the
834 Exponent_Mask : constant Float_Word :=
835 Float_Word (IEEE_Emax - IEEE_Emin + 2) *
837 -- Value needed to mask out the exponent field. This assumes that the
838 -- range IEEE_Emin - 1 .. IEEE_Emax + contains 2**N values, for some N
841 function To_Float is new Ada.Unchecked_Conversion (Float_Rep, T);
843 type Float_Access is access all T;
844 function To_Address is
845 new Ada.Unchecked_Conversion (Float_Access, System.Address);
847 XA : constant System.Address := To_Address (Float_Access (X));
850 pragma Import (Ada, R);
851 for R'Address use XA;
852 -- R is a view of the input floating-point parameter. Note that we
853 -- must avoid copying the actual bits of this parameter in float
854 -- form (since it may be a signalling NaN.
856 E : constant IEEE_Exponent_Range :=
857 Integer ((R (Most_Significant_Word) and Exponent_Mask) /
860 -- Mask/Shift T to only get bits from the exponent. Then convert biased
861 -- value to integer value.
864 -- Float_Rep representation of significant of X.all
869 -- All denormalized numbers are valid, so the only invalid numbers
870 -- are overflows and NaNs, both with exponent = Emax + 1.
872 return E /= IEEE_Emax + 1;
876 -- All denormalized numbers except 0.0 are invalid
878 -- Set exponent of X to zero, so we end up with the significand, which
879 -- definitely is a valid number and can be converted back to a float.
882 SR (Most_Significant_Word) :=
883 (SR (Most_Significant_Word)
884 and not Exponent_Mask) + Float_Word (IEEE_Bias) * Exponent_Factor;
886 return (E in IEEE_Emin .. IEEE_Emax) or else
887 ((E = IEEE_Emin - 1) and then abs To_Float (SR) = 1.0);
890 ---------------------
891 -- Unaligned_Valid --
892 ---------------------
894 function Unaligned_Valid (A : System.Address) return Boolean is
895 subtype FS is String (1 .. T'Size / Character'Size);
896 type FSP is access FS;
898 function To_FSP is new Ada.Unchecked_Conversion (Address, FSP);
903 -- Note that we have to be sure that we do not load the value into a
904 -- floating-point register, since a signalling NaN may cause a trap.
905 -- The following assignment is what does the actual alignment, since
906 -- we know that the target Local_T is aligned.
908 To_FSP (Local_T'Address).all := To_FSP (A).all;
910 -- Now that we have an aligned value, we can use the normal aligned
911 -- version of Valid to obtain the required result.
913 return Valid (Local_T'Access);