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