1------------------------------------------------------------------------------ 2-- -- 3-- GNAT RUN-TIME COMPONENTS -- 4-- -- 5-- S Y S T E M . A R I T H _ D O U B L E -- 6-- -- 7-- B o d y -- 8-- -- 9-- Copyright (C) 1992-2020, 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. -- 17-- -- 18-- As a special exception under Section 7 of GPL version 3, you are granted -- 19-- additional permissions described in the GCC Runtime Library Exception, -- 20-- version 3.1, as published by the Free Software Foundation. -- 21-- -- 22-- You should have received a copy of the GNU General Public License and -- 23-- a copy of the GCC Runtime Library Exception along with this program; -- 24-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- 25-- <http://www.gnu.org/licenses/>. -- 26-- -- 27-- GNAT was originally developed by the GNAT team at New York University. -- 28-- Extensive contributions were provided by Ada Core Technologies Inc. -- 29-- -- 30------------------------------------------------------------------------------ 31 32with Ada.Unchecked_Conversion; 33 34package body System.Arith_Double is 35 36 pragma Suppress (Overflow_Check); 37 pragma Suppress (Range_Check); 38 39 function To_Uns is new Ada.Unchecked_Conversion (Double_Int, Double_Uns); 40 function To_Int is new Ada.Unchecked_Conversion (Double_Uns, Double_Int); 41 42 Double_Size : constant Natural := Double_Int'Size; 43 Single_Size : constant Natural := Double_Int'Size / 2; 44 45 ----------------------- 46 -- Local Subprograms -- 47 ----------------------- 48 49 function "+" (A, B : Single_Uns) return Double_Uns is 50 (Double_Uns (A) + Double_Uns (B)); 51 function "+" (A : Double_Uns; B : Single_Uns) return Double_Uns is 52 (A + Double_Uns (B)); 53 -- Length doubling additions 54 55 function "*" (A, B : Single_Uns) return Double_Uns is 56 (Double_Uns (A) * Double_Uns (B)); 57 -- Length doubling multiplication 58 59 function "/" (A : Double_Uns; B : Single_Uns) return Double_Uns is 60 (A / Double_Uns (B)); 61 -- Length doubling division 62 63 function "&" (Hi, Lo : Single_Uns) return Double_Uns is 64 (Shift_Left (Double_Uns (Hi), Single_Size) or Double_Uns (Lo)); 65 -- Concatenate hi, lo values to form double result 66 67 function "abs" (X : Double_Int) return Double_Uns is 68 (if X = Double_Int'First 69 then 2 ** (Double_Size - 1) 70 else Double_Uns (Double_Int'(abs X))); 71 -- Convert absolute value of X to unsigned. Note that we can't just use 72 -- the expression of the Else since it overflows for X = Double_Int'First. 73 74 function "rem" (A : Double_Uns; B : Single_Uns) return Double_Uns is 75 (A rem Double_Uns (B)); 76 -- Length doubling remainder 77 78 function Le3 (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) return Boolean; 79 -- Determines if (3 * Single_Size)-bit value X1&X2&X3 <= Y1&Y2&Y3 80 81 function Lo (A : Double_Uns) return Single_Uns is 82 (Single_Uns (A and (2 ** Single_Size - 1))); 83 -- Low order half of double value 84 85 function Hi (A : Double_Uns) return Single_Uns is 86 (Single_Uns (Shift_Right (A, Single_Size))); 87 -- High order half of double value 88 89 procedure Sub3 (X1, X2, X3 : in out Single_Uns; Y1, Y2, Y3 : Single_Uns); 90 -- Computes X1&X2&X3 := X1&X2&X3 - Y1&Y1&Y3 mod 2 ** (3 * Single_Size) 91 92 function To_Neg_Int (A : Double_Uns) return Double_Int; 93 -- Convert to negative integer equivalent. If the input is in the range 94 -- 0 .. 2 ** (Double_Size - 1), then the corresponding nonpositive signed 95 -- integer (obtained by negating the given value) is returned, otherwise 96 -- constraint error is raised. 97 98 function To_Pos_Int (A : Double_Uns) return Double_Int; 99 -- Convert to positive integer equivalent. If the input is in the range 100 -- 0 .. 2 ** (Double_Size - 1) - 1, then the corresponding non-negative 101 -- signed integer is returned, otherwise constraint error is raised. 102 103 procedure Raise_Error; 104 pragma No_Return (Raise_Error); 105 -- Raise constraint error with appropriate message 106 107 -------------------------- 108 -- Add_With_Ovflo_Check -- 109 -------------------------- 110 111 function Add_With_Ovflo_Check (X, Y : Double_Int) return Double_Int is 112 R : constant Double_Int := To_Int (To_Uns (X) + To_Uns (Y)); 113 114 begin 115 if X >= 0 then 116 if Y < 0 or else R >= 0 then 117 return R; 118 end if; 119 120 else -- X < 0 121 if Y > 0 or else R < 0 then 122 return R; 123 end if; 124 end if; 125 126 Raise_Error; 127 end Add_With_Ovflo_Check; 128 129 ------------------- 130 -- Double_Divide -- 131 ------------------- 132 133 procedure Double_Divide 134 (X, Y, Z : Double_Int; 135 Q, R : out Double_Int; 136 Round : Boolean) 137 is 138 Xu : constant Double_Uns := abs X; 139 Yu : constant Double_Uns := abs Y; 140 141 Yhi : constant Single_Uns := Hi (Yu); 142 Ylo : constant Single_Uns := Lo (Yu); 143 144 Zu : constant Double_Uns := abs Z; 145 Zhi : constant Single_Uns := Hi (Zu); 146 Zlo : constant Single_Uns := Lo (Zu); 147 148 T1, T2 : Double_Uns; 149 Du, Qu, Ru : Double_Uns; 150 Den_Pos : Boolean; 151 152 begin 153 if Yu = 0 or else Zu = 0 then 154 Raise_Error; 155 end if; 156 157 -- Set final signs (RM 4.5.5(27-30)) 158 159 Den_Pos := (Y < 0) = (Z < 0); 160 161 -- Compute Y * Z. Note that if the result overflows Double_Uns, then 162 -- the rounded result is zero, except for the very special case where 163 -- X = -2 ** (Double_Size - 1) and abs(Y*Z) = 2 ** Double_Size, when 164 -- Round is True. 165 166 if Yhi /= 0 then 167 if Zhi /= 0 then 168 169 -- Handle the special case when Round is True 170 171 if Yhi = 1 172 and then Zhi = 1 173 and then Ylo = 0 174 and then Zlo = 0 175 and then X = Double_Int'First 176 and then Round 177 then 178 Q := (if Den_Pos then -1 else 1); 179 else 180 Q := 0; 181 end if; 182 183 R := X; 184 return; 185 else 186 T2 := Yhi * Zlo; 187 end if; 188 189 else 190 T2 := Ylo * Zhi; 191 end if; 192 193 T1 := Ylo * Zlo; 194 T2 := T2 + Hi (T1); 195 196 if Hi (T2) /= 0 then 197 198 -- Handle the special case when Round is True 199 200 if Hi (T2) = 1 201 and then Lo (T2) = 0 202 and then Lo (T1) = 0 203 and then X = Double_Int'First 204 and then Round 205 then 206 Q := (if Den_Pos then -1 else 1); 207 else 208 Q := 0; 209 end if; 210 211 R := X; 212 return; 213 end if; 214 215 Du := Lo (T2) & Lo (T1); 216 217 -- Check overflow case of largest negative number divided by -1 218 219 if X = Double_Int'First and then Du = 1 and then not Den_Pos then 220 Raise_Error; 221 end if; 222 223 -- Perform the actual division 224 225 pragma Assert (Du /= 0); 226 -- Multiplication of 2-limb arguments Yu and Zu leads to 4-limb result 227 -- (where each limb is a single value). Cases where 4 limbs are needed 228 -- require Yhi/=0 and Zhi/=0 and lead to early exit. Remaining cases 229 -- where 3 limbs are needed correspond to Hi(T2)/=0 and lead to early 230 -- exit. Thus, at this point, the result fits in 2 limbs which are 231 -- exactly Lo(T2) and Lo(T1), which corresponds to the value of Du. 232 -- As the case where one of Yu or Zu is null also led to early exit, 233 -- we have Du/=0 here. 234 Qu := Xu / Du; 235 Ru := Xu rem Du; 236 237 -- Deal with rounding case 238 239 if Round and then Ru > (Du - Double_Uns'(1)) / Double_Uns'(2) then 240 Qu := Qu + Double_Uns'(1); 241 end if; 242 243 -- Case of dividend (X) sign positive 244 245 if X >= 0 then 246 R := To_Int (Ru); 247 Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu)); 248 249 -- Case of dividend (X) sign negative 250 251 -- We perform the unary minus operation on the unsigned value 252 -- before conversion to signed, to avoid a possible overflow 253 -- for value -2 ** (Double_Size - 1), both for computing R and Q. 254 255 else 256 R := To_Int (-Ru); 257 Q := (if Den_Pos then To_Int (-Qu) else To_Int (Qu)); 258 end if; 259 end Double_Divide; 260 261 --------- 262 -- Le3 -- 263 --------- 264 265 function Le3 (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) return Boolean is 266 begin 267 if X1 < Y1 then 268 return True; 269 elsif X1 > Y1 then 270 return False; 271 elsif X2 < Y2 then 272 return True; 273 elsif X2 > Y2 then 274 return False; 275 else 276 return X3 <= Y3; 277 end if; 278 end Le3; 279 280 ------------------------------- 281 -- Multiply_With_Ovflo_Check -- 282 ------------------------------- 283 284 function Multiply_With_Ovflo_Check (X, Y : Double_Int) return Double_Int is 285 Xu : constant Double_Uns := abs X; 286 Xhi : constant Single_Uns := Hi (Xu); 287 Xlo : constant Single_Uns := Lo (Xu); 288 289 Yu : constant Double_Uns := abs Y; 290 Yhi : constant Single_Uns := Hi (Yu); 291 Ylo : constant Single_Uns := Lo (Yu); 292 293 T1, T2 : Double_Uns; 294 295 begin 296 if Xhi /= 0 then 297 if Yhi /= 0 then 298 Raise_Error; 299 else 300 T2 := Xhi * Ylo; 301 end if; 302 303 elsif Yhi /= 0 then 304 T2 := Xlo * Yhi; 305 306 else -- Yhi = Xhi = 0 307 T2 := 0; 308 end if; 309 310 -- Here we have T2 set to the contribution to the upper half of the 311 -- result from the upper halves of the input values. 312 313 T1 := Xlo * Ylo; 314 T2 := T2 + Hi (T1); 315 316 if Hi (T2) /= 0 then 317 Raise_Error; 318 end if; 319 320 T2 := Lo (T2) & Lo (T1); 321 322 if X >= 0 then 323 if Y >= 0 then 324 return To_Pos_Int (T2); 325 pragma Annotate (CodePeer, Intentional, "precondition", 326 "Intentional Unsigned->Signed conversion"); 327 else 328 return To_Neg_Int (T2); 329 end if; 330 else -- X < 0 331 if Y < 0 then 332 return To_Pos_Int (T2); 333 pragma Annotate (CodePeer, Intentional, "precondition", 334 "Intentional Unsigned->Signed conversion"); 335 else 336 return To_Neg_Int (T2); 337 end if; 338 end if; 339 340 end Multiply_With_Ovflo_Check; 341 342 ----------------- 343 -- Raise_Error -- 344 ----------------- 345 346 procedure Raise_Error is 347 begin 348 raise Constraint_Error with "Double arithmetic overflow"; 349 end Raise_Error; 350 351 ------------------- 352 -- Scaled_Divide -- 353 ------------------- 354 355 procedure Scaled_Divide 356 (X, Y, Z : Double_Int; 357 Q, R : out Double_Int; 358 Round : Boolean) 359 is 360 Xu : constant Double_Uns := abs X; 361 Xhi : constant Single_Uns := Hi (Xu); 362 Xlo : constant Single_Uns := Lo (Xu); 363 364 Yu : constant Double_Uns := abs Y; 365 Yhi : constant Single_Uns := Hi (Yu); 366 Ylo : constant Single_Uns := Lo (Yu); 367 368 Zu : Double_Uns := abs Z; 369 Zhi : Single_Uns := Hi (Zu); 370 Zlo : Single_Uns := Lo (Zu); 371 372 D : array (1 .. 4) of Single_Uns; 373 -- The dividend, four digits (D(1) is high order) 374 375 Qd : array (1 .. 2) of Single_Uns; 376 -- The quotient digits, two digits (Qd(1) is high order) 377 378 S1, S2, S3 : Single_Uns; 379 -- Value to subtract, three digits (S1 is high order) 380 381 Qu : Double_Uns; 382 Ru : Double_Uns; 383 -- Unsigned quotient and remainder 384 385 Mask : Single_Uns; 386 -- Mask of bits used to compute the scaling factor below 387 388 Scale : Natural; 389 -- Scaling factor used for multiple-precision divide. Dividend and 390 -- Divisor are multiplied by 2 ** Scale, and the final remainder is 391 -- divided by the scaling factor. The reason for this scaling is to 392 -- allow more accurate estimation of quotient digits. 393 394 Shift : Natural; 395 -- Shift factor used to compute the scaling factor above 396 397 T1, T2, T3 : Double_Uns; 398 -- Temporary values 399 400 begin 401 -- First do the multiplication, giving the four digit dividend 402 403 T1 := Xlo * Ylo; 404 D (4) := Lo (T1); 405 D (3) := Hi (T1); 406 407 if Yhi /= 0 then 408 T1 := Xlo * Yhi; 409 T2 := D (3) + Lo (T1); 410 D (3) := Lo (T2); 411 D (2) := Hi (T1) + Hi (T2); 412 413 if Xhi /= 0 then 414 T1 := Xhi * Ylo; 415 T2 := D (3) + Lo (T1); 416 D (3) := Lo (T2); 417 T3 := D (2) + Hi (T1); 418 T3 := T3 + Hi (T2); 419 D (2) := Lo (T3); 420 D (1) := Hi (T3); 421 422 T1 := (D (1) & D (2)) + Double_Uns'(Xhi * Yhi); 423 D (1) := Hi (T1); 424 D (2) := Lo (T1); 425 426 else 427 D (1) := 0; 428 end if; 429 430 else 431 if Xhi /= 0 then 432 T1 := Xhi * Ylo; 433 T2 := D (3) + Lo (T1); 434 D (3) := Lo (T2); 435 D (2) := Hi (T1) + Hi (T2); 436 437 else 438 D (2) := 0; 439 end if; 440 441 D (1) := 0; 442 end if; 443 444 -- Now it is time for the dreaded multiple precision division. First an 445 -- easy case, check for the simple case of a one digit divisor. 446 447 if Zhi = 0 then 448 if D (1) /= 0 or else D (2) >= Zlo then 449 Raise_Error; 450 451 -- Here we are dividing at most three digits by one digit 452 453 else 454 T1 := D (2) & D (3); 455 T2 := Lo (T1 rem Zlo) & D (4); 456 457 Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo); 458 Ru := T2 rem Zlo; 459 end if; 460 461 -- If divisor is double digit and dividend is too large, raise error 462 463 elsif (D (1) & D (2)) >= Zu then 464 Raise_Error; 465 466 -- This is the complex case where we definitely have a double digit 467 -- divisor and a dividend of at least three digits. We use the classical 468 -- multiple-precision division algorithm (see section (4.3.1) of Knuth's 469 -- "The Art of Computer Programming", Vol. 2 for a description 470 -- (algorithm D). 471 472 else 473 -- First normalize the divisor so that it has the leading bit on. 474 -- We do this by finding the appropriate left shift amount. 475 476 Shift := Single_Size / 2; 477 Mask := Shift_Left (2 ** (Single_Size / 2) - 1, Shift); 478 Scale := 0; 479 480 while Shift /= 0 loop 481 if (Hi (Zu) and Mask) = 0 then 482 Scale := Scale + Shift; 483 Zu := Shift_Left (Zu, Shift); 484 end if; 485 486 Shift := Shift / 2; 487 Mask := Shift_Left (Mask, Shift); 488 end loop; 489 490 Zhi := Hi (Zu); 491 Zlo := Lo (Zu); 492 493 pragma Assert (Zhi /= 0); 494 -- We have Hi(Zu)/=0 before normalization. The sequence of Shift_Left 495 -- operations results in the leading bit of Zu being 1 by moving the 496 -- leftmost 1-bit in Zu to leading position, thus Zhi=Hi(Zu)/=0 here. 497 498 -- Note that when we scale up the dividend, it still fits in four 499 -- digits, since we already tested for overflow, and scaling does 500 -- not change the invariant that (D (1) & D (2)) < Zu. 501 502 T1 := Shift_Left (D (1) & D (2), Scale); 503 D (1) := Hi (T1); 504 T2 := Shift_Left (0 & D (3), Scale); 505 D (2) := Lo (T1) or Hi (T2); 506 T3 := Shift_Left (0 & D (4), Scale); 507 D (3) := Lo (T2) or Hi (T3); 508 D (4) := Lo (T3); 509 510 -- Loop to compute quotient digits, runs twice for Qd(1) and Qd(2) 511 512 for J in 0 .. 1 loop 513 514 -- Compute next quotient digit. We have to divide three digits by 515 -- two digits. We estimate the quotient by dividing the leading 516 -- two digits by the leading digit. Given the scaling we did above 517 -- which ensured the first bit of the divisor is set, this gives 518 -- an estimate of the quotient that is at most two too high. 519 520 Qd (J + 1) := (if D (J + 1) = Zhi 521 then 2 ** Single_Size - 1 522 else Lo ((D (J + 1) & D (J + 2)) / Zhi)); 523 524 -- Compute amount to subtract 525 526 T1 := Qd (J + 1) * Zlo; 527 T2 := Qd (J + 1) * Zhi; 528 S3 := Lo (T1); 529 T1 := Hi (T1) + Lo (T2); 530 S2 := Lo (T1); 531 S1 := Hi (T1) + Hi (T2); 532 533 -- Adjust quotient digit if it was too high 534 535 -- We use the version of the algorithm in the 2nd Edition of 536 -- "The Art of Computer Programming". This had a bug not 537 -- discovered till 1995, see Vol 2 errata: 538 -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. 539 -- Under rare circumstances the expression in the test could 540 -- overflow. This version was further corrected in 2005, see 541 -- Vol 2 errata: 542 -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz. 543 -- This implementation is not impacted by these bugs, due to the 544 -- use of a word-size comparison done in function Le3 instead of 545 -- a comparison on two-word integer quantities in the original 546 -- algorithm. 547 548 loop 549 exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3)); 550 Qd (J + 1) := Qd (J + 1) - 1; 551 Sub3 (S1, S2, S3, 0, Zhi, Zlo); 552 end loop; 553 554 -- Now subtract S1&S2&S3 from D1&D2&D3 ready for next step 555 556 Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3); 557 end loop; 558 559 -- The two quotient digits are now set, and the remainder of the 560 -- scaled division is in D3&D4. To get the remainder for the 561 -- original unscaled division, we rescale this dividend. 562 563 -- We rescale the divisor as well, to make the proper comparison 564 -- for rounding below. 565 566 Qu := Qd (1) & Qd (2); 567 Ru := Shift_Right (D (3) & D (4), Scale); 568 Zu := Shift_Right (Zu, Scale); 569 end if; 570 571 -- Deal with rounding case 572 573 if Round and then Ru > (Zu - Double_Uns'(1)) / Double_Uns'(2) then 574 575 -- Protect against wrapping around when rounding, by signaling 576 -- an overflow when the quotient is too large. 577 578 if Qu = Double_Uns'Last then 579 Raise_Error; 580 end if; 581 582 Qu := Qu + Double_Uns'(1); 583 end if; 584 585 -- Set final signs (RM 4.5.5(27-30)) 586 587 -- Case of dividend (X * Y) sign positive 588 589 if (X >= 0 and then Y >= 0) or else (X < 0 and then Y < 0) then 590 R := To_Pos_Int (Ru); 591 Q := (if Z > 0 then To_Pos_Int (Qu) else To_Neg_Int (Qu)); 592 593 -- Case of dividend (X * Y) sign negative 594 595 else 596 R := To_Neg_Int (Ru); 597 Q := (if Z > 0 then To_Neg_Int (Qu) else To_Pos_Int (Qu)); 598 end if; 599 end Scaled_Divide; 600 601 ---------- 602 -- Sub3 -- 603 ---------- 604 605 procedure Sub3 (X1, X2, X3 : in out Single_Uns; Y1, Y2, Y3 : Single_Uns) is 606 begin 607 if Y3 > X3 then 608 if X2 = 0 then 609 X1 := X1 - 1; 610 end if; 611 612 X2 := X2 - 1; 613 end if; 614 615 X3 := X3 - Y3; 616 617 if Y2 > X2 then 618 X1 := X1 - 1; 619 end if; 620 621 X2 := X2 - Y2; 622 X1 := X1 - Y1; 623 end Sub3; 624 625 ------------------------------- 626 -- Subtract_With_Ovflo_Check -- 627 ------------------------------- 628 629 function Subtract_With_Ovflo_Check (X, Y : Double_Int) return Double_Int is 630 R : constant Double_Int := To_Int (To_Uns (X) - To_Uns (Y)); 631 632 begin 633 if X >= 0 then 634 if Y > 0 or else R >= 0 then 635 return R; 636 end if; 637 638 else -- X < 0 639 if Y <= 0 or else R < 0 then 640 return R; 641 end if; 642 end if; 643 644 Raise_Error; 645 end Subtract_With_Ovflo_Check; 646 647 ---------------- 648 -- To_Neg_Int -- 649 ---------------- 650 651 function To_Neg_Int (A : Double_Uns) return Double_Int is 652 R : constant Double_Int := 653 (if A = 2 ** (Double_Size - 1) then Double_Int'First else -To_Int (A)); 654 -- Note that we can't just use the expression of the Else, because it 655 -- overflows for A = 2 ** (Double_Size - 1). 656 begin 657 if R <= 0 then 658 return R; 659 else 660 Raise_Error; 661 end if; 662 end To_Neg_Int; 663 664 ---------------- 665 -- To_Pos_Int -- 666 ---------------- 667 668 function To_Pos_Int (A : Double_Uns) return Double_Int is 669 R : constant Double_Int := To_Int (A); 670 begin 671 if R >= 0 then 672 return R; 673 else 674 Raise_Error; 675 end if; 676 end To_Pos_Int; 677 678end System.Arith_Double; 679