1------------------------------------------------------------------------------
2--                                                                          --
3--                         GNAT COMPILER COMPONENTS                         --
4--                                                                          --
5--                       S Y S T E M . F A T _ G E N                        --
6--                                                                          --
7--                                 B o d y                                  --
8--                                                                          --
9--          Copyright (C) 1992-2018, Free Software Foundation, Inc.         --
10--                                                                          --
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.                                     --
17--                                                                          --
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.               --
21--                                                                          --
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/>.                                          --
26--                                                                          --
27-- GNAT was originally developed  by the GNAT team at  New York University. --
28-- Extensive contributions were provided by Ada Core Technologies Inc.      --
29--                                                                          --
30------------------------------------------------------------------------------
31
32--  The implementation here is portable to any IEEE implementation. It does
33--  not handle nonbinary 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.
38
39with Ada.Unchecked_Conversion;
40with System;
41package body System.Fat_Gen is
42
43   Float_Radix        : constant T := T (T'Machine_Radix);
44   Radix_To_M_Minus_1 : constant T := Float_Radix ** (T'Machine_Mantissa - 1);
45
46   pragma Assert (T'Machine_Radix = 2);
47   --  This version does not handle radix 16
48
49   --  Constants for Decompose and Scaling
50
51   Rad    : constant T := T (T'Machine_Radix);
52   Invrad : constant T := 1.0 / Rad;
53
54   subtype Expbits is Integer range 0 .. 6;
55   --  2 ** (2 ** 7) might overflow.  How big can radix-16 exponents get?
56
57   Log_Power : constant array (Expbits) of Integer := (1, 2, 4, 8, 16, 32, 64);
58
59   R_Power : constant array (Expbits) of T :=
60     (Rad **  1,
61      Rad **  2,
62      Rad **  4,
63      Rad **  8,
64      Rad ** 16,
65      Rad ** 32,
66      Rad ** 64);
67
68   R_Neg_Power : constant array (Expbits) of T :=
69     (Invrad **  1,
70      Invrad **  2,
71      Invrad **  4,
72      Invrad **  8,
73      Invrad ** 16,
74      Invrad ** 32,
75      Invrad ** 64);
76
77   -----------------------
78   -- Local Subprograms --
79   -----------------------
80
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.
86
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.
91
92   --------------
93   -- Adjacent --
94   --------------
95
96   function Adjacent (X, Towards : T) return T is
97   begin
98      if Towards = X then
99         return X;
100      elsif Towards > X then
101         return Succ (X);
102      else
103         return Pred (X);
104      end if;
105   end Adjacent;
106
107   -------------
108   -- Ceiling --
109   -------------
110
111   function Ceiling (X : T) return T is
112      XT : constant T := Truncation (X);
113   begin
114      if X <= 0.0 then
115         return XT;
116      elsif X = XT then
117         return X;
118      else
119         return XT + 1.0;
120      end if;
121   end Ceiling;
122
123   -------------
124   -- Compose --
125   -------------
126
127   function Compose (Fraction : T; Exponent : UI) return T is
128      Arg_Frac : T;
129      Arg_Exp  : UI;
130      pragma Unreferenced (Arg_Exp);
131   begin
132      Decompose (Fraction, Arg_Frac, Arg_Exp);
133      return Scaling (Arg_Frac, Exponent);
134   end Compose;
135
136   ---------------
137   -- Copy_Sign --
138   ---------------
139
140   function Copy_Sign (Value, Sign : T) return T is
141      Result : T;
142
143      function Is_Negative (V : T) return Boolean;
144      pragma Import (Intrinsic, Is_Negative);
145
146   begin
147      Result := abs Value;
148
149      if Is_Negative (Sign) then
150         return -Result;
151      else
152         return Result;
153      end if;
154   end Copy_Sign;
155
156   ---------------
157   -- Decompose --
158   ---------------
159
160   procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
161      X : constant T := T'Machine (XX);
162
163   begin
164      if X = 0.0 then
165
166         --  The normalized exponent of zero is zero, see RM A.5.2(15)
167
168         Frac := X;
169         Expo := 0;
170
171      --  Check for infinities, transfinites, whatnot
172
173      elsif X > T'Safe_Last then
174         Frac := Invrad;
175         Expo := T'Machine_Emax + 1;
176
177      elsif X < T'Safe_First then
178         Frac := -Invrad;
179         Expo := T'Machine_Emax + 2;    -- how many extra negative values?
180
181      else
182         --  Case of nonzero finite x. Essentially, we just multiply
183         --  by Rad ** (+-2**N) to reduce the range.
184
185         declare
186            Ax : T  := abs X;
187            Ex : UI := 0;
188
189         --  Ax * Rad ** Ex is invariant
190
191         begin
192            if Ax >= 1.0 then
193               while Ax >= R_Power (Expbits'Last) loop
194                  Ax := Ax * R_Neg_Power (Expbits'Last);
195                  Ex := Ex + Log_Power (Expbits'Last);
196               end loop;
197
198               --  Ax < Rad ** 64
199
200               for N in reverse Expbits'First .. Expbits'Last - 1 loop
201                  if Ax >= R_Power (N) then
202                     Ax := Ax * R_Neg_Power (N);
203                     Ex := Ex + Log_Power (N);
204                  end if;
205
206                  --  Ax < R_Power (N)
207
208               end loop;
209
210               --  1 <= Ax < Rad
211
212               Ax := Ax * Invrad;
213               Ex := Ex + 1;
214
215            else
216               --  0 < ax < 1
217
218               while Ax < R_Neg_Power (Expbits'Last) loop
219                  Ax := Ax * R_Power (Expbits'Last);
220                  Ex := Ex - Log_Power (Expbits'Last);
221               end loop;
222
223               --  Rad ** -64 <= Ax < 1
224
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);
229                  end if;
230
231                  --  R_Neg_Power (N) <= Ax < 1
232
233               end loop;
234            end if;
235
236            Frac := (if X > 0.0 then Ax else -Ax);
237            Expo := Ex;
238         end;
239      end if;
240   end Decompose;
241
242   --------------
243   -- Exponent --
244   --------------
245
246   function Exponent (X : T) return UI is
247      X_Frac : T;
248      X_Exp  : UI;
249      pragma Unreferenced (X_Frac);
250   begin
251      Decompose (X, X_Frac, X_Exp);
252      return X_Exp;
253   end Exponent;
254
255   -----------
256   -- Floor --
257   -----------
258
259   function Floor (X : T) return T is
260      XT : constant T := Truncation (X);
261   begin
262      if X >= 0.0 then
263         return XT;
264      elsif XT = X then
265         return X;
266      else
267         return XT - 1.0;
268      end if;
269   end Floor;
270
271   --------------
272   -- Fraction --
273   --------------
274
275   function Fraction (X : T) return T is
276      X_Frac : T;
277      X_Exp  : UI;
278      pragma Unreferenced (X_Exp);
279   begin
280      Decompose (X, X_Frac, X_Exp);
281      return X_Frac;
282   end Fraction;
283
284   ---------------------
285   -- Gradual_Scaling --
286   ---------------------
287
288   function Gradual_Scaling  (Adjustment : UI) return T is
289      Y  : T;
290      Y1 : T;
291      Ex : UI := Adjustment;
292
293   begin
294      if Adjustment < T'Machine_Emin - 1 then
295         Y  := 2.0 ** T'Machine_Emin;
296         Y1 := Y;
297         Ex := Ex - T'Machine_Emin;
298         while Ex < 0 loop
299            Y := T'Machine (Y / 2.0);
300
301            if Y = 0.0 then
302               return Y1;
303            end if;
304
305            Ex := Ex + 1;
306            Y1 := Y;
307         end loop;
308
309         return Y1;
310
311      else
312         return Scaling (1.0, Adjustment);
313      end if;
314   end Gradual_Scaling;
315
316   ------------------
317   -- Leading_Part --
318   ------------------
319
320   function Leading_Part (X : T; Radix_Digits : UI) return T is
321      L    : UI;
322      Y, Z : T;
323
324   begin
325      if Radix_Digits >= T'Machine_Mantissa then
326         return X;
327
328      elsif Radix_Digits <= 0 then
329         raise Constraint_Error;
330
331      else
332         L := Exponent (X) - Radix_Digits;
333         Y := Truncation (Scaling (X, -L));
334         Z := Scaling (Y, L);
335         return Z;
336      end if;
337   end Leading_Part;
338
339   -------------
340   -- Machine --
341   -------------
342
343   --  The trick with Machine is to force the compiler to store the result
344   --  in memory so that we do not have extra precision used. The compiler
345   --  is clever, so we have to outwit its possible optimizations. We do
346   --  this by using an intermediate pragma Volatile location.
347
348   function Machine (X : T) return T is
349      Temp : T;
350      pragma Volatile (Temp);
351   begin
352      Temp := X;
353      return Temp;
354   end Machine;
355
356   ----------------------
357   -- Machine_Rounding --
358   ----------------------
359
360   --  For now, the implementation is identical to that of Rounding, which is
361   --  a permissible behavior, but is not the most efficient possible approach.
362
363   function Machine_Rounding (X : T) return T is
364      Result : T;
365      Tail   : T;
366
367   begin
368      Result := Truncation (abs X);
369      Tail   := abs X - Result;
370
371      if Tail >= 0.5 then
372         Result := Result + 1.0;
373      end if;
374
375      if X > 0.0 then
376         return Result;
377
378      elsif X < 0.0 then
379         return -Result;
380
381      --  For zero case, make sure sign of zero is preserved
382
383      else
384         return X;
385      end if;
386   end Machine_Rounding;
387
388   -----------
389   -- Model --
390   -----------
391
392   --  We treat Model as identical to Machine. This is true of IEEE and other
393   --  nice floating-point systems, but not necessarily true of all systems.
394
395   function Model (X : T) return T is
396   begin
397      return T'Machine (X);
398   end Model;
399
400   ----------
401   -- Pred --
402   ----------
403
404   function Pred (X : T) return T is
405      X_Frac : T;
406      X_Exp  : UI;
407
408   begin
409      --  Zero has to be treated specially, since its exponent is zero
410
411      if X = 0.0 then
412         return -Succ (X);
413
414      --  Special treatment for most negative number
415
416      elsif X = T'First then
417
418         --  If not generating infinities, we raise a constraint error
419
420         if T'Machine_Overflows then
421            raise Constraint_Error with "Pred of largest negative number";
422
423         --  Otherwise generate a negative infinity
424
425         else
426            return X / (X - X);
427         end if;
428
429      --  For infinities, return unchanged
430
431      elsif X < T'First or else X > T'Last then
432         return X;
433
434      --  Subtract from the given number a number equivalent to the value
435      --  of its least significant bit. Given that the most significant bit
436      --  represents a value of 1.0 * radix ** (exp - 1), the value we want
437      --  is obtained by shifting this by (mantissa-1) bits to the right,
438      --  i.e. decreasing the exponent by that amount.
439
440      else
441         Decompose (X, X_Frac, X_Exp);
442
443         --  A special case, if the number we had was a positive power of
444         --  two, then we want to subtract half of what we would otherwise
445         --  subtract, since the exponent is going to be reduced.
446
447         --  Note that X_Frac has the same sign as X, so if X_Frac is 0.5,
448         --  then we know that we have a positive number (and hence a
449         --  positive power of 2).
450
451         if X_Frac = 0.5 then
452            return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
453
454         --  Otherwise the exponent is unchanged
455
456         else
457            return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa);
458         end if;
459      end if;
460   end Pred;
461
462   ---------------
463   -- Remainder --
464   ---------------
465
466   function Remainder (X, Y : T) return T is
467      A        : T;
468      B        : T;
469      Arg      : T;
470      P        : T;
471      P_Frac   : T;
472      Sign_X   : T;
473      IEEE_Rem : T;
474      Arg_Exp  : UI;
475      P_Exp    : UI;
476      K        : UI;
477      P_Even   : Boolean;
478
479      Arg_Frac : T;
480      pragma Unreferenced (Arg_Frac);
481
482   begin
483      if Y = 0.0 then
484         raise Constraint_Error;
485      end if;
486
487      if X > 0.0 then
488         Sign_X :=  1.0;
489         Arg := X;
490      else
491         Sign_X := -1.0;
492         Arg := -X;
493      end if;
494
495      P := abs Y;
496
497      if Arg < P then
498         P_Even := True;
499         IEEE_Rem := Arg;
500         P_Exp := Exponent (P);
501
502      else
503         Decompose (Arg, Arg_Frac, Arg_Exp);
504         Decompose (P,   P_Frac,   P_Exp);
505
506         P := Compose (P_Frac, Arg_Exp);
507         K := Arg_Exp - P_Exp;
508         P_Even := True;
509         IEEE_Rem := Arg;
510
511         for Cnt in reverse 0 .. K loop
512            if IEEE_Rem >= P then
513               P_Even := False;
514               IEEE_Rem := IEEE_Rem - P;
515            else
516               P_Even := True;
517            end if;
518
519            P := P * 0.5;
520         end loop;
521      end if;
522
523      --  That completes the calculation of modulus remainder. The final
524      --  step is get the IEEE remainder. Here we need to compare Rem with
525      --  (abs Y) / 2. We must be careful of unrepresentable Y/2 value
526      --  caused by subnormal numbers
527
528      if P_Exp >= 0 then
529         A := IEEE_Rem;
530         B := abs Y * 0.5;
531
532      else
533         A := IEEE_Rem * 2.0;
534         B := abs Y;
535      end if;
536
537      if A > B or else (A = B and then not P_Even) then
538         IEEE_Rem := IEEE_Rem - abs Y;
539      end if;
540
541      return Sign_X * IEEE_Rem;
542   end Remainder;
543
544   --------------
545   -- Rounding --
546   --------------
547
548   function Rounding (X : T) return T is
549      Result : T;
550      Tail   : T;
551
552   begin
553      Result := Truncation (abs X);
554      Tail   := abs X - Result;
555
556      if Tail >= 0.5 then
557         Result := Result + 1.0;
558      end if;
559
560      if X > 0.0 then
561         return Result;
562
563      elsif X < 0.0 then
564         return -Result;
565
566      --  For zero case, make sure sign of zero is preserved
567
568      else
569         return X;
570      end if;
571   end Rounding;
572
573   -------------
574   -- Scaling --
575   -------------
576
577   --  Return x * rad ** adjustment quickly, or quietly underflow to zero,
578   --  or overflow naturally.
579
580   function Scaling (X : T; Adjustment : UI) return T is
581   begin
582      if X = 0.0 or else Adjustment = 0 then
583         return X;
584      end if;
585
586      --  Nonzero x essentially, just multiply repeatedly by Rad ** (+-2**n)
587
588      declare
589         Y  : T  := X;
590         Ex : UI := Adjustment;
591
592      --  Y * Rad ** Ex is invariant
593
594      begin
595         if Ex < 0 then
596            while Ex <= -Log_Power (Expbits'Last) loop
597               Y := Y * R_Neg_Power (Expbits'Last);
598               Ex := Ex + Log_Power (Expbits'Last);
599            end loop;
600
601            --  -64 < Ex <= 0
602
603            for N in reverse Expbits'First .. Expbits'Last - 1 loop
604               if Ex <= -Log_Power (N) then
605                  Y := Y * R_Neg_Power (N);
606                  Ex := Ex + Log_Power (N);
607               end if;
608
609               --  -Log_Power (N) < Ex <= 0
610
611            end loop;
612
613            --  Ex = 0
614
615         else
616            --  Ex >= 0
617
618            while Ex >= Log_Power (Expbits'Last) loop
619               Y := Y * R_Power (Expbits'Last);
620               Ex := Ex - Log_Power (Expbits'Last);
621            end loop;
622
623            --  0 <= Ex < 64
624
625            for N in reverse Expbits'First .. Expbits'Last - 1 loop
626               if Ex >= Log_Power (N) then
627                  Y := Y * R_Power (N);
628                  Ex := Ex - Log_Power (N);
629               end if;
630
631               --  0 <= Ex < Log_Power (N)
632
633            end loop;
634
635            --  Ex = 0
636
637         end if;
638
639         return Y;
640      end;
641   end Scaling;
642
643   ----------
644   -- Succ --
645   ----------
646
647   function Succ (X : T) return T is
648      X_Frac : T;
649      X_Exp  : UI;
650      X1, X2 : T;
651
652   begin
653      --  Treat zero specially since it has a zero exponent
654
655      if X = 0.0 then
656         X1 := 2.0 ** T'Machine_Emin;
657
658         --  Following loop generates smallest denormal
659
660         loop
661            X2 := T'Machine (X1 / 2.0);
662            exit when X2 = 0.0;
663            X1 := X2;
664         end loop;
665
666         return X1;
667
668      --  Special treatment for largest positive number
669
670      elsif X = T'Last then
671
672         --  If not generating infinities, we raise a constraint error
673
674         if T'Machine_Overflows then
675            raise Constraint_Error with "Succ of largest negative number";
676
677         --  Otherwise generate a positive infinity
678
679         else
680            return X / (X - X);
681         end if;
682
683      --  For infinities, return unchanged
684
685      elsif X < T'First or else X > T'Last then
686         return X;
687
688      --  Add to the given number a number equivalent to the value
689      --  of its least significant bit. Given that the most significant bit
690      --  represents a value of 1.0 * radix ** (exp - 1), the value we want
691      --  is obtained by shifting this by (mantissa-1) bits to the right,
692      --  i.e. decreasing the exponent by that amount.
693
694      else
695         Decompose (X, X_Frac, X_Exp);
696
697         --  A special case, if the number we had was a negative power of two,
698         --  then we want to add half of what we would otherwise add, since the
699         --  exponent is going to be reduced.
700
701         --  Note that X_Frac has the same sign as X, so if X_Frac is -0.5,
702         --  then we know that we have a negative number (and hence a negative
703         --  power of 2).
704
705         if X_Frac = -0.5 then
706            return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
707
708         --  Otherwise the exponent is unchanged
709
710         else
711            return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa);
712         end if;
713      end if;
714   end Succ;
715
716   ----------------
717   -- Truncation --
718   ----------------
719
720   --  The basic approach is to compute
721
722   --    T'Machine (RM1 + N) - RM1
723
724   --  where N >= 0.0 and RM1 = radix ** (mantissa - 1)
725
726   --  This works provided that the intermediate result (RM1 + N) does not
727   --  have extra precision (which is why we call Machine). When we compute
728   --  RM1 + N, the exponent of N will be normalized and the mantissa shifted
729   --  appropriately so the lower order bits, which cannot contribute to the
730   --  integer part of N, fall off on the right. When we subtract RM1 again,
731   --  the significant bits of N are shifted to the left, and what we have is
732   --  an integer, because only the first e bits are different from zero
733   --  (assuming binary radix here).
734
735   function Truncation (X : T) return T is
736      Result : T;
737
738   begin
739      Result := abs X;
740
741      if Result >= Radix_To_M_Minus_1 then
742         return T'Machine (X);
743
744      else
745         Result :=
746           T'Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1;
747
748         if Result > abs X then
749            Result := Result - 1.0;
750         end if;
751
752         if X > 0.0 then
753            return  Result;
754
755         elsif X < 0.0 then
756            return -Result;
757
758         --  For zero case, make sure sign of zero is preserved
759
760         else
761            return X;
762         end if;
763      end if;
764   end Truncation;
765
766   -----------------------
767   -- Unbiased_Rounding --
768   -----------------------
769
770   function Unbiased_Rounding (X : T) return T is
771      Abs_X  : constant T := abs X;
772      Result : T;
773      Tail   : T;
774
775   begin
776      Result := Truncation (Abs_X);
777      Tail   := Abs_X - Result;
778
779      if Tail > 0.5 then
780         Result := Result + 1.0;
781
782      elsif Tail = 0.5 then
783         Result := 2.0 * Truncation ((Result / 2.0) + 0.5);
784      end if;
785
786      if X > 0.0 then
787         return Result;
788
789      elsif X < 0.0 then
790         return -Result;
791
792      --  For zero case, make sure sign of zero is preserved
793
794      else
795         return X;
796      end if;
797   end Unbiased_Rounding;
798
799   -----------
800   -- Valid --
801   -----------
802
803   function Valid (X : not null access T) return Boolean is
804      IEEE_Emin : constant Integer := T'Machine_Emin - 1;
805      IEEE_Emax : constant Integer := T'Machine_Emax - 1;
806
807      IEEE_Bias : constant Integer := -(IEEE_Emin - 1);
808
809      subtype IEEE_Exponent_Range is
810        Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
811
812      --  The implementation of this floating point attribute uses a
813      --  representation type Float_Rep that allows direct access to the
814      --  exponent and mantissa parts of a floating point number.
815
816      --  The Float_Rep type is an array of Float_Word elements. This
817      --  representation is chosen to make it possible to size the type based
818      --  on a generic parameter. Since the array size is known at compile
819      --  time, efficient code can still be generated. The size of Float_Word
820      --  elements should be large enough to allow accessing the exponent in
821      --  one read, but small enough so that all floating point object sizes
822      --  are a multiple of the Float_Word'Size.
823
824      --  The following conditions must be met for all possible instantiations
825      --  of the attributes package:
826
827      --    - T'Size is an integral multiple of Float_Word'Size
828
829      --    - The exponent and sign are completely contained in a single
830      --      component of Float_Rep, named Most_Significant_Word (MSW).
831
832      --    - The sign occupies the most significant bit of the MSW and the
833      --      exponent is in the following bits. Unused bits (if any) are in
834      --      the least significant part.
835
836      type Float_Word is mod 2**Positive'Min (System.Word_Size, 32);
837      type Rep_Index is range 0 .. 7;
838
839      Rep_Words : constant Positive :=
840                    (T'Size + Float_Word'Size - 1) / Float_Word'Size;
841      Rep_Last  : constant Rep_Index :=
842                    Rep_Index'Min
843                      (Rep_Index (Rep_Words - 1),
844                       (T'Mantissa + 16) / Float_Word'Size);
845      --  Determine the number of Float_Words needed for representing the
846      --  entire floating-point value. Do not take into account excessive
847      --  padding, as occurs on IA-64 where 80 bits floats get padded to 128
848      --  bits. In general, the exponent field cannot be larger than 15 bits,
849      --  even for 128-bit floating-point types, so the final format size
850      --  won't be larger than T'Mantissa + 16.
851
852      type Float_Rep is
853         array (Rep_Index range 0 .. Rep_Index (Rep_Words - 1)) of Float_Word;
854
855      pragma Suppress_Initialization (Float_Rep);
856      --  This pragma suppresses the generation of an initialization procedure
857      --  for type Float_Rep when operating in Initialize/Normalize_Scalars
858      --  mode. This is not just a matter of efficiency, but of functionality,
859      --  since Valid has a pragma Inline_Always, which is not permitted if
860      --  there are nested subprograms present.
861
862      Most_Significant_Word : constant Rep_Index :=
863                                Rep_Last * Standard'Default_Bit_Order;
864      --  Finding the location of the Exponent_Word is a bit tricky. In general
865      --  we assume Word_Order = Bit_Order.
866
867      Exponent_Factor : constant Float_Word :=
868                          2**(Float_Word'Size - 1) /
869                            Float_Word (IEEE_Emax - IEEE_Emin + 3) *
870                              Boolean'Pos (Most_Significant_Word /= 2) +
871                                Boolean'Pos (Most_Significant_Word = 2);
872      --  Factor that the extracted exponent needs to be divided by to be in
873      --  range 0 .. IEEE_Emax - IEEE_Emin + 2. Special case: Exponent_Factor
874      --  is 1 for x86/IA64 double extended (GCC adds unused bits to the type).
875
876      Exponent_Mask : constant Float_Word :=
877                        Float_Word (IEEE_Emax - IEEE_Emin + 2) *
878                          Exponent_Factor;
879      --  Value needed to mask out the exponent field. This assumes that the
880      --  range IEEE_Emin - 1 .. IEEE_Emax + contains 2**N values, for some N
881      --  in Natural.
882
883      function To_Float is new Ada.Unchecked_Conversion (Float_Rep, T);
884
885      type Float_Access is access all T;
886      function To_Address is
887         new Ada.Unchecked_Conversion (Float_Access, System.Address);
888
889      XA : constant System.Address := To_Address (Float_Access (X));
890
891      R : Float_Rep;
892      pragma Import (Ada, R);
893      for R'Address use XA;
894      --  R is a view of the input floating-point parameter. Note that we
895      --  must avoid copying the actual bits of this parameter in float
896      --  form (since it may be a signalling NaN).
897
898      E  : constant IEEE_Exponent_Range :=
899             Integer ((R (Most_Significant_Word) and Exponent_Mask) /
900                                                        Exponent_Factor)
901               - IEEE_Bias;
902      --  Mask/Shift T to only get bits from the exponent. Then convert biased
903      --  value to integer value.
904
905      SR : Float_Rep;
906      --  Float_Rep representation of significant of X.all
907
908   begin
909      if T'Denorm then
910
911         --  All denormalized numbers are valid, so the only invalid numbers
912         --  are overflows and NaNs, both with exponent = Emax + 1.
913
914         return E /= IEEE_Emax + 1;
915
916      end if;
917
918      --  All denormalized numbers except 0.0 are invalid
919
920      --  Set exponent of X to zero, so we end up with the significand, which
921      --  definitely is a valid number and can be converted back to a float.
922
923      SR := R;
924      SR (Most_Significant_Word) :=
925           (SR (Most_Significant_Word)
926             and not Exponent_Mask) + Float_Word (IEEE_Bias) * Exponent_Factor;
927
928      return (E in IEEE_Emin .. IEEE_Emax) or else
929         ((E = IEEE_Emin - 1) and then abs To_Float (SR) = 1.0);
930   end Valid;
931
932end System.Fat_Gen;
933