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