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