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