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