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