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