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