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