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