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 --                            $Revision: 1.1 $
  10 --                                                                          --
  11 --          Copyright (C) 1992-2001 Free Software Foundation, Inc.          --
  12 --                                                                          --
  13 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  14 -- terms of the  GNU General Public License as published  by the Free Soft- --
  15 -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  16 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  17 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  18 -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  19 -- for  more details.  You should have  received  a copy of the GNU General --
  20 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  21 -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  22 -- MA 02111-1307, USA.                                                      --
  23 --                                                                          --
  24 -- As a special exception,  if other files  instantiate  generics from this --
  25 -- unit, or you link  this unit with other files  to produce an executable, --
  26 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  27 -- covered  by the  GNU  General  Public  License.  This exception does not --
  28 -- however invalidate  any other reasons why  the executable file  might be --
  29 -- covered by the  GNU Public License.                                      --
  30 --                                                                          --
  31 -- GNAT was originally developed  by the GNAT team at  New York University. --
  32 -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
  33 --                                                                          --
  34 ------------------------------------------------------------------------------
  35
  36 --  The implementation here is portable to any IEEE implementation. It does
  37 --  not handle non-binary radix, and also assumes that model numbers and
  38 --  machine numbers are basically identical, which is not true of all possible
  39 --  floating-point implementations. On a non-IEEE machine, this body must be
  40 --  specialized appropriately, or better still, its generic instantiations
  41 --  should be replaced by efficient machine-specific code.
  42
  43 with Ada.Unchecked_Conversion; use Ada;
  44 with System;
  45 package body System.Fat_Gen is
  46
  47    Float_Radix        : constant T := T (T'Machine_Radix);
  48    Float_Radix_Inv    : constant T := 1.0 / Float_Radix;
  49    Radix_To_M_Minus_1 : constant T := Float_Radix ** (T'Machine_Mantissa - 1);
  50
  51    pragma Assert (T'Machine_Radix = 2);
  52    --  This version does not handle radix 16
  53
  54    --  Constants for Decompose and Scaling
  55
  56    Rad    : constant T := T (T'Machine_Radix);
  57    Invrad : constant T := 1.0 / Rad;
  58
  59    subtype Expbits is Integer range 0 .. 6;
  60    --  2 ** (2 ** 7) might overflow.  how big can radix-16 exponents get?
  61
  62    Log_Power : constant array (Expbits) of Integer := (1, 2, 4, 8, 16, 32, 64);
  63
  64    R_Power : constant array (Expbits) of T :=
  65      (Rad **  1,
  66       Rad **  2,
  67       Rad **  4,
  68       Rad **  8,
  69       Rad ** 16,
  70       Rad ** 32,
  71       Rad ** 64);
  72
  73    R_Neg_Power : constant array (Expbits) of T :=
  74      (Invrad **  1,
  75       Invrad **  2,
  76       Invrad **  4,
  77       Invrad **  8,
  78       Invrad ** 16,
  79       Invrad ** 32,
  80       Invrad ** 64);
  81
  82    -----------------------
  83    -- Local Subprograms --
  84    -----------------------
  85
  86    procedure Decompose (XX : T; Frac : out T; Expo : out UI);
  87    --  Decomposes a floating-point number into fraction and exponent parts
  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
 103       elsif Towards > X then
 104          return Succ (X);
 105
 106       else
 107          return Pred (X);
 108       end if;
 109    end Adjacent;
 110
 111    -------------
 112    -- Ceiling --
 113    -------------
 114
 115    function Ceiling (X : T) return T is
 116       XT : constant T := Truncation (X);
 117
 118    begin
 119       if X <= 0.0 then
 120          return XT;
 121
 122       elsif X = XT then
 123          return X;
 124
 125       else
 126          return XT + 1.0;
 127       end if;
 128    end Ceiling;
 129
 130    -------------
 131    -- Compose --
 132    -------------
 133
 134    function Compose (Fraction : T; Exponent : UI) return T is
 135       Arg_Frac : T;
 136       Arg_Exp  : UI;
 137
 138    begin
 139       Decompose (Fraction, Arg_Frac, Arg_Exp);
 140       return Scaling (Arg_Frac, Exponent);
 141    end Compose;
 142
 143    ---------------
 144    -- Copy_Sign --
 145    ---------------
 146
 147    function Copy_Sign (Value, Sign : T) return T is
 148       Result : T;
 149
 150       function Is_Negative (V : T) return Boolean;
 151       pragma Import (Intrinsic, Is_Negative);
 152
 153    begin
 154       Result := abs Value;
 155
 156       if Is_Negative (Sign) then
 157          return -Result;
 158       else
 159          return Result;
 160       end if;
 161    end Copy_Sign;
 162
 163    ---------------
 164    -- Decompose --
 165    ---------------
 166
 167    procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
 168       X : T := T'Machine (XX);
 169
 170    begin
 171       if X = 0.0 then
 172          Frac := X;
 173          Expo := 0;
 174
 175          --  More useful would be defining Expo to be T'Machine_Emin - 1 or
 176          --  T'Machine_Emin - T'Machine_Mantissa, which would preserve
 177          --  monotonicity of the exponent function ???
 178
 179       --  Check for infinities, transfinites, whatnot.
 180
 181       elsif X > T'Safe_Last then
 182          Frac := Invrad;
 183          Expo := T'Machine_Emax + 1;
 184
 185       elsif X < T'Safe_First then
 186          Frac := -Invrad;
 187          Expo := T'Machine_Emax + 2;    -- how many extra negative values?
 188
 189       else
 190          --  Case of nonzero finite x. Essentially, we just multiply
 191          --  by Rad ** (+-2**N) to reduce the range.
 192
 193          declare
 194             Ax : T  := abs X;
 195             Ex : UI := 0;
 196
 197          --  Ax * Rad ** Ex is invariant.
 198
 199          begin
 200             if Ax >= 1.0 then
 201                while Ax >= R_Power (Expbits'Last) loop
 202                   Ax := Ax * R_Neg_Power (Expbits'Last);
 203                   Ex := Ex + Log_Power (Expbits'Last);
 204                end loop;
 205
 206                --  Ax < Rad ** 64
 207
 208                for N in reverse Expbits'First .. Expbits'Last - 1 loop
 209                   if Ax >= R_Power (N) then
 210                      Ax := Ax * R_Neg_Power (N);
 211                      Ex := Ex + Log_Power (N);
 212                   end if;
 213
 214                   --  Ax < R_Power (N)
 215                end loop;
 216
 217                --  1 <= Ax < Rad
 218
 219                Ax := Ax * Invrad;
 220                Ex := Ex + 1;
 221
 222             else
 223                --  0 < ax < 1
 224
 225                while Ax < R_Neg_Power (Expbits'Last) loop
 226                   Ax := Ax * R_Power (Expbits'Last);
 227                   Ex := Ex - Log_Power (Expbits'Last);
 228                end loop;
 229
 230                --  Rad ** -64 <= Ax < 1
 231
 232                for N in reverse Expbits'First .. Expbits'Last - 1 loop
 233                   if Ax < R_Neg_Power (N) then
 234                      Ax := Ax * R_Power (N);
 235                      Ex := Ex - Log_Power (N);
 236                   end if;
 237
 238                   --  R_Neg_Power (N) <= Ax < 1
 239                end loop;
 240             end if;
 241
 242             if X > 0.0 then
 243                Frac := Ax;
 244             else
 245                Frac := -Ax;
 246             end if;
 247
 248             Expo := Ex;
 249          end;
 250       end if;
 251    end Decompose;
 252
 253    --------------
 254    -- Exponent --
 255    --------------
 256
 257    function Exponent (X : T) return UI is
 258       X_Frac : T;
 259       X_Exp  : UI;
 260
 261    begin
 262       Decompose (X, X_Frac, X_Exp);
 263       return X_Exp;
 264    end Exponent;
 265
 266    -----------
 267    -- Floor --
 268    -----------
 269
 270    function Floor (X : T) return T is
 271       XT : constant T := Truncation (X);
 272
 273    begin
 274       if X >= 0.0 then
 275          return XT;
 276
 277       elsif XT = X then
 278          return X;
 279
 280       else
 281          return XT - 1.0;
 282       end if;
 283    end Floor;
 284
 285    --------------
 286    -- Fraction --
 287    --------------
 288
 289    function Fraction (X : T) return T is
 290       X_Frac : T;
 291       X_Exp  : UI;
 292
 293    begin
 294       Decompose (X, X_Frac, X_Exp);
 295       return X_Frac;
 296    end Fraction;
 297
 298    ---------------------
 299    -- Gradual_Scaling --
 300    ---------------------
 301
 302    function Gradual_Scaling  (Adjustment : UI) return T is
 303       Y  : T;
 304       Y1 : T;
 305       Ex : UI := Adjustment;
 306
 307    begin
 308       if Adjustment < T'Machine_Emin then
 309          Y  := 2.0 ** T'Machine_Emin;
 310          Y1 := Y;
 311          Ex := Ex - T'Machine_Emin;
 312
 313          while Ex <= 0 loop
 314             Y := T'Machine (Y / 2.0);
 315
 316             if Y = 0.0 then
 317                return Y1;
 318             end if;
 319
 320             Ex := Ex + 1;
 321             Y1 := Y;
 322          end loop;
 323
 324          return Y1;
 325
 326       else
 327          return Scaling (1.0, Adjustment);
 328       end if;
 329    end Gradual_Scaling;
 330
 331    ------------------
 332    -- Leading_Part --
 333    ------------------
 334
 335    function Leading_Part (X : T; Radix_Digits : UI) return T is
 336       L    : UI;
 337       Y, Z : T;
 338
 339    begin
 340       if Radix_Digits >= T'Machine_Mantissa then
 341          return X;
 342
 343       else
 344          L := Exponent (X) - Radix_Digits;
 345          Y := Truncation (Scaling (X, -L));
 346          Z := Scaling (Y, L);
 347          return Z;
 348       end if;
 349
 350    end Leading_Part;
 351
 352    -------------
 353    -- Machine --
 354    -------------
 355
 356    --  The trick with Machine is to force the compiler to store the result
 357    --  in memory so that we do not have extra precision used. The compiler
 358    --  is clever, so we have to outwit its possible optimizations! We do
 359    --  this by using an intermediate pragma Volatile location.
 360
 361    function Machine (X : T) return T is
 362       Temp : T;
 363       pragma Volatile (Temp);
 364
 365    begin
 366       Temp := X;
 367       return Temp;
 368    end Machine;
 369
 370    -----------
 371    -- Model --
 372    -----------
 373
 374    --  We treat Model as identical to Machine. This is true of IEEE and other
 375    --  nice floating-point systems, but not necessarily true of all systems.
 376
 377    function Model (X : T) return T is
 378    begin
 379       return Machine (X);
 380    end Model;
 381
 382    ----------
 383    -- Pred --
 384    ----------
 385
 386    --  Subtract from the given number a number equivalent to the value of its
 387    --  least significant bit. Given that the most significant bit represents
 388    --  a value of 1.0 * radix ** (exp - 1), the value we want is obtained by
 389    --  shifting this by (mantissa-1) bits to the right, i.e. decreasing the
 390    --  exponent by that amount.
 391
 392    --  Zero has to be treated specially, since its exponent is zero
 393
 394    function Pred (X : T) return T is
 395       X_Frac : T;
 396       X_Exp  : UI;
 397
 398    begin
 399       if X = 0.0 then
 400          return -Succ (X);
 401
 402       else
 403          Decompose (X, X_Frac, X_Exp);
 404
 405          --  A special case, if the number we had was a positive power of
 406          --  two, then we want to subtract half of what we would otherwise
 407          --  subtract, since the exponent is going to be reduced.
 408
 409          if X_Frac = 0.5 and then X > 0.0 then
 410             return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
 411
 412          --  Otherwise the exponent stays the same
 413
 414          else
 415             return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa);
 416          end if;
 417       end if;
 418    end Pred;
 419
 420    ---------------
 421    -- Remainder --
 422    ---------------
 423
 424    function Remainder (X, Y : T) return T is
 425       A        : T;
 426       B        : T;
 427       Arg      : T;
 428       P        : T;
 429       Arg_Frac : T;
 430       P_Frac   : T;
 431       Sign_X   : T;
 432       IEEE_Rem : T;
 433       Arg_Exp  : UI;
 434       P_Exp    : UI;
 435       K        : UI;
 436       P_Even   : Boolean;
 437
 438    begin
 439       if X > 0.0 then
 440          Sign_X :=  1.0;
 441          Arg := X;
 442       else
 443          Sign_X := -1.0;
 444          Arg := -X;
 445       end if;
 446
 447       P := abs Y;
 448
 449       if Arg < P then
 450          P_Even := True;
 451          IEEE_Rem := Arg;
 452          P_Exp := Exponent (P);
 453
 454       else
 455          Decompose (Arg, Arg_Frac, Arg_Exp);
 456          Decompose (P,   P_Frac,   P_Exp);
 457
 458          P := Compose (P_Frac, Arg_Exp);
 459          K := Arg_Exp - P_Exp;
 460          P_Even := True;
 461          IEEE_Rem := Arg;
 462
 463          for Cnt in reverse 0 .. K loop
 464             if IEEE_Rem >= P then
 465                P_Even := False;
 466                IEEE_Rem := IEEE_Rem - P;
 467             else
 468                P_Even := True;
 469             end if;
 470
 471             P := P * 0.5;
 472          end loop;
 473       end if;
 474
 475       --  That completes the calculation of modulus remainder. The final
 476       --  step is get the IEEE remainder. Here we need to compare Rem with
 477       --  (abs Y) / 2. We must be careful of unrepresentable Y/2 value
 478       --  caused by subnormal numbers
 479
 480       if P_Exp >= 0 then
 481          A := IEEE_Rem;
 482          B := abs Y * 0.5;
 483
 484       else
 485          A := IEEE_Rem * 2.0;
 486          B := abs Y;
 487       end if;
 488
 489       if A > B or else (A = B and then not P_Even) then
 490          IEEE_Rem := IEEE_Rem - abs Y;
 491       end if;
 492
 493       return Sign_X * IEEE_Rem;
 494
 495    end Remainder;
 496
 497    --------------
 498    -- Rounding --
 499    --------------
 500
 501    function Rounding (X : T) return T is
 502       Result : T;
 503       Tail   : T;
 504
 505    begin
 506       Result := Truncation (abs X);
 507       Tail   := abs X - Result;
 508
 509       if Tail >= 0.5  then
 510          Result := Result + 1.0;
 511       end if;
 512
 513       if X > 0.0 then
 514          return Result;
 515
 516       elsif X < 0.0 then
 517          return -Result;
 518
 519       --  For zero case, make sure sign of zero is preserved
 520
 521       else
 522          return X;
 523       end if;
 524
 525    end Rounding;
 526
 527    -------------
 528    -- Scaling --
 529    -------------
 530
 531    --  Return x * rad ** adjustment quickly,
 532    --  or quietly underflow to zero, or overflow naturally.
 533
 534    function Scaling (X : T; Adjustment : UI) return T is
 535    begin
 536       if X = 0.0 or else Adjustment = 0 then
 537          return X;
 538       end if;
 539
 540       --  Nonzero x. essentially, just multiply repeatedly by Rad ** (+-2**n).
 541
 542       declare
 543          Y  : T  := X;
 544          Ex : UI := Adjustment;
 545
 546       --  Y * Rad ** Ex is invariant
 547
 548       begin
 549          if Ex < 0 then
 550             while Ex <= -Log_Power (Expbits'Last) loop
 551                Y := Y * R_Neg_Power (Expbits'Last);
 552                Ex := Ex + Log_Power (Expbits'Last);
 553             end loop;
 554
 555             --  -64 < Ex <= 0
 556
 557             for N in reverse Expbits'First .. Expbits'Last - 1 loop
 558                if Ex <= -Log_Power (N) then
 559                   Y := Y * R_Neg_Power (N);
 560                   Ex := Ex + Log_Power (N);
 561                end if;
 562
 563                --  -Log_Power (N) < Ex <= 0
 564             end loop;
 565
 566             --  Ex = 0
 567
 568          else
 569             --  Ex >= 0
 570
 571             while Ex >= Log_Power (Expbits'Last) loop
 572                Y := Y * R_Power (Expbits'Last);
 573                Ex := Ex - Log_Power (Expbits'Last);
 574             end loop;
 575
 576             --  0 <= Ex < 64
 577
 578             for N in reverse Expbits'First .. Expbits'Last - 1 loop
 579                if Ex >= Log_Power (N) then
 580                   Y := Y * R_Power (N);
 581                   Ex := Ex - Log_Power (N);
 582                end if;
 583
 584                --  0 <= Ex < Log_Power (N)
 585             end loop;
 586
 587             --  Ex = 0
 588          end if;
 589          return Y;
 590       end;
 591    end Scaling;
 592
 593    ----------
 594    -- Succ --
 595    ----------
 596
 597    --  Similar computation to that of Pred: find value of least significant
 598    --  bit of given number, and add. Zero has to be treated specially since
 599    --  the exponent can be zero, and also we want the smallest denormal if
 600    --  denormals are supported.
 601
 602    function Succ (X : T) return T is
 603       X_Frac : T;
 604       X_Exp  : UI;
 605       X1, X2 : T;
 606
 607    begin
 608       if X = 0.0 then
 609          X1 := 2.0 ** T'Machine_Emin;
 610
 611          --  Following loop generates smallest denormal
 612
 613          loop
 614             X2 := T'Machine (X1 / 2.0);
 615             exit when X2 = 0.0;
 616             X1 := X2;
 617          end loop;
 618
 619          return X1;
 620
 621       else
 622          Decompose (X, X_Frac, X_Exp);
 623
 624          --  A special case, if the number we had was a negative power of
 625          --  two, then we want to add half of what we would otherwise add,
 626          --  since the exponent is going to be reduced.
 627
 628          if X_Frac = 0.5 and then X < 0.0 then
 629             return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
 630
 631          --  Otherwise the exponent stays the same
 632
 633          else
 634             return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa);
 635          end if;
 636       end if;
 637    end Succ;
 638
 639    ----------------
 640    -- Truncation --
 641    ----------------
 642
 643    --  The basic approach is to compute
 644
 645    --    T'Machine (RM1 + N) - RM1.
 646
 647    --  where N >= 0.0 and RM1 = radix ** (mantissa - 1)
 648
 649    --  This works provided that the intermediate result (RM1 + N) does not
 650    --  have extra precision (which is why we call Machine). When we compute
 651    --  RM1 + N, the exponent of N will be normalized and the mantissa shifted
 652    --  shifted appropriately so the lower order bits, which cannot contribute
 653    --  to the integer part of N, fall off on the right. When we subtract RM1
 654    --  again, the significant bits of N are shifted to the left, and what we
 655    --  have is an integer, because only the first e bits are different from
 656    --  zero (assuming binary radix here).
 657
 658    function Truncation (X : T) return T is
 659       Result : T;
 660
 661    begin
 662       Result := abs X;
 663
 664       if Result >= Radix_To_M_Minus_1 then
 665          return Machine (X);
 666
 667       else
 668          Result := Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1;
 669
 670          if Result > abs X  then
 671             Result := Result - 1.0;
 672          end if;
 673
 674          if X > 0.0 then
 675             return  Result;
 676
 677          elsif X < 0.0 then
 678             return -Result;
 679
 680          --  For zero case, make sure sign of zero is preserved
 681
 682          else
 683             return X;
 684          end if;
 685       end if;
 686
 687    end Truncation;
 688
 689    -----------------------
 690    -- Unbiased_Rounding --
 691    -----------------------
 692
 693    function Unbiased_Rounding (X : T) return T is
 694       Abs_X  : constant T := abs X;
 695       Result : T;
 696       Tail   : T;
 697
 698    begin
 699       Result := Truncation (Abs_X);
 700       Tail   := Abs_X - Result;
 701
 702       if Tail > 0.5  then
 703          Result := Result + 1.0;
 704
 705       elsif Tail = 0.5 then
 706          Result := 2.0 * Truncation ((Result / 2.0) + 0.5);
 707       end if;
 708
 709       if X > 0.0 then
 710          return Result;
 711
 712       elsif X < 0.0 then
 713          return -Result;
 714
 715       --  For zero case, make sure sign of zero is preserved
 716
 717       else
 718          return X;
 719       end if;
 720
 721    end Unbiased_Rounding;
 722
 723    -----------
 724    -- Valid --
 725    -----------
 726
 727    function Valid (X : access T) return Boolean is
 728
 729       IEEE_Emin : constant Integer := T'Machine_Emin - 1;
 730       IEEE_Emax : constant Integer := T'Machine_Emax - 1;
 731
 732       IEEE_Bias : constant Integer := -(IEEE_Emin - 1);
 733
 734       subtype IEEE_Exponent_Range is
 735         Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
 736
 737       --  The implementation of this floating point attribute uses
 738       --  a representation type Float_Rep that allows direct access to
 739       --  the exponent and mantissa parts of a floating point number.
 740
 741       --  The Float_Rep type is an array of Float_Word elements. This
 742       --  representation is chosen to make it possible to size the
 743       --  type based on a generic parameter.
 744
 745       --  The following conditions must be met for all possible
 746       --  instantiations of the attributes package:
 747
 748       --    - T'Size is an integral multiple of Float_Word'Size
 749
 750       --    - The exponent and sign are completely contained in a single
 751       --      component of Float_Rep, named Most_Significant_Word (MSW).
 752
 753       --    - The sign occupies the most significant bit of the MSW
 754       --      and the exponent is in the following bits.
 755       --      Unused bits (if any) are in the least significant part.
 756
 757       type Float_Word is mod 2**32;
 758       type Rep_Index is range 0 .. 7;
 759
 760       Rep_Last : constant Rep_Index := (T'Size - 1) / Float_Word'Size;
 761
 762       type Float_Rep is array (Rep_Index range 0 .. Rep_Last) of Float_Word;
 763
 764       Most_Significant_Word : constant Rep_Index :=
 765                                 Rep_Last * Standard'Default_Bit_Order;
 766       --  Finding the location of the Exponent_Word is a bit tricky.
 767       --  In general we assume Word_Order = Bit_Order.
 768       --  This expression needs to be refined for VMS.
 769
 770       Exponent_Factor : constant Float_Word :=
 771                           2**(Float_Word'Size - 1) /
 772                             Float_Word (IEEE_Emax - IEEE_Emin + 3) *
 773                               Boolean'Pos (T'Size /= 96) +
 774                                 Boolean'Pos (T'Size = 96);
 775       --  Factor that the extracted exponent needs to be divided by
 776       --  to be in range 0 .. IEEE_Emax - IEEE_Emin + 2.
 777       --  Special kludge: Exponent_Factor is 0 for x86 double extended
 778       --  as GCC adds 16 unused bits to the type.
 779
 780       Exponent_Mask : constant Float_Word :=
 781                         Float_Word (IEEE_Emax - IEEE_Emin + 2) *
 782                           Exponent_Factor;
 783       --  Value needed to mask out the exponent field.
 784       --  This assumes that the range IEEE_Emin - 1 .. IEEE_Emax + 1
 785       --  contains 2**N values, for some N in Natural.
 786
 787       function To_Float is new Unchecked_Conversion (Float_Rep, T);
 788
 789       type Float_Access is access all T;
 790       function To_Address is
 791          new Unchecked_Conversion (Float_Access, System.Address);
 792
 793       XA : constant System.Address := To_Address (Float_Access (X));
 794
 795       R : Float_Rep;
 796       pragma Import (Ada, R);
 797       for R'Address use XA;
 798       --  R is a view of the input floating-point parameter. Note that we
 799       --  must avoid copying the actual bits of this parameter in float
 800       --  form (since it may be a signalling NaN.
 801
 802       E  : constant IEEE_Exponent_Range :=
 803              Integer ((R (Most_Significant_Word) and Exponent_Mask) /
 804                                                         Exponent_Factor)
 805                - IEEE_Bias;
 806       --  Mask/Shift T to only get bits from the exponent
 807       --  Then convert biased value to integer value.
 808
 809       SR : Float_Rep;
 810       --  Float_Rep representation of significant of X.all
 811
 812    begin
 813       if T'Denorm then
 814
 815          --  All denormalized numbers are valid, so only invalid numbers
 816          --  are overflows and NaN's, both with exponent = Emax + 1.
 817
 818          return E /= IEEE_Emax + 1;
 819
 820       end if;
 821
 822       --  All denormalized numbers except 0.0 are invalid
 823
 824       --  Set exponent of X to zero, so we end up with the significand, which
 825       --  definitely is a valid number and can be converted back to a float.
 826
 827       SR := R;
 828       SR (Most_Significant_Word) :=
 829            (SR (Most_Significant_Word)
 830              and not Exponent_Mask) + Float_Word (IEEE_Bias) * Exponent_Factor;
 831
 832       return (E in IEEE_Emin .. IEEE_Emax) or else
 833          ((E = IEEE_Emin - 1) and then abs To_Float (SR) = 1.0);
 834    end Valid;
 835
 836 end System.Fat_Gen;