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