1 /*==============================================================================
   2 |
   3 |  NAME
   4 |
   5 |      ppBigInt.cpp
   6 |
   7 |  DESCRIPTION
   8 |
   9 |      Classes for non-negative multiple precision integer arithmetic.
  10 |
  11 |      User manual and technical documentation are described in detail in my web page at
  12 |      http://seanerikoconnor.freeservers.com/Mathematics/AbstractAlgebra/PrimitivePolynomials/overview.html
  13 |
  14 |  NOTES
  15 |
  16 |      The algorithms are based on
  17 |
  18 |          D. E. Knuth, THE ART OF COMPUTER PROGRAMMING, VOL. 2, 3rd ed.,
  19 |          Addison-Wesley, 1981, pgs. 250-265.
  20 |
  21 |          Errata for Volume 2:
  22 |          http://www-cs-faculty.stanford.edu/~knuth/taocp.html
  23 |
  24 |      Member functions are exception-safe:  operations
  25 |      leave arithmetic objects in a valid state when exceptions
  26 |      are thrown and no memory is leaked.  When possible, the
  27 |      operations either succeed or else leave the object unchanged after
  28 |      throwing an exception.  In general, we do this by constructing
  29 |      a new representation completely without risk of exceptions
  30 |      before disposing of the old one, releasing any resources
  31 |      before we throw, making sure all objects are in a valid
  32 |      (free-able) state before throwing.  Constructors who throw
  33 |      don't leave any object behind and destructors shouldn't
  34 |      throw at all.  See "Standard Library Exception Safety" in
  35 |
  36 |          The C++ Programming Language, Special Edition, Bjarne
  37 |          Stroustrup 2000 AT&T, Addison-Wesley, Appendix E.
  38 |
  39 |  LEGAL
  40 |
  41 |     Primpoly Version 16.1 - A Program for Computing Primitive Polynomials.
  42 |     Copyright (C) 1999-2021 by Sean Erik O'Connor.  All Rights Reserved.
  43 |
  44 |     This program is free software: you can redistribute it and/or modify
  45 |     it under the terms of the GNU General Public License as published by
  46 |     the Free Software Foundation, either version 3 of the License, or
  47 |     (at your option) any later version.
  48 |
  49 |     This program is distributed in the hope that it will be useful,
  50 |     but WITHOUT ANY WARRANTY; without even the implied warranty of
  51 |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  52 |     GNU General Public License for more details.
  53 |
  54 |     You should have received a copy of the GNU General Public License
  55 |     along with this program.  If not, see http://www.gnu.org/licenses/.
  56 |     
  57 |     The author's address is seanerikoconnor!AT!gmail!DOT!com
  58 |     with the !DOT! replaced by . and the !AT! replaced by @
  59 |
  60 ==============================================================================*/
  61 
  62 
  63 /*------------------------------------------------------------------------------
  64 |                                Includes                                      |
  65 ------------------------------------------------------------------------------*/
  66 
  67 #include <cstdlib>      // abort()
  68 #include <iostream>     // Basic stream I/O.
  69 #include <new>          // set_new_handler()
  70 #include <limits>       // Numeric limits.
  71 #include <cmath>        // Basic math functions e.g. sqrt()
  72 #include <complex>      // Complex data type and operations.
  73 #include <fstream>      // File stream I/O.
  74 #include <sstream>      // String stream I/O.
  75 #include <vector>       // STL vector class.
  76 #include <string>       // STL string class.
  77 #include <algorithm>    // Iterators.
  78 #include <stdexcept>    // Exceptions.
  79 #include <cassert>      // assert()
  80 
  81 using namespace std ;   // So we don't need to say std::vector everywhere.
  82 
  83 
  84 
  85 /*------------------------------------------------------------------------------
  86 |                                PP Include Files                              |
  87 ------------------------------------------------------------------------------*/
  88 
  89 #include "Primpoly.h"          // Global functions.
  90 #include "ppArith.h"           // Basic arithmetic functions.
  91 #include "ppBigInt.h"          // Arbitrary precision integer arithmetic.
  92 #include "ppOperationCount.h"  // OperationCount collection for factoring and poly finding.
  93 #include "ppFactor.h"          // Prime factorization and Euler Phi.
  94 #include "ppPolynomial.h"      // Polynomial operations and mod polynomial operations.
  95 #include "ppParser.h"          // Parsing of polynomials and I/O services.
  96 #include "ppUnitTest.h"        // Complete unit test.
  97 
  98 
  99 
 100 /*------------------------------------------------------------------------------
 101 |                     Initialize Class Variables                               |
 102 ------------------------------------------------------------------------------*/
 103 
 104 // Pointer to number system base.  Used for unit testing only.
 105 ppuint * BigInt::pbase = 0 ;
 106 
 107 
 108 
 109 /*=============================================================================
 110 |
 111 | NAME
 112 |
 113 |     BigInt::BigInt()
 114 |
 115 | DESCRIPTION
 116 |
 117 |    Default construct a multiple precision integer u with no digits.
 118 |
 119 | EXAMPLE
 120 |
 121 |    Called during default construction e.g.
 122 |        BigInt n ;
 123 |        foo( BigInt x )
 124 |
 125 +============================================================================*/
 126 
 127 BigInt::BigInt()
 128        : digit_( 0 )                // Construct a vector with zero length.
 129 {
 130     // Make certain we have zero digits.
 131     digit_.clear() ;
 132 }
 133 
 134 
 135 
 136 
 137 /*=============================================================================
 138 |
 139 | NAME
 140 |
 141 |     ~BigInt
 142 |
 143 | DESCRIPTION
 144 |
 145 |    Default destructor for a BigInt type.
 146 |
 147 | EXAMPLE
 148 |
 149 |    Called upon leaving scope:
 150 |    {
 151 |         BigInt n ;
 152 |    } <--- called here.
 153 |
 154 +============================================================================*/
 155 
 156 BigInt::~BigInt()
 157 {
 158     // Digits are in a vector which frees itself when it goes out of scope.
 159     // Other class variables are primitive types.
 160     // We obviously don't throw any exceptions from this destructor so
 161     // terminate() won't ever be called.
 162 }
 163 
 164 
 165 
 166 /*=============================================================================
 167 |
 168 | NAME
 169 |
 170 |     BigInt::BigInt( ppuint d )
 171 |
 172 | DESCRIPTION
 173 |
 174 |     Construct a BigInt from an unsigned 64-bit integer.
 175 |
 176 | EXAMPLE
 177 |
 178 |    try
 179 |    {
 180 |        ppuint d{ 123 } ;
 181 |        BigInt u( d ) ;
 182 |    }
 183 |    catch( BigIntMathError & e )
 184 |    {
 185 |        cerr << e.what() << endl ;
 186 |    }
 187 |
 188 | METHOD
 189 |
 190 |    The unsigned integer d written in BigInt format using
 191 |    base b is
 192 |                                         n-1
 193 |       d  =  (u     . . . u )   = u    b    + ... + u  b +  u
 194 |               n-1         0  b    n-1               1       0
 195 |   
 196 |       The first digit of u is
 197 |   
 198 |       d mod b = u
 199 |                  0
 200 |
 201 |       after which the remainder is
 202 |
 203 |                                n-2
 204 |           | d / b |  =  u    b    + ... + u
 205 |           --     --      n-1               1
 206 |   
 207 |       then we set d <-- | d / b |  and repeat to get u , etc.
 208 |                         --     --                     1
 209 |   
 210 |       So the algorithm for extracting the BigInt digits is
 211 |   
 212 |       do until d == 0:
 213 |   
 214 |           U = d mod b
 215 |   
 216 |           u = |  d / b |
 217 |               --      --
 218 |
 219 +============================================================================*/
 220 
 221 BigInt::BigInt( const ppuint d )
 222        : digit_( 0 )                // Construct a vector with zero length.
 223 {
 224     ppuint b{ base_() } ;
 225 
 226     try
 227     {
 228         // Early return if d = 0;  d will never be negative.
 229         if (ppuint d2{ d } ; d2 == 0)
 230         {
 231             digit_.push_back( 0 ) ;
 232         }
 233         else
 234             while (d2 != 0)
 235             {
 236                 ppuint digit{ d2 % b } ;
 237 
 238                 // Can throw from argument operator=().
 239                 digit_.push_back( digit ) ;
 240 
 241                 d2 = d2 / b ;
 242             }
 243     }
 244     catch( bad_alloc & e )
 245     {
 246         ostringstream os ;
 247         os << "BigInt::BigInt( const ppuint d ) constructor had a bad_alloc memory allocation exception"
 248            << "while adding digits to BigInt" ;
 249         throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
 250     }
 251 }
 252 
 253 
 254 /*=============================================================================
 255  |
 256  | NAME
 257  |
 258  |     BigInt::BigInt( ppuint32 d )
 259  |
 260  | DESCRIPTION
 261  |
 262  |     Construct a BigInt from an unsigned 32-bit integer.
 263  |
 264  | EXAMPLE
 265  |
 266  |    try
 267  |    {
 268  |        ppuint32 d{ 123 } ;
 269  |        BigInt u( d ) ;
 270  |    }
 271  |    catch( BigIntMathError & e )
 272  |    {
 273  |        cerr << e.what() << endl ;
 274  |    }
 275  |
 276  +============================================================================*/
 277 
 278 BigInt::BigInt( const ppuint32 d )
 279 : digit_( 0 )                // Construct a vector with zero length.
 280 {
 281     // Cast the 32-bit number into a 64-bit number without any loss of precision,
 282     // and then into a BigInt.
 283     BigInt w { static_cast<ppuint>( d ) } ;
 284 
 285     // Swap the digits of w into this number.
 286     swap( w.digit_, digit_ ) ;
 287 }
 288 
 289 
 290 /*=============================================================================
 291 |
 292 | NAME
 293 |
 294 |     BigInt::BigInt( string )
 295 |
 296 | DESCRIPTION
 297 |
 298 |  Construct a BigInt
 299 |
 300 |   (u    ... u  )
 301 |     n-1      0  b
 302 |
 303 |   from a decimal string
 304 |
 305 |   (U    ... U  )
 306 |     m-1      0  10
 307 |
 308 |
 309 |  EXAMPLE
 310 |
 311 |    BigInt v( "1234567890" ) ;
 312 |
 313 +============================================================================*/
 314 
 315 BigInt::BigInt( const string & s )
 316        : digit_( 0 )                  // Construct a vector with zero length.
 317 {
 318     // Construct temporary big integer w = 0 or else throw exception upwards.
 319     // If this fails, memory for w is automatically released during stack unwind,
 320     // this BigInt object isn't constructed.
 321     BigInt w( 0u ) ;
 322 
 323     #ifdef DEBUG_PP_BIGINT
 324     cout << "BigInt( string ) digits " ;
 325     #endif
 326 
 327     //  Use Horner's rule on base 10 numerals to convert decimal string
 328     //  of digits to base b.  e.g. 123 = 10 * (10 * ((10 * 0) + 1) + 2) + 3
 329     for (auto & c : s)
 330     {
 331         w *= static_cast<ppuint>( 10u ) ;
 332 
 333         // This only works for ASCII characters.
 334         int digit = 0 ;
 335         if (isdigit( c ))
 336         {
 337             char asciiDigit[ 2 ] ;
 338             asciiDigit[ 0 ] = c ;
 339             asciiDigit[ 1 ] = '\0' ;
 340             digit = atoi( asciiDigit ) ;
 341         }
 342         else
 343         {
 344             ostringstream os ;
 345             os << "BigInt::BigInt( string ):  range error from character = " << c ;
 346             throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
 347         }
 348 
 349         #ifdef DEBUG_PP_BIGINT
 350         cout << digit << " " ;
 351         #endif
 352 
 353         w += static_cast<ppuint>( digit ) ;
 354     }
 355 
 356     #ifdef DEBUG_PP_BIGINT
 357     cout << endl ;
 358     #endif
 359 
 360     // Swap the digits of w into this number.
 361     swap( w.digit_, digit_ ) ;
 362 
 363     // As we go out of scope, w's destructor will be called.
 364 }
 365 
 366 
 367 
 368 /*=============================================================================
 369 |
 370 | NAME
 371 |
 372 |     BigInt::BigInt
 373 |
 374 | DESCRIPTION
 375 |
 376 |     Copy constructor.
 377 |
 378 | EXAMPLE
 379 |
 380 |        BigInt u, v ;
 381 |        BigInt u( v ) ;
 382 |
 383 +============================================================================*/
 384 
 385 BigInt::BigInt( const BigInt & u )
 386        : digit_( u.digit_ )
 387 {
 388 
 389 }
 390 
 391 
 392 
 393 /*=============================================================================
 394 |
 395 | NAME
 396 |
 397 |     BigInt::operator=
 398 |
 399 | DESCRIPTION
 400 |
 401 |    Assignment operator.
 402 |
 403 | EXAMPLE
 404 |
 405 |    BigInt m = n ;
 406 |
 407 +============================================================================*/
 408 
 409 BigInt & BigInt::operator=( const BigInt & n )
 410 {
 411     // Check for assigning to oneself:  just pass back a reference to the unchanged object.
 412     if (this == &n)
 413         return *this ;
 414 
 415     try
 416     {
 417         vector<ppuint> tempDigits( n.digit_ ) ;
 418 
 419         // Move the old values into the temporary, and the new values into the object.
 420         // The temporary containing the old values will be destroyed when we leave scope.
 421         swap( digit_, tempDigits ) ;
 422     }
 423     catch( bad_alloc & e )
 424     {
 425         ostringstream os ;
 426         os << "BigInt::operator=( const BigInt & n ) had a bad_alloc memory allocation exception" ;
 427         throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
 428     }
 429 
 430     // Return a reference to assigned object to make chaining possible.
 431     return *this ;
 432 }
 433 
 434 
 435 
 436 /*=============================================================================
 437 |
 438 | NAME
 439 |
 440 |     operator ppuint
 441 |
 442 | DESCRIPTION
 443 |
 444 |    Convert an n-place number
 445 |
 446 |    (u    ... u  )
 447 |      n-1      0  b
 448 |
 449 |    as an unsigned integer.  If the number is too
 450 |    large, throw a BigIntOverflow exception.
 451 |
 452 +============================================================================*/
 453 
 454 BigInt::operator ppuint() const
 455 {
 456     ppuint result{ 0 } ;
 457     ppuint b{ base_() } ;
 458 
 459     for (int i = static_cast<unsigned int>( digit_.size()) - 1 ;  i >= 0 ;  --i)
 460     {
 461         // Check for overflow before it happens:  Will result * base + digit > (maximum unsigned integer)?
 462         if (result > (numeric_limits<ppuint>::max() - digit_[ i ]) / b )
 463         {
 464             ostringstream os ;
 465             os << "BigInt::operator ppuint() is about to overflow "
 466                << " partial result [" << result << "] * base [" << b << "] + digit [" << digit_[ i ]
 467                << "] > numeric limit for ppuint [" << numeric_limits<ppuint>::max() << "]" ;
 468             throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
 469         }
 470         else
 471             result = result * b + digit_[ i ] ;
 472     }
 473 
 474     return result ;
 475 }
 476 
 477 
 478 
 479 /*=============================================================================
 480 |
 481 | NAME
 482 |
 483 |     Bigint::toString
 484 |
 485 | DESCRIPTION
 486 |
 487 |    Print n-place number
 488 |
 489 |    (u    ... u  )
 490 |      n-1      0  b
 491 |
 492 |    as a decimal string,
 493 |
 494 |    (U    ... U  )
 495 |      m-1      0  10
 496 |
 497 | EXAMPLE
 498 |
 499 |     // Set up multiprecision calculations up to the size p^n.
 500 |     Bigint n = "31415926535897932" ;
 501 |     string s = n.toString() ;
 502 |
 503 +============================================================================*/
 504 
 505 string BigInt::toString() const
 506 {
 507     const ppuint decimal_base{ static_cast<ppuint>( 10u ) } ;
 508 
 509     // Output stream to hold the digits.
 510     ostringstream os ;
 511     string s = "" ;
 512 
 513     // Pull out the decimal digits in reverse:
 514     //
 515     //    do until u == 0:
 516     //        U = u mod 10
 517     //        u = |  u / 10 |
 518     //            --       --
 519     BigInt u( *this ) ;
 520 
 521     // Special cases.
 522     if (u == BigInt( static_cast<ppuint>( 0u ) ))
 523         return "0" ;
 524     else if (u == BigInt( static_cast<ppuint>( 1u ) ))
 525         return "1" ;
 526 
 527     #ifdef DEBUG_PP_BIGINT
 528     cout << "toString digits = "  ;
 529     #endif
 530 
 531     while (u != BigInt( static_cast<ppuint>( 0u ) ))
 532     {
 533         ppuint digit{ u % decimal_base } ;
 534 
 535         #ifdef DEBUG_PP_BIGINT
 536         // This line was recursing in and out of BigInt::toString and BigInt::printNumber, ending in a memory violation.
 537         // cout << "toString:  number to convert u = " ; printNumber( u, cout ) ;
 538         cout << digit << " "  ;
 539         #endif
 540 
 541         if (!(os << digit))
 542         {
 543             ostringstream os ;
 544             os << "BigInt::toString can't convert digit = " << digit ;
 545             throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
 546         }
 547 
 548         u /= decimal_base ;
 549     }
 550 
 551     #ifdef DEBUG_PP_BIGINT
 552     cout << endl ;
 553     #endif
 554 
 555     s = os.str() ;
 556     reverse( s.begin(), s.end() ) ;
 557 
 558     return s ;
 559 }
 560 
 561 
 562 
 563 /*=============================================================================
 564 |
 565 | NAME
 566 |
 567 |     Operator << for BigInt
 568 |
 569 | DESCRIPTION
 570 |
 571 |     Print a BigInt to the output stream.
 572 |
 573 | EXAMPLE
 574 |
 575 |     try
 576 |     {
 577 |         BigInt x( "123" ) ;
 578 |         cout << x << endl ;
 579 |     }
 580 |     catch( BigIntRangeError & e )
 581 |     {
 582 |         cout << "Error in outputting BigInt x" << endl ;
 583 |     }
 584 |
 585 +============================================================================*/
 586 
 587 ostream & operator<<( ostream & out, const BigInt & u )
 588 {
 589     // Convert to a string first, then output as a string.
 590     // Can throw an exception, which we'll catch at higher level.
 591     out << u.toString() ;
 592 
 593     return out ;
 594 }
 595 
 596 
 597 
 598 /*=============================================================================
 599 |
 600 | NAME
 601 |
 602 |     Operator >> for BigInt
 603 |
 604 | DESCRIPTION
 605 |
 606 |     BigInt stream input.
 607 |
 608 | EXAMPLE
 609 |
 610 |     try
 611 |     {
 612 |         BigInt x ;  // Declare a BigInt object.
 613 |         cin >> x ;  // Read a long decimal string from the standard input.
 614 |                     // e.g. type in 31415926535897932
 615 |     }
 616 |     catch( BigIntRangeError & e )
 617 |     {
 618 |         cout << "Error in inputting BigInt x = 31415926535897932" << endl ;
 619 |     }
 620 |
 621 +============================================================================*/
 622 
 623 istream & operator>>( istream & in, BigInt & u )
 624 {
 625     // Input the number as a string first.
 626     string s ;
 627     in >> s ;
 628 
 629     // Convert to BigInt and copy into argument.
 630     // We may throw a BigIntRangeError at this point if s is corrupted.
 631     u = BigInt( s ) ;
 632 
 633     return in ;
 634 }
 635 
 636 
 637 
 638 /*=============================================================================
 639 |
 640 | NAME
 641 |
 642 |    +
 643 |
 644 | DESCRIPTION
 645 |
 646 |    Add n-place numbers,
 647 |
 648 |    (u    ... u  )  +
 649 |      n-1      0  b
 650 |
 651 |    (v    ... v  )  =
 652 |      n-1      0  b
 653 |
 654 |    (w    ... w  )  +
 655 |      n-1      0  b
 656 |
 657 |
 658 | EXAMPLE
 659 |
 660 |    try
 661 |    {
 662 |        BigInt u, v, w ;
 663 |        u = v + w ;
 664 |    }
 665 |    catch( BigIntMathError & e )
 666 |    {
 667 |        cerr << e.what() << endl ;
 668 |    }
 669 |
 670 +============================================================================*/
 671 
 672 const BigInt operator+( const BigInt & u, const BigInt & v )
 673 {
 674     // Do + in terms of += to maintain consistency.
 675     // Copy construct temporary copy and add to it.
 676     // Return value optimization compiles away the copy constructor.
 677     // const on return type disallows doing (u+v) = w ;
 678     return BigInt( u ) += v ;
 679 }
 680 
 681 
 682 
 683 /*=============================================================================
 684 |
 685 | NAME
 686 |
 687 |     +
 688 |
 689 | DESCRIPTION
 690 |
 691 |    Add n-place number and a digit,
 692 |
 693 |    (u    ... u  )  + d
 694 |      n-1      0  b
 695 |
 696 | EXAMPLE
 697 |
 698 |    try
 699 |    {
 700 |        BigInt u, w ;
 701 |        ppuint d ;
 702 |        u = v + d ;
 703 |    }
 704 |    catch( BigIntMathError & e )
 705 |    {
 706 |        cerr << e.what() << endl ;
 707 |    }
 708 |
 709 +============================================================================*/
 710 
 711 const BigInt operator+( const BigInt & u, ppuint d )
 712 {
 713     // Do + in terms of += to maintain consistency.
 714     // Return value optimization compiles away the copy constructor.
 715     // const on return type disallows doing (u+v) = w ;
 716     return BigInt( u ) += d ;
 717 }
 718 
 719 
 720 
 721 /*=============================================================================
 722 |
 723 | NAME
 724 |
 725 |     +=
 726 |
 727 | DESCRIPTION
 728 |
 729 |    Add n-place number and an m-place number,
 730 |
 731 |    (u    ... u  )  +=
 732 |      n-1      0  b
 733 |
 734 |    (0 ... v    ... v  ) 
 735 |            m-1      0  b
 736 |
 737 | EXAMPLE
 738 |
 739 |    try
 740 |    {
 741 |        BigInt u, v ;
 742 |        u += v ;
 743 |    }
 744 |    catch( BigIntMathError & e )
 745 |    {
 746 |        cerr << e.what() << endl ;
 747 |    }
 748 |
 749 +============================================================================*/
 750 
 751 BigInt & BigInt::operator+=( const BigInt & v )
 752 {
 753     // Pull out the number of digits and base.
 754     int  n{ static_cast<int>( digit_.size() ) } ;
 755     int  m{ static_cast<int>( v.digit_.size() ) } ;
 756     ppuint b{ base_() } ;
 757 
 758 #ifdef DEBUG_PP_BIGINT
 759     printNumber( *this, cout ) ;
 760     printNumber( v,  cout ) ;
 761 #endif
 762 
 763     // Allocate temporary space for the sum.
 764     BigInt w ;
 765 
 766     ppuint carry{ 0 } ;  // Always < 2b
 767     ppuint sum{ 0 };  // Always in [0, 2b)
 768 
 769     //  Add and carry, starting with the least significant
 770     //  digits in each number.
 771     for (int i = 0 ;  i < (m < n ? m : n) ;  ++i)
 772     {
 773         sum         = digit_[ i ] + v.digit_[ i ] + carry ;
 774         // Can throw from argument operator=().
 775         w.digit_.push_back( sum % b ) ;
 776         carry       = sum / b ;
 777     }
 778 
 779     if (n < m)
 780         for (int i = n ;  i < m ;  ++i)
 781         {
 782             sum         = 0 + v.digit_[ i ] + carry ;
 783             // Can throw from argument operator=().
 784             w.digit_.push_back( sum % b ) ;
 785             carry       = sum / b ;
 786         }
 787     else if (m < n)
 788         for (int i = m ;  i < n ;  ++i)
 789         {
 790             sum         = digit_[ i ] + 0 + carry ;
 791             // Can throw from argument operator=().
 792             w.digit_.push_back( sum % b ) ;
 793             carry       = sum / b ;
 794         }
 795 
 796 
 797     // Add the carry digit as the most significant digit.
 798     if (carry != 0)
 799         // Can throw from argument operator=().
 800         w.digit_.push_back( carry ) ;
 801 
 802     // Swap the digits of w into this number.
 803     swap( w.digit_, digit_ ) ;
 804 
 805     // Return current object now containing the sum.
 806     return *this ;
 807 
 808     // As we go out of scope, w's destructor will be called.
 809 }
 810 
 811 
 812 
 813 /*=============================================================================
 814 |
 815 | NAME
 816 |
 817 |     +=
 818 |
 819 | DESCRIPTION
 820 |
 821 |    Add an unsigned integer d to an n-place number u
 822 |
 823 |    (u    ... u  )   + d
 824 |      n-1      0  b
 825 |
 826 | EXAMPLE
 827 |
 828 |    try
 829 |    {
 830 |        BigInt u ;
 831 |        ppuint d ;
 832 |        u += d ;
 833 |    }
 834 |    catch( BigIntMathError & e )
 835 |    {
 836 |        cerr << e.what() << endl ;
 837 |    }
 838 |
 839 +============================================================================*/
 840 
 841 BigInt & BigInt::operator+=( const ppuint d )
 842 {
 843     // Pull out the number of digits and base.
 844     int  n = static_cast<int>( digit_.size() ) ;
 845     ppuint b{ base_() } ;
 846 
 847     if (d >= b)
 848     {
 849         // TODO:   We should do BigInt + here instead of aborting.
 850         ostringstream os ;
 851         os << "BigInt::operator+=  Adding BigInt = " << *this << " + d = " << d << " but d >= base = " << b
 852            << " which is bigger than the largest value a digit of BigInt can have.  Please email me to request a feature enhancement" ;
 853         throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
 854     }
 855 
 856     // Allocate temporary space for the sum.
 857     BigInt w ;
 858 
 859     // Add the first digit.
 860     ppuint sum{ digit_[ 0 ] + d } ; // Always in [0, 2b)
 861     w.digit_.push_back( sum % b ) ;
 862 
 863     ppuint carry{ sum / b } ; // Always < 2b
 864 
 865     //  Add and carry in the other digits if needed.
 866     for (int i = 1 ;  i < n ;  ++i)
 867     {
 868         sum         = digit_[ i ] + carry ;
 869         // Can throw from argument operator=().
 870         w.digit_.push_back( sum % b ) ;
 871         carry       = sum / b ;
 872     }
 873 
 874     // Add the carry digit as the most significant digit.
 875     if (carry != 0)
 876         // Can throw from argument operator=().
 877         w.digit_.push_back( carry ) ;
 878 
 879     // Swap the digits of w into this number.
 880     swap( w.digit_, digit_ ) ;
 881 
 882     // Return current object now containing the sum.
 883     return *this ;
 884 
 885     // As we go out of scope, w's destructor will be called.
 886 }
 887 
 888 
 889 
 890 /*=============================================================================
 891 |
 892 | NAME
 893 |
 894 |     Prefix ++
 895 |
 896 | DESCRIPTION
 897 |
 898 |     Prefix increment operator.
 899 |
 900 | EXAMPLE
 901 |
 902 |     BigInt i ;
 903 |     ++i ;
 904 |     i.operator++() ; // Equivalent
 905 |
 906 +============================================================================*/
 907 
 908 BigInt & BigInt::operator++()
 909 {
 910     // Add 1.
 911     (*this) += 1 ;
 912 
 913     // Return reference to incremented BigInt object so we
 914     // can modify the object, e.g.  ++i++ ;
 915     return *this ;
 916 }
 917 
 918 
 919 
 920 /*=============================================================================
 921 |
 922 | NAME
 923 |
 924 |     Postfix ++
 925 |
 926 | DESCRIPTION
 927 |
 928 |     Postfix
 929 |
 930 | EXAMPLE
 931 |
 932 |     BigInt i ;
 933 |     ++i ;
 934 |     i.operator++(0) ; // Equivalent
 935 |
 936 +============================================================================*/
 937 
 938 const BigInt BigInt::operator++( int ) // Dummy argument to distinguish postfix form.
 939 {
 940     BigInt save = *this ;
 941 
 942     // Add 1 using prefix operator.
 943     ++(*this) ;
 944 
 945     // Return const copy of original
 946     // to prevent doing i++++ which is the same as i.operator++(0).operator++(0)
 947     // which does i += 1, returns i, then does i += 1 which would increment
 948     // by 1;  not what we wanted.
 949     return save ;
 950 }
 951 
 952 
 953 
 954 /*=============================================================================
 955 |
 956 | NAME
 957 |
 958 |     -
 959 |
 960 | DESCRIPTION
 961 |
 962 |    BigInt - BigInt
 963 |
 964 | EXAMPLE
 965 |
 966 |     BigInt u = 3 ;
 967 |     BigInt v = 4 ;
 968 |     BigInt w = u - v ;
 969 |
 970 +============================================================================*/
 971 
 972 const BigInt operator-( const BigInt & u, const BigInt & v )
 973 {
 974     // Do - in terms of -= to maintain consistency.
 975     // Return value optimization compiles away the copy constructor.
 976     // const on return type disallows doing (u-v) = w ;
 977     return BigInt( u ) -= v ;
 978 }
 979 
 980 
 981 
 982 /*=============================================================================
 983 |
 984 | NAME
 985 |
 986 |     -
 987 |
 988 | DESCRIPTION
 989 |
 990 |     BigInt - integer
 991 |
 992 | EXAMPLE
 993 |
 994 |     BigInt u{ 3 } ;
 995 |     ppuint v{ 4 } ;
 996 |     BigInt w = u - v ;
 997 |
 998 +============================================================================*/
 999 
1000 const BigInt operator-( const BigInt & u, const ppuint d )
1001 {
1002     // Do - in terms of -= to maintain consistency.
1003     // Return value optimization compiles away the copy constructor.
1004     // const on return type disallows doing (u-v) = w ;
1005     return BigInt( u ) -= d ;
1006 }
1007 
1008 
1009 
1010 /*=============================================================================
1011 |
1012 | NAME
1013 |
1014 |     -=
1015 |
1016 | DESCRIPTION
1017 |
1018 |     u -= v
1019 |
1020 |    Subtract u - v = (u1 ... un) - (v1 ... vm) = (w1 ... wn) = w.
1021 |    Assume u >= v.  If u < v, the carry is negative, and we abort.
1022 |
1023 +============================================================================*/
1024 
1025 BigInt & BigInt::operator-=( const BigInt & v )
1026 {
1027     // Pull out the number of digits and the base.
1028     int  n = static_cast<int>( digit_.size()) ;
1029     int  m = static_cast<int>( v.digit_.size()) ;
1030     ppuint b{ base_() } ;
1031 
1032     if (m > n)
1033     {
1034         ostringstream os ;
1035         os << "BigInt::operator-= " << " negative result for u - v = " << *this << " - " << v ;
1036         throw BigIntUnderflow( os.str(), __FILE__, __LINE__ ) ;
1037     }
1038 
1039     // Allocate temporary space for the difference.
1040     BigInt w ;
1041 
1042     // Subtract u - v starting with nth digits assuming n >= m.
1043     ppuint borrow{ 0 } ;
1044     for (int i = 0 ;  i < n ;  ++i)
1045     {
1046         // t in [-b, b)
1047         ppsint t = digit_[ i ] - (i >= m ? 0 : v.digit_[ i ]) + borrow ;
1048 
1049         // Subtract, allowing for where u < v and we must borrow.
1050         // Can throw from argument operator=().
1051         w.digit_.push_back( (t > 0) ? t % b : (t + b) % b ) ;
1052 
1053         borrow = (t < 0) ? -1 : 0 ; // -1 if u - v + borrow < 0, else 0
1054     }
1055 
1056     if (borrow == -1)
1057     {
1058         ostringstream os ;
1059         os << "BigInt::operator-= " << " negative result for u - v = " << *this << " - " << v ;
1060         throw BigIntUnderflow( os.str(), __FILE__, __LINE__ ) ;
1061     }
1062 
1063     // Trim leading zeros, if any, but stop if q = 0.
1064     while (w.digit_.size() > 1 && w.digit_.back() == 0)
1065         w.digit_.pop_back() ;
1066 
1067     // Swap the digits of w into this number.
1068     swap( w.digit_, digit_ ) ;
1069 
1070     // Return (reference to) the sum.
1071     return *this ;
1072 
1073     // As we go out of scope, w's destructor will be called.
1074 }
1075 
1076 
1077 
1078 /*=============================================================================
1079 |
1080 | NAME
1081 |
1082 |     -=
1083 |
1084 | DESCRIPTION
1085 |
1086 +============================================================================*/
1087 
1088 BigInt & BigInt::operator-=( const ppuint u )
1089 {
1090     // Get the base and number of digits.
1091     ppuint b{ base_() } ;
1092     int  n = static_cast<int>( digit_.size()) ;
1093 
1094     if (u >= b)
1095     {
1096         // TODO:   We should do BigInt - here instead of aborting.
1097         ostringstream os ;
1098         os << "BigInt::operator-=( u ) BigInt = " << *this << " - u = " << u << " >= base = " << b
1099            << " please email me to request a feature enhancement" ;
1100         throw BigIntUnderflow( os.str(), __FILE__, __LINE__ ) ;
1101     }
1102 
1103     // Subtract 1 from the least significant digit.
1104     ppsint t = digit_[ 0 ] - u ;
1105 
1106     // Subtract and allow for borrow.
1107     digit_[ 0 ] = (t > 0) ? t % b : (t + b) % b ;
1108     ppsint borrow  = (t < 0) ? -1 : 0 ; // -1 if u - 1 + borrow < 0, else 0
1109 
1110     // Propagate the subtraction up to the (n-1)st digit.
1111     for (int i = 1 ;  i < n ;  ++i)
1112     {
1113         t = digit_[ i ] + borrow ;
1114 
1115         // Subtract and allow for borrow.
1116         digit_[ i ] = (t > 0) ? t % b : (t + b) % b ;
1117         borrow      = (t < 0) ? -1 : 0 ; // -1 if u - 1 + borrow < 0, else 0
1118     }
1119 
1120     if (borrow == -1)
1121     {
1122         ostringstream os ;
1123         os << "BigInt::operator-= " << " underflow, borrow = -1" ;
1124         throw BigIntUnderflow( os.str(), __FILE__, __LINE__ ) ;
1125     }
1126 
1127     // Trim leading zeros, if any, but stop if n = 0.
1128     while (digit_.size() > 1 && digit_.back() == 0)
1129         digit_.pop_back() ;
1130 
1131     // Return (reference to) the sum.
1132     return *this ;
1133 }
1134 
1135 
1136 
1137 /*=============================================================================
1138 |
1139 | NAME
1140 |
1141 |     --
1142 |
1143 | DESCRIPTION
1144 |
1145 |     Subtract u - 1.
1146 |
1147 |     ( u  u     . . . u )
1148 |        n  n-1         0
1149 |    -                1
1150 |    --------------------
1151 |     ( w  w     . . . w )
1152 |        n  n-1         0
1153 |
1154 |    Subtract u - 1 = (u1 ... un) - (0 ... 1) = (w1 ... wn) = w.
1155 |    Assume u >= 1.  If u < 0, the carry is negative, and we abort.
1156 |
1157 +============================================================================*/
1158 
1159 BigInt & BigInt::operator--()
1160 {
1161     ppuint b{ base_() } ;
1162 
1163     // Subtract 1 from the least significant digit.
1164     ppsint t = digit_[ 0 ] - 1 ;
1165 
1166     // Subtract and allow for borrow.
1167     digit_[ 0 ] = (t > 0) ? t % b : (t + b) % b ;
1168     ppsint borrow        = (t < 0) ? -1 : 0 ; // -1 if u - 1 + borrow < 0, else 0
1169 
1170     // Subtract u - 1 starting with nth digits.
1171     for (unsigned int i = 1 ;  i < digit_.size() ;  ++i)
1172     {
1173         t = digit_[ i ] + borrow ;
1174 
1175         // Subtract and allow for borrow.
1176         digit_[ i ] = (t > 0) ? t % b : (t + b) % b ;
1177         borrow        = (t < 0) ? -1 : 0 ; // -1 if u - 1 + borrow < 0, else 0
1178     }
1179 
1180     if (borrow == -1)
1181     {
1182         ostringstream os ;
1183         os << "BigInt::operator-- " << " underflow, borrow = -1" ;
1184         throw BigIntUnderflow( os.str(), __FILE__, __LINE__ ) ;
1185     }
1186 
1187     return *this ;
1188 }
1189 
1190 
1191 
1192 /*=============================================================================
1193 |
1194 | NAME
1195 |
1196 |     --
1197 |
1198 | DESCRIPTION
1199 |
1200 +============================================================================*/
1201 
1202 BigInt BigInt::operator--( int )
1203 {
1204     BigInt save = *this ;
1205 
1206     // Subtract 1.
1207     --(*this) ;
1208 
1209     return save ;
1210 }
1211 
1212 
1213 
1214 /*=============================================================================
1215 |
1216 | NAME
1217 |
1218 |     *
1219 |
1220 | DESCRIPTION
1221 |
1222 +============================================================================*/
1223 
1224 //
1225 //  Multiply (u1 ... un) * digit = (w1 ... wn).
1226 //  Numbers are right justified, so that un and wn are in array
1227 //  location n.
1228 const BigInt operator*( const BigInt & u, const BigInt & v )
1229 {
1230     // Do * in terms of *= to maintain consistency.
1231     // Return value optimization compiles away the copy constructor.
1232     // const on return type disallows doing (u*v) = w ;
1233     return BigInt( u ) *= v ;
1234 }
1235 
1236 
1237 
1238 /*=============================================================================
1239 |
1240 | NAME
1241 |
1242 |     *
1243 |
1244 | DESCRIPTION
1245 |
1246 +============================================================================*/
1247 
1248 //
1249 //  Multiply (u1 ... un) * digit = (w1 ... wn).
1250 //  Numbers are right justified, so that un and wn are in array
1251 //  location n.
1252 const BigInt operator*( const BigInt & u, const ppuint digit )
1253 {
1254     // Do * in terms of *= to maintain consistency.
1255     // Return value optimization compiles away the copy constructor.
1256     // const on return type disallows doing (u*v) = w ;
1257     return BigInt( u ) *= digit ;
1258 }
1259 
1260 
1261 
1262 /*=============================================================================
1263 |
1264 | NAME
1265 |
1266 |     *=
1267 |
1268 | DESCRIPTION
1269 |
1270 |    Multiply
1271 |
1272 |    (u    ...  u  )
1273 |      m-1       0  b
1274 |
1275 |      x
1276 |
1277 |    (v    ...  v  )
1278 |      n-1       0  b
1279 |
1280 |     =
1281 |
1282 |    (u    ...  u  )
1283 |      m+n-1     0  b
1284 |
1285 +============================================================================*/
1286 
1287 BigInt & BigInt::operator*=( const BigInt & v )
1288 {
1289     // Pull out the number of digits and base.
1290     int m = static_cast<int>( digit_.size() ) ;
1291     int n = static_cast<int>( v.digit_.size() ) ;
1292 
1293     ppuint b{ base_() } ;
1294 
1295     if (m  == 0 || n == 0)
1296     {
1297         ostringstream os ;
1298         os << "BigInt::operator*= Internal error.  Num digits in u = " << m << " and v = " << n ;
1299         throw BigIntMathError( os.str(), __FILE__, __LINE__ ) ;
1300     }
1301 
1302     // Allocate temporary space for the product.
1303     BigInt w ;
1304     try
1305     {
1306         w.digit_.resize( m + n ) ;
1307     }
1308     catch( length_error & e )
1309     {
1310         ostringstream os ;
1311         os << "BigInt::operator*=() had a length_error exception during w.digit_resize(m+n = " << m+n << ")" ;
1312         throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
1313     }
1314 
1315     ppuint carry{ 0 } ;
1316 
1317     // Zero out leading digits of product.
1318     for (int j = m - 1 ;  j >= 0 ;  --j)
1319         w.digit_[ j ] = 0 ;
1320 
1321     // Multiply u by each digit of v from right to left.
1322     for (int j = 0 ;  j < n ;  ++j)
1323     {
1324         // Skip if digit of v is zero.
1325         if (v.digit_[ j ] == 0)
1326         {
1327             w.digit_[ j + m ] = 0 ;
1328             continue ;
1329         }
1330 
1331         // Multiply u by the jth digit of v.
1332         carry = 0 ;
1333         for (int i = 0 ;  i < m ;  ++i)
1334         {
1335             ppuint t{ digit_[ i ] * v.digit_[ j ] + w.digit_[ i + j ] + carry } ;
1336             w.digit_[ i + j ] = t % b ;
1337             carry             = t / b ;
1338         }
1339 
1340         w.digit_[ j + m ] = carry ;
1341     }
1342 
1343     // Trim a leading zero digit.
1344     if (w.digit_[ m + n - 1 ] == 0)
1345     {
1346         w.digit_.pop_back() ;
1347     }
1348 
1349     // Swap the digits of w into this number.
1350     swap( w.digit_, digit_ ) ;
1351 
1352     // Trim off leading zero digits.
1353 
1354     // Return (reference to) the sum.
1355     return *this ;
1356 
1357     // As we go out of scope, w's destructor will be called.
1358 }
1359 
1360 
1361 /*=============================================================================
1362 |
1363 | NAME
1364 |
1365 |     *=
1366 |
1367 | DESCRIPTION
1368 |
1369 |
1370 |  Multiply u * d =
1371 |
1372 |    (u    ...  u  )
1373 |      m-1       0  b
1374 |
1375 |      x
1376 |
1377 |      d
1378 |
1379 +============================================================================*/
1380 
1381 BigInt & BigInt::operator*=( ppuint d )
1382 {
1383     // Pull out the number of digits.
1384     int  n{ static_cast<int>( digit_.size() ) } ;
1385     ppuint b{ base_() } ;
1386 
1387     // Allocate temporary space for the product.
1388     BigInt w ;
1389 
1390     if (d > b)
1391     {
1392         // TODO:   We should do BigInt * here instead of aborting.
1393         //      w *= static_cast<BigInt>( d ) ;
1394         ostringstream os ;
1395         os << "BigInt::operator*=   digit d = " << d << " > base b = " << b
1396            << " please email me to request a feature enhancement" ;
1397         throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
1398     }
1399     // In this special case, we just shift digits left and zero fill.
1400     // But do nothing if the number is zero.
1401     else if (d == b && *this != static_cast<BigInt>( static_cast<ppuint>( 0u )))
1402     {
1403         w.digit_.push_back( 0 ) ;
1404 
1405         for (int i = 0 ;  i < n ;  ++i)
1406             w.digit_.push_back( digit_[ i ] ) ;
1407     }
1408     else
1409     {
1410         ppuint carry = 0 ;
1411 
1412         // Multiply digits and carry.
1413         for (int i = 0 ;  i < n ;  ++i)
1414         {
1415             ppuint t = digit_[ i ] * d + carry ;
1416             ppuint r = t % b ;
1417 
1418             w.digit_.push_back( r ) ;
1419             carry = t / b ;
1420 
1421             #ifdef DEBUG_PP_BIGINT
1422             cout << "BigInt::operator*=(ppuint)" << endl ;
1423             cout << "d        = " << d << endl ;
1424             cout << "digit[i] = " << digit_[ i ] << " i = " << i << endl ;
1425             cout << "t        = " << t << endl ;
1426             cout << "r        = " << r << endl ;
1427             cout << "carry    = " << carry << endl ;
1428             #endif
1429         }
1430 
1431         // Additional carry beyond the nth digit.
1432         if (carry)
1433             w.digit_.push_back( carry ) ;
1434     }
1435 
1436     // Swap the digits of w into this number.
1437     swap( w.digit_, digit_ ) ;
1438 
1439     // Return (reference to) the product.
1440     return *this ;
1441 
1442     // As we go out of scope, w's destructor will be called.
1443 }
1444 
1445 
1446 
1447 /*=============================================================================
1448 |
1449 | NAME
1450 |
1451 |     /
1452 |
1453 | DESCRIPTION
1454 |
1455 |    Integer division
1456 |
1457 | EXAMPLE
1458 |
1459 |  u / v =
1460 |
1461 +============================================================================*/
1462 
1463 const BigInt operator/( const BigInt & u, const BigInt & v )
1464 {
1465     // Do / in terms of /= to maintain consistency.
1466     // Return value optimization compiles away the copy constructor.
1467     // const on return type disallows doing (u/v) = w ;
1468     return BigInt( u ) /= v ;
1469 }
1470 
1471 
1472 
1473 BigInt & BigInt::operator/=( const BigInt & v )
1474 {
1475     BigInt r ; // Thrown away.
1476     BigInt q ;
1477 
1478     divMod( *this, v, q, r ) ;
1479 
1480     // Swap the digits of q and current object.
1481     swap( q.digit_, digit_ ) ;
1482 
1483     return *this ;
1484 }
1485 
1486 
1487 
1488 /*=============================================================================
1489 |
1490 | NAME
1491 |
1492 |     /
1493 |
1494 | DESCRIPTION
1495 |
1496 |     Divide by a digit.
1497 |
1498 +============================================================================*/
1499 
1500 const BigInt operator/( const BigInt & u, ppuint d )
1501 {
1502     // Do / in terms of /= to maintain consistency.
1503     // Return value optimization compiles away the copy constructor.
1504     // const on return type disallows doing (u/v) = w ;
1505     return BigInt( u ) /= d ;
1506 }
1507 
1508 /*=============================================================================
1509  |
1510  | NAME
1511  |
1512  |     /=
1513  |
1514  | DESCRIPTION
1515  |
1516  |     Divide by a digit.
1517  |
1518  +============================================================================*/
1519 
1520 BigInt & BigInt::operator/=( ppuint d )
1521 {
1522 
1523     BigInt q ;
1524     ppuint   r ; // Thrown away.
1525 
1526     divMod( *this, d, q, r ) ;
1527 
1528     // Swap the digits of q and current object.
1529     swap( q.digit_, digit_ ) ;
1530 
1531     return *this ;
1532 }
1533 
1534 
1535 
1536 /*=============================================================================
1537 |
1538 | NAME
1539 |
1540 |     %
1541 |
1542 | DESCRIPTION
1543 |
1544 |     u mod v  for BigInts
1545 |
1546 +============================================================================*/
1547 
1548 BigInt operator%( const BigInt & u, const BigInt & v )
1549 {
1550     BigInt q ; // Thrown away.
1551     BigInt r ;
1552 
1553     divMod( u, v, q, r ) ;
1554 
1555     // Create a copy of r upon return and destruct r.
1556     return r ;
1557 }
1558 
1559 
1560 
1561 /*=============================================================================
1562 |
1563 | NAME
1564 |
1565 |     %
1566 |
1567 | DESCRIPTION
1568 |
1569 |     u mod d  for BigInt u and integer d
1570 |
1571 +============================================================================*/
1572 
1573 ppuint operator%( const BigInt & u, const ppuint d )
1574 {
1575     BigInt q ; // Thrown away.
1576     ppuint r ;
1577 
1578     divMod( u, d, q, r ) ;
1579 
1580     return r ;
1581 }
1582 
1583 
1584 
1585 /*=============================================================================
1586 |
1587 | NAME
1588 |
1589 |     divMod
1590 |
1591 | DESCRIPTION
1592 |
1593 |    Compute the quotient and remainder for
1594 |    u / d where u is a BigInt and d is an integer.
1595 |
1596 |        (u    ... u  )
1597 |          n-1      0  b
1598 |    q = ------------------ = ( w    ... w  )
1599 |                d               n-1      0  b
1600 |
1601 |    r = (u    ... u  )   mod d
1602 |          n-1      0  b
1603 |
1604 +============================================================================*/
1605 
1606 void divMod( const BigInt & u,
1607              const ppuint   d,
1608              BigInt &       q,
1609              ppuint &       r )
1610 {
1611     // Get the number of digits and the base.
1612     ppuint b{ u.base_() } ;
1613     int  n{ static_cast<int>( u.digit_.size() ) } ;
1614 
1615     // Don't allow zero divide.
1616     if (d == static_cast<ppuint>( 0u ))
1617     {
1618         ostringstream os ;
1619         os << "BigInt::divMod  Divide by 0" ;
1620         throw BigIntZeroDivide( os.str(), __FILE__, __LINE__ ) ;
1621     }
1622 
1623     // u has no digits.
1624     if (n <= 0)
1625     {
1626         ostringstream os ;
1627         os << "BigInt::divMod  u/d = qd+r u has zero digits" ;
1628         throw BigIntMathError( os.str(), __FILE__, __LINE__ ) ;
1629     }
1630 
1631     // If q is not empty, clear it.
1632     if (q.digit_.size() > 0)
1633         q.digit_.clear() ;
1634 
1635     // Call multiprecision divide.
1636     if (d > b)
1637     {
1638         // TODO:   We should do BigInt divMod here instead of aborting.
1639         //    BigInt rr ;
1640         //    divMod( u, static_cast<BigInt>( d ), q, rr ) ;
1641         //    r = static_cast<ppuint>( rr ) ;
1642         ostringstream os ;
1643         os << "BigInt::divMod " << " d = " << d << " > base b = " << b
1644            << " please email me to request a feature enhancement" ;
1645         throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
1646     }
1647     // In this special case, we just shift digits right.
1648     else if (d == b)
1649     {
1650         r = 0 ;
1651 
1652         // u is zero;  leave it alone.
1653         if (u == static_cast<BigInt>( 0u ))
1654             ;
1655         // The single digit of u is shifted out.  Only zero remains.
1656         else if (n == 1)
1657         {
1658             r = u.digit_[ 0 ] ;
1659             q.digit_.push_back( 0 ) ;
1660         }
1661         // Remainder r is the zeroth digit.  Shift all other digits into the quotient.
1662         else
1663         {
1664             r = u.digit_[ 0 ] ;
1665 
1666             for (int j = 1 ; j <= n-1 ;  ++j)
1667                 q.digit_.push_back( u.digit_[ j ] ) ;
1668 
1669             #ifdef DEBUG_PP_BIGINT
1670             cout << "BigInt::divMod: the special case of divisor d = " << d << " equals " << " base b = " << b << endl ;
1671             #endif
1672         }
1673     }
1674     else
1675     {
1676         r = 0 ;
1677 
1678         // Long division from most to least significant digit.
1679         for (int j = n - 1 ; j >= 0 ;  --j)
1680         {
1681             ppuint t         = r * b + u.digit_[ j ] ;
1682             q.digit_.push_back( t / d ) ;
1683             r             = t % d ;
1684         }
1685 
1686         // Put the digits into correct order.
1687         reverse( q.digit_.begin(), q.digit_.end() ) ;
1688 
1689         // Trim leading zeros, if any, but stop if q = 0.
1690         while (q.digit_.size() > 1 && q.digit_.back() == 0)
1691             q.digit_.pop_back() ;
1692     }
1693 
1694     return ;
1695 }
1696 
1697 
1698 
1699 /*=============================================================================
1700 |
1701 | NAME
1702 |
1703 |     divMod
1704 |
1705 | DESCRIPTION
1706 |
1707 |     Compute u / v = q, r where all quantities are BigInts.
1708 |
1709 |    u = (u      ... u  u  )
1710 |          m+n-1      1  0  b
1711 |
1712 |    v = (v    ... v  v  )
1713 |          n-1      1  0  b
1714 |
1715 |    v   != 0,   n > 1
1716 |     n-1
1717 |
1718 |    q =  | u / v | = (q  q... q  )
1719 |         --     --     m  m-1  0  b
1720 |
1721 |    r = u mod v = (r   ... r  )
1722 |                    n-1     0  b
1723 |
1724 |    u = q v + r
1725 |
1726 |    Throw BigIntZeroDivide if v = 0.
1727 |
1728 | EXAMPLE
1729 |
1730 |   b = 10
1731 |
1732 |   n = 4 
1733 |   v = (v    v    v   v ) = 3457
1734 |         n-1  n-2  1   0
1735 |
1736 |   m+n = 6 
1737 |   u = (u      u       u  u  u  u ) = 398765
1738 |         m+n-1  m+n-2   3  2  1  0
1739 |
1740 |   m = 2 
1741 |
1742 |           v3 v2 v1 v0
1743 |                        --------------------
1744 |                        ) u5 u4 u3 u2 u1 u0
1745 |       v = 3  4  5  7   )  3  9  8  7  6  5 = u
1746 |
1747 |
1748 |   Now copy v2 = v and u2 = u
1749 |
1750 |                            |       b       |   |   10      |
1751 |   Normalizing factor = d = |   ----------  | = |   -----   | = 2
1752 |                            --  vn-1 + 1   --   --  3 + 1  --
1753 |
1754 |   v2 = v2 * d
1755 |   u2 = u2 * d
1756 |
1757 |   But let's use u and v instead of u2 and v2 for simplicity of notation.
1758 |   The normalized division is then.  Note extra 0 appended at the left of u.
1759 |
1760 |           v3 v2 v1 v0                      1   1   5
1761 |                         ----------------------------
1762 |                         ) u6  u5  u4  u3  u2  u1  u0
1763 |           6   9  1  4   )  0   7   9   7   5   3   0
1764 |          vn-1            um+n
1765 |                         -  0   6   9   1   4
1766 |                            -----------------
1767 |                            0   1   0   6   1   3
1768 |                            -   0   6   9   1   4
1769 |                                -----------------
1770 |                                0   3   6   9   9   0
1771 |                             -  0   3   4   5   7   0
1772 |                                ---------------------
1773 |                                    0   2   4   2   0
1774 |
1775 |   Recall n = 4  m = 2   m+n = 6 
1776 |
1777 |                                        *** Correction needed? ***
1778 |
1779 |        j      q           r            q >= b?    q >= b r + u  ?
1780 |                                                               j+n-2
1781 |
1782 |        2   | 0 7 |      0 10 + 7       no         no
1783 |            | --- |= 1   mod 6 = 1
1784 |            -  6 -
1785 |
1786 |    subtraction (uj+n...uj) = 01061
1787 |
1788 |    trial quotient q2 = 1 r2 = 4
1789 |    corrected trial quotient q2 = 1 r2 = 4
1790 |    final trial quotient q2 = 1 r2 = 4
1791 |    i = 0 j+i = 1 q2 = 1 borrow = 0
1792 |    v2[ i ] = 4 u2[j+i] = 3 t2 = -1
1793 |    corrected borrow = -1u2.[j+i] = 9
1794 |    i = 1 j+i = 2 q2 = 1 borrow = -1
1795 |    v2[ i ] = 1 u2[j+i] = 1 t2 = -1
1796 |    corrected borrow = -1u2.[j+i] = 9
1797 |    i = 2 j+i = 3 q2 = 1 borrow = -1
1798 |    v2[ i ] = 9 u2[j+i] = 6 t2 = -4
1799 |    corrected borrow = -1u2.[j+i] = 6
1800 |    i = 3 j+i = 4 q2 = 1 borrow = -1
1801 |    v2[ i ] = 6 u2[j+i] = 0 t2 = -7
1802 |    corrected borrow = -1u2.[j+i] = 3
1803 |    i = 4 j+i = 5 q2 = 1 borrow = -1
1804 |    v2[ i ] = 0 u2[j+i] = 1 t2 = 0
1805 |    corrected borrow = 0u2.[j+i] = 0
1806 |    j = 1
1807 |    After subtraction (uj+n...uj) = 03699
1808 |    trial quotient q2 = 6 r2 = 0
1809 |    corrected trial quotient q2 = 5 r2 = 6
1810 |    final trial quotient q2 = 5 r2 = 6
1811 |    i = 0 j+i = 0 q2 = 5 borrow = 0
1812 |    v2[ i ] = 4 u2[j+i] = 0 t2 = -20
1813 |    corrected borrow = -2u2.[j+i] = 0
1814 |    i = 1 j+i = 1 q2 = 5 borrow = -2
1815 |    v2[ i ] = 1 u2[j+i] = 9 t2 = 2
1816 |    corrected borrow = 0u2.[j+i] = 2
1817 |    i = 2 j+i = 2 q2 = 5 borrow = 0
1818 |    v2[ i ] = 9 u2[j+i] = 9 t2 = -36
1819 |    corrected borrow = -4u2.[j+i] = 4
1820 |    i = 3 j+i = 3 q2 = 5 borrow = -4
1821 |    v2[ i ] = 6 u2[j+i] = 6 t2 = -28
1822 |    corrected borrow = -3u2.[j+i] = 2
1823 |    i = 4 j+i = 4 q2 = 5 borrow = -3
1824 |    v2[ i ] = 0 u2[j+i] = 3 t2 = 0
1825 |    corrected borrow = 0u2.[j+i] = 0
1826 |    j = 0
1827 |    After subtraction (uj+n...uj) = 02420
1828 |    m = 2n = 4
1829 |    Quotient q = 115
1830 |    Normalized remainder r = 1210
1831 |
1832 +============================================================================*/
1833 
1834 void divMod( const BigInt & u,
1835              const BigInt & v,
1836              BigInt &       q,
1837              BigInt &       r )
1838 {
1839     ppuint b{ u.base_() } ;
1840     int  m_n{ static_cast<int>( u.digit_.size()) } ;
1841     int    n{ static_cast<int>( v.digit_.size()) } ;
1842     int    m{ m_n - n } ; //  Compute m from the fact that u has m + n digits, and v has n digits.
1843 
1844     #ifdef DEBUG_PP_BIGINT
1845     cout << "\t------------divMod( " << u << " , " << v << " BigInt )---------" << endl ;
1846     cout << "\tm = " << m_n << endl ;
1847     cout << "\tn = " << n << endl ;
1848     cout << "\tu = " ; printNumber( u, cout ) ;
1849     cout << "\tv = " ; printNumber( v, cout ) ;
1850     #endif
1851 
1852     if (n < 1 || m_n < 1)
1853     {
1854         ostringstream os ;
1855         os << "BigInt::divMod() zero digits for either "
1856            << "num digits of u = " << m_n << " or num digits of v = " << n ;
1857         throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
1858     }
1859 
1860     // One digit division.
1861     if (n == 1)
1862     {
1863         // Make r have at least one digit.
1864         try
1865         {
1866             r.digit_.resize( 1 ) ;
1867         }
1868         catch( length_error & e )
1869         {
1870             ostringstream os ;
1871             os << "BigInt::divMod() had a length_error exception during r.digit_.resize( 1 )" ;
1872             throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
1873         }
1874 
1875         divMod( u, v.digit_[0], q, r.digit_[0] ) ;
1876         return ;
1877     }
1878 
1879     #ifdef DEBUG_PP_BIGINT
1880     cout << "\tm = " << m << endl ;
1881     #endif
1882     // If u < v, return q = 0, r = u.
1883     if (m < 0)
1884      {
1885         q = BigInt( 0u ) ;
1886         r = u ;
1887         return ;
1888     }
1889 
1890     // Allocate temporary space for normalized copies of u and v.  
1891     // Destructors for these temporaries will be called when we go out 
1892     // of scope at the end of this function.
1893     BigInt u2( u ) ;
1894     BigInt v2( v ) ;
1895 
1896     // Add an extra zero digit to the beginning of u:  u[ m+n ] = 0.
1897     u2.digit_.push_back( 0 ) ;
1898 
1899     // Add an extra zero digit to the beginning of v:  v[ n ] = 0.
1900     v2.digit_.push_back( 0 ) ;
1901 
1902     #ifdef DEBUG_PP_BIGINT
1903     cout << "\tu2 = " ; printNumber( u2, cout ) ;
1904     cout << "\tv2 = " ; printNumber( v2, cout ) ;
1905     #endif
1906 
1907 
1908     // Make certain the number of digits in q and r are sufficient.
1909     try
1910     {
1911         q.digit_.resize( m+1 ) ;
1912         r.digit_.resize( n ) ;
1913     }
1914     catch( length_error & e )
1915     {
1916         ostringstream os ;
1917         os << "BigInt::divMod() had a length_error exception during q.digit_resize( m = " << m+1 << " )"
1918            << "or r.digit_.resize( n = " << n << " )" ;
1919         throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
1920     }
1921 
1922     //  Assume v    != 0, and n > 1.  If not, skip leading zero digits.
1923     //          n-1
1924     //
1925     int k = static_cast<int>( v2.digit_.size() ) - 1 ;
1926     while( v2.digit_[ k ] == 0 && k >= 0)
1927         --k ;
1928 
1929     if (k < 0)
1930     {
1931         ostringstream os ;
1932         os << "BigInt::divMod " << " zero divide, k = " << k ;
1933         throw BigIntZeroDivide( os.str(), __FILE__, __LINE__ ) ;
1934     }
1935 
1936     if (n <= 1 || v2.digit_[ n-1 ] == 0)
1937     {
1938         ostringstream os ;
1939         os << "BigInt::divMod  Zero divide n = " << n << "v2.digit_[ " << n-1 << " ]" << v2.digit_[ n-1 ] ;
1940         throw BigIntZeroDivide( os.str(), __FILE__, __LINE__ ) ;
1941     }
1942 
1943     //  Normalizer.
1944     ppuint d = b / (v2.digit_[ n-1 ] + 1) ;
1945 
1946     #ifdef DEBUG_PP_BIGINT
1947     cout << "\tnormalizer d = " << d << endl ;
1948     #endif
1949 
1950     //  Skip normalizing step if d = 1.
1951     if (d == 1)
1952     {
1953         u2.digit_[ m+n ] = 0 ;
1954     }
1955     else
1956     {
1957         // Normalize (u     ... u   ) = (u      ... u )  d
1958         //             m+n-1     0        m+n-1      0
1959 
1960         ppuint carry = 0 ;
1961         for (int j = 0 ;  j <= m + n - 1 ;  ++j)
1962         {
1963             ppuint t{ u2.digit_[ j ] * d + carry } ;
1964 
1965             u2.digit_[ j ] = t % b ;
1966             carry          = t / b ;
1967         }
1968         u2.digit_[ m+n ] = carry ;
1969 
1970         // Normalize (v   ... v  ) = (v   ... v  )  d
1971         //             n-1     0       n-1     0
1972         //
1973         //
1974         carry = 0 ;
1975         for (int j = 0 ;  j <= n-1 ;  ++j)
1976         {
1977             ppuint t = v2.digit_[ j ] * d + carry ;
1978 
1979             v2.digit_[ j ] = t % b ;
1980             carry          = t / b ;
1981         }
1982 
1983         // No carry can occur, so flag an error if it happens.
1984         if (carry != 0)
1985         {
1986             ostringstream os ;
1987             os << "BigInt::divMod " << " carry = " << carry ;
1988             throw BigIntMathError( os.str(), __FILE__, __LINE__ ) ;
1989         }
1990     }
1991 
1992     #ifdef DEBUG_PP_BIGINT
1993     cout << "\tnormalized u2 = " ; printNumber( u2, cout ) ;
1994     cout << "\tnormalized v2 = " ; printNumber( v2, cout ) ;
1995     #endif
1996 
1997 
1998     //  Find the quotient digit,
1999     //
2000     //         |  (u    ... u  )   |
2001     //         |    j+n      j  b  |
2002     //   q  =  | ----------------- |
2003     //    j    |  (v   ... v  )    |
2004     //         |    n-1     0  b   |
2005     //         --                 --
2006     for (int j = m ;  j >= 0 ;  --j)
2007     {
2008         // First determine a trial quotient digit
2009         //
2010         //       |  (u    u     )    |
2011         //   ^   |    j+n  j+n-1  b  |
2012         //   q = | ----------------- |
2013         //       |      v            |
2014         //       |       n-1         |
2015         //       --                 --
2016         //
2017         ppuint temp = u2.digit_[ j+n ] * b + u2.digit_[ j+n-1 ] ;
2018         ppuint q2   = temp / v2.digit_[ n-1 ] ;
2019         ppuint r2   = temp % v2.digit_[ n-1 ] ;
2020 
2021         #ifdef DEBUG_PP_BIGINT
2022         cout << "\ttrial quotient q2 = " << q2 << " r2 = " << r2 << endl ;
2023         #endif
2024 
2025         // Correction if necessary.
2026         if (q2 >= b || (q2 * v2.digit_[ n-2 ] > b * r2 + u2.digit_[ j+n-2 ]))
2027         {
2028             --q2 ;
2029             r2 += v2.digit_[ n-1 ] ;
2030         }
2031 
2032         #ifdef DEBUG_PP_BIGINT
2033         cout << "\tcorrected trial quotient q2 = " << q2 << " r2 = " << r2 << endl ;
2034         #endif
2035 
2036         // Low probability repeat correction if necessary.
2037         if (r2 < b)
2038         {
2039             if (q2 >= b || (q2 * v2.digit_[ n-2 ] > b * r2 + u2.digit_[ j+n-2 ]))
2040             {
2041                 --q2 ;
2042                 // We don't use the remainder since this is the last correction.
2043                 ///r2 += v2.digit_[ n-1 ] ;
2044             }
2045         }
2046 
2047         #ifdef DEBUG_PP_BIGINT
2048         cout << "\trepeat corrected trial quotient q2 = " << q2 << " r2 = " << r2 << endl ;
2049         #endif
2050 
2051         #ifdef DEBUG_PP_BIGINT
2052         cout << "\tfinal trial quotient q2 = " << q2 << " r2 = " << r2 << endl ;
2053         #endif
2054 
2055         // Multiply and subtract:
2056         //
2057         //         (u    u      ... u  )
2058         //           j+n  j+n-1      j
2059         //
2060         //  - q2 * (0    v    ...   v  )
2061         //                n-1        0
2062         //
2063         //  =
2064         //         (u    u      ... u  )
2065         //           j+n  j+n-1      j
2066 
2067         ppsint borrow{ 0 } ;
2068         for (int i = 0 ;  i <= n ;  ++i)
2069         {
2070             ppsint t2  = borrow + u2.digit_[ j + i ] - q2 * v2.digit_[ i ] ;
2071 
2072             #ifdef DEBUG_PP_BIGINT
2073             cout << "\t\ti = " << i << " j+i = " << j+i << " q2 = " << q2 << " borrow = " << borrow << endl ;
2074             cout << "\t\tv2[ i ] = " << getDigit( v2, i ) << " u2[j+i] = " << getDigit( u2, j+i ) << " t2 = " << t2 << endl ;
2075             #endif
2076 
2077             // Digit subtraction is positive:  don't need to borrow.
2078             if (t2 >= 0)
2079             {
2080                 u2.digit_[ j + i ] = t2 ;
2081                 borrow = 0 ;
2082             }
2083             //  If t2 is negative, we need to borrow.
2084             else
2085             {
2086                     // e.g. t = -1 -> -10,
2087                     // borrow = floor( (t+1)/b ) - 1 = 0/10 - 1 -> -9/10 - 1 => -1;
2088                     //  -11 -> -20:  borrow -2
2089                     borrow = -1 + (t2 + 1) / (ppsint) b ;
2090 
2091                     // Add the positive borrow to the digit to force it positive.
2092                     u2.digit_[ j + i ] = -borrow * (ppsint) b + t2 ;
2093             }
2094             #ifdef DEBUG_PP_BIGINT
2095             cout << "\t\tcorrected borrow = " << borrow << "u2.[j+i] = "<< getDigit( u2,j+i ) << endl ;
2096             #endif
2097         }
2098 
2099         #ifdef DEBUG_PP_BIGINT
2100         cout << "\tj = " << j << endl ;
2101         cout << "\tAfter subtraction (uj+n...uj) = " ;
2102         for (int k = j+n ;  k >= j ;  --k)
2103            cout << getDigit( u2, k ) ;
2104         cout << endl ;
2105         #endif
2106 
2107         // Save the quotient.
2108         q.digit_[ j ] = q2 ;
2109 
2110 
2111         // Decrease q2 and add back correction if q2 is too big.
2112         // Probability of this step given random digits is about 2 / b.
2113         //
2114         //      (u    u      ... u )
2115         //        n+j  n+j-1      j
2116         //
2117         //  +   (0    v    ...   v )
2118         //             n-1        0
2119         //  =
2120         //      (u    u      ... u )
2121         //        n+j  n+j-1      j
2122         //
2123         // Ignore the carry to the left of u    since it cancels with the earlier borrow.
2124         //                                  n+j
2125         if (borrow < 0)
2126         {
2127             --q.digit_[ j ] ;
2128 
2129             ppuint carry = 0 ;
2130             for (int i = 0 ;  i <= n ;  ++i)
2131             {
2132                ppsint t = u2.digit_[ j + i ] + v2.digit_[ i ] + carry ;
2133 
2134                 u2.digit_[ j + i ] = t % b ;
2135                 carry              = t / b ;
2136             }
2137         }
2138     } // end for j
2139 
2140 
2141     // Quotient is (q    ... q  )
2142     //               m        0  b
2143 
2144     // Trim leading zeros, if any, but stop if q = 0.
2145     while (q.digit_.size() > 1 && q.digit_.back() == 0)
2146         q.digit_.pop_back() ;
2147 
2148     #ifdef DEBUG_PP_BIGINT
2149     cout << "\tm = " << m << " n = " << n << endl ;
2150     cout << "\tQuotient q = " ; printNumber( q, cout ) ;
2151     #endif
2152 
2153 
2154     // Remainder:  (r    ... r  )   = (u    ... u  )  / d
2155     //               n-1      0  b      n-1      0  b
2156     if (d != 1)
2157     {
2158         ppuint remainder = 0 ;
2159         for (int j = n-1 ;  j >= 0 ;  --j)
2160         {
2161             ppuint t{ u2.digit_[ j ] + remainder * b } ;
2162 
2163             r.digit_[ j ] = t / d ;
2164             remainder     = t % d ;
2165 
2166             #ifdef DEBUG_PP_BIGINT
2167             cout << "\tRemainder normalization:  j = " << j << " t = " << t << endl ;
2168             cout << "\tr[ j ] = " << r.digit_[ j ] << " remainder = " << remainder << endl ;
2169             #endif
2170         }
2171     }
2172     // No need to normalize r, just copy digits of u.
2173     else
2174        for (int j = n-1 ;  j >= 0 ;  --j)
2175             r.digit_[ j ] = u2.digit_[ j ] ;
2176 
2177 
2178     #ifdef DEBUG_PP_BIGINT
2179     cout << "\tRemainder r = " ; printNumber( r, cout ) ;
2180     #endif
2181 
2182     // Trim leading zeros, if any, but stop if r = 0.
2183     while (r.digit_.size() > 1 && r.digit_.back() == 0)
2184         r.digit_.pop_back() ;
2185 
2186     #ifdef DEBUG_PP_BIGINT
2187     cout << "\tNormalized remainder r = " ;  printNumber( r, cout ) ;
2188     #endif
2189 
2190     return ;
2191 }
2192 
2193 
2194 
2195 /*=============================================================================
2196 |
2197 | NAME
2198 |
2199 |    power
2200 |
2201 | DESCRIPTION
2202 |
2203 | METHOD
2204 |
2205 |     Speed up by doing multiplication by repeated squaring.
2206 |     Expand exponent n in binary.  Discard the leading 1 bit.
2207 |     Thereafter, square for every 0 bit;  square and multiply
2208 |     by p for every 1 bit.
2209 |
2210 |    100 = 0  . . . 0  1  1  0  0  1  0  0 (binary representation)
2211 |          |<- ignore ->| S  S SX SX SX  S (exponentiation rule,
2212 |          S = square, X = multiply by x)
2213 |
2214 +============================================================================*/
2215 
2216 BigInt power( ppuint p, ppuint n )
2217 {
2218     // Special cases first.
2219     if (p == static_cast<ppuint>( 0u ))
2220     {
2221         //   0
2222         //  0  = undefined
2223         if (n == static_cast<ppuint>( 0u ))
2224         {
2225             ostringstream os ;
2226             os << "BigInt::power  Attempt to do 0 ^ 0 in pow" ;
2227             throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
2228         }
2229         //   nonzero
2230         //  0        = 0
2231         else
2232             return BigInt( static_cast<ppuint>( 0u ) ) ;
2233     }
2234    else if (n == static_cast<ppuint>( 0u ))
2235    {
2236        //        0
2237        // nonzero  = 1
2238        return BigInt( static_cast<ppuint>( 1u ) ) ;
2239    }
2240 
2241     // Find the number of the leading bit.
2242     int bitNum = sizeof( n ) * 8u - 1u ; // Number of highest possible bit integer n.
2243 
2244     #ifdef DEBUG_PP_BIGINT
2245     cout << "\tn = " << n << " p = " << p << endl ;
2246     cout << "initial max bitNum = " << bitNum << endl ;
2247     #endif // NDEBUG
2248 
2249     while (!testBit( n, bitNum ))
2250         --bitNum ;
2251 
2252     #ifdef DEBUG_PP_BIGINT
2253     cout << "after skipping leading 0 bits, bitNum = " << bitNum << endl ;
2254     #endif // NDEBUG
2255 
2256     if (bitNum == -1)
2257     {
2258         ostringstream os ;
2259         os << "BigInt::power bitNum == -1 internal error in BigInt" ;
2260         throw BigIntMathError( os.str(), __FILE__, __LINE__ ) ;
2261     }
2262 
2263     // Exponentiation by repeated squaring.  Skip over the leading 1 bit.
2264     // Thereafter, square for every 0 bit;  square and multiply by p for
2265     // every 1 bit.
2266 
2267     BigInt w( p ) ;
2268     while( --bitNum >= 0 )
2269     {
2270         // Square.
2271         w *= w ;
2272 
2273         // Times p.
2274         if (testBit( n, bitNum ))
2275             w *= p ;
2276 
2277             #ifdef DEBUG_PP_BIGINT
2278             cout << "bitNum = " << bitNum << endl ;
2279             cout << "testBit( n ) = " << testBit(n,bitNum) << endl ;
2280             cout << "intermediate w = " << w << endl ;
2281             #endif // NDEBUG
2282     }
2283 
2284     // Return a copy of power.
2285     // As we go out of scope, power's destructor will be called.
2286     return w ;
2287 }
2288 
2289 
2290 
2291 /*=============================================================================
2292 |
2293 | NAME
2294 |
2295 |     ceilLg
2296 |
2297 | DESCRIPTION
2298 |
2299 |     --            --
2300 |     |  log  ( n )  |
2301 |            2
2302 |
2303 | EXAMPLE
2304 |
2305 |      BigInt n = 6 ;
2306 |      n.ceilLg() => 3
2307 |
2308 | NOTES
2309 |
2310 |      Let n be written in binary, m bits long with a leading 1 bit:
2311 |
2312 |      n = 1 *  ...  *
2313 |
2314 |                 m       m
2315 |      Then n <= 2 - 1 < 2
2316 | 
2317 |      Since log  = lg is monotonically increasing,
2318 |               2
2319 |
2320 |      lg( n ) < m.
2321 |      We define ceilLg( n ) = m
2322 |
2323 |      n = 6 = 110
2324 |      There are m = 3 bits counting from leading 1 bit to lsb.
2325 |      lg( n ) = 2.58496... = < 3 = ceilLg( n )
2326 |
2327 +============================================================================*/
2328 
2329 int BigInt::ceilLg()
2330 {
2331     int bitNum ;
2332 
2333     for (bitNum = maxBitNumber();  testBit( bitNum ) == false && bitNum >= 0 ;  --bitNum)
2334     #ifdef DEBUG_PP_BIGINT
2335     cout << "bitNum = " << bitNum << " testBit = " << testBit( bitNum ) << endl
2336     #endif // NDEBUG
2337      ;
2338 
2339     return bitNum + 1 ;
2340 }
2341 
2342 
2343 
2344 /*=============================================================================
2345 |
2346 | NAME
2347 |
2348 |     ==
2349 |
2350 | DESCRIPTION
2351 |
2352 +============================================================================*/
2353 
2354 bool operator==( const BigInt & u, const BigInt & v )
2355 {
2356     //  Get the number of digits.
2357     int  m     = static_cast<int>( u.digit_.size()) ;
2358     int  n     = static_cast<int>( v.digit_.size()) ;
2359     int max_n = (n > m) ? n : m ;
2360 
2361     // Different number of non-zero digits.
2362     if (m != n)
2363         return false ;
2364 
2365     for (int i = 0 ;  i < max_n ;   ++i)
2366     {
2367         // Two digits are not equal.
2368         if (u.digit_[ i ] != v.digit_[ i ])
2369             return false ;
2370     }
2371 
2372     // All digits equal.
2373     return true ;
2374 }
2375 
2376 
2377 
2378 /*=============================================================================
2379 |
2380 | NAME
2381 |
2382 |     ==
2383 |
2384 | DESCRIPTION
2385 |
2386 +============================================================================*/
2387 
2388 bool operator==( const BigInt & u, const ppuint d )
2389 {
2390     ppuint b{ u.base_() } ;
2391 
2392     // TODO:   We should do BigInt == here instead of aborting.
2393     if (d > b)
2394     {
2395         ostringstream os ;
2396         os << "BigInt::operator== " << " digit d = " << d << " > base b = " << u.base_()
2397            << " please email me to request a feature enhancement" ;
2398         throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
2399     }
2400     // Special case check to see if u = 10 in base b.
2401     else if (d == b)
2402     {
2403         if (u.digit_.size() != 2 || u.digit_[ 0 ] != 0 || u.digit_[ 1 ] != 1)
2404             return false ;
2405     }
2406     // More than one base b digit or no digits.
2407     else if (u.digit_.size() != 1)
2408         return false ;
2409     // One digit which equals d.
2410     else if (u.digit_[ 0 ] == d)
2411         return true ;
2412 
2413     return false ;
2414 }
2415 
2416 
2417 /*=============================================================================
2418 |
2419 | NAME
2420 |
2421 |     !=
2422 |
2423 | DESCRIPTION
2424 |
2425 +============================================================================*/
2426 
2427 bool operator!=( const BigInt & u, const BigInt & v )
2428 {
2429 return !(u == v) ;
2430 }
2431 
2432 
2433 
2434 /*=============================================================================
2435 |
2436 | NAME
2437 |
2438 |     !=
2439 |
2440 | DESCRIPTION
2441 |
2442 +============================================================================*/
2443 
2444 bool operator!=( const BigInt & u, const ppuint d )
2445 {
2446     return !(u == d) ;
2447 }
2448 
2449 
2450 
2451 /*=============================================================================
2452 |
2453 | NAME
2454 |
2455 |    >
2456 |
2457 | DESCRIPTION
2458 |
2459 +============================================================================*/
2460 
2461 bool operator>( const BigInt & u, const BigInt & v )
2462 {
2463     //  Get the number of digits.
2464     int m = static_cast<int>( u.digit_.size()) ;
2465     int n = static_cast<int>( v.digit_.size()) ;
2466 
2467     // u has more non-zero digits than v.
2468     if (m > n)
2469         return true ;
2470 
2471     // u has fewer nonzero digits than v.
2472     else if (m < n)
2473        return false ;
2474 
2475     // Same number of digits.  Is a leading digit of u bigger?
2476     for (int i = n-1 ;  i >= 0 ;   --i)
2477     {
2478         if (u.digit_[ i ] > v.digit_[ i ])
2479             return true ;
2480         else if (u.digit_[ i ] < v.digit_[ i ])
2481             return false ;
2482     }
2483 
2484     // All digits equal.
2485     return false ;
2486 }
2487 
2488 
2489 
2490 /*=============================================================================
2491 |
2492 | NAME
2493 |
2494 |    >
2495 |
2496 | DESCRIPTION
2497 |
2498 +============================================================================*/
2499 
2500 bool operator>( const BigInt & u, const ppuint d )
2501 {
2502     ppuint b{ u.base_() } ;
2503 
2504     if (d >= b)
2505     {
2506         ostringstream os ;
2507         os << "BigInt::operator> " << "d = " << d << " > base = " << b ;
2508         throw BigIntOverflow( os.str(), __FILE__, __LINE__ ) ;
2509     }
2510 
2511     // More digits:  definitely greater.
2512     if (u.digit_.size() != 1)
2513         return true ;
2514 
2515     // One digit and greater.
2516     if (u.digit_[ 0 ] > d)
2517         return true ;
2518 
2519     return false ;
2520 }
2521 
2522 
2523 
2524 /*=============================================================================
2525 |
2526 | NAME
2527 |
2528 |     <
2529 |
2530 | DESCRIPTION
2531 |
2532 +============================================================================*/
2533 
2534 bool operator<( const BigInt & u, const BigInt & v )
2535 {
2536     return !(u > v || u == v) ;
2537 }
2538 
2539 
2540 
2541 /*=============================================================================
2542 |
2543 | NAME
2544 |
2545 |    <
2546 |
2547 | DESCRIPTION
2548 |
2549 +============================================================================*/
2550 
2551 bool operator<( const BigInt & u, const ppuint d )
2552 {
2553     return !(u > d || u == d) ;
2554 }
2555 
2556 
2557 
2558 /*=============================================================================
2559 |
2560 | NAME
2561 |
2562 | DESCRIPTION
2563 |
2564 +============================================================================*/
2565 
2566 bool operator>=( const BigInt & u, const BigInt & v )
2567 {
2568     return (u > v) || u == v ;
2569 }
2570 
2571 
2572 
2573 /*=============================================================================
2574 |
2575 | NAME
2576 |
2577 |     >=
2578 |
2579 | DESCRIPTION
2580 |
2581 +============================================================================*/
2582 
2583 bool operator>=( const BigInt & u, ppuint d )
2584 {
2585     return (u > d) || u == d ;
2586 }
2587 
2588 
2589 
2590 /*=============================================================================
2591 |
2592 | NAME
2593 |
2594 | DESCRIPTION
2595 |
2596 +============================================================================*/
2597 
2598 bool operator<=( const BigInt & u, const BigInt & v )
2599 {
2600     return u < v || u == v ;
2601 }
2602 
2603 
2604 
2605 /*=============================================================================
2606 |
2607 | NAME
2608 |
2609 |     <=
2610 |
2611 | DESCRIPTION
2612 |
2613 +============================================================================*/
2614 
2615 bool operator<=( const BigInt & u, const ppuint d )
2616 {
2617     return u < d || u == d ;
2618 }
2619 
2620 
2621 
2622 /*=============================================================================
2623 |
2624 | NAME
2625 |
2626 | DESCRIPTION
2627 |
2628 +============================================================================*/
2629 
2630 BigInt operator&( const BigInt & u, const BigInt & n )
2631 {
2632     // Allocate temporary space for the result which will
2633     // be destructed as this function goes out of scope.
2634     BigInt w ;
2635 
2636     // Return a copy of the result.
2637     return w ;
2638 }
2639 
2640 
2641 
2642 /*=============================================================================
2643 |
2644 | NAME
2645 |
2646 |     Left shift operator << for BigInt
2647 |
2648 | DESCRIPTION
2649 |
2650 |     Bit shift left.
2651 |
2652 +============================================================================*/
2653 
2654 BigInt operator<<( const BigInt & u, ppuint n )
2655 {
2656     // Allocate temporary space for the result which will
2657     // be destructed as this function goes out of scope.
2658     BigInt w ;
2659 
2660     // Return a copy of the result.
2661     return w ;
2662 }
2663 
2664 
2665 
2666 /*=============================================================================
2667 |
2668 | NAME
2669 |
2670 |     testBit
2671 |
2672 | DESCRIPTION
2673 |
2674 |     Bit test for BigInt precision integers.
2675 |
2676 +============================================================================*/
2677 
2678 const bool BigInt::testBit( const int bitNum ) const
2679 {
2680     int n = static_cast<int>( digit_.size()) ;
2681 
2682     // Compute the digit in which the bit lies.
2683     //
2684     // e.g. if there are 7 bits per digit, bitNum = 8 lies in digit = 1.
2685     // and is subbit number 1.
2686     //
2687     //      87 6543210    Test the 8th bit.
2688     // 0000010 0000000    bitNum    = 8
2689     //       1       0    digitNum  = 1
2690     // 6543210 6543210    subBitNum = 8 - 1 * 7 = 1
2691     int digitNum = bitNum / numBitsPerDigit_() ;
2692 
2693     if (digitNum > n - 1)
2694     {
2695         ostringstream os ;
2696         os << "BigInt::testBit( " << bitNum << " ) is out of range" << endl
2697            << "The number has " << numBitsPerDigit_() * n << " bits" << endl ;
2698                 throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
2699     }
2700 
2701     int subBitNum = bitNum - digitNum * numBitsPerDigit_() ;
2702 
2703     return (digit_[ digitNum ] & (static_cast<ppuint>( 1u ) << subBitNum)) ? true : false ;
2704 }
2705 
2706 //  Bit test for low precision integers.
2707 const bool testBit( const ppuint n, const int bitNum )
2708 {
2709     return (n & (static_cast<ppuint>( 1u ) << bitNum)) ? true : false ;
2710 }
2711 
2712 
2713 
2714 /*=============================================================================
2715 |
2716 | NAME
2717 |
2718 |     maxBitNumber
2719 |
2720 | DESCRIPTION
2721 |
2722 |     Highest bit number in a BigInt number.
2723 |
2724 +============================================================================*/
2725 
2726 int BigInt::maxBitNumber() const
2727 {
2728     // Highest bit number = number of digits * bits per digit - 1
2729     return numBitsPerDigit_() * static_cast<int>( digit_.size()) - 1 ;
2730 }
2731 
2732 /*=============================================================================
2733 |
2734 | NAME
2735 |
2736 |     getBase
2737 |
2738 | DESCRIPTION
2739 |
2740 |     Return the BigInt base.
2741 |
2742 +============================================================================*/
2743 
2744 const ppuint BigInt::getBase()
2745 {
2746     return base_() ;
2747 }
2748 
2749 
2750 /*=============================================================================
2751 |
2752 | NAME
2753 |
2754 |     getDigit
2755 |
2756 | DESCRIPTION
2757 |
2758 |     Return the nth digit of the BigInt.  For testing only.
2759 |
2760 +============================================================================*/
2761 
2762 const ppuint getDigit( const BigInt & u, const int n )
2763 {
2764     if (n >= 0 && n < (int)u.digit_.size())
2765         return u.digit_[ n ] ;
2766     else
2767     {
2768         ostringstream os ;
2769         os << "BigInt::getDigit( " << n << " ) is out of range"
2770            << "The number has " << u.digit_.size() << " digits" ;
2771                 throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
2772     }
2773 }
2774 
2775 
2776 /*=============================================================================
2777 |
2778 | NAME
2779 |
2780 |     getNumDigits
2781 |
2782 | DESCRIPTION
2783 |
2784 |     Return the number BigInt digits n.
2785 |
2786 +============================================================================*/
2787 
2788 const int getNumDigits(const BigInt & u )
2789 {
2790     return static_cast<int>( u.digit_.size() ) ;
2791 }
2792 
2793 
2794 /*=============================================================================
2795 |
2796 | NAME
2797 |
2798 |     setBase
2799 |
2800 | DESCRIPTION
2801 |
2802 |     Forcibly reset the base for all BigInt numbers.  Used for testing only.
2803 |
2804 +============================================================================*/
2805 
2806 void setBase( const BigInt & u, const ppuint base )
2807 {
2808    *u.pbase = base ;
2809 }
2810 
2811 
2812 /*=============================================================================
2813 |
2814 | NAME
2815 |
2816 |     printNumber
2817 |
2818 | DESCRIPTION
2819 |
2820 |     Print out a BigInt quantity to an output stream.  If the BigInt is not
2821 |     defined or initialized, print out "undefined".
2822 |
2823 | EXAMPLE
2824 |
2825 |     For standard output to the console,
2826 |         BigInt w( "68719476735", cout ) ;
2827 |        printNumber(r)
2828 |     gives
2829 |        68719476735 [digits = 31 2147483647 base b = 2147483648 number of digits = 2]
2830 |
2831 +============================================================================*/
2832 
2833 void printNumber( const BigInt & u, ostream & out = cout )
2834 {
2835     // If the BigInt was not initialized. Don't throw an exception since this is
2836     // only a debug print routine.
2837     if (u.base_() == 0 || u.digit_.size() == 0)
2838     {
2839         out << "undefined" << endl ;
2840         return ;
2841     }
2842 
2843     // The number itself in decimal.
2844     out << u.toString() ;
2845 
2846     // Base b digits.
2847     out << "    [ " << getNumDigits( u ) << " digit" << (getNumDigits(u) > 1 ? "s" : "") << " = " ;
2848 
2849     for (int i = getNumDigits( u ) - 1 ;  i >= 0 ; --i)
2850         out << getDigit( u, i ) << " " ;
2851 
2852     out << " base b = " << BigInt::getBase() << " ]" << endl ;
2853 }
2854 
2855 
2856 /*=============================================================================
2857 |
2858 | NAME
2859 |
2860 |     printNumber
2861 |
2862 | DESCRIPTION
2863 |
2864 |     Print out a BigInt quantity to the standard output stream.
2865 |     You can call it in the lldb debugger from the command line.
2866 |
2867 | EXAMPLE
2868 |
2869 |    (lldb) expr printNumber(r)
2870 |    68719476735 [digits = 31 2147483647 base b = 2147483648 number of digits = 2]
2871 |
2872 +============================================================================*/
2873 
2874 void printNumber( const BigInt & u )
2875 {
2876     printNumber( u, cout ) ;
2877 }
2878 
2879 /*=============================================================================
2880 |
2881 | NAME
2882 |
2883 |     Bigint::base_
2884 |
2885 | DESCRIPTION
2886 |
2887 |     Initialize the base of the multiprecision arithmetic system.
2888 |     Throw BigIntRangeError.
2889 |
2890 | EXAMPLE
2891 |
2892 |     BigInt::useBase()
2893 |     {
2894 |         ppuint b{ base_() } ;
2895 |         ...
2896 |     }
2897 |
2898 +============================================================================*/
2899 
2900 ppuint & BigInt::base_()
2901 {
2902     // Base of the number system used for each digit.  If a digit has can hold N bits,
2903     // the signed integer maximum is
2904     //          N-1
2905     //         2    - 1
2906     //
2907     // If we let the base be
2908     //
2909     //          N/2 - 1
2910     //     b = 2
2911     //
2912     //       2    N-2
2913     // then b  = 2    will not only fit into a signed integer but will also be a power
2914     //
2915     // of 2, making bit shifting and testing much easier.
2916     //
2917     // e.g. for N = 8, the largest positive number = 01111111 = 127, base b = 8 whose
2918     // square, 64 is less than the signed maximum.
2919     //
2920     static ppuint base = static_cast<ppuint>(1u) << numBitsPerDigit_() ;
2921 
2922     // Point to it for testing.  Use reference?
2923     pbase = &base ;
2924 
2925     // Detect if the base is way out of range, and build up an exception.
2926     if (base <= 0)
2927     {
2928         ostringstream os ;
2929         os << "BigInt::base_() base  = " << base << " <= 0 " ;
2930         throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
2931     }
2932 
2933     return base ;
2934 }
2935 
2936 
2937 
2938 /*=============================================================================
2939 |
2940 | NAME
2941 |
2942 |     Bigint::numBitsPerDigit_
2943 |
2944 | DESCRIPTION
2945 |
2946 |     Initialize the base of the multiprecision arithmetic system.
2947 |     Throw BigIntRangeError.
2948 |
2949 | EXAMPLE
2950 |
2951 |     BigInt::useBase()
2952 |     {
2953 |         ppuint b = base_() ;
2954 |         ...
2955 |     }
2956 |
2957 +============================================================================*/
2958 
2959 int & BigInt::numBitsPerDigit_()
2960 {
2961     // Base of the number system used for each digit.  If a digit has can hold N bits,
2962     // we let
2963     //          N/2 - 1
2964     //     b = 2
2965     //
2966     //          2
2967     // so that b  fits into a digit and b is an integer number of bits in
2968     // length, making bit shifting and testing easier.
2969     //
2970     // Use half the number of bits in an unsigned integer.
2971     static int num_bits_per_digit = (sizeof( ppuint ) * 8) / 2 - 1 ;
2972 
2973     #ifdef DEBUG_PP_BIGINT
2974     cout << "numBitsPerDigit_():" << endl ;
2975     cout << "sizeof( ppuint ) = " << sizeof( ppuint ) << " bytes" << endl ;
2976     cout << "num_bits_per_digit = " << num_bits_per_digit << endl ;
2977     cout << "base = 1 << num_bits_per_digit = " << (static_cast<ppuint>(1u) << num_bits_per_digit) << endl ;
2978     #endif
2979 
2980     if (num_bits_per_digit <= 0)
2981     {
2982         ostringstream os ;
2983         os << "BigInt::numBitsPerDigit()_  num_bits_per_digit  = " << num_bits_per_digit << " <= 0" ;
2984         throw BigIntRangeError( os.str(), __FILE__, __LINE__ ) ;
2985     }
2986 
2987     return num_bits_per_digit ;
2988 }