1/*==============================================================================
2|
3| NAME
4|
5| ppArith.cpp
6|
7| DESCRIPTION
8|
9| Miscellaneous integer multiple precision math routines.
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| LEGAL
15|
16| Primpoly Version 16.3 - A Program for Computing Primitive Polynomials.
17| Copyright (C) 1999-2024 by Sean Erik O'Connor. All Rights Reserved.
18|
19| This program is free software: you can redistribute it and/or modify
20| it under the terms of the GNU General Public License as published by
21| the Free Software Foundation, either version 3 of the License, or
22| (at your option) any later version.
23|
24| This program is distributed in the hope that it will be useful,
25| but WITHOUT ANY WARRANTY; without even the implied warranty of
26| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27| GNU General Public License for more details.
28|
29| You should have received a copy of the GNU General Public License
30| along with this program. If not, see http://www.gnu.org/licenses/.
31|
32| The author's address is seanerikoconnor!AT!gmail!DOT!com
33| with the !DOT! replaced by . and the !AT! replaced by @
34|
35==============================================================================*/
36
37/*------------------------------------------------------------------------------
38| Include Files |
39------------------------------------------------------------------------------*/
40
41#include <cstdlib> // abort()
42#include <iostream> // Basic stream I/O.
43#include <new> // set_new_handler()
44#include <cmath> // Basic math functions e.g. sqrt()
45#include <complex> // Complex data type and operations.
46#include <fstream> // File stream I/O.
47#include <sstream> // String stream I/O.
48#include <vector> // STL vector class.
49#include <string> // STL string class.
50#include <algorithm> // Iterators.
51#include <stdexcept> // Exceptions.
52#include <cassert> // assert()
53
54using namespace std ;
55
56/*------------------------------------------------------------------------------
57| PP Include Files |
58------------------------------------------------------------------------------*/
59
60#include "Primpoly.hpp" // Global functions.
61#include "ppArith.hpp" // Basic arithmetic functions.
62#include "ppBigInt.hpp" // Arbitrary precision integer arithmetic.
63#include "ppOperationCount.hpp" // OperationCount collection for factoring and poly finding.
64#include "ppFactor.hpp" // Prime factorization and Euler Phi.
65#include "ppPolynomial.hpp" // Polynomial operations and mod polynomial operations.
66#include "ppParser.hpp" // Parsing of polynomials and I/O services.
67#include "ppUnitTest.hpp" // Complete unit test.
68
69/*=============================================================================
70 |
71 | NAME
72 |
73 | ModP::operator()
74 |
75 | DESCRIPTION
76 |
77 | Computes k = n mod p where 0 <= k < p for both positive and negative
78 | arguments n for p >= 1.
79 | Takes a surprising proportion of the total compute time according to the
80 | profiler, so worth optimizing.
81 |
82 | EXAMPLE
83 |
84 | The C language gives -5 % 3 = - ( -(-5) mod 3 ) = -( 5 mod 3 ) = -2.
85 | The result is nonzero, so we add p=3 to give 1.
86 |
87 | C computes -3 % 3 = - ( -(-3) mod 3 ) = 0, which we leave unchanged.
88 |
89 | METHOD
90 |
91 | For n >= 0, the C language defines r = n mod p by the equation
92 |
93 | n = kp + r 0 <= r < p
94 |
95 | but when n < 0, C returns the quantity
96 |
97 | r = - ( (-n) mod p )
98 |
99 | To put the result into the correct range 0 to p-1, add p to r if
100 | r is non-zero.
101 |
102 | By the way, dear old FORTRAN's MOD function does the same thing.
103 |
104 +============================================================================*/
105
106template <typename UIntType, typename SIntType>
107inline UIntType ModP<UIntType,SIntType>::operator()( SIntType n )
108{
109 if (p_ <= 0)
110 {
111 ostringstream os ;
112 os << "ModP::operator() The modulus p = " << p_ << " <= 0" ;
113 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
114 }
115
116 if (n >= 0)
117 {
118 // Avoid integer division in this special case.
119 if (p_ == 2)
120 return n - ((n >> 1) << 1) ;
121 else
122 return n % p_ ;
123 }
124 // Both types must be signed if we want n % p to compute correctly for n < 0.
125 else
126 return( (n % static_cast<SIntType>(p_)) + p_ ) ;
127}
128
129/*=============================================================================
130|
131| NAME
132|
133| PowerMod::operator()
134|
135| DESCRIPTION
136|
137| n
138| Compute a (modulo p) for generic integer types.
139| 0
140| where a >= 0, p >= 2, n >= 0. Other cases including 0
141| will throw an ArithModPError.
142|
143| We have two versions: one is generic for any integer type, and
144| the other is specialized for ppuint.
145|
146| EXAMPLE
147|
148| BigInt a = 2 ;
149| BigInt n = 3 ;
150| BigInt p = 7 ;
151| try
152| {
153| PowerMod<BigInt> powermod( p ) ;
154| BigInt pwr = powermod( a, n ) ;
155| }
156| catch( ArithModPError e )
157| {
158| cout << "should have gotten pwr = 1\n" ;
159| }
160|
161| METHOD
162|
163| Multiplication by repeated squaring.
164|
165+============================================================================*/
166
167// Generic.
168template <typename IntType>
169IntType PowerMod<IntType>::operator()( const IntType & a, const IntType & n )
170{
171 IntType a1 = a ;
172
173 // Out of range conditions.
174 if (a < static_cast<IntType>( 0u ) ||
175 n < static_cast<IntType>( 0u ) ||
176 p_ <= static_cast<IntType>( 1u ) ||
177 (a == static_cast<IntType>( 0u ) && n == static_cast<IntType>( 0u )))
178 {
179 ostringstream os ;
180 os << "PowerMod::operator() One or more parameters out of range: a = "
181 << a << " n = " << n << " p_ = " << p_ ;
182 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
183 }
184
185 // Quick return for 0 ^ n, a^0 and a^1.
186 if (a == static_cast<IntType>( 0u ))
187 return static_cast<IntType>( 0u ) ;
188
189 if (n == static_cast<IntType>( 0u ))
190 return static_cast<IntType>( 1u ) ;
191
192 if (n == static_cast<IntType>( 1u ))
193 return a % static_cast<IntType>( p_ ) ;
194
195 int bitNum = n.maxBitNumber() ; // Index to highest bit.
196
197 #ifdef DEBUG_PP_ARITH
198 cout << "initial max bitNum = " << bitNum << endl ;
199 cout << "a = " << a << endl ;
200 #endif
201
202 // Find the index of the leading 1 bit.
203 while (!n.testBit( bitNum ))
204 --bitNum ;
205
206 #ifdef DEBUG_PP_ARITH
207 cout << "after skipping leading 0 bits, bitNum = " << bitNum << endl ;
208 #endif
209
210 if (bitNum == -1)
211 {
212 ostringstream os ;
213 os << "PowerMod::operator() " << "bitNum == -1" ;
214 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
215 }
216
217 #ifdef DEBUG_PP_ARITH
218 cout << "\nAfter skipping zero bits, bitNum = " << bitNum << endl ;
219 #endif
220
221 // Exponentiation by repeated squaring. Discard the leading 1 bit.
222 // Thereafter, square for every 0 bit; square and multiply by a for
223 // every 1 bit.
224 while ( --bitNum >= 0 )
225 {
226 a1 = (a1 * a1) % static_cast<IntType>( p_ ) ; // Square mod p.
227
228 if (n.testBit( bitNum ))
229 a1 = (a1 * a) % static_cast<IntType>( p_ ) ; // Times a mod p.
230
231 #ifdef DEBUG_PP_ARITH
232 cout << "S " ;
233 if (n.testBit( bitNum ))
234 cout << "X " ;
235 cout << "Bit num = " << bitNum << " a1 = " << a1 << endl ;
236 #endif
237 }
238
239 #ifdef DEBUG_PP_ARITH
240 cout << "Out of the loop bitNum = " << bitNum << " a1 = " << a1 << endl ;
241 #endif
242
243 return a1 ;
244}
245
246
247
248/*=============================================================================
249 |
250 | NAME
251 |
252 | addMod()
253 |
254 | DESCRIPTION
255 |
256 | Computes (x + y) mod p for generic integer types, up to the maximum
257 | unsigned type. We extend the range to the limit by testing and compensating
258 | for carries internally.
259 |
260 | Example:
261 | a = 4294967290
262 | n = 4294967294
263 | (a + b) % n = 4294967286
264 |
265 +============================================================================*/
266
267template <typename IntType>
268IntType addMod( IntType a, IntType b, IntType n )
269{
270 #ifdef DEBUG_PP_ARITH
271 cout << "addMod" << endl ;
272 cout << " sizeof IntType = " << 8 * sizeof( IntType ) << " bits" << endl ;
273 cout << " a = " << a << " b = " << b << " n = " << n << endl ;
274 #endif
275
276 // Assume the positive arguments might not already be in the range [0, n).
277 a %= n ;
278 b %= n ;
279
280 // Add with carry. Carry bit thrown away by the language.
281 IntType c = a + b ;
282
283 #ifdef DEBUG_PP_ARITH
284 cout << " a mod n = " << a << " b mod n = " << b << endl ;
285 cout << " c = a + b (discarding carry bit) = " << c << endl ;
286 #endif
287
288 // Test for a carry. Most computer languages don't give us access to
289 // the carry bit, so we must infer whether a carry occurred during addition.
290 // Logically,
291 // (no carry for a+b) => a+b <= a && a+b <= b
292 // The contrapositive is
293 // a+b < a || a+b < b => carry
294 if (c < a || c < b)
295 {
296 // Lemma 1. If a carry occured, n < c < 2n.
297 // Pf.
298 // m
299 // When a carry occurs, c = a + b > 2 > n. But 0 <= a, b < n. So n < c and a + b < n + n = 2n.
300 //
301 // Lemma 2. If a carry occured, we can simply compute c - n as usual, discarding the carry.
302 // Pf.
303 // Since a carry occurred, we an write
304 // m
305 // c = a + b = 2 + c'
306 //
307 // Let's assume the computer hardware does subtraction in two's complement.
308 // Now use two's complement arithmetic on c' - n and see the results.
309 //
310 // c' + TwosComplement( n ) <discarding carry bit> =
311 //
312 // c' + (~n + 1) <discarding carry bit> =
313 // m m
314 // (c' + (2 - n)) mod 2 =
315 // m m m m m m
316 // (c - 2 + 2 + (2 - n)) mod 2 = (c + (2 - n)) mod 2 =
317 //
318 // c + TwosComplement( n ) <discarding carry bit> =
319 //
320 // c - n done the usual way in computer hardware.
321 //
322 // Th. If a carry occurs, c mod n = c - n.
323 // Pf. See Lemma 1 and 2.
324 //
325 // e.g. Let m = 4 bits and
326 //
327 // a = 1 1 0 1 = 13
328 // b = 1 1 0 1 = 13
329 // n = 1 1 1 0 = 14
330 //
331 // c = 1 1 0 1 0
332 // c' = 1 0 1 0
333 // ~n = 0 0 0 1 + 1 =
334 // 0 0 1 0
335 // c' - n = c' + ~n =
336 // 1 0 1 0
337 // + 0 0 1 0
338 // = 1 1 0 0 = 12 = 26 mod 14
339
340 c -= n ;
341
342 #ifdef DEBUG_PP_ARITH
343 cout << " Carry!" << " c < a = " << (c < a ) << "|| c < b = " << (c < b) << endl ;
344 #endif
345 }
346 else
347 {
348 c %= n ;
349
350 #ifdef DEBUG_PP_ARITH
351 cout << " No carry" << endl ;
352 #endif
353 }
354
355 return c ;
356}
357
358
359
360/*=============================================================================
361 |
362 | NAME
363 |
364 | timesTwoMod()
365 |
366 | DESCRIPTION
367 |
368 | Computes (2 * x) mod p for generic integer types, up to the maximum
369 | unsigned type. We extend the range to the limit by testing and compensating
370 | for carries internally.
371 |
372 | Example:
373 | a = 10376293541461622782
374 | n = 18446744073709551610
375 | (2 * a) % n = 2305843009213693954
376 |
377 +============================================================================*/
378
379template <typename IntType>
380IntType timesTwoMod( IntType a, IntType n )
381{
382 #ifdef DEBUG_PP_ARITH
383 cout << "timesTwoMod" << endl ;
384 cout << " sizeof IntType = " << 8 * sizeof( IntType ) << " bits" << endl ;
385 cout << " a = " << a << " n = " << n << endl ;
386 #endif
387
388 // Assume the positive argument might not already be in the range [0, n).
389 a %= n ;
390
391 #ifdef DEBUG_PP_ARITH
392 cout << " a mod n = " << a << endl ;
393 #endif
394
395 // High bit mask.
396 IntType mask = ((IntType)1 << (8 * sizeof( IntType ) - 1)) ;
397
398 IntType c = (a << 1) ; // Lose the carry bit.
399
400 #ifdef DEBUG_PP_ARITH
401 cout << " mask = " << hex << mask << dec << endl ;
402 cout << " c = 2 a (losing the carry bit) = " << c << endl ;
403 #endif
404
405 // High bit is set so (2 a) will have a carry.
406 if (mask & a)
407 {
408 // See notes in function addMod() above for why this works correctly.
409 c -= n ;
410
411 #ifdef DEBUG_PP_ARITH
412 cout << " Carry!" << endl ;
413 #endif
414 }
415 else
416 {
417 c %= n ;
418
419 #ifdef DEBUG_PP_ARITH
420 cout << " No carry" << endl ;
421 #endif
422 }
423
424 return c ;
425}
426
427
428
429/*=============================================================================
430 |
431 | NAME
432 |
433 | multiplyMod()
434 |
435 | DESCRIPTION
436 |
437 | Computes (x * y) mod p for generic integer types, up to the maximum
438 | unsigned type. We extend the range this far by testing and compensating
439 | for carries internally.
440 |
441 | Multiply a * b mod n
442 |
443 | m-1 m-2
444 | a b = a (b 2 + b 2 + ... + b 2 + b )
445 | m-1 m-2 1 0
446 |
447 | Expanding using Horner's rule,
448 |
449 | a b = (0 2 + a b ) 2 + a b ) 2 + ... + a b ) 2 + a b
450 | m-1 m-2 1 0
451 |
452 | where b = 0 or 1
453 | i
454 | At each step call specialized functions to do r = 2 a mod n and r = r + a mod n
455 | which will work even when carries occur.
456 |
457 | The process takes O( ln n ) operations. An example,
458 |
459 | a = 4294967290
460 | b = 4294967290
461 | n = 4294967294
462 | (a * b) % n = 16
463 |
464 | Reference: "Comments on 'A Computer Algorithm for Calculating the Product
465 | A B Modulo M'" K. R. Sloane, IEEE Transactions on Computers, Vol C-34, No. 3,
466 | March 1985.
467 |
468 +============================================================================*/
469
470template<typename IntType>
471IntType multiplyMod( const IntType a, const IntType b, const IntType n )
472{
473 // High bit mask.
474 const int numBits = 8 * sizeof( IntType ) ;
475 IntType mask = ((IntType)1 << (numBits - 1)) ;
476
477 IntType r { 0 } ;
478
479 for (int i = numBits-1 ; i >= 0 ; --i, mask >>= 1)
480 {
481 r = timesTwoMod( r, n ) ;
482
483 if (mask & b)
484 r = addMod( r, a, n ) ;
485 }
486
487 return r ;
488}
489
490
491
492/*=============================================================================
493 |
494 | NAME
495 |
496 | PowerMod::operator()
497 |
498 | DESCRIPTION
499 |
500 | Specialized for ppuint type.
501 |
502 +============================================================================*/
503
504template<>
505ppuint PowerMod<ppuint>::operator()( const ppuint & a, const ppuint & n )
506{
507 // mask = 1000 ... 000. That is, all bits of mask are zero except for the
508 // most significant bit of the computer word holding its value.
509 ppuint mask{ static_cast<ppuint>(1) << (8 * sizeof( ppuint ) - 1) } ;
510 int bit_count{ 0 } ; // Number of bits in exponent to the right of the leading bit.
511 ppuint n1{ n } ;
512 ppuint product{ a } ;
513
514 // Out of range conditions.
515 if (p_ <= 1 || (a == 0 && n1 == 0))
516 {
517 ostringstream os ;
518 os << "PowerMod<ppuint>::operator() parameters are out of range: a = " << a << " n = " << n << " p_ = " << p_ ;
519 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
520 }
521
522 // n 0 1
523 // Quick return for 0 , a and a
524 if (a == 0)
525 return 0 ;
526
527 if (n == 0)
528 return 1 ;
529
530 if (n == 1)
531 return a % p_ ;
532
533 #ifdef DEBUG_PP_ARITH
534 cout << "a = " << a << endl ;
535 cout << "n = " << n << endl ;
536 cout << "p = " << p_ << endl ;
537 cout << "n1 (before shifting) = " << n1 << endl ;
538 #endif
539
540 // Advance the leading bit of the exponent up to the word's left hand boundary.
541 // Count how many bits were to the right of the leading bit.
542 while (! (n1 & mask))
543 {
544 n1 <<= 1 ;
545 ++bit_count ;
546 }
547
548 bit_count = (8 * sizeof( ppuint )) - bit_count ;
549
550 #ifdef DEBUG_PP_ARITH
551 cout << "n1 = " << n1 << endl ;
552 cout << "mask = " << mask << endl ;
553 cout << "bit_count = " << bit_count << endl ;
554 #endif
555
556 // Exponentiation by repeated squaring. Discard the leading 1 bit.
557 // Thereafter, square for every 0 bit; square and multiply by x for every 1 bit.
558 while ( --bit_count > 0 )
559 {
560 #ifdef DEBUG_PP_ARITH
561 cout << "product (before squaring) = " << product << " n1 = " << n1 << endl ;
562 #endif
563
564 // Expose the next bit.
565 n1 <<= 1 ;
566
567 // Out of range conditions.
568 if (product > BigInt::getBase() || a > BigInt::getBase())
569 {
570 // Square modulo p.
571 product = multiplyMod( product, product, p_ ) ;
572
573 // Leading bit is 1: multiply by a modulo p.
574 if (n1 & mask)
575 product = multiplyMod( a, product, p_ ) ;
576 }
577 else
578 {
579 // Square modulo p.
580 product = (product * product) % p_ ;
581
582 // Leading bit is 1: multiply by a modulo p.
583 if (n1 & mask)
584 product = (a * product) % p_ ;
585 }
586
587 #ifdef DEBUG_PP_ARITH
588 cout << "S " ;
589 if (n1 & mask)
590 cout << "X " ;
591 cout << "product = " << product << " n1 = " << n1 << endl ;
592 #endif
593 }
594
595 return product ;
596}
597
598
599
600/*=============================================================================
601|
602| NAME
603|
604| IsPrimitiveRoot::operator()
605|
606| DESCRIPTION
607|
608| Returns true if a is a primitive root of prime p, and false otherwise.
609| Throws ArithModPError if a < 1 or p < 2 or p is even.
610| p must be a prime number, but we only test that p is 2 or odd.
611|
612| EXAMPLE
613|
614| ppuint p{ 7 } ; ppuint a{ 3 } ;
615| IsPrimitiveRoot isRoot( p ) ;
616| try
617| {
618| bool isPrim = isRoot( a ) ;
619| }
620| catch( ArithModPError e )
621| {
622| cout << "out of range" << endl ;
623| }
624| 1 2 3
625| 3 is a primitive root of the prime p = 7 since 3 = 3, 3 = 2, 3 = 6,
626| 4 5 p-1 6
627| 3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3 = 3 = 1 (mod 7).
628|
629| So 3 is a primitive root of 7 because it has maximal period. a = 2
630| 3
631| isn't a primitive root of p = 7 because already 2 = 1 (mod 7).
632|
633| METHOD
634|
635| From number theory, a is a primitive root of the prime p iff
636| (p-1)/q
637| a != 1 (mod p) for all prime divisors q of (p-1).
638| (p-1)
639| Don't need to check if a = 1 (mod p): Fermat's little
640| theorem guarantees it.
641|
642+============================================================================*/
643
644bool IsPrimitiveRoot::operator()( ppuint a )
645{
646 PowerMod<ppuint> powermod( p_ ) ;
647
648 // a = 0 (mod p); Zero can't be a primitive root of p.
649 if (a == 0)
650 return false ;
651
652
653 // Error for out of range inputs, including p
654 // being an even number greater than 2.
655 if (p_ < 2 || a < 1 || (p_ > 2 && (p_ % 2 == 0)))
656 {
657 ostringstream os ;
658 os << "IsPrimitiveRoot::operator() Inputs are out of range: p = " << p_ << " a = " << a ;
659 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
660 }
661
662 // Special cases:
663 // 1 is the only primitive root of 2; i.e. 1 generates the unit elements
664 // of GF( 2 ); 2 is the only primitive root of 3, and 2 and 3 are the only
665 // primitive roots of 5.
666 //
667 if ( (p_ == 2 && a == 1) ||
668 (p_ == 3 && a == 2) ||
669 (p_ == 5 && (a == 2 || a == 3)) ||
670 (p_ == 7 && (a == 3 || a == 5)) ||
671 (p_ == 11 && (a == 2 || a == 6 || a == 7 || a == 8)) ||
672 (p_ == 13 && (a == 2 || a == 6 || a == 7 || a == 11)))
673 {
674 return true ;
675 }
676
677 // Reduce a down modulo p.
678 a = a % p_ ;
679
680 // a = 0 (mod p); Zero can't be a primitive root of p.
681 if (a == 0)
682 return false ;
683
684 // Factor p-1 into distinct primes.
685 Factorization<ppuint> factorization( p_ - 1 ) ;
686
687 // p-1
688 // a
689 // Test --- != 1 (mod p) for all primes q | (p-1).
690 // q
691 // If so, we have a primitive root of p, otherwise, we exit early.
692 //
693 for (unsigned int i = 0 ; i < factorization.numDistinctFactors() ; ++i)
694 {
695 if (powermod( a, (p_ - 1) / factorization.primeFactor( i )) == 1)
696 {
697 return false ;
698 }
699 }
700
701 return true ;
702}
703
704
705
706/*=============================================================================
707 |
708 | NAME
709 |
710 | inverse_mod_p
711 |
712 | DESCRIPTION
713 |
714 | -1
715 | Find the inverse of u modulo p. In other words, find u (mod p)
716 | p >= 2.
717 |
718 | EXAMPLE
719 | -1
720 | For p=7, and u = 2, u = 4 because 2 * 4 = 8 (mod 7) = 1.
721 |
722 | METHOD
723 |
724 | Do extended Euclid's algorithm on u and v to find u1, u2, and u3 such that
725 |
726 | u u1 + v u2 = u3 = gcd( u, v ).
727 |
728 | Now let v = p = prime number, so gcd = u3 = 1. Then we get
729 |
730 | u u1 + p ? = 1
731 |
732 | or u u1 = 0 (mod p) which makes u (mod p) the unique multiplicative
733 | inverse of u.
734 |
735 +============================================================================*/
736
737ppsint InverseModP::operator()( ppsint u )
738{
739 ModP<ppsint,ppsint> mod( p_ ) ;
740
741 // Do Euclid's algorithm to find the quantities.
742 ppsint t1 = 0 ;
743 ppsint t3 = 0 ;
744 ppsint q = 0 ;
745
746 ppsint u1 = 1 ;
747 ppsint u3 = u ;
748 ppsint v1 = 0 ;
749 ppsint v3 = p_ ;
750
751 ppsint inv_v = 0 ;
752
753 while( v3 != 0)
754 {
755 // Casting to do a floor function.
756 q = static_cast<ppsint>(u3 / v3) ;
757
758 t1 = u1 - v1 * q ;
759 t3 = u3 - v3 * q ;
760
761 u1 = v1 ;
762 u3 = v3 ;
763
764 v1 = t1 ;
765 v3 = t3 ;
766 }
767
768 inv_v = mod( u1 ) ;
769
770 // -1
771 // Self check: does u u = 1 (mod p)?
772 //
773 if ( mod( u * inv_v ) != 1)
774 {
775 ostringstream os ;
776 os << "InverseModP::operator(): inverse self check failed: u = " << u << " * u^(-1) = " << inv_v << " != 1" ;
777 throw ArithModPError( os.str(), __FILE__, __LINE__ ) ;
778 }
779
780 return inv_v ;
781}
782
783
784/*=============================================================================
785 |
786 | NAME
787 |
788 | constCoeffTest
789 |
790 | DESCRIPTION
791 |
792 | n
793 | Test if a = (-1) a (mod p) where a is the constant coefficient
794 | 0 0
795 | of polynomial f(x) and a is the number from the orderR() test.
796 | Return true if we pass the test, false otherwise.
797 |
798 | EXAMPLE
799 | 11 3
800 | Let p = 5, n = 11, f( x ) = x + 2 x + 1, and a = 4.
801 | Thus a = 1.
802 | 0
803 | Then return true since
804 |
805 | 11
806 | (-1) * 1 (mod 5) = -1 (mod 5) = 4 (mod 5).
807 |
808 | METHOD
809 | n
810 | We test if (a - (-1) a ) mod p = 0.
811 | 0
812 |
813 +============================================================================*/
814
815bool ArithModP::constCoeffTest( ppsint a, ppsint a0, int n )
816{
817 ppsint constant_coeff = a0 ;
818
819 ModP<ppuint,ppsint> mod( p_ ) ; // Initialize the functionoid.
820
821 if (n % 2 != 0)
822 constant_coeff = -constant_coeff ; // (-1)^n < 0 for odd n.
823
824 return( (mod( a - constant_coeff ) == 0) ? true : false ) ;
825}
826
827
828
829/*=============================================================================
830 |
831 | NAME
832 |
833 | constCoeffIsPrimitiveRoot
834 |
835 | DESCRIPTION
836 |
837 | n
838 | Test if (-1) a (mod p) is a primitive root of the prime p where
839 | 0
840 | a is the constant term of the polynomial f(x).
841 | 0
842 |
843 | EXAMPLE
844 | 11 3
845 | Let p = 7, n = 11, f( x ) = x + 2 x + 4. Thus a = 4.
846 | 0
847 | Then return true since
848 |
849 | 11
850 | (-1) * 4 (mod 7) = -4 (mod 7) = 3 (mod 7), and 3 is a
851 | 1 2 3
852 | primitive root of the prime 7 since 3 = 3, 3 = 2, 3 = 6,
853 | 4 5 7-1
854 | 3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3 = 1 (mod 7).
855 |
856 | METHOD
857 |
858 | n
859 | We test if (-1) a mod p is a primitive root mod p.
860 | 0
861 |
862 +============================================================================*/
863
864bool ArithModP::constCoeffIsPrimitiveRoot( ppuint a0, int n )
865{
866 ppsint constant_coeff = a0 ;
867
868 ModP<ppuint,ppsint> mod( p_ ) ;
869 IsPrimitiveRoot isroot( p_ ) ;
870
871 // (-1)^n < 0 for odd n.
872 if (n % 2 != 0)
873 constant_coeff = -constant_coeff ;
874
875 return isroot( mod( constant_coeff ) ) ;
876}
877
878
879
880/*=============================================================================
881|
882| NAME
883|
884| gcd
885|
886| DESCRIPTION
887|
888| Euclid's method for computing greatest common divisor.
889|
890| METHOD
891|
892| Iterate the equation
893|
894| gcd( u, v ) = gcd( v, r )
895|
896| where r = | u / v | = u - q v
897| -- --
898| until v is zero, terminating in gcd( u, 0 ) = u.
899|
900| Let's use Mathematica to illustrate the method for large integers,
901|
902| gcd[ u_, v_] :=
903| (* Compute GCD using Euclid's method. *)
904| Module[ {u2, v2, r}, u2 = u ; v2 = v ;
905| Print[ "Computing GCD( u = ", u, " v = ", v , " )" ] ;
906| While[ v2 != 0,
907| r = Mod[u2, v2]; u2 = v2; v2 = r;
908| Print[ "r = ", r, " u2 = ", u2, " v2 = ", v2]];
909| Return[ u2 ]]
910|
911| gcd[ 779953197883173551166308319545, 1282866356929526866866376009397 ]
912|
913| Computing GCD( u = 779953197883173551166308319545 v = 1282866356929526866866376009397 )
914| u = 779953197883173551166308319545 v = 1282866356929526866866376009397
915| r = 779953197883173551166308319545 u2 = 1282866356929526866866376009397 v2 = 779953197883173551166308319545
916| r = 502913159046353315700067689852 u2 = 779953197883173551166308319545 v2 = 502913159046353315700067689852
917| r = 277040038836820235466240629693 u2 = 502913159046353315700067689852 v2 = 277040038836820235466240629693
918| r = 225873120209533080233827060159 u2 = 277040038836820235466240629693 v2 = 225873120209533080233827060159
919| r = 51166918627287155232413569534 u2 = 225873120209533080233827060159 v2 = 51166918627287155232413569534
920| r = 21205445700384459304172782023 u2 = 51166918627287155232413569534 v2 = 21205445700384459304172782023
921| r = 8756027226518236624068005488 u2 = 21205445700384459304172782023 v2 = 8756027226518236624068005488
922| r = 3693391247347986056036771047 u2 = 8756027226518236624068005488 v2 = 3693391247347986056036771047
923| r = 1369244731822264511994463394 u2 = 3693391247347986056036771047 v2 = 1369244731822264511994463394
924| r = 954901783703457032047844259 u2 = 1369244731822264511994463394 v2 = 954901783703457032047844259
925| r = 414342948118807479946619135 u2 = 954901783703457032047844259 v2 = 414342948118807479946619135
926| r = 126215887465842072154605989 u2 = 414342948118807479946619135 v2 = 126215887465842072154605989
927| r = 35695285721281263482801168 u2 = 126215887465842072154605989 v2 = 35695285721281263482801168
928| r = 19130030301998281706202485 u2 = 35695285721281263482801168 v2 = 19130030301998281706202485
929| r = 16565255419282981776598683 u2 = 19130030301998281706202485 v2 = 16565255419282981776598683
930| r = 2564774882715299929603802 u2 = 16565255419282981776598683 v2 = 2564774882715299929603802
931| r = 1176606122991182198975871 u2 = 2564774882715299929603802 v2 = 1176606122991182198975871
932| r = 211562636732935531652060 u2 = 1176606122991182198975871 v2 = 211562636732935531652060
933| r = 118792939326504540715571 u2 = 211562636732935531652060 v2 = 118792939326504540715571
934| r = 92769697406430990936489 u2 = 118792939326504540715571 v2 = 92769697406430990936489
935| r = 26023241920073549779082 u2 = 92769697406430990936489 v2 = 26023241920073549779082
936| r = 14699971646210341599243 u2 = 26023241920073549779082 v2 = 14699971646210341599243
937| r = 11323270273863208179839 u2 = 14699971646210341599243 v2 = 11323270273863208179839
938| r = 3376701372347133419404 u2 = 11323270273863208179839 v2 = 3376701372347133419404
939| r = 1193166156821807921627 u2 = 3376701372347133419404 v2 = 1193166156821807921627
940| r = 990369058703517576150 u2 = 1193166156821807921627 v2 = 990369058703517576150
941| r = 202797098118290345477 u2 = 990369058703517576150 v2 = 202797098118290345477
942| r = 179180666230356194242 u2 = 202797098118290345477 v2 = 179180666230356194242
943| r = 23616431887934151235 u2 = 179180666230356194242 v2 = 23616431887934151235
944| r = 13865643014817135597 u2 = 23616431887934151235 v2 = 13865643014817135597
945| r = 9750788873117015638 u2 = 13865643014817135597 v2 = 9750788873117015638
946| r = 4114854141700119959 u2 = 9750788873117015638 v2 = 4114854141700119959
947| r = 1521080589716775720 u2 = 4114854141700119959 v2 = 1521080589716775720
948| r = 1072692962266568519 u2 = 1521080589716775720 v2 = 1072692962266568519
949| r = 448387627450207201 u2 = 1072692962266568519 v2 = 448387627450207201
950| r = 175917707366154117 u2 = 448387627450207201 v2 = 175917707366154117
951| r = 96552212717898967 u2 = 175917707366154117 v2 = 96552212717898967
952| r = 79365494648255150 u2 = 96552212717898967 v2 = 79365494648255150
953| r = 17186718069643817 u2 = 79365494648255150 v2 = 17186718069643817
954| r = 10618622369679882 u2 = 17186718069643817 v2 = 10618622369679882
955| r = 6568095699963935 u2 = 10618622369679882 v2 = 6568095699963935
956| r = 4050526669715947 u2 = 6568095699963935 v2 = 4050526669715947
957| r = 2517569030247988 u2 = 4050526669715947 v2 = 2517569030247988
958| r = 1532957639467959 u2 = 2517569030247988 v2 = 1532957639467959
959| r = 984611390780029 u2 = 1532957639467959 v2 = 984611390780029
960| r = 548346248687930 u2 = 984611390780029 v2 = 548346248687930
961| r = 436265142092099 u2 = 548346248687930 v2 = 436265142092099
962| r = 112081106595831 u2 = 436265142092099 v2 = 112081106595831
963| r = 100021822304606 u2 = 112081106595831 v2 = 100021822304606
964| r = 12059284291225 u2 = 100021822304606 v2 = 12059284291225
965| r = 3547547974806 u2 = 12059284291225 v2 = 3547547974806
966| r = 1416640366807 u2 = 3547547974806 v2 = 1416640366807
967| r = 714267241192 u2 = 1416640366807 v2 = 714267241192
968| r = 702373125615 u2 = 714267241192 v2 = 702373125615
969| r = 11894115577 u2 = 702373125615 v2 = 11894115577
970| r = 620306572 u2 = 11894115577 v2 = 620306572
971| r = 108290709 u2 = 620306572 v2 = 108290709
972| r = 78853027 u2 = 108290709 v2 = 78853027
973| r = 29437682 u2 = 78853027 v2 = 29437682
974| r = 19977663 u2 = 29437682 v2 = 19977663
975| r = 9460019 u2 = 19977663 v2 = 9460019
976| r = 1057625 u2 = 9460019 v2 = 1057625
977| r = 999019 u2 = 1057625 v2 = 999019
978| r = 58606 u2 = 999019 v2 = 58606
979| r = 2717 u2 = 58606 v2 = 2717
980| r = 1549 u2 = 2717 v2 = 1549
981| r = 1168 u2 = 1549 v2 = 1168
982| r = 381 u2 = 1168 v2 = 381
983| r = 25 u2 = 381 v2 = 25
984| r = 6 u2 = 25 v2 = 6
985| r = 1 u2 = 6 v2 = 1
986| r = 0 u2 = 1 v2 = 0
987| ---> 1
988|
989+============================================================================*/
990
991template <typename IntType>
992IntType gcd( const IntType & u, const IntType & v )
993{
994 IntType r ;
995 IntType u2 { u } ;
996 IntType v2 { v } ;
997
998 #ifdef DEBUG_PP_ARITH
999 cout << "gcd: u = " << u << " v = " << v << endl ;
1000 #endif
1001
1002 while (v2 != static_cast<IntType>(0u))
1003 {
1004 r = u2 % v2 ;
1005 u2 = v2 ;
1006 v2 = r ;
1007
1008 #ifdef DEBUG_PP_ARITH
1009 cout << " r = " << r << " u2 = " << u2 << " v2 = " << v2 << endl ;
1010 #endif
1011 }
1012
1013 return u2 ;
1014}
1015
1016
1017/*==============================================================================
1018| Forced Template Instantiations |
1019==============================================================================*/
1020
1021// C++ doesn't automatically generate templated classes or functions until they
1022// are first used. So depending on the order of compilation, the linker may say
1023// the templated functions are undefined.
1024//
1025// Therefore, explicitly instantiate ALL templates here:
1026
1027template ppuint gcd( const ppuint &, const ppuint & ) ;
1028template BigInt gcd( const BigInt &, const BigInt & ) ;
1029
1030template ppuint addMod( ppuint a, ppuint b, ppuint n ) ;
1031template ppuint timesTwoMod( ppuint, ppuint ) ;
1032template ppuint multiplyMod( const ppuint, const ppuint, const ppuint ) ;
1033
1034// We already specialized this function for ppuint in the source code implementation above, so we can omit
1035// the template instantiation: template ppuint PowerMod<ppuint>::operator()( const ppuint & a, const ppuint & n ) ;
1036template BigInt PowerMod<BigInt>::operator()( const BigInt &, const BigInt & ) ;
1037
1038template ModP<ppuint,ppsint>::ModP( ppuint ) ;
1039template ModP<ppuint,ppsint>::ModP( const ModP & ) ;
1040template ppuint ModP<ppuint,ppsint>::operator()( ppsint ) ;