OSDN Git Service

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