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