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