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