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