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