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