gcc/ada/eval_fat.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT COMPILER COMPONENTS                         --
   4 --                                                                          --
   5 --                             E V A L _ F A T                              --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --          Copyright (C) 1992-2007, Free Software Foundation, Inc.         --
  10 --                                                                          --
  11 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  12 -- terms of the  GNU General Public License as published  by the Free Soft- --
  13 -- ware  Foundation;  either version 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.  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 COPYING3.  If not, go to --
  19 -- http://www.gnu.org/licenses for a complete copy of the license.          --
  20 --                                                                          --
  21 -- GNAT was originally developed  by the GNAT team at  New York University. --
  22 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  23 --                                                                          --
  24 ------------------------------------------------------------------------------
  25
  26 with Einfo;    use Einfo;
  27 with Errout;   use Errout;
  28 with Sem_Util; use Sem_Util;
  29 with Ttypef;   use Ttypef;
  30 with Targparm; use Targparm;
  31
  32 package body Eval_Fat is
  33
  34    Radix : constant Int := 2;
  35    --  This code is currently only correct for the radix 2 case. We use the
  36    --  symbolic value Radix where possible to help in the unlikely case of
  37    --  anyone ever having to adjust this code for another value, and for
  38    --  documentation purposes.
  39
  40    --  Another assumption is that the range of the floating-point type is
  41    --  symmetric around zero.
  42
  43    type Radix_Power_Table is array (Int range 1 .. 4) of Int;
  44
  45    Radix_Powers : constant Radix_Power_Table :=
  46                     (Radix ** 1, Radix ** 2, Radix ** 3, Radix ** 4);
  47
  48    -----------------------
  49    -- Local Subprograms --
  50    -----------------------
  51
  52    procedure Decompose
  53      (RT       : R;
  54       X        : T;
  55       Fraction : out T;
  56       Exponent : out UI;
  57       Mode     : Rounding_Mode := Round);
  58    --  Decomposes a non-zero floating-point number into fraction and exponent
  59    --  parts. The fraction is in the interval 1.0 / Radix .. T'Pred (1.0) and
  60    --  uses Rbase = Radix. The result is rounded to a nearest machine number.
  61
  62    procedure Decompose_Int
  63      (RT       : R;
  64       X        : T;
  65       Fraction : out UI;
  66       Exponent : out UI;
  67       Mode     : Rounding_Mode);
  68    --  This is similar to Decompose, except that the Fraction value returned
  69    --  is an integer representing the value Fraction * Scale, where Scale is
  70    --  the value (Radix ** Machine_Mantissa (RT)). The value is obtained by
  71    --  using biased rounding (halfway cases round away from zero), round to
  72    --  even, a floor operation or a ceiling operation depending on the setting
  73    --  of Mode (see corresponding descriptions in Urealp).
  74
  75    function Machine_Emin (RT : R) return Int;
  76    --  Return value of the Machine_Emin attribute
  77
  78    --------------
  79    -- Adjacent --
  80    --------------
  81
  82    function Adjacent (RT : R; X, Towards : T) return T is
  83    begin
  84       if Towards = X then
  85          return X;
  86       elsif Towards > X then
  87          return Succ (RT, X);
  88       else
  89          return Pred (RT, X);
  90       end if;
  91    end Adjacent;
  92
  93    -------------
  94    -- Ceiling --
  95    -------------
  96
  97    function Ceiling (RT : R; X : T) return T is
  98       XT : constant T := Truncation (RT, X);
  99    begin
 100       if UR_Is_Negative (X) then
 101          return XT;
 102       elsif X = XT then
 103          return X;
 104       else
 105          return XT + Ureal_1;
 106       end if;
 107    end Ceiling;
 108
 109    -------------
 110    -- Compose --
 111    -------------
 112
 113    function Compose (RT : R; Fraction : T; Exponent : UI) return T is
 114       Arg_Frac : T;
 115       Arg_Exp  : UI;
 116       pragma Warnings (Off, Arg_Exp);
 117    begin
 118       Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
 119       return Scaling (RT, Arg_Frac, Exponent);
 120    end Compose;
 121
 122    ---------------
 123    -- Copy_Sign --
 124    ---------------
 125
 126    function Copy_Sign (RT : R; Value, Sign : T) return T is
 127       pragma Warnings (Off, RT);
 128       Result : T;
 129
 130    begin
 131       Result := abs Value;
 132
 133       if UR_Is_Negative (Sign) then
 134          return -Result;
 135       else
 136          return Result;
 137       end if;
 138    end Copy_Sign;
 139
 140    ---------------
 141    -- Decompose --
 142    ---------------
 143
 144    procedure Decompose
 145      (RT       : R;
 146       X        : T;
 147       Fraction : out T;
 148       Exponent : out UI;
 149       Mode     : Rounding_Mode := Round)
 150    is
 151       Int_F : UI;
 152
 153    begin
 154       Decompose_Int (RT, abs X, Int_F, Exponent, Mode);
 155
 156       Fraction := UR_From_Components
 157        (Num      => Int_F,
 158         Den      => UI_From_Int (Machine_Mantissa (RT)),
 159         Rbase    => Radix,
 160         Negative => False);
 161
 162       if UR_Is_Negative (X) then
 163          Fraction := -Fraction;
 164       end if;
 165
 166       return;
 167    end Decompose;
 168
 169    -------------------
 170    -- Decompose_Int --
 171    -------------------
 172
 173    --  This procedure should be modified with care, as there are many non-
 174    --  obvious details that may cause problems that are hard to detect. For
 175    --  zero arguments, Fraction and Exponent are set to zero. Note that sign
 176    --  of zero cannot be preserved.
 177
 178    procedure Decompose_Int
 179      (RT       : R;
 180       X        : T;
 181       Fraction : out UI;
 182       Exponent : out UI;
 183       Mode     : Rounding_Mode)
 184    is
 185       Base : Int := Rbase (X);
 186       N    : UI  := abs Numerator (X);
 187       D    : UI  := Denominator (X);
 188
 189       N_Times_Radix : UI;
 190
 191       Even : Boolean;
 192       --  True iff Fraction is even
 193
 194       Most_Significant_Digit : constant UI :=
 195                                  Radix ** (Machine_Mantissa (RT) - 1);
 196
 197       Uintp_Mark : Uintp.Save_Mark;
 198       --  The code is divided into blocks that systematically release
 199       --  intermediate values (this routine generates lots of junk!)
 200
 201    begin
 202       if N = Uint_0 then
 203          Fraction := Uint_0;
 204          Exponent := Uint_0;
 205          return;
 206       end if;
 207
 208       Calculate_D_And_Exponent_1 : begin
 209          Uintp_Mark := Mark;
 210          Exponent := Uint_0;
 211
 212          --  In cases where Base > 1, the actual denominator is Base**D. For
 213          --  cases where Base is a power of Radix, use the value 1 for the
 214          --  Denominator and adjust the exponent.
 215
 216          --  Note: Exponent has different sign from D, because D is a divisor
 217
 218          for Power in 1 .. Radix_Powers'Last loop
 219             if Base = Radix_Powers (Power) then
 220                Exponent := -D * Power;
 221                Base := 0;
 222                D := Uint_1;
 223                exit;
 224             end if;
 225          end loop;
 226
 227          Release_And_Save (Uintp_Mark, D, Exponent);
 228       end Calculate_D_And_Exponent_1;
 229
 230       if Base > 0 then
 231          Calculate_Exponent : begin
 232             Uintp_Mark := Mark;
 233
 234             --  For bases that are a multiple of the Radix, divide the base by
 235             --  Radix and adjust the Exponent. This will help because D will be
 236             --  much smaller and faster to process.
 237
 238             --  This occurs for decimal bases on machines with binary floating-
 239             --  point for example. When calculating 1E40, with Radix = 2, N
 240             --  will be 93 bits instead of 133.
 241
 242             --        N            E
 243             --      ------  * Radix
 244             --           D
 245             --       Base
 246
 247             --                  N                        E
 248             --    =  --------------------------  *  Radix
 249             --                     D        D
 250             --         (Base/Radix)  * Radix
 251
 252             --             N                  E-D
 253             --    =  ---------------  *  Radix
 254             --                    D
 255             --        (Base/Radix)
 256
 257             --  This code is commented out, because it causes numerous
 258             --  failures in the regression suite. To be studied ???
 259
 260             while False and then Base > 0 and then Base mod Radix = 0 loop
 261                Base := Base / Radix;
 262                Exponent := Exponent + D;
 263             end loop;
 264
 265             Release_And_Save (Uintp_Mark, Exponent);
 266          end Calculate_Exponent;
 267
 268          --  For remaining bases we must actually compute the exponentiation
 269
 270          --  Because the exponentiation can be negative, and D must be integer,
 271          --  the numerator is corrected instead.
 272
 273          Calculate_N_And_D : begin
 274             Uintp_Mark := Mark;
 275
 276             if D < 0 then
 277                N := N * Base ** (-D);
 278                D := Uint_1;
 279             else
 280                D := Base ** D;
 281             end if;
 282
 283             Release_And_Save (Uintp_Mark, N, D);
 284          end Calculate_N_And_D;
 285
 286          Base := 0;
 287       end if;
 288
 289       --  Now scale N and D so that N / D is a value in the interval [1.0 /
 290       --  Radix, 1.0) and adjust Exponent accordingly, so the value N / D *
 291       --  Radix ** Exponent remains unchanged.
 292
 293       --  Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0
 294
 295       --  N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.
 296       --  As this scaling is not possible for N is Uint_0, zero is handled
 297       --  explicitly at the start of this subprogram.
 298
 299       Calculate_N_And_Exponent : begin
 300          Uintp_Mark := Mark;
 301
 302          N_Times_Radix := N * Radix;
 303          while not (N_Times_Radix >= D) loop
 304             N := N_Times_Radix;
 305             Exponent := Exponent - 1;
 306             N_Times_Radix := N * Radix;
 307          end loop;
 308
 309          Release_And_Save (Uintp_Mark, N, Exponent);
 310       end Calculate_N_And_Exponent;
 311
 312       --  Step 2 - Adjust D so N / D < 1
 313
 314       --  Scale up D so N / D < 1, so N < D
 315
 316       Calculate_D_And_Exponent_2 : begin
 317          Uintp_Mark := Mark;
 318
 319          while not (N < D) loop
 320
 321             --  As N / D >= 1, N / (D * Radix) will be at least 1 / Radix, so
 322             --  the result of Step 1 stays valid
 323
 324             D := D * Radix;
 325             Exponent := Exponent + 1;
 326          end loop;
 327
 328          Release_And_Save (Uintp_Mark, D, Exponent);
 329       end Calculate_D_And_Exponent_2;
 330
 331       --  Here the value N / D is in the range [1.0 / Radix .. 1.0)
 332
 333       --  Now find the fraction by doing a very simple-minded division until
 334       --  enough digits have been computed.
 335
 336       --  This division works for all radices, but is only efficient for a
 337       --  binary radix. It is just like a manual division algorithm, but
 338       --  instead of moving the denominator one digit right, we move the
 339       --  numerator one digit left so the numerator and denominator remain
 340       --  integral.
 341
 342       Fraction := Uint_0;
 343       Even := True;
 344
 345       Calculate_Fraction_And_N : begin
 346          Uintp_Mark := Mark;
 347
 348          loop
 349             while N >= D loop
 350                N := N - D;
 351                Fraction := Fraction + 1;
 352                Even := not Even;
 353             end loop;
 354
 355             --  Stop when the result is in [1.0 / Radix, 1.0)
 356
 357             exit when Fraction >= Most_Significant_Digit;
 358
 359             N := N * Radix;
 360             Fraction := Fraction * Radix;
 361             Even := True;
 362          end loop;
 363
 364          Release_And_Save (Uintp_Mark, Fraction, N);
 365       end Calculate_Fraction_And_N;
 366
 367       Calculate_Fraction_And_Exponent : begin
 368          Uintp_Mark := Mark;
 369
 370          --  Determine correct rounding based on the remainder which is in
 371          --  N and the divisor D. The rounding is performed on the absolute
 372          --  value of X, so Ceiling and Floor need to check for the sign of
 373          --  X explicitly.
 374
 375          case Mode is
 376             when Round_Even =>
 377
 378                --  This rounding mode should not be used for static
 379                --  expressions, but only for compile-time evaluation of
 380                --  non-static expressions.
 381
 382                if (Even and then N * 2 > D)
 383                      or else
 384                   (not Even and then N * 2 >= D)
 385                then
 386                   Fraction := Fraction + 1;
 387                end if;
 388
 389             when Round   =>
 390
 391                --  Do not round to even as is done with IEEE arithmetic, but
 392                --  instead round away from zero when the result is exactly
 393                --  between two machine numbers. See RM 4.9(38).
 394
 395                if N * 2 >= D then
 396                   Fraction := Fraction + 1;
 397                end if;
 398
 399             when Ceiling =>
 400                if N > Uint_0 and then not UR_Is_Negative (X) then
 401                   Fraction := Fraction + 1;
 402                end if;
 403
 404             when Floor   =>
 405                if N > Uint_0 and then UR_Is_Negative (X) then
 406                   Fraction := Fraction + 1;
 407                end if;
 408          end case;
 409
 410          --  The result must be normalized to [1.0/Radix, 1.0), so adjust if
 411          --  the result is 1.0 because of rounding.
 412
 413          if Fraction = Most_Significant_Digit * Radix then
 414             Fraction := Most_Significant_Digit;
 415             Exponent := Exponent + 1;
 416          end if;
 417
 418          --  Put back sign after applying the rounding
 419
 420          if UR_Is_Negative (X) then
 421             Fraction := -Fraction;
 422          end if;
 423
 424          Release_And_Save (Uintp_Mark, Fraction, Exponent);
 425       end Calculate_Fraction_And_Exponent;
 426    end Decompose_Int;
 427
 428    --------------
 429    -- Exponent --
 430    --------------
 431
 432    function Exponent (RT : R; X : T) return UI is
 433       X_Frac : UI;
 434       X_Exp  : UI;
 435       pragma Warnings (Off, X_Frac);
 436    begin
 437       Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);
 438       return X_Exp;
 439    end Exponent;
 440
 441    -----------
 442    -- Floor --
 443    -----------
 444
 445    function Floor (RT : R; X : T) return T is
 446       XT : constant T := Truncation (RT, X);
 447
 448    begin
 449       if UR_Is_Positive (X) then
 450          return XT;
 451
 452       elsif XT = X then
 453          return X;
 454
 455       else
 456          return XT - Ureal_1;
 457       end if;
 458    end Floor;
 459
 460    --------------
 461    -- Fraction --
 462    --------------
 463
 464    function Fraction (RT : R; X : T) return T is
 465       X_Frac : T;
 466       X_Exp  : UI;
 467       pragma Warnings (Off, X_Exp);
 468    begin
 469       Decompose (RT, X, X_Frac, X_Exp);
 470       return X_Frac;
 471    end Fraction;
 472
 473    ------------------
 474    -- Leading_Part --
 475    ------------------
 476
 477    function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
 478       RD : constant UI := UI_Min (Radix_Digits, Machine_Mantissa (RT));
 479       L  : UI;
 480       Y  : T;
 481    begin
 482       L := Exponent (RT, X) - RD;
 483       Y := UR_From_Uint (UR_Trunc (Scaling (RT, X, -L)));
 484       return Scaling (RT, Y, L);
 485    end Leading_Part;
 486
 487    -------------
 488    -- Machine --
 489    -------------
 490
 491    function Machine
 492      (RT    : R;
 493       X     : T;
 494       Mode  : Rounding_Mode;
 495       Enode : Node_Id) return T
 496    is
 497       X_Frac : T;
 498       X_Exp  : UI;
 499       Emin   : constant UI := UI_From_Int (Machine_Emin (RT));
 500
 501    begin
 502       Decompose (RT, X, X_Frac, X_Exp, Mode);
 503
 504       --  Case of denormalized number or (gradual) underflow
 505
 506       --  A denormalized number is one with the minimum exponent Emin, but that
 507       --  breaks the assumption that the first digit of the mantissa is a one.
 508       --  This allows the first non-zero digit to be in any of the remaining
 509       --  Mant - 1 spots. The gap between subsequent denormalized numbers is
 510       --  the same as for the smallest normalized numbers. However, the number
 511       --  of significant digits left decreases as a result of the mantissa now
 512       --  having leading seros.
 513
 514       if X_Exp < Emin then
 515          declare
 516             Emin_Den : constant UI :=
 517                          UI_From_Int
 518                            (Machine_Emin (RT) - Machine_Mantissa (RT) + 1);
 519          begin
 520             if X_Exp < Emin_Den or not Denorm_On_Target then
 521                if UR_Is_Negative (X) then
 522                   Error_Msg_N
 523                     ("floating-point value underflows to -0.0?", Enode);
 524                   return Ureal_M_0;
 525
 526                else
 527                   Error_Msg_N
 528                     ("floating-point value underflows to 0.0?", Enode);
 529                   return Ureal_0;
 530                end if;
 531
 532             elsif Denorm_On_Target then
 533
 534                --  Emin - Mant <= X_Exp < Emin, so result is denormal. Handle
 535                --  gradual underflow by first computing the number of
 536                --  significant bits still available for the mantissa and
 537                --  then truncating the fraction to this number of bits.
 538
 539                --  If this value is different from the original fraction,
 540                --  precision is lost due to gradual underflow.
 541
 542                --  We probably should round here and prevent double rounding as
 543                --  a result of first rounding to a model number and then to a
 544                --  machine number. However, this is an extremely rare case that
 545                --  is not worth the extra complexity. In any case, a warning is
 546                --  issued in cases where gradual underflow occurs.
 547
 548                declare
 549                   Denorm_Sig_Bits : constant UI := X_Exp - Emin_Den + 1;
 550
 551                   X_Frac_Denorm   : constant T := UR_From_Components
 552                     (UR_Trunc (Scaling (RT, abs X_Frac, Denorm_Sig_Bits)),
 553                      Denorm_Sig_Bits,
 554                      Radix,
 555                      UR_Is_Negative (X));
 556
 557                begin
 558                   if X_Frac_Denorm /= X_Frac then
 559                      Error_Msg_N
 560                        ("gradual underflow causes loss of precision?",
 561                         Enode);
 562                      X_Frac := X_Frac_Denorm;
 563                   end if;
 564                end;
 565             end if;
 566          end;
 567       end if;
 568
 569       return Scaling (RT, X_Frac, X_Exp);
 570    end Machine;
 571
 572    ------------------
 573    -- Machine_Emin --
 574    ------------------
 575
 576    function Machine_Emin (RT : R) return Int is
 577       Digs : constant UI := Digits_Value (RT);
 578       Emin : Int;
 579
 580    begin
 581       if Vax_Float (RT) then
 582          if Digs = VAXFF_Digits then
 583             Emin := VAXFF_Machine_Emin;
 584
 585          elsif Digs = VAXDF_Digits then
 586             Emin := VAXDF_Machine_Emin;
 587
 588          else
 589             pragma Assert (Digs = VAXGF_Digits);
 590             Emin := VAXGF_Machine_Emin;
 591          end if;
 592
 593       elsif Is_AAMP_Float (RT) then
 594          if Digs = AAMPS_Digits then
 595             Emin := AAMPS_Machine_Emin;
 596
 597          else
 598             pragma Assert (Digs = AAMPL_Digits);
 599             Emin := AAMPL_Machine_Emin;
 600          end if;
 601
 602       else
 603          if Digs = IEEES_Digits then
 604             Emin := IEEES_Machine_Emin;
 605
 606          elsif Digs = IEEEL_Digits then
 607             Emin := IEEEL_Machine_Emin;
 608
 609          else
 610             pragma Assert (Digs = IEEEX_Digits);
 611             Emin := IEEEX_Machine_Emin;
 612          end if;
 613       end if;
 614
 615       return Emin;
 616    end Machine_Emin;
 617
 618    ----------------------
 619    -- Machine_Mantissa --
 620    ----------------------
 621
 622    function Machine_Mantissa (RT : R) return Nat is
 623       Digs : constant UI := Digits_Value (RT);
 624       Mant : Nat;
 625
 626    begin
 627       if Vax_Float (RT) then
 628          if Digs = VAXFF_Digits then
 629             Mant := VAXFF_Machine_Mantissa;
 630
 631          elsif Digs = VAXDF_Digits then
 632             Mant := VAXDF_Machine_Mantissa;
 633
 634          else
 635             pragma Assert (Digs = VAXGF_Digits);
 636             Mant := VAXGF_Machine_Mantissa;
 637          end if;
 638
 639       elsif Is_AAMP_Float (RT) then
 640          if Digs = AAMPS_Digits then
 641             Mant := AAMPS_Machine_Mantissa;
 642
 643          else
 644             pragma Assert (Digs = AAMPL_Digits);
 645             Mant := AAMPL_Machine_Mantissa;
 646          end if;
 647
 648       else
 649          if Digs = IEEES_Digits then
 650             Mant := IEEES_Machine_Mantissa;
 651
 652          elsif Digs = IEEEL_Digits then
 653             Mant := IEEEL_Machine_Mantissa;
 654
 655          else
 656             pragma Assert (Digs = IEEEX_Digits);
 657             Mant := IEEEX_Machine_Mantissa;
 658          end if;
 659       end if;
 660
 661       return Mant;
 662    end Machine_Mantissa;
 663
 664    -------------------
 665    -- Machine_Radix --
 666    -------------------
 667
 668    function Machine_Radix (RT : R) return Nat is
 669       pragma Warnings (Off, RT);
 670    begin
 671       return Radix;
 672    end Machine_Radix;
 673
 674    -----------
 675    -- Model --
 676    -----------
 677
 678    function Model (RT : R; X : T) return T is
 679       X_Frac : T;
 680       X_Exp  : UI;
 681    begin
 682       Decompose (RT, X, X_Frac, X_Exp);
 683       return Compose (RT, X_Frac, X_Exp);
 684    end Model;
 685
 686    ----------
 687    -- Pred --
 688    ----------
 689
 690    function Pred (RT : R; X : T) return T is
 691    begin
 692       return -Succ (RT, -X);
 693    end Pred;
 694
 695    ---------------
 696    -- Remainder --
 697    ---------------
 698
 699    function Remainder (RT : R; X, Y : T) return T is
 700       A        : T;
 701       B        : T;
 702       Arg      : T;
 703       P        : T;
 704       Arg_Frac : T;
 705       P_Frac   : T;
 706       Sign_X   : T;
 707       IEEE_Rem : T;
 708       Arg_Exp  : UI;
 709       P_Exp    : UI;
 710       K        : UI;
 711       P_Even   : Boolean;
 712
 713       pragma Warnings (Off, Arg_Frac);
 714
 715    begin
 716       if UR_Is_Positive (X) then
 717          Sign_X :=  Ureal_1;
 718       else
 719          Sign_X := -Ureal_1;
 720       end if;
 721
 722       Arg := abs X;
 723       P   := abs Y;
 724
 725       if Arg < P then
 726          P_Even := True;
 727          IEEE_Rem := Arg;
 728          P_Exp := Exponent (RT, P);
 729
 730       else
 731          --  ??? what about zero cases?
 732          Decompose (RT, Arg, Arg_Frac, Arg_Exp);
 733          Decompose (RT, P,   P_Frac,   P_Exp);
 734
 735          P := Compose (RT, P_Frac, Arg_Exp);
 736          K := Arg_Exp - P_Exp;
 737          P_Even := True;
 738          IEEE_Rem := Arg;
 739
 740          for Cnt in reverse 0 .. UI_To_Int (K) loop
 741             if IEEE_Rem >= P then
 742                P_Even := False;
 743                IEEE_Rem := IEEE_Rem - P;
 744             else
 745                P_Even := True;
 746             end if;
 747
 748             P := P * Ureal_Half;
 749          end loop;
 750       end if;
 751
 752       --  That completes the calculation of modulus remainder. The final step
 753       --  is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
 754
 755       if P_Exp >= 0 then
 756          A := IEEE_Rem;
 757          B := abs Y * Ureal_Half;
 758
 759       else
 760          A := IEEE_Rem * Ureal_2;
 761          B := abs Y;
 762       end if;
 763
 764       if A > B or else (A = B and then not P_Even) then
 765          IEEE_Rem := IEEE_Rem - abs Y;
 766       end if;
 767
 768       return Sign_X * IEEE_Rem;
 769    end Remainder;
 770
 771    --------------
 772    -- Rounding --
 773    --------------
 774
 775    function Rounding (RT : R; X : T) return T is
 776       Result : T;
 777       Tail   : T;
 778
 779    begin
 780       Result := Truncation (RT, abs X);
 781       Tail   := abs X - Result;
 782
 783       if Tail >= Ureal_Half  then
 784          Result := Result + Ureal_1;
 785       end if;
 786
 787       if UR_Is_Negative (X) then
 788          return -Result;
 789       else
 790          return Result;
 791       end if;
 792    end Rounding;
 793
 794    -------------
 795    -- Scaling --
 796    -------------
 797
 798    function Scaling (RT : R; X : T; Adjustment : UI) return T is
 799       pragma Warnings (Off, RT);
 800
 801    begin
 802       if Rbase (X) = Radix then
 803          return UR_From_Components
 804            (Num      => Numerator (X),
 805             Den      => Denominator (X) - Adjustment,
 806             Rbase    => Radix,
 807             Negative => UR_Is_Negative (X));
 808
 809       elsif Adjustment >= 0 then
 810          return X * Radix ** Adjustment;
 811       else
 812          return X / Radix ** (-Adjustment);
 813       end if;
 814    end Scaling;
 815
 816    ----------
 817    -- Succ --
 818    ----------
 819
 820    function Succ (RT : R; X : T) return T is
 821       Emin     : constant UI := UI_From_Int (Machine_Emin (RT));
 822       Mantissa : constant UI := UI_From_Int (Machine_Mantissa (RT));
 823       Exp      : UI := UI_Max (Emin, Exponent (RT, X));
 824       Frac     : T;
 825       New_Frac : T;
 826
 827    begin
 828       if UR_Is_Zero (X) then
 829          Exp := Emin;
 830       end if;
 831
 832       --  Set exponent such that the radix point will be directly following the
 833       --  mantissa after scaling.
 834
 835       if Denorm_On_Target or Exp /= Emin then
 836          Exp := Exp - Mantissa;
 837       else
 838          Exp := Exp - 1;
 839       end if;
 840
 841       Frac := Scaling (RT, X, -Exp);
 842       New_Frac := Ceiling (RT, Frac);
 843
 844       if New_Frac = Frac then
 845          if New_Frac = Scaling (RT, -Ureal_1, Mantissa - 1) then
 846             New_Frac := New_Frac + Scaling (RT, Ureal_1, Uint_Minus_1);
 847          else
 848             New_Frac := New_Frac + Ureal_1;
 849          end if;
 850       end if;
 851
 852       return Scaling (RT, New_Frac, Exp);
 853    end Succ;
 854
 855    ----------------
 856    -- Truncation --
 857    ----------------
 858
 859    function Truncation (RT : R; X : T) return T is
 860       pragma Warnings (Off, RT);
 861    begin
 862       return UR_From_Uint (UR_Trunc (X));
 863    end Truncation;
 864
 865    -----------------------
 866    -- Unbiased_Rounding --
 867    -----------------------
 868
 869    function Unbiased_Rounding (RT : R; X : T) return T is
 870       Abs_X  : constant T := abs X;
 871       Result : T;
 872       Tail   : T;
 873
 874    begin
 875       Result := Truncation (RT, Abs_X);
 876       Tail   := Abs_X - Result;
 877
 878       if Tail > Ureal_Half  then
 879          Result := Result + Ureal_1;
 880
 881       elsif Tail = Ureal_Half then
 882          Result := Ureal_2 *
 883                      Truncation (RT, (Result / Ureal_2) + Ureal_Half);
 884       end if;
 885
 886       if UR_Is_Negative (X) then
 887          return -Result;
 888       elsif UR_Is_Positive (X) then
 889          return Result;
 890
 891       --  For zero case, make sure sign of zero is preserved
 892
 893       else
 894          return X;
 895       end if;
 896    end Unbiased_Rounding;
 897
 898 end Eval_Fat;