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