1------------------------------------------------------------------------------ 2-- -- 3-- GNAT COMPILER COMPONENTS -- 4-- -- 5-- S Y S T E M . B I G N U M S -- 6-- -- 7-- B o d y -- 8-- -- 9-- Copyright (C) 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. -- 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 32-- This package provides arbitrary precision signed integer arithmetic for 33-- use in computing intermediate values in expressions for the case where 34-- pragma Overflow_Check (Eliminate) is in effect. 35 36with System; use System; 37with System.Secondary_Stack; use System.Secondary_Stack; 38with System.Storage_Elements; use System.Storage_Elements; 39 40package body System.Bignums is 41 42 use Interfaces; 43 -- So that operations on Unsigned_32 are available 44 45 type DD is mod Base ** 2; 46 -- Double length digit used for intermediate computations 47 48 function MSD (X : DD) return SD is (SD (X / Base)); 49 function LSD (X : DD) return SD is (SD (X mod Base)); 50 -- Most significant and least significant digit of double digit value 51 52 function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y)); 53 -- Compose double digit value from two single digit values 54 55 subtype LLI is Long_Long_Integer; 56 57 One_Data : constant Digit_Vector (1 .. 1) := (1 => 1); 58 -- Constant one 59 60 Zero_Data : constant Digit_Vector (1 .. 0) := (1 .. 0 => 0); 61 -- Constant zero 62 63 ----------------------- 64 -- Local Subprograms -- 65 ----------------------- 66 67 function Add (X, Y : Digit_Vector; X_Neg, Y_Neg : Boolean) return Bignum 68 with Pre => X'First = 1 and then Y'First = 1; 69 -- This procedure adds two signed numbers returning the Sum, it is used 70 -- for both addition and subtraction. The value computed is X + Y, with 71 -- X_Neg and Y_Neg giving the signs of the operands. 72 73 function Allocate_Bignum (Len : Length) return Bignum 74 with Post => Allocate_Bignum'Result.Len = Len; 75 -- Allocate Bignum value of indicated length on secondary stack. On return 76 -- the Neg and D fields are left uninitialized. 77 78 type Compare_Result is (LT, EQ, GT); 79 -- Indicates result of comparison in following call 80 81 function Compare 82 (X, Y : Digit_Vector; 83 X_Neg, Y_Neg : Boolean) return Compare_Result 84 with Pre => X'First = 1 and then Y'First = 1; 85 -- Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the 86 -- result of the signed comparison. 87 88 procedure Div_Rem 89 (X, Y : Bignum; 90 Quotient : out Bignum; 91 Remainder : out Bignum; 92 Discard_Quotient : Boolean := False; 93 Discard_Remainder : Boolean := False); 94 -- Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The 95 -- values of X and Y are not modified. If Discard_Quotient is True, then 96 -- Quotient is undefined on return, and if Discard_Remainder is True, then 97 -- Remainder is undefined on return. Service routine for Big_Div/Rem/Mod. 98 99 procedure Free_Bignum (X : Bignum) is null; 100 -- Called to free a Bignum value used in intermediate computations. In 101 -- this implementation using the secondary stack, it does nothing at all, 102 -- because we rely on Mark/Release, but it may be of use for some 103 -- alternative implementation. 104 105 function Normalize 106 (X : Digit_Vector; 107 Neg : Boolean := False) return Bignum; 108 -- Given a digit vector and sign, allocate and construct a Bignum value. 109 -- Note that X may have leading zeroes which must be removed, and if the 110 -- result is zero, the sign is forced positive. 111 112 --------- 113 -- Add -- 114 --------- 115 116 function Add (X, Y : Digit_Vector; X_Neg, Y_Neg : Boolean) return Bignum is 117 begin 118 -- If signs are the same, we are doing an addition, it is convenient to 119 -- ensure that the first operand is the longer of the two. 120 121 if X_Neg = Y_Neg then 122 if X'Last < Y'Last then 123 return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg); 124 125 -- Here signs are the same, and the first operand is the longer 126 127 else 128 pragma Assert (X_Neg = Y_Neg and then X'Last >= Y'Last); 129 130 -- Do addition, putting result in Sum (allowing for carry) 131 132 declare 133 Sum : Digit_Vector (0 .. X'Last); 134 RD : DD; 135 136 begin 137 RD := 0; 138 for J in reverse 1 .. X'Last loop 139 RD := RD + DD (X (J)); 140 141 if J >= 1 + (X'Last - Y'Last) then 142 RD := RD + DD (Y (J - (X'Last - Y'Last))); 143 end if; 144 145 Sum (J) := LSD (RD); 146 RD := RD / Base; 147 end loop; 148 149 Sum (0) := SD (RD); 150 return Normalize (Sum, X_Neg); 151 end; 152 end if; 153 154 -- Signs are different so really this is a subtraction, we want to make 155 -- sure that the largest magnitude operand is the first one, and then 156 -- the result will have the sign of the first operand. 157 158 else 159 declare 160 CR : constant Compare_Result := Compare (X, Y, False, False); 161 162 begin 163 if CR = EQ then 164 return Normalize (Zero_Data); 165 166 elsif CR = LT then 167 return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg); 168 169 else 170 pragma Assert (X_Neg /= Y_Neg and then CR = GT); 171 172 -- Do subtraction, putting result in Diff 173 174 declare 175 Diff : Digit_Vector (1 .. X'Length); 176 RD : DD; 177 178 begin 179 RD := 0; 180 for J in reverse 1 .. X'Last loop 181 RD := RD + DD (X (J)); 182 183 if J >= 1 + (X'Last - Y'Last) then 184 RD := RD - DD (Y (J - (X'Last - Y'Last))); 185 end if; 186 187 Diff (J) := LSD (RD); 188 RD := (if RD < Base then 0 else -1); 189 end loop; 190 191 return Normalize (Diff, X_Neg); 192 end; 193 end if; 194 end; 195 end if; 196 end Add; 197 198 --------------------- 199 -- Allocate_Bignum -- 200 --------------------- 201 202 function Allocate_Bignum (Len : Length) return Bignum is 203 Addr : Address; 204 205 begin 206 -- Change the if False here to if True to get allocation on the heap 207 -- instead of the secondary stack, which is convenient for debugging 208 -- System.Bignum itself. 209 210 if False then 211 declare 212 B : Bignum; 213 begin 214 B := new Bignum_Data'(Len, False, (others => 0)); 215 return B; 216 end; 217 218 -- Normal case of allocation on the secondary stack 219 220 else 221 -- Note: The approach used here is designed to avoid strict aliasing 222 -- warnings that appeared previously using unchecked conversion. 223 224 SS_Allocate (Addr, Storage_Offset (4 + 4 * Len)); 225 226 declare 227 B : Bignum; 228 for B'Address use Addr'Address; 229 pragma Import (Ada, B); 230 231 BD : Bignum_Data (Len); 232 for BD'Address use Addr; 233 pragma Import (Ada, BD); 234 235 -- Expose a writable view of discriminant BD.Len so that we can 236 -- initialize it. We need to use the exact layout of the record 237 -- to ensure that the Length field has 24 bits as expected. 238 239 type Bignum_Data_Header is record 240 Len : Length; 241 Neg : Boolean; 242 end record; 243 244 for Bignum_Data_Header use record 245 Len at 0 range 0 .. 23; 246 Neg at 3 range 0 .. 7; 247 end record; 248 249 BDH : Bignum_Data_Header; 250 for BDH'Address use BD'Address; 251 pragma Import (Ada, BDH); 252 253 pragma Assert (BDH.Len'Size = BD.Len'Size); 254 255 begin 256 BDH.Len := Len; 257 return B; 258 end; 259 end if; 260 end Allocate_Bignum; 261 262 ------------- 263 -- Big_Abs -- 264 ------------- 265 266 function Big_Abs (X : Bignum) return Bignum is 267 begin 268 return Normalize (X.D); 269 end Big_Abs; 270 271 ------------- 272 -- Big_Add -- 273 ------------- 274 275 function Big_Add (X, Y : Bignum) return Bignum is 276 begin 277 return Add (X.D, Y.D, X.Neg, Y.Neg); 278 end Big_Add; 279 280 ------------- 281 -- Big_Div -- 282 ------------- 283 284 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result 285 -- varies with the signs of the operands. 286 287 -- A B A/B A B A/B 288 -- 289 -- 10 5 2 -10 5 -2 290 -- 11 5 2 -11 5 -2 291 -- 12 5 2 -12 5 -2 292 -- 13 5 2 -13 5 -2 293 -- 14 5 2 -14 5 -2 294 -- 295 -- A B A/B A B A/B 296 -- 297 -- 10 -5 -2 -10 -5 2 298 -- 11 -5 -2 -11 -5 2 299 -- 12 -5 -2 -12 -5 2 300 -- 13 -5 -2 -13 -5 2 301 -- 14 -5 -2 -14 -5 2 302 303 function Big_Div (X, Y : Bignum) return Bignum is 304 Q, R : Bignum; 305 begin 306 Div_Rem (X, Y, Q, R, Discard_Remainder => True); 307 Q.Neg := Q.Len > 0 and then (X.Neg xor Y.Neg); 308 return Q; 309 end Big_Div; 310 311 ------------- 312 -- Big_Exp -- 313 ------------- 314 315 function Big_Exp (X, Y : Bignum) return Bignum is 316 317 function "**" (X : Bignum; Y : SD) return Bignum; 318 -- Internal routine where we know right operand is one word 319 320 ---------- 321 -- "**" -- 322 ---------- 323 324 function "**" (X : Bignum; Y : SD) return Bignum is 325 begin 326 case Y is 327 328 -- X ** 0 is 1 329 330 when 0 => 331 return Normalize (One_Data); 332 333 -- X ** 1 is X 334 335 when 1 => 336 return Normalize (X.D); 337 338 -- X ** 2 is X * X 339 340 when 2 => 341 return Big_Mul (X, X); 342 343 -- For X greater than 2, use the recursion 344 345 -- X even, X ** Y = (X ** (Y/2)) ** 2; 346 -- X odd, X ** Y = (X ** (Y/2)) ** 2 * X; 347 348 when others => 349 declare 350 XY2 : constant Bignum := X ** (Y / 2); 351 XY2S : constant Bignum := Big_Mul (XY2, XY2); 352 Res : Bignum; 353 354 begin 355 Free_Bignum (XY2); 356 357 -- Raise storage error if intermediate value is getting too 358 -- large, which we arbitrarily define as 200 words for now! 359 360 if XY2S.Len > 200 then 361 Free_Bignum (XY2S); 362 raise Storage_Error with 363 "exponentiation result is too large"; 364 end if; 365 366 -- Otherwise take care of even/odd cases 367 368 if (Y and 1) = 0 then 369 return XY2S; 370 371 else 372 Res := Big_Mul (XY2S, X); 373 Free_Bignum (XY2S); 374 return Res; 375 end if; 376 end; 377 end case; 378 end "**"; 379 380 -- Start of processing for Big_Exp 381 382 begin 383 -- Error if right operand negative 384 385 if Y.Neg then 386 raise Constraint_Error with "exponentiation to negative power"; 387 388 -- X ** 0 is always 1 (including 0 ** 0, so do this test first) 389 390 elsif Y.Len = 0 then 391 return Normalize (One_Data); 392 393 -- 0 ** X is always 0 (for X non-zero) 394 395 elsif X.Len = 0 then 396 return Normalize (Zero_Data); 397 398 -- (+1) ** Y = 1 399 -- (-1) ** Y = +/-1 depending on whether Y is even or odd 400 401 elsif X.Len = 1 and then X.D (1) = 1 then 402 return Normalize 403 (X.D, Neg => X.Neg and then ((Y.D (Y.Len) and 1) = 1)); 404 405 -- If the absolute value of the base is greater than 1, then the 406 -- exponent must not be bigger than one word, otherwise the result 407 -- is ludicrously large, and we just signal Storage_Error right away. 408 409 elsif Y.Len > 1 then 410 raise Storage_Error with "exponentiation result is too large"; 411 412 -- Special case (+/-)2 ** K, where K is 1 .. 31 using a shift 413 414 elsif X.Len = 1 and then X.D (1) = 2 and then Y.D (1) < 32 then 415 declare 416 D : constant Digit_Vector (1 .. 1) := 417 (1 => Shift_Left (SD'(1), Natural (Y.D (1)))); 418 begin 419 return Normalize (D, X.Neg); 420 end; 421 422 -- Remaining cases have right operand of one word 423 424 else 425 return X ** Y.D (1); 426 end if; 427 end Big_Exp; 428 429 ------------ 430 -- Big_EQ -- 431 ------------ 432 433 function Big_EQ (X, Y : Bignum) return Boolean is 434 begin 435 return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ; 436 end Big_EQ; 437 438 ------------ 439 -- Big_GE -- 440 ------------ 441 442 function Big_GE (X, Y : Bignum) return Boolean is 443 begin 444 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT; 445 end Big_GE; 446 447 ------------ 448 -- Big_GT -- 449 ------------ 450 451 function Big_GT (X, Y : Bignum) return Boolean is 452 begin 453 return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT; 454 end Big_GT; 455 456 ------------ 457 -- Big_LE -- 458 ------------ 459 460 function Big_LE (X, Y : Bignum) return Boolean is 461 begin 462 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT; 463 end Big_LE; 464 465 ------------ 466 -- Big_LT -- 467 ------------ 468 469 function Big_LT (X, Y : Bignum) return Boolean is 470 begin 471 return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT; 472 end Big_LT; 473 474 ------------- 475 -- Big_Mod -- 476 ------------- 477 478 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result 479 -- of Rem and Mod vary with the signs of the operands. 480 481 -- A B A mod B A rem B A B A mod B A rem B 482 483 -- 10 5 0 0 -10 5 0 0 484 -- 11 5 1 1 -11 5 4 -1 485 -- 12 5 2 2 -12 5 3 -2 486 -- 13 5 3 3 -13 5 2 -3 487 -- 14 5 4 4 -14 5 1 -4 488 489 -- A B A mod B A rem B A B A mod B A rem B 490 491 -- 10 -5 0 0 -10 -5 0 0 492 -- 11 -5 -4 1 -11 -5 -1 -1 493 -- 12 -5 -3 2 -12 -5 -2 -2 494 -- 13 -5 -2 3 -13 -5 -3 -3 495 -- 14 -5 -1 4 -14 -5 -4 -4 496 497 function Big_Mod (X, Y : Bignum) return Bignum is 498 Q, R : Bignum; 499 500 begin 501 -- If signs are same, result is same as Rem 502 503 if X.Neg = Y.Neg then 504 return Big_Rem (X, Y); 505 506 -- Case where Mod is different 507 508 else 509 -- Do division 510 511 Div_Rem (X, Y, Q, R, Discard_Quotient => True); 512 513 -- Zero result is unchanged 514 515 if R.Len = 0 then 516 return R; 517 518 -- Otherwise adjust result 519 520 else 521 declare 522 T1 : constant Bignum := Big_Sub (Y, R); 523 begin 524 T1.Neg := Y.Neg; 525 Free_Bignum (R); 526 return T1; 527 end; 528 end if; 529 end if; 530 end Big_Mod; 531 532 ------------- 533 -- Big_Mul -- 534 ------------- 535 536 function Big_Mul (X, Y : Bignum) return Bignum is 537 Result : Digit_Vector (1 .. X.Len + Y.Len) := (others => 0); 538 -- Accumulate result (max length of result is sum of operand lengths) 539 540 L : Length; 541 -- Current result digit 542 543 D : DD; 544 -- Result digit 545 546 begin 547 for J in 1 .. X.Len loop 548 for K in 1 .. Y.Len loop 549 L := Result'Last - (X.Len - J) - (Y.Len - K); 550 D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L)); 551 Result (L) := LSD (D); 552 D := D / Base; 553 554 -- D is carry which must be propagated 555 556 while D /= 0 and then L >= 1 loop 557 L := L - 1; 558 D := D + DD (Result (L)); 559 Result (L) := LSD (D); 560 D := D / Base; 561 end loop; 562 563 -- Must not have a carry trying to extend max length 564 565 pragma Assert (D = 0); 566 end loop; 567 end loop; 568 569 -- Return result 570 571 return Normalize (Result, X.Neg xor Y.Neg); 572 end Big_Mul; 573 574 ------------ 575 -- Big_NE -- 576 ------------ 577 578 function Big_NE (X, Y : Bignum) return Boolean is 579 begin 580 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ; 581 end Big_NE; 582 583 ------------- 584 -- Big_Neg -- 585 ------------- 586 587 function Big_Neg (X : Bignum) return Bignum is 588 begin 589 return Normalize (X.D, not X.Neg); 590 end Big_Neg; 591 592 ------------- 593 -- Big_Rem -- 594 ------------- 595 596 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result 597 -- varies with the signs of the operands. 598 599 -- A B A rem B A B A rem B 600 601 -- 10 5 0 -10 5 0 602 -- 11 5 1 -11 5 -1 603 -- 12 5 2 -12 5 -2 604 -- 13 5 3 -13 5 -3 605 -- 14 5 4 -14 5 -4 606 607 -- A B A rem B A B A rem B 608 609 -- 10 -5 0 -10 -5 0 610 -- 11 -5 1 -11 -5 -1 611 -- 12 -5 2 -12 -5 -2 612 -- 13 -5 3 -13 -5 -3 613 -- 14 -5 4 -14 -5 -4 614 615 function Big_Rem (X, Y : Bignum) return Bignum is 616 Q, R : Bignum; 617 begin 618 Div_Rem (X, Y, Q, R, Discard_Quotient => True); 619 R.Neg := R.Len > 0 and then X.Neg; 620 return R; 621 end Big_Rem; 622 623 ------------- 624 -- Big_Sub -- 625 ------------- 626 627 function Big_Sub (X, Y : Bignum) return Bignum is 628 begin 629 -- If right operand zero, return left operand (avoiding sharing) 630 631 if Y.Len = 0 then 632 return Normalize (X.D, X.Neg); 633 634 -- Otherwise add negative of right operand 635 636 else 637 return Add (X.D, Y.D, X.Neg, not Y.Neg); 638 end if; 639 end Big_Sub; 640 641 ------------- 642 -- Compare -- 643 ------------- 644 645 function Compare 646 (X, Y : Digit_Vector; 647 X_Neg, Y_Neg : Boolean) return Compare_Result 648 is 649 begin 650 -- Signs are different, that's decisive, since 0 is always plus 651 652 if X_Neg /= Y_Neg then 653 return (if X_Neg then LT else GT); 654 655 -- Lengths are different, that's decisive since no leading zeroes 656 657 elsif X'Last /= Y'Last then 658 return (if (X'Last > Y'Last) xor X_Neg then GT else LT); 659 660 -- Need to compare data 661 662 else 663 for J in X'Range loop 664 if X (J) /= Y (J) then 665 return (if (X (J) > Y (J)) xor X_Neg then GT else LT); 666 end if; 667 end loop; 668 669 return EQ; 670 end if; 671 end Compare; 672 673 ------------- 674 -- Div_Rem -- 675 ------------- 676 677 procedure Div_Rem 678 (X, Y : Bignum; 679 Quotient : out Bignum; 680 Remainder : out Bignum; 681 Discard_Quotient : Boolean := False; 682 Discard_Remainder : Boolean := False) 683 is 684 begin 685 -- Error if division by zero 686 687 if Y.Len = 0 then 688 raise Constraint_Error with "division by zero"; 689 end if; 690 691 -- Handle simple cases with special tests 692 693 -- If X < Y then quotient is zero and remainder is X 694 695 if Compare (X.D, Y.D, False, False) = LT then 696 Remainder := Normalize (X.D); 697 Quotient := Normalize (Zero_Data); 698 return; 699 700 -- If both X and Y are less than 2**63-1, we can use Long_Long_Integer 701 -- arithmetic. Note it is good not to do an accurate range check against 702 -- Long_Long_Integer since -2**63 / -1 overflows! 703 704 elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31)) 705 and then 706 (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31)) 707 then 708 declare 709 A : constant LLI := abs (From_Bignum (X)); 710 B : constant LLI := abs (From_Bignum (Y)); 711 begin 712 Quotient := To_Bignum (A / B); 713 Remainder := To_Bignum (A rem B); 714 return; 715 end; 716 717 -- Easy case if divisor is one digit 718 719 elsif Y.Len = 1 then 720 declare 721 ND : DD; 722 Div : constant DD := DD (Y.D (1)); 723 724 Result : Digit_Vector (1 .. X.Len); 725 Remdr : Digit_Vector (1 .. 1); 726 727 begin 728 ND := 0; 729 for J in 1 .. X.Len loop 730 ND := Base * ND + DD (X.D (J)); 731 Result (J) := SD (ND / Div); 732 ND := ND rem Div; 733 end loop; 734 735 Quotient := Normalize (Result); 736 Remdr (1) := SD (ND); 737 Remainder := Normalize (Remdr); 738 return; 739 end; 740 end if; 741 742 -- The complex full multi-precision case. We will employ algorithm 743 -- D defined in the section "The Classical Algorithms" (sec. 4.3.1) 744 -- of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd 745 -- edition. The terminology is adjusted for this section to match that 746 -- reference. 747 748 -- We are dividing X.Len digits of X (called u here) by Y.Len digits 749 -- of Y (called v here), developing the quotient and remainder. The 750 -- numbers are represented using Base, which was chosen so that we have 751 -- the operations of multiplying to single digits (SD) to form a double 752 -- digit (DD), and dividing a double digit (DD) by a single digit (SD) 753 -- to give a single digit quotient and a single digit remainder. 754 755 -- Algorithm D from Knuth 756 757 -- Comments here with square brackets are directly from Knuth 758 759 Algorithm_D : declare 760 761 -- The following lower case variables correspond exactly to the 762 -- terminology used in algorithm D. 763 764 m : constant Length := X.Len - Y.Len; 765 n : constant Length := Y.Len; 766 b : constant DD := Base; 767 768 u : Digit_Vector (0 .. m + n); 769 v : Digit_Vector (1 .. n); 770 q : Digit_Vector (0 .. m); 771 r : Digit_Vector (1 .. n); 772 773 u0 : SD renames u (0); 774 v1 : SD renames v (1); 775 v2 : SD renames v (2); 776 777 d : DD; 778 j : Length; 779 qhat : DD; 780 rhat : DD; 781 temp : DD; 782 783 begin 784 -- Initialize data of left and right operands 785 786 for J in 1 .. m + n loop 787 u (J) := X.D (J); 788 end loop; 789 790 for J in 1 .. n loop 791 v (J) := Y.D (J); 792 end loop; 793 794 -- [Division of nonnegative integers.] Given nonnegative integers u 795 -- = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we 796 -- form the quotient u / v = (q0,ql..qm) and the remainder u mod v = 797 -- (r1,r2..rn). 798 799 pragma Assert (v1 /= 0); 800 pragma Assert (n > 1); 801 802 -- Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n) 803 -- equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to 804 -- (v1,v2..vn) times d. Note the introduction of a new digit position 805 -- u0 at the left of u1; if d = 1 all we need to do in this step is 806 -- to set u0 = 0. 807 808 d := b / (DD (v1) + 1); 809 810 if d = 1 then 811 u0 := 0; 812 813 else 814 declare 815 Carry : DD; 816 Tmp : DD; 817 818 begin 819 -- Multiply Dividend (u) by d 820 821 Carry := 0; 822 for J in reverse 1 .. m + n loop 823 Tmp := DD (u (J)) * d + Carry; 824 u (J) := LSD (Tmp); 825 Carry := Tmp / Base; 826 end loop; 827 828 u0 := SD (Carry); 829 830 -- Multiply Divisor (v) by d 831 832 Carry := 0; 833 for J in reverse 1 .. n loop 834 Tmp := DD (v (J)) * d + Carry; 835 v (J) := LSD (Tmp); 836 Carry := Tmp / Base; 837 end loop; 838 839 pragma Assert (Carry = 0); 840 end; 841 end if; 842 843 -- D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7, 844 -- will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn) 845 -- to get a single quotient digit qj. 846 847 j := 0; 848 849 -- Loop through digits 850 851 loop 852 -- Note: In the original printing, step D3 was as follows: 853 854 -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise 855 -- set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than 856 -- (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and 857 -- repeat this test 858 859 -- This had a bug not discovered till 1995, see Vol 2 errata: 860 -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under 861 -- rare circumstances the expression in the test could overflow. 862 -- This version was further corrected in 2005, see Vol 2 errata: 863 -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz. 864 -- The code below is the fixed version of this step. 865 866 -- D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to 867 -- to (uj,uj+1) mod v1. 868 869 temp := u (j) & u (j + 1); 870 qhat := temp / DD (v1); 871 rhat := temp mod DD (v1); 872 873 -- D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2): 874 -- if so, decrease qhat by 1, increase rhat by v1, and repeat this 875 -- test if rhat < b. [The test on v2 determines at at high speed 876 -- most of the cases in which the trial value qhat is one too 877 -- large, and eliminates all cases where qhat is two too large.] 878 879 while qhat >= b 880 or else DD (v2) * qhat > LSD (rhat) & u (j + 2) 881 loop 882 qhat := qhat - 1; 883 rhat := rhat + DD (v1); 884 exit when rhat >= b; 885 end loop; 886 887 -- D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by 888 -- (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step 889 -- consists of a simple multiplication by a one-place number, 890 -- combined with a subtraction. 891 892 -- The digits (uj,uj+1..uj+n) are always kept positive; if the 893 -- result of this step is actually negative then (uj,uj+1..uj+n) 894 -- is left as the true value plus b**(n+1), i.e. as the b's 895 -- complement of the true value, and a "borrow" to the left is 896 -- remembered. 897 898 declare 899 Borrow : SD; 900 Carry : DD; 901 Temp : DD; 902 903 Negative : Boolean; 904 -- Records if subtraction causes a negative result, requiring 905 -- an add back (case where qhat turned out to be 1 too large). 906 907 begin 908 Borrow := 0; 909 for K in reverse 1 .. n loop 910 Temp := qhat * DD (v (K)) + DD (Borrow); 911 Borrow := MSD (Temp); 912 913 if LSD (Temp) > u (j + K) then 914 Borrow := Borrow + 1; 915 end if; 916 917 u (j + K) := u (j + K) - LSD (Temp); 918 end loop; 919 920 Negative := u (j) < Borrow; 921 u (j) := u (j) - Borrow; 922 923 -- D5. [Test remainder.] Set qj = qhat. If the result of step 924 -- D4 was negative, we will do the add back step (step D6). 925 926 q (j) := LSD (qhat); 927 928 if Negative then 929 930 -- D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn) 931 -- to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left 932 -- of uj, and it is be ignored since it cancels with the 933 -- borrow that occurred in D4.) 934 935 q (j) := q (j) - 1; 936 937 Carry := 0; 938 for K in reverse 1 .. n loop 939 Temp := DD (v (K)) + DD (u (j + K)) + Carry; 940 u (j + K) := LSD (Temp); 941 Carry := Temp / Base; 942 end loop; 943 944 u (j) := u (j) + SD (Carry); 945 end if; 946 end; 947 948 -- D7. [Loop on j.] Increase j by one. Now if j <= m, go back to 949 -- D3 (the start of the loop on j). 950 951 j := j + 1; 952 exit when not (j <= m); 953 end loop; 954 955 -- D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and 956 -- the desired remainder may be obtained by dividing (um+1..um+n) 957 -- by d. 958 959 if not Discard_Quotient then 960 Quotient := Normalize (q); 961 end if; 962 963 if not Discard_Remainder then 964 declare 965 Remdr : DD; 966 967 begin 968 Remdr := 0; 969 for K in 1 .. n loop 970 Remdr := Base * Remdr + DD (u (m + K)); 971 r (K) := SD (Remdr / d); 972 Remdr := Remdr rem d; 973 end loop; 974 975 pragma Assert (Remdr = 0); 976 end; 977 978 Remainder := Normalize (r); 979 end if; 980 end Algorithm_D; 981 end Div_Rem; 982 983 ----------------- 984 -- From_Bignum -- 985 ----------------- 986 987 function From_Bignum (X : Bignum) return Long_Long_Integer is 988 begin 989 if X.Len = 0 then 990 return 0; 991 992 elsif X.Len = 1 then 993 return (if X.Neg then -LLI (X.D (1)) else LLI (X.D (1))); 994 995 elsif X.Len = 2 then 996 declare 997 Mag : constant DD := X.D (1) & X.D (2); 998 begin 999 if X.Neg and then Mag <= 2 ** 63 then 1000 return -LLI (Mag); 1001 elsif Mag < 2 ** 63 then 1002 return LLI (Mag); 1003 end if; 1004 end; 1005 end if; 1006 1007 raise Constraint_Error with "expression value out of range"; 1008 end From_Bignum; 1009 1010 ------------------------- 1011 -- Bignum_In_LLI_Range -- 1012 ------------------------- 1013 1014 function Bignum_In_LLI_Range (X : Bignum) return Boolean is 1015 begin 1016 -- If length is 0 or 1, definitely fits 1017 1018 if X.Len <= 1 then 1019 return True; 1020 1021 -- If length is greater than 2, definitely does not fit 1022 1023 elsif X.Len > 2 then 1024 return False; 1025 1026 -- Length is 2, more tests needed 1027 1028 else 1029 declare 1030 Mag : constant DD := X.D (1) & X.D (2); 1031 begin 1032 return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63); 1033 end; 1034 end if; 1035 end Bignum_In_LLI_Range; 1036 1037 --------------- 1038 -- Normalize -- 1039 --------------- 1040 1041 function Normalize 1042 (X : Digit_Vector; 1043 Neg : Boolean := False) return Bignum 1044 is 1045 B : Bignum; 1046 J : Length; 1047 1048 begin 1049 J := X'First; 1050 while J <= X'Last and then X (J) = 0 loop 1051 J := J + 1; 1052 end loop; 1053 1054 B := Allocate_Bignum (X'Last - J + 1); 1055 B.Neg := B.Len > 0 and then Neg; 1056 B.D := X (J .. X'Last); 1057 return B; 1058 end Normalize; 1059 1060 --------------- 1061 -- To_Bignum -- 1062 --------------- 1063 1064 function To_Bignum (X : Long_Long_Integer) return Bignum is 1065 R : Bignum; 1066 1067 begin 1068 if X = 0 then 1069 R := Allocate_Bignum (0); 1070 1071 -- One word result 1072 1073 elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then 1074 R := Allocate_Bignum (1); 1075 R.D (1) := SD (abs (X)); 1076 1077 -- Largest negative number annoyance 1078 1079 elsif X = Long_Long_Integer'First then 1080 R := Allocate_Bignum (2); 1081 R.D (1) := 2 ** 31; 1082 R.D (2) := 0; 1083 1084 -- Normal two word case 1085 1086 else 1087 R := Allocate_Bignum (2); 1088 R.D (2) := SD (abs (X) mod Base); 1089 R.D (1) := SD (abs (X) / Base); 1090 end if; 1091 1092 R.Neg := X < 0; 1093 return R; 1094 end To_Bignum; 1095 1096end System.Bignums; 1097