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