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