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