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 27with Einfo; use Einfo; 28with Errout; use Errout; 29with Sem_Util; use Sem_Util; 30with Ttypef; use Ttypef; 31with Targparm; use Targparm; 32 33package 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 972end Eval_Fat; 973