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