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