gcc/ada/libgnat/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-2023, 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 --  This implementation is portable to any IEEE implementation. It does not
  33 --  handle nonbinary radix, and also assumes that model numbers and machine
  34 --  numbers are basically identical, which is not true of all possible
  35 --  floating-point implementations.
  36
  37 with Ada.Unchecked_Conversion;
  38 with System.Unsigned_Types;
  39
  40 pragma Warnings (Off, "non-static constant in preelaborated unit");
  41 --  Every constant is static given our instantiation model
  42
  43 package body System.Fat_Gen is
  44    pragma Assert (T'Machine_Radix = 2);
  45    --  This version does not handle radix 16
  46
  47    Rad : constant T := T (T'Machine_Radix);
  48    --  Renaming for the machine radix
  49
  50    Mantissa : constant Integer := T'Machine_Mantissa;
  51    --  Renaming for the machine mantissa
  52
  53    Invrad : constant T := 1.0 / Rad;
  54    --  Smallest positive mantissa in the canonical form (RM A.5.3(4))
  55
  56    --  Small : constant T := Rad ** (T'Machine_Emin - 1);
  57    --  Smallest positive normalized number
  58
  59    --  Tiny : constant T := Rad ** (T'Machine_Emin - Mantissa);
  60    --  Smallest positive denormalized number
  61
  62    RM1 : constant T := Rad ** (Mantissa - 1);
  63    --  Smallest positive member of the large consecutive integers. It is equal
  64    --  to the ratio Small / Tiny, which means that multiplying by it normalizes
  65    --  any nonzero denormalized number.
  66
  67    IEEE_Emin : constant Integer := T'Machine_Emin - 1;
  68    IEEE_Emax : constant Integer := T'Machine_Emax - 1;
  69    --  The mantissa is a fraction with first digit set in Ada whereas it is
  70    --  shifted by 1 digit to the left in the IEEE floating-point format.
  71
  72    subtype IEEE_Erange is Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
  73    --  The IEEE floating-point format extends the machine range by 1 to the
  74    --  left for denormalized numbers and 1 to the right for infinities/NaNs.
  75
  76    IEEE_Ebias : constant Integer := -(IEEE_Emin - 1);
  77    --  The exponent is biased such that denormalized numbers have it zero
  78
  79    --  The implementation uses a representation type Float_Rep that allows
  80    --  direct access to exponent and mantissa of the floating point number.
  81
  82    --  The Float_Rep type is a simple array of Float_Word elements. This
  83    --  representation is chosen to make it possible to size the type based
  84    --  on a generic parameter. Since the array size is known at compile
  85    --  time, efficient code can still be generated. The size of Float_Word
  86    --  elements should be large enough to allow accessing the exponent in
  87    --  one read, but small enough so that all floating-point object sizes
  88    --  are a multiple of Float_Word'Size.
  89
  90    --  The following conditions must be met for all possible instantiations
  91    --  of the attribute package:
  92
  93    --    - T'Size is an integral multiple of Float_Word'Size
  94
  95    --    - The exponent and sign are completely contained in a single
  96    --      component of Float_Rep, named Most Significant Word (MSW).
  97
  98    --    - The sign occupies the most significant bit of the MSW and the
  99    --      exponent is in the following bits.
 100
 101    --  The low-level primitives Copy_Sign, Decompose, Finite_Succ, Scaling and
 102    --  Valid are implemented by accessing the bit pattern of the floating-point
 103    --  number. Only the normalization of denormalized numbers, if any, and the
 104    --  gradual underflow are left to the hardware, mainly because there is some
 105    --  leeway for the hardware implementation in this area: for example the MSB
 106    --  of the mantissa, that is 1 for normalized numbers and 0 for denormalized
 107    --  numbers, may or may not be stored by the hardware.
 108
 109    Siz : constant := 16;
 110    type Float_Word is mod 2**Siz;
 111    --  We use the GCD of the size of all the supported floating-point formats
 112
 113    N  : constant Natural := (T'Size + Siz - 1) / Siz;
 114    NR : constant Natural := (Mantissa + 16 + Siz - 1) / Siz;
 115    Rep_Last : constant Natural := Natural'Min (N, NR) - 1;
 116    --  Determine the number of Float_Words needed for representing the
 117    --  entire floating-point value. Do not take into account excessive
 118    --  padding, as occurs on IA-64 where 80 bits floats get padded to 128
 119    --  bits. In general, the exponent field cannot be larger than 15 bits,
 120    --  even for 128-bit floating-point types, so the final format size
 121    --  won't be larger than Mantissa + 16.
 122
 123    type Float_Rep is array (Natural range 0 .. N - 1) of Float_Word;
 124    pragma Suppress_Initialization (Float_Rep);
 125    --  This pragma suppresses the generation of an initialization procedure
 126    --  for type Float_Rep when operating in Initialize/Normalize_Scalars mode.
 127
 128    MSW : constant Natural := Rep_Last * Standard'Default_Bit_Order;
 129    --  Finding the location of the Exponent_Word is a bit tricky. In general
 130    --  we assume Word_Order = Bit_Order.
 131
 132    Exp_Factor : constant Float_Word :=
 133                   2**(Siz - 1) / Float_Word (IEEE_Emax - IEEE_Emin + 3);
 134    --  Factor that the extracted exponent needs to be divided by to be in
 135    --  range 0 .. IEEE_Emax - IEEE_Emin + 2
 136
 137    Exp_Mask : constant Float_Word :=
 138                 Float_Word (IEEE_Emax - IEEE_Emin + 2) * Exp_Factor;
 139    --  Value needed to mask out the exponent field. This assumes that the
 140    --  range 0 .. IEEE_Emax - IEEE_Emin + 2 contains 2**N values, for some
 141    --  N in Natural.
 142
 143    Sign_Mask : constant Float_Word := 2**(Siz - 1);
 144    --  Value needed to mask out the sign field
 145
 146    -----------------------
 147    -- Local Subprograms --
 148    -----------------------
 149
 150    procedure Decompose (XX : T; Frac : out T; Expo : out UI);
 151    --  Decomposes a floating-point number into fraction and exponent parts.
 152    --  Both results are signed, with Frac having the sign of XX, and UI has
 153    --  the sign of the exponent. The absolute value of Frac is in the range
 154    --  0.0 <= Frac < 1.0. If Frac = 0.0 or -0.0, then Expo is always zero.
 155
 156    function Finite_Succ (X : T) return T;
 157    --  Return the successor of X, a finite number not equal to T'Last
 158
 159    --------------
 160    -- Adjacent --
 161    --------------
 162
 163    function Adjacent (X, Towards : T) return T is
 164    begin
 165       if Towards = X then
 166          return X;
 167       elsif Towards > X then
 168          return Succ (X);
 169       else
 170          return Pred (X);
 171       end if;
 172    end Adjacent;
 173
 174    -------------
 175    -- Ceiling --
 176    -------------
 177
 178    function Ceiling (X : T) return T is
 179       XT : constant T := Truncation (X);
 180    begin
 181       if X <= 0.0 then
 182          return XT;
 183       elsif X = XT then
 184          return X;
 185       else
 186          return XT + 1.0;
 187       end if;
 188    end Ceiling;
 189
 190    -------------
 191    -- Compose --
 192    -------------
 193
 194    function Compose (Fraction : T; Exponent : UI) return T is
 195       Arg_Frac : T;
 196       Arg_Exp  : UI;
 197    begin
 198       Decompose (Fraction, Arg_Frac, Arg_Exp);
 199       return Scaling (Arg_Frac, Exponent);
 200    end Compose;
 201
 202    ---------------
 203    -- Copy_Sign --
 204    ---------------
 205
 206    function Copy_Sign (Value, Sign : T) return T is
 207       S : constant T := T'Machine (Sign);
 208
 209       Rep_S : Float_Rep;
 210       for Rep_S'Address use S'Address;
 211       --  Rep_S is a view of the Sign parameter
 212
 213       V : T := T'Machine (Value);
 214
 215       Rep_V : Float_Rep;
 216       for Rep_V'Address use V'Address;
 217       --  Rep_V is a view of the Value parameter
 218
 219    begin
 220       Rep_V (MSW) :=
 221         (Rep_V (MSW) and not Sign_Mask) or (Rep_S (MSW) and Sign_Mask);
 222       return V;
 223    end Copy_Sign;
 224
 225    ---------------
 226    -- Decompose --
 227    ---------------
 228
 229    procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
 230       X : T := T'Machine (XX);
 231
 232       Rep : Float_Rep;
 233       for Rep'Address use X'Address;
 234       --  Rep is a view of the input floating-point parameter
 235
 236       Exp : constant IEEE_Erange :=
 237               Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias;
 238       --  Mask/Shift X to only get bits from the exponent. Then convert biased
 239       --  value to final value.
 240
 241       Minus : constant Boolean := (Rep (MSW) and Sign_Mask) /= 0;
 242       --  Mask/Shift X to only get bit from the sign
 243
 244    begin
 245       --  The normalized exponent of zero is zero, see RM A.5.3(15)
 246
 247       if X = 0.0 then
 248          Expo := 0;
 249          Frac := X;
 250
 251       --  Check for infinities and NaNs
 252
 253       elsif Exp = IEEE_Emax + 1 then
 254          Expo := T'Machine_Emax + 1;
 255          Frac := (if Minus then -Invrad else Invrad);
 256
 257       --  Check for nonzero denormalized numbers
 258
 259       elsif Exp = IEEE_Emin - 1 then
 260          --  Normalize by multiplying by Radix ** (Mantissa - 1)
 261
 262          Decompose (X * RM1, Frac, Expo);
 263          Expo := Expo - (Mantissa - 1);
 264
 265       --  Case of normalized numbers
 266
 267       else
 268          --  The Ada exponent is the IEEE exponent plus 1, see above
 269
 270          Expo := Exp + 1;
 271
 272          --  Set Ada exponent of X to zero, so we end up with the fraction
 273
 274          Rep (MSW) := (Rep (MSW) and not Exp_Mask) +
 275                         Float_Word (IEEE_Ebias - 1) * Exp_Factor;
 276          Frac := X;
 277       end if;
 278    end Decompose;
 279
 280    --------------
 281    -- Exponent --
 282    --------------
 283
 284    function Exponent (X : T) return UI is
 285       X_Frac : T;
 286       X_Exp  : UI;
 287    begin
 288       Decompose (X, X_Frac, X_Exp);
 289       return X_Exp;
 290    end Exponent;
 291
 292    -----------------
 293    -- Finite_Succ --
 294    -----------------
 295
 296    function Finite_Succ (X : T) return T is
 297       XX : T := T'Machine (X);
 298
 299       Rep : Float_Rep;
 300       for Rep'Address use XX'Address;
 301       --  Rep is a view of the input floating-point parameter
 302
 303    begin
 304       --  If the floating-point type does not support denormalized numbers,
 305       --  there is a couple of problematic values, namely -Small and Zero,
 306       --  because the increment is equal to Small in these cases.
 307
 308       if not T'Denorm then
 309          declare
 310             Small : constant T := Rad ** (T'Machine_Emin - 1);
 311             --  Smallest positive normalized number declared here and not at
 312             --  library level for the sake of the CCG compiler, which cannot
 313             --  currently compile the constant because the target is C90.
 314
 315          begin
 316             if X = -Small then
 317                XX := 0.0;
 318                return -XX;
 319             elsif X = 0.0 then
 320                return Small;
 321             end if;
 322          end;
 323       end if;
 324
 325       --  In all the other cases, the increment is equal to 1 in the binary
 326       --  integer representation of the number if X is nonnegative and equal
 327       --  to -1 if X is negative.
 328
 329       if XX >= 0.0 then
 330          --  First clear the sign of negative Zero
 331
 332          Rep (MSW) := Rep (MSW) and not Sign_Mask;
 333
 334          --  Deal with big endian
 335
 336          if MSW = 0 then
 337             for J in reverse 0 .. Rep_Last loop
 338                Rep (J) := Rep (J) + 1;
 339
 340                --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
 341                --  so, when it has been flipped, its status must be reanalyzed.
 342
 343                if Mantissa = 64 and then J = 1 then
 344
 345                   --  If the MSB changed from denormalized to normalized, then
 346                   --  keep it normalized since the exponent will be bumped.
 347
 348                   if Rep (J) = 2**(Siz - 1) then
 349                      null;
 350
 351                   --  If the MSB changed from normalized, restore it since we
 352                   --  cannot denormalize in this context.
 353
 354                   elsif Rep (J) = 0 then
 355                      Rep (J) := 2**(Siz - 1);
 356
 357                   else
 358                      exit;
 359                   end if;
 360
 361                --  In other cases, stop if there is no carry
 362
 363                else
 364                   exit when Rep (J) > 0;
 365                end if;
 366             end loop;
 367
 368          --  Deal with little endian
 369
 370          else
 371             for J in 0 .. Rep_Last loop
 372                Rep (J) := Rep (J) + 1;
 373
 374                --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
 375                --  so, when it has been flipped, its status must be reanalyzed.
 376
 377                if Mantissa = 64 and then J = Rep_Last - 1 then
 378
 379                   --  If the MSB changed from denormalized to normalized, then
 380                   --  keep it normalized since the exponent will be bumped.
 381
 382                   if Rep (J) = 2**(Siz - 1) then
 383                      null;
 384
 385                   --  If the MSB changed from normalized, restore it since we
 386                   --  cannot denormalize in this context.
 387
 388                   elsif Rep (J) = 0 then
 389                      Rep (J) := 2**(Siz - 1);
 390
 391                   else
 392                      exit;
 393                   end if;
 394
 395                --  In other cases, stop if there is no carry
 396
 397                else
 398                   exit when Rep (J) > 0;
 399                end if;
 400             end loop;
 401          end if;
 402
 403       else
 404          if MSW = 0 then
 405             for J in reverse 0 .. Rep_Last loop
 406                Rep (J) := Rep (J) - 1;
 407
 408                --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
 409                --  so, when it has been flipped, its status must be reanalyzed.
 410
 411                if Mantissa = 64 and then J = 1 then
 412
 413                   --  If the MSB changed from normalized to denormalized, then
 414                   --  keep it normalized if the exponent is not 1.
 415
 416                   if Rep (J) = 2**(Siz - 1) - 1 then
 417                      if Rep (0) /= 2**(Siz - 1) + 1 then
 418                         Rep (J) := 2**Siz - 1;
 419                      end if;
 420
 421                   else
 422                      exit;
 423                   end if;
 424
 425                --  In other cases, stop if there is no borrow
 426
 427                else
 428                   exit when Rep (J) < 2**Siz - 1;
 429                end if;
 430             end loop;
 431
 432          else
 433             for J in 0 .. Rep_Last loop
 434                Rep (J) := Rep (J) - 1;
 435
 436                --  For 80-bit IEEE Extended, the MSB of the mantissa is stored
 437                --  so, when it has been flipped, its status must be reanalyzed.
 438
 439                if Mantissa = 64 and then J = Rep_Last - 1 then
 440
 441                   --  If the MSB changed from normalized to denormalized, then
 442                   --  keep it normalized if the exponent is not 1.
 443
 444                   if Rep (J) = 2**(Siz - 1) - 1 then
 445                      if Rep (Rep_Last) /= 2**(Siz - 1) + 1 then
 446                         Rep (J) := 2**Siz - 1;
 447                      end if;
 448
 449                   else
 450                      exit;
 451                   end if;
 452
 453                --  In other cases, stop if there is no borrow
 454
 455                else
 456                   exit when Rep (J) < 2**Siz - 1;
 457                end if;
 458             end loop;
 459          end if;
 460       end if;
 461
 462       return XX;
 463    end Finite_Succ;
 464
 465    -----------
 466    -- Floor --
 467    -----------
 468
 469    function Floor (X : T) return T is
 470       XT : constant T := Truncation (X);
 471    begin
 472       if X >= 0.0 then
 473          return XT;
 474       elsif XT = X then
 475          return X;
 476       else
 477          return XT - 1.0;
 478       end if;
 479    end Floor;
 480
 481    --------------
 482    -- Fraction --
 483    --------------
 484
 485    function Fraction (X : T) return T is
 486       X_Frac : T;
 487       X_Exp  : UI;
 488    begin
 489       Decompose (X, X_Frac, X_Exp);
 490       return X_Frac;
 491    end Fraction;
 492
 493    ------------------
 494    -- Leading_Part --
 495    ------------------
 496
 497    function Leading_Part (X : T; Radix_Digits : UI) return T is
 498       L    : UI;
 499       Y, Z : T;
 500
 501    begin
 502       if Radix_Digits >= Mantissa then
 503          return X;
 504
 505       elsif Radix_Digits <= 0 then
 506          raise Constraint_Error;
 507
 508       else
 509          L := Exponent (X) - Radix_Digits;
 510          Y := Truncation (Scaling (X, -L));
 511          Z := Scaling (Y, L);
 512          return Z;
 513       end if;
 514    end Leading_Part;
 515
 516    -------------
 517    -- Machine --
 518    -------------
 519
 520    --  The trick with Machine is to force the compiler to store the result
 521    --  in memory so that we do not have extra precision used. The compiler
 522    --  is clever, so we have to outwit its possible optimizations. We do
 523    --  this by using an intermediate pragma Volatile location.
 524
 525    function Machine (X : T) return T is
 526       Temp : T;
 527       pragma Volatile (Temp);
 528    begin
 529       Temp := X;
 530       return Temp;
 531    end Machine;
 532
 533    ----------------------
 534    -- Machine_Rounding --
 535    ----------------------
 536
 537    --  For now, the implementation is identical to that of Rounding, which is
 538    --  a permissible behavior, but is not the most efficient possible approach.
 539
 540    function Machine_Rounding (X : T) return T is
 541       Result : T;
 542       Tail   : T;
 543
 544    begin
 545       Result := Truncation (abs X);
 546       Tail   := abs X - Result;
 547
 548       if Tail >= 0.5 then
 549          Result := Result + 1.0;
 550       end if;
 551
 552       if X > 0.0 then
 553          return Result;
 554
 555       elsif X < 0.0 then
 556          return -Result;
 557
 558       --  For zero case, make sure sign of zero is preserved
 559
 560       else
 561          return X;
 562       end if;
 563    end Machine_Rounding;
 564
 565    -----------
 566    -- Model --
 567    -----------
 568
 569    --  We treat Model as identical to Machine. This is true of IEEE and other
 570    --  nice floating-point systems, but not necessarily true of all systems.
 571
 572    function Model (X : T) return T is
 573    begin
 574       return T'Machine (X);
 575    end Model;
 576
 577    ----------
 578    -- Pred --
 579    ----------
 580
 581    function Pred (X : T) return T is
 582    begin
 583       --  Special treatment for largest negative number: raise Constraint_Error
 584
 585       if X = T'First then
 586          raise Constraint_Error with "Pred of largest negative number";
 587
 588       --  For finite numbers, use the symmetry around zero of floating point
 589
 590       elsif X > T'First and then X <= T'Last then
 591          pragma Annotate (CodePeer, Intentional, "test always true",
 592                           "Check for invalid float");
 593          pragma Annotate (CodePeer, Intentional, "condition predetermined",
 594                           "Check for invalid float");
 595          return -Finite_Succ (-X);
 596
 597       --  For infinities and NaNs, return unchanged
 598
 599       else
 600          return X;
 601          pragma Annotate (CodePeer, Intentional, "dead code",
 602                           "Check float range.");
 603       end if;
 604    end Pred;
 605
 606    ---------------
 607    -- Remainder --
 608    ---------------
 609
 610    function Remainder (X, Y : T) return T is
 611       A        : T;
 612       B        : T;
 613       Arg      : T;
 614       P        : T;
 615       P_Frac   : T;
 616       Sign_X   : T;
 617       IEEE_Rem : T;
 618       Arg_Exp  : UI;
 619       P_Exp    : UI;
 620       K        : UI;
 621       P_Even   : Boolean;
 622
 623       Arg_Frac : T;
 624
 625    begin
 626       if Y = 0.0 then
 627          raise Constraint_Error;
 628       end if;
 629
 630       if X > 0.0 then
 631          Sign_X :=  1.0;
 632          Arg := X;
 633       else
 634          Sign_X := -1.0;
 635          Arg := -X;
 636       end if;
 637
 638       P := abs Y;
 639
 640       if Arg < P then
 641          P_Even := True;
 642          IEEE_Rem := Arg;
 643          P_Exp := Exponent (P);
 644
 645       else
 646          Decompose (Arg, Arg_Frac, Arg_Exp);
 647          Decompose (P,   P_Frac,   P_Exp);
 648
 649          P := Compose (P_Frac, Arg_Exp);
 650          K := Arg_Exp - P_Exp;
 651          P_Even := True;
 652          IEEE_Rem := Arg;
 653
 654          for Cnt in reverse 0 .. K loop
 655             if IEEE_Rem >= P then
 656                P_Even := False;
 657                IEEE_Rem := IEEE_Rem - P;
 658             else
 659                P_Even := True;
 660             end if;
 661
 662             P := P * 0.5;
 663          end loop;
 664       end if;
 665
 666       --  That completes the calculation of modulus remainder. The final
 667       --  step is get the IEEE remainder. Here we need to compare Rem with
 668       --  (abs Y) / 2. We must be careful of unrepresentable Y/2 value
 669       --  caused by subnormal numbers
 670
 671       if P_Exp >= 0 then
 672          A := IEEE_Rem;
 673          B := abs Y * 0.5;
 674
 675       else
 676          A := IEEE_Rem * 2.0;
 677          B := abs Y;
 678       end if;
 679
 680       if A > B or else (A = B and then not P_Even) then
 681          IEEE_Rem := IEEE_Rem - abs Y;
 682       end if;
 683
 684       return Sign_X * IEEE_Rem;
 685    end Remainder;
 686
 687    --------------
 688    -- Rounding --
 689    --------------
 690
 691    function Rounding (X : T) return T is
 692       Result : T;
 693       Tail   : T;
 694
 695    begin
 696       Result := Truncation (abs X);
 697       Tail   := abs X - Result;
 698
 699       if Tail >= 0.5 then
 700          Result := Result + 1.0;
 701       end if;
 702
 703       if X > 0.0 then
 704          return Result;
 705
 706       elsif X < 0.0 then
 707          return -Result;
 708
 709       --  For zero case, make sure sign of zero is preserved
 710
 711       else
 712          return X;
 713       end if;
 714    end Rounding;
 715
 716    -------------
 717    -- Scaling --
 718    -------------
 719
 720    function Scaling (X : T; Adjustment : UI) return T is
 721       pragma Assert (Mantissa <= 64);
 722       --  This implementation handles only 80-bit IEEE Extended or smaller
 723
 724       package UST renames System.Unsigned_Types;
 725       use type UST.Long_Long_Unsigned;
 726
 727       XX : T := T'Machine (X);
 728
 729       Rep : Float_Rep;
 730       for Rep'Address use XX'Address;
 731       --  Rep is a view of the input floating-point parameter
 732
 733       Exp : constant IEEE_Erange :=
 734               Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias;
 735       --  Mask/Shift X to only get bits from the exponent. Then convert biased
 736       --  value to final value.
 737
 738       Minus : constant Boolean := (Rep (MSW) and Sign_Mask) /= 0;
 739       --  Mask/Shift X to only get bit from the sign
 740
 741       Expi, Expf : IEEE_Erange;
 742
 743    begin
 744       --  Check for zero, infinities, NaNs as well as no adjustment
 745
 746       if X = 0.0 or else Exp = IEEE_Emax + 1 or else Adjustment = 0 then
 747          return X;
 748
 749       --  Check for nonzero denormalized numbers
 750
 751       elsif Exp = IEEE_Emin - 1 then
 752          --  Check for zero result to protect the subtraction below
 753
 754          if Adjustment < -(Mantissa - 1) then
 755             XX := 0.0;
 756             return (if Minus then -XX else XX);
 757
 758          --  Normalize by multiplying by Radix ** (Mantissa - 1)
 759
 760          else
 761             return Scaling (XX * RM1, Adjustment - (Mantissa - 1));
 762          end if;
 763
 764       --  Case of normalized numbers
 765
 766       else
 767          --  Check for overflow
 768
 769          if Adjustment > IEEE_Emax - Exp then
 770             --  Optionally raise Constraint_Error as per RM A.5.3(29)
 771
 772             if T'Machine_Overflows then
 773                raise Constraint_Error with "Too large exponent";
 774
 775             else
 776                XX := 0.0;
 777                return (if Minus then -1.0 / XX else 1.0 / XX);
 778                pragma Annotate (CodePeer, Intentional, "overflow check",
 779                                 "Infinity produced");
 780                pragma Annotate (CodePeer, Intentional, "divide by zero",
 781                                 "Infinity produced");
 782             end if;
 783
 784          --  Check for underflow
 785
 786          elsif Adjustment < IEEE_Emin - Exp then
 787             --  Check for possibly gradual underflow (up to the hardware)
 788
 789             if Adjustment >= IEEE_Emin - Mantissa - Exp then
 790                Expf := IEEE_Emin;
 791                Expi := Exp + Adjustment - Expf;
 792
 793             --  Case of zero result
 794
 795             else
 796                XX := 0.0;
 797                return (if Minus then -XX else XX);
 798             end if;
 799
 800          --  Case of normalized results
 801
 802          else
 803             Expf := Exp + Adjustment;
 804             Expi := 0;
 805          end if;
 806
 807          Rep (MSW) := (Rep (MSW) and not Exp_Mask) +
 808                         Float_Word (IEEE_Ebias + Expf) * Exp_Factor;
 809
 810          if Expi < 0 then
 811             --  Given that Expi >= -Mantissa, only -64 is problematic
 812
 813             if Expi = -64 then
 814                pragma Annotate
 815                  (CodePeer, Intentional, "test always false",
 816                   "test always false in some instantiations");
 817                XX := XX / 2.0;
 818                Expi := -63;
 819             end if;
 820
 821             XX := XX / T (UST.Long_Long_Unsigned (2) ** (-Expi));
 822          end if;
 823
 824          return XX;
 825       end if;
 826    end Scaling;
 827
 828    ----------
 829    -- Succ --
 830    ----------
 831
 832    function Succ (X : T) return T is
 833    begin
 834       --  Special treatment for largest positive number: raise Constraint_Error
 835
 836       if X = T'Last then
 837          raise Constraint_Error with "Succ of largest positive number";
 838
 839       --  For finite numbers, call the specific routine
 840
 841       elsif X >= T'First and then X < T'Last then
 842          pragma Annotate (CodePeer, Intentional, "test always true",
 843                           "Check for invalid float");
 844          pragma Annotate (CodePeer, Intentional, "condition predetermined",
 845                           "Check for invalid float");
 846          return Finite_Succ (X);
 847
 848       --  For infinities and NaNs, return unchanged
 849
 850       else
 851          return X;
 852          pragma Annotate (CodePeer, Intentional, "dead code",
 853                           "Check float range.");
 854       end if;
 855    end Succ;
 856
 857    ----------------
 858    -- Truncation --
 859    ----------------
 860
 861    --  The basic approach is to compute
 862
 863    --    T'Machine (RM1 + N) - RM1
 864
 865    --  where N >= 0.0 and RM1 = Radix ** (Mantissa - 1)
 866
 867    --  This works provided that the intermediate result (RM1 + N) does not
 868    --  have extra precision (which is why we call Machine). When we compute
 869    --  RM1 + N, the exponent of N will be normalized and the mantissa shifted
 870    --  appropriately so the lower order bits, which cannot contribute to the
 871    --  integer part of N, fall off on the right. When we subtract RM1 again,
 872    --  the significant bits of N are shifted to the left, and what we have is
 873    --  an integer, because only the first e bits are different from zero
 874    --  (assuming binary radix here).
 875
 876    function Truncation (X : T) return T is
 877       Result : T;
 878
 879    begin
 880       Result := abs X;
 881
 882       if Result >= RM1 then
 883          return T'Machine (X);
 884
 885       else
 886          Result := T'Machine (RM1 + Result) - RM1;
 887
 888          if Result > abs X then
 889             Result := Result - 1.0;
 890          end if;
 891
 892          if X > 0.0 then
 893             return Result;
 894
 895          elsif X < 0.0 then
 896             return -Result;
 897
 898          --  For zero case, make sure sign of zero is preserved
 899
 900          else
 901             return X;
 902          end if;
 903       end if;
 904    end Truncation;
 905
 906    -----------------------
 907    -- Unbiased_Rounding --
 908    -----------------------
 909
 910    function Unbiased_Rounding (X : T) return T is
 911       Abs_X  : constant T := abs X;
 912       Result : T;
 913       Tail   : T;
 914
 915    begin
 916       Result := Truncation (Abs_X);
 917       Tail   := Abs_X - Result;
 918
 919       if Tail > 0.5 then
 920          Result := Result + 1.0;
 921
 922       elsif Tail = 0.5 then
 923          Result := 2.0 * Truncation ((Result / 2.0) + 0.5);
 924       end if;
 925
 926       if X > 0.0 then
 927          return Result;
 928
 929       elsif X < 0.0 then
 930          return -Result;
 931
 932       --  For zero case, make sure sign of zero is preserved
 933
 934       else
 935          return X;
 936       end if;
 937    end Unbiased_Rounding;
 938
 939    -----------
 940    -- Valid --
 941    -----------
 942
 943    function Valid (X : not null access T) return Boolean is
 944       type Access_T is access all T;
 945       function To_Address is
 946         new Ada.Unchecked_Conversion (Access_T, System.Address);
 947
 948       Rep : Float_Rep;
 949       for Rep'Address use To_Address (Access_T (X));
 950       --  Rep is a view of the input floating-point parameter. Note that we
 951       --  must avoid reading the actual bits of this parameter in float form
 952       --  since it may be a signalling NaN.
 953
 954       Exp : constant IEEE_Erange :=
 955               Integer ((Rep (MSW) and Exp_Mask) / Exp_Factor) - IEEE_Ebias;
 956       --  Mask/Shift X to only get bits from the exponent. Then convert biased
 957       --  value to final value.
 958
 959    begin
 960       if Exp = IEEE_Emax + 1 then
 961          --  This is an infinity or a NaN, i.e. always invalid
 962
 963          return False;
 964
 965       elsif Exp in IEEE_Emin .. IEEE_Emax then
 966          --  This is a normalized number, i.e. always valid
 967
 968          return True;
 969
 970       else pragma Assert (Exp = IEEE_Emin - 1);
 971          --  This is a denormalized number, valid if T'Denorm is True or 0.0
 972
 973          if T'Denorm then
 974             return True;
 975
 976          --  Note that we cannot do a direct comparison with 0.0 because the
 977          --  hardware may evaluate it to True for all denormalized numbers.
 978
 979          else
 980             --  First clear the sign bit (the exponent is already zero)
 981
 982             Rep (MSW) := Rep (MSW) and not Sign_Mask;
 983
 984             return (for all J in 0 .. Rep_Last => Rep (J) = 0);
 985          end if;
 986       end if;
 987    end Valid;
 988
 989 end System.Fat_Gen;