1/*==============================================================================
2|
3| NAME
4|
5| ppPolynomial.cpp
6|
7| DESCRIPTION
8|
9| Polynomial arithmetic and polynomial exponentiation classes.
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.4 - A Program for Computing Primitive Polynomials.
17| Copyright (C) 1999-2025 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 <iomanip> // I/O manipulators.
44#include <new> // set_new_handler()
45#include <cmath> // Basic math functions e.g. sqrt()
46#include <complex> // Complex data type and operations.
47#include <fstream> // File stream I/O.
48#include <sstream> // String stream I/O.
49#include <vector> // STL vector class.
50#include <string> // STL string class.
51#include <algorithm> // Iterators.
52#include <stdexcept> // Exceptions.
53#include <cassert> // assert()
54
55using namespace std ;
56
57
58/*------------------------------------------------------------------------------
59| PP Include Files |
60------------------------------------------------------------------------------*/
61
62#include "Primpoly.hpp" // Global functions.
63#include "ppArith.hpp" // Basic arithmetic functions.
64#include "ppBigInt.hpp" // Arbitrary precision integer arithmetic.
65#include "ppOperationCount.hpp" // OperationCount collection for factoring and poly finding.
66#include "ppFactor.hpp" // Prime factorization and Euler Phi.
67#include "ppPolynomial.hpp" // Polynomial operations and mod polynomial operations.
68#include "ppParser.hpp" // Parsing of polynomials and I/O services.
69#include "ppUnitTest.hpp" // Complete unit test.
70
71
72/*=============================================================================
73|
74| NAME
75|
76| Polynomial()
77|
78| DESCRIPTION
79|
80| Default constructor for Polynomial class. Constructs the zero degree
81| polynomial p(x) = 0 (mod 2)
82|
83| EXAMPLE
84|
85| Polynomial p ;
86|
87+============================================================================*/
88
89Polynomial::Polynomial()
90 : f_()
91 , n_( 0 )
92 , p_( 2 )
93 , mod( p_ )
94{
95 // f(x) = 0
96 f_.push_back( 0 ) ;
97}
98
99
100/*=============================================================================
101 |
102 | NAME
103 |
104 | Polynomial()
105 |
106 | DESCRIPTION
107 |
108 | Constructor for a polynomial from a vector of integers.
109 |
110 | EXAMPLE
111 |
112 | vector<ppunit> v { 1, 2, 3 } ;
113 | Polynomial p{ v } ;
114 | Polynomial p( v ) ;
115 |
116 +============================================================================*/
117
118Polynomial::Polynomial( const vector<ppuint> v )
119: n_{ static_cast<int>( v.size() - 1) }
120, p_( 2 )
121, mod( p_ )
122{
123 // Copy over the polynomial coefficients.
124 f_ = v ;
125}
126
127
128/*=============================================================================
129|
130| NAME
131|
132| ~Polynomial
133|
134| DESCRIPTION
135|
136| Destructor.
137|
138| EXAMPLE
139|
140| void add1( const Polynomial & f )
141| {
142| const Polynomial g{ 1u } ;
143| return f + g ;
144| // Destructor for g automatically called as f goes out of scope.
145| }
146|
147+============================================================================*/
148
149Polynomial::~Polynomial()
150{
151 // vector f_ frees itself and mod_ calls its own destructor.
152}
153
154
155/*=============================================================================
156|
157| NAME
158|
159| Polynomial
160|
161| DESCRIPTION
162|
163| Copy constructor.
164|
165| EXAMPLE
166|
167| try
168| {
169| Polynomial f ;
170| Polynomial f( g ) ;
171| }
172| catch( PolynomialError & e )
173| {
174| cout << "Error in constructing polynomial f(x) or g(x)" << endl ;
175| }
176|
177+============================================================================*/
178
179Polynomial::Polynomial( const Polynomial & g )
180 : f_( g.f_ )
181 , n_( g.n_ )
182 , p_( g.p_ )
183 , mod( p_ )
184{
185 // The classes in the initializer above all copy themselves.
186}
187
188
189/*=============================================================================
190|
191| NAME
192|
193| Polynomial::operator=
194|
195| DESCRIPTION
196|
197| Safe assigment operator for polynomials, f( x ) = g( x )
198| which leaves the polynomial f( x ) unchanged if an exception is thrown.
199|
200| EXAMPLE
201|
202| try
203| {
204| Polynomial f ;
205| Polynomial g ;
206| f = g ;
207| }
208| catch( PolynomialError & e )
209| {
210| cout << "Error in constructing polynomial f(x)" << endl ;
211| }
212|
213+============================================================================*/
214
215Polynomial & Polynomial::operator=( const Polynomial & g )
216{
217 // Check for assigning to oneself by comparing the unique pointers to the classes for speed.
218 // If the user does f = f, just pass back a reference to the unchanged polynomial f.
219 if (this == &g)
220 return *this ;
221
222 // Assign the scalars first.
223 n_ = g.n_ ;
224 p_ = g.p_ ;
225
226 // And the modulus functionoid.
227 mod( g.p_ ) ;
228
229 // Overwrite the old polynomial coefficients in f_ with the new coefficients in g.f_:
230 // 1) Delete the old polynomial coefficients, i.e. destruct the vector valued member variable f_.
231 // 2) Construct a new vector f_.
232 // 3) Copy the coefficients from g to f_.
233 // But here's the problem: if we fail to construct the new f_ and throw an exception,
234 // e.g by requesting a bad vector size, we've destroyed the existing polynomial coefficients f_.
235 //
236 // The following solution guarantees that if f = g throws an exception, the value of f is unchanged.
237 try
238 {
239 // Create a temporary copy of the polynomial coefficients.
240 vector<ppuint> tempCoeff{ g.f_ } ;
241
242 // Move the old values into the temporary, and the new values into the object.
243 // The library function swap() exchanges the values in the two containers,
244 // guarantees no exceptions will be thrown.
245 // The temporary containing the old values will be destroyed when we leave scope.
246 swap( f_, tempCoeff ) ;
247 }
248 // Failed to construct tempPoly.
249 catch( bad_alloc & e )
250 {
251 throw PolynomialError( "Polynomial.operator=() had a bad_alloc memory allocation failure", __FILE__, __LINE__ ) ;
252 }
253
254 // Return a reference to the altered object.
255 return *this ;
256}
257
258
259/*=============================================================================
260|
261| NAME
262|
263| Polynomial::operator=( string )
264|
265| DESCRIPTION
266|
267| Assigment operator for string to polynomial, f( x ) = "polynomial"
268|
269| EXAMPLE
270|
271| try
272| {
273| Polynomial f ;
274| f = "x^2 + 1" ;
275| }
276| catch( PolynomialError & e )
277| {
278| cout << "Error in constructing polynomial f(x)" << endl ;
279| }
280|
281+============================================================================*/
282
283Polynomial & Polynomial::operator=( string s )
284{
285 try
286 {
287 // Construct a new polynomial from the string.
288 Polynomial g( s ) ;
289
290 n_ = g.n_ ;
291 p_ = g.p_ ;
292 mod( g.p_ ) ;
293
294 // Right, no exceptions were thrown from the constructor, so
295 // we've got a new polynomial object now.
296 // Trash the existing polynomial.
297 f_.clear() ;
298 f_ = g.f_ ;
299 }
300 catch( bad_alloc & e )
301 {
302 throw PolynomialRangeError( "Polynomial.operator=(string) had a bad_alloc memory allocation failure", __FILE__, __LINE__ ) ;
303 }
304
305 // Return a reference to assigned object to make chaining possible.
306 return *this ;
307}
308
309
310/*=============================================================================
311|
312| NAME
313|
314| Polynomial, construct from string.
315|
316| DESCRIPTION
317|
318| Construct a polynomial from a string.
319|
320| EXAMPLE
321|
322| try
323| {
324| Polynomial f( "x^2 + 2x + 1, 3" ) ;
325| }
326| catch( PolynomialRangeError & e )
327| {
328| cout << "Error in construction polynomial f(x) from string" << endl ;
329| }
330|
331 , n_( 0 )
332 , p_( 2 )
333+============================================================================*/
334
335Polynomial::Polynomial( string s, ppuint p )
336 : f_()
337 , mod( 0 )
338 , n_( 0 )
339 , p_( 2 )
340{
341 // The polynomial string must have at least one character in it.
342 // Not really an exception but more of a user input error.
343 if (s.empty())
344 throw PolynomialRangeError( "polynomial string is empty" ) ;
345
346 try
347 {
348 // Initialize an LALR(1) parser for polynomials.
349 PolyParser< PolySymbol, PolyValue > pp ;
350
351 PolyValue v = pp.parse( s ) ;
352
353 // Get the modulus specified by the polynomial.
354 p_ = v.scalar_ ;
355
356 // If the modulus is explicitly input, use that instead of the polynomial's modulus.
357 if (p > 0)
358 p_ = p ;
359
360 // Sanity check the modulus.
361 if (p_ <= 0)
362 {
363 ostringstream os ;
364 os << "Error. Polynomial modulus p must be > 0 but p = " << p_ << endl ;
365 throw PolynomialRangeError( os.str() ) ;
366 };
367 // TODO: check upper range.
368
369 mod.set( p_ ) ;
370
371 // Sanity check the degree of the polynomial.
372 n_ = static_cast<int>( v.f_.size() ) - 1 ;
373 if (n_ < 0)
374 {
375 ostringstream os ;
376 os << "Error. Polynomial degree n must be >= 0 but n = " << n_ << endl ;
377 throw PolynomialRangeError( os.str() ) ;
378 }
379
380 // Reduce all the (positive) polynomial coefficients modulo p.
381 vector<ppuint>::iterator i ;
382 for (i = v.f_.begin() ; i != v.f_.end() ; ++i)
383 *i = mod( *i ) ;
384
385 // Copy over the polynomial coefficients.
386 f_ = v.f_ ;
387 }
388 catch( ParserError & e )
389 {
390 ostringstream os ;
391 os << "Error in parser converting polynomial from string " << s << " for p = " << p_ << " " << e.what() ;
392 throw PolynomialRangeError( os.str() ) ;
393 }
394}
395
396
397/*=============================================================================
398|
399| NAME
400|
401| string
402|
403| DESCRIPTION
404|
405| Convert a polynomial to a string.
406|
407| EXAMPLE
408|
409| try
410| {
411| Polynomial f( "x^2 + 1,3" ) ;
412| string poly = f ;
413| }
414| catch( PolynomialRangeError & e )
415| cout << "Error in construction polynomial f(x) from string" << endl ;
416|
417+============================================================================*/
418
419// Operator casting to string type.
420Polynomial::operator string() const
421{
422 // Set up a string stream for convenience.
423 ostringstream os ;
424
425 // Spin out the polynomial coefficients from high to low degree.
426 // Special case of f(x) = const.
427 if (f_.size() == 1)
428 {
429 if (!(os << f_[ 0 ]))
430 {
431 ostringstream os ;
432 os << "Error in converting polynomial to string: "
433 << " with n = " << n_ << " and p = " << p_ ;
434 throw PolynomialRangeError( os.str() ) ;
435 }
436 }
437 else
438 {
439 int lowestDegreeTerm = -1 ;
440 for (int deg = n_ ; deg >= 0 ; --deg)
441 if (f_[ deg ] != 0)
442 lowestDegreeTerm = deg ;
443
444 for (int deg = n_ ; deg >= 0 ; --deg)
445 {
446 if (f_[ deg ] != 0)
447 {
448 // x^n has a nonzero coefficient.
449 ppuint coeff = f_[ deg ] ;
450
451 // Print coeff of x^n unless it is 1.
452 // But print the constant term regardless.
453 if (coeff != 1 || deg == 0)
454 {
455 string extraBlank = deg == 0 ? "" : " " ;
456
457 if (!(os << coeff << extraBlank))
458 {
459 ostringstream os ;
460 os << "Error in converting polynomial to string: "
461 << " with n = " << n_ << " and p = " << p_ ;
462 throw PolynomialRangeError( os.str() ) ;
463 }
464 }
465
466 // Print x (no exponent).
467 if (deg == 1)
468 {
469 if ( !(os << "x") )
470 {
471 ostringstream os ;
472 os << "Error in converting polynomial to string: "
473 << " with n = " << n_ << " and p = " << p_ ;
474 throw PolynomialRangeError( os.str() ) ;
475 }
476 }
477 // Print x^n ... x^2
478 else if (deg != 0)
479 {
480 if (!(os << "x ^ " << deg))
481 {
482 ostringstream os ;
483 os << "Error in converting polynomial to string: "
484 << " with n = " << n_ << " and p = " << p_
485 << " in file " << __FILE__ << " line " << __LINE__ ;
486 throw PolynomialRangeError( os.str() ) ;
487 }
488 }
489
490 // Print +, but only when followed by a lower degree term or constant.
491 // e.g. x^2 + 2 x, x + 3
492 if (lowestDegreeTerm != -1 && deg > lowestDegreeTerm)
493 if (!(os << " + "))
494 {
495 ostringstream os ;
496 os << "Error in converting polynomial to string with n = " << n_ << " and p = " << p_ ;
497 throw PolynomialRangeError( os.str() ) ;
498 }
499 } // end coeff != 0.
500 } // end for deg
501 } // end f(x) = const.
502
503 // Print out the modulus.
504 if (!(os << ", " << p_ ))
505 {
506 ostringstream os ;
507 os << "Error in converting polynomial to string: "
508 << " with n = " << n_ << " and p = " << p_ ;
509 throw PolynomialRangeError( os.str() ) ;
510 }
511
512 // Return the string from the stream.
513 return os.str() ;
514}
515
516
517/*=============================================================================
518|
519| NAME
520|
521| operator << for Polynomial
522|
523| DESCRIPTION
524|
525| Print a polynomial to the output stream.
526|
527| EXAMPLE
528|
529| try
530| {
531| Polynomial f( "x^2 + 1,3" ) ;
532| cout << f << endl ;
533| }
534| catch( PolynomialRangeError & e )
535| {
536| cout << "Error in outputting polynomial f(x)" << endl ;
537| }
538|
539+============================================================================*/
540
541ostream & operator<<( ostream & out, const Polynomial & p )
542{
543 // Cast to polynomial to a string first, then output as a string.
544 // May throw a PolynomialRangeError.
545 out << static_cast<string>( p ) ;
546
547 return out ;
548}
549
550
551/*=============================================================================
552|
553| NAME
554|
555| Operator >> for Polynomial
556|
557| DESCRIPTION
558|
559| Polynomial stream input.
560|
561| EXAMPLE
562|
563| try
564| {
565|
566| Polynomial f ;
567| cin >> f ;
568| }
569| catch( PolynomialRangeError & e )
570| {
571| cout << "Error in inputting polynomial f(x)" << endl ;
572| }
573|
574+============================================================================*/
575
576istream & operator>>( istream & in, Polynomial & p )
577{
578 // Input as a string.
579 string s ;
580 in >> s ;
581
582 // Copy into argument polynomial. Can throw an exception.
583 p = Polynomial( s ) ;
584
585 return in ;
586}
587
588
589/*=============================================================================
590 |
591 | NAME
592 |
593 | Polynomial operator==
594 |
595 | DESCRIPTION
596 |
597 | Polynomial equality test operator.
598 |
599 | EXAMPLE
600 |
601 | try
602 | {
603 | Polynomial f1( "2x^2 + 1, 3" ) ;
604 | Polynomial f2( "2x^2 + 1", 3 ) ;
605 |
606 | if (f1 == f2)
607 | cout << "f1 = " << f1 << " == " << f2 << endl ;
608 | else
609 | cout << "f1 = " << f1 << " != " << f2 << endl ;
610 |
611 +============================================================================*/
612
613bool operator==( const Polynomial & p1, const Polynomial & p2 )
614{
615 // The degrees and moduli have to match.
616 if ((p1.n_ != p2.n_) || (p1.p_ != p2.p_))
617 return false ;
618
619 // Test coefficients for equality.
620 for (int i = 0 ; i <= p1.n_ ; ++i)
621 if (p1.f_[ i ] != p2.f_[ i ])
622 return false ;
623
624 return true ;
625}
626
627bool operator!=( const Polynomial & p1, const Polynomial & p2 )
628{
629 return !( p1 == p2) ;
630}
631
632
633/*=============================================================================
634|
635| NAME
636|
637| Polynomial operator[]
638|
639| DESCRIPTION
640|
641| Polynomial indexing operator which allows an lvalue: f[ 33 ] = 2 ;
642| If we don't have a coefficient of this degree, create it and backfill
643| earlier coefficients with zero.
644|
645| Throws an exception if out of bounds.
646|
647| EXAMPLE
648|
649| try
650| {
651| Polynomial f( "2x^2 + 1, 3" ) ;
652|
653| f[ 5 ] = 2 ;
654|
655| // Now f(x) = 2 x^5 + 0 x^4 + 0 x^3 + 2 x^2 + 0 x + 1
656| // f_.size() => 5 + 1 = 6
657|
658| }
659| catch( PolynomialRangeError & e )
660| {
661| cout << "Error in assigning polynomial f(x) coefficient" << endl ;
662| }
663|
664+============================================================================*/
665
666ppuint & Polynomial::operator[]( int i )
667{
668 // We attempt to access beyond the current degree.
669 if (i > n_)
670 {
671 try
672 {
673 // Extend the vector size with zeros.
674 f_.resize( i+1, 0 ) ;
675 n_ = i ;
676 }
677 // Failed to resize the polynomial.
678 catch( length_error & e )
679 {
680 throw PolynomialError( "Polynomial.operator[]: failed to resize", __FILE__, __LINE__ ) ;
681 }
682 }
683
684 // Return a reference to the coefficient.
685 return f_[ i ] ;
686}
687
688
689/*=============================================================================
690|
691| NAME
692|
693| Polynomial operator[]
694|
695| DESCRIPTION
696|
697| Polynomial indexing operator for read only access: int coeff = f[ 33 ] ;
698| Throws an exception if out of bounds.
699|
700| EXAMPLE
701|
702| try
703| {
704| Polynomial f ;
705| int value = f[ 33 ] ;
706| }
707| catch( PolynomialRangeError & e )
708| {
709| cout << "Error in getting polynomial f(x) coefficient" << endl ;
710| }
711|
712+============================================================================*/
713
714const ppuint Polynomial::operator[]( int i ) const
715{
716 // We throw our own exception here.
717 if (i > n_)
718 {
719 ostringstream os ;
720 os << "Error accessing polynomial with coefficients p[0]...p[n] = (" ;
721 for (int j = 0 ; j <= n_ ; ++j)
722 os << f_[ j ] << " " ;
723 os << ")" << endl
724 << " at index i = " << i
725 << " of degree n = " << n_ << " modulo p = " << p_ ;
726 throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
727 }
728
729 return f_[ i ] ;
730}
731
732
733/*=============================================================================
734|
735| NAME
736|
737| deg
738| modulus
739| setModulus
740|
741| DESCRIPTION
742|
743| Return the degree of f(x).
744|
745| EXAMPLE
746|
747| try
748| {
749| Polynomial f ;
750| cout << "Degree of f(x) = " << f.deg() << endl ;
751| }
752|
753+============================================================================*/
754
755int Polynomial::deg() const
756{
757 return n_ ;
758}
759
760// Return the modulus p of f(x).
761ppuint Polynomial::modulus() const
762{
763 return p_ ;
764}
765
766void Polynomial::setModulus( ppuint p )
767{
768 p_ = p ;
769
770 // And the modulus functionoid.
771 mod( p_ ) ;
772}
773
774
775/*=============================================================================
776|
777| NAME
778|
779| Polynomial operator+=
780|
781| DESCRIPTION
782|
783| Polynomial sum f(x) += g(x)
784|
785|
786| EXAMPLE
787|
788| try
789| {
790| Polynomial f ;
791| Polynomial g ;
792| f += g ;
793| }
794| catch( PolynomialRangeError & e )
795| {
796| cout << "Error in polynomial sum f(x) += g(x)" << endl ;
797| }
798|
799+============================================================================*/
800
801Polynomial & Polynomial::operator+=( const Polynomial & g )
802{
803 // f(x) = x^2 + + 1
804 // g(x) = x + 3
805 //
806 // f(x) = x + 3
807 // g(x) = x^2 + + 1
808 //
809 int minDeg = (n_ < g.n_) ? n_ : g.n_ ;
810
811 // Add coefficients modulo p for smaller degree terms.
812 for (int i = 0 ; i <= minDeg ; ++i)
813 f_[ i ] = mod( f_[ i ] + g.f_[ i ] ) ;
814
815 // If g(x) has larger degree, extend f(x) and copy the coefficients of g(x).
816 if (g.n_ > n_)
817 {
818 // Extend the vector size with zeros.
819 try
820 {
821 f_.resize( g.n_ + 1, 0 ) ;
822 }
823 // Failed to resize the polynomial.
824 catch( length_error & e )
825 {
826 throw PolynomialError( "Polynomial::operator+= had a length_error exception while resizing the polynomial", __FILE__, __LINE__ ) ;
827 }
828
829 for (int i = n_ + 1 ; i <= g.n_ ; ++i)
830 f_[ i ] = g.f_[ i ] ;
831 }
832
833 // Trim leading zero coefficients, but leave a constant term of zero.
834 while( f_.back() == 0 && f_.size() > 1)
835 {
836 f_.pop_back() ;
837 --n_ ;
838 }
839
840 // Return current object now containing the sum.
841 return *this ;
842}
843
844
845/*=============================================================================
846|
847| NAME
848|
849| Polynomial operator+()
850|
851| DESCRIPTION
852|
853| Add polynomials.
854|
855+============================================================================*/
856
857const Polynomial operator+( const Polynomial & f, const Polynomial &g )
858{
859 // Do + in terms of += to maintain consistency.
860 // Copy construct temporary copy, then add to it.
861 // Return value optimization compiles away the copy constructor.
862 // const on return type disallows doing (u+v) = w ;
863 return Polynomial( f ) += g ;
864}
865
866/*=============================================================================
867|
868| NAME
869|
870| Polynomial operator*=()
871|
872| DESCRIPTION
873|
874| Scalar multiply polynomials.
875|
876+============================================================================*/
877
878Polynomial & Polynomial::operator*=( const ppuint k )
879{
880 // Multiply coefficients modulo p.
881 for (int i = 0 ; i <= n_ ; ++i)
882 f_[ i ] = mod( f_[ i ] * k ) ;
883
884 // Return current object now containing the scalar product.
885 return *this ;
886}
887
888
889/*=============================================================================
890|
891| NAME
892|
893| Polynomial operator*()
894|
895| DESCRIPTION
896|
897| Scalar multiply polynomials.
898|
899+============================================================================*/
900
901const Polynomial operator*( const Polynomial & f, const ppuint k )
902{
903 // Do * in terms of *= to maintain consistency.
904 // Copy construct temporary copy, then add to it.
905 // Return value optimization compiles away the copy constructor.
906 // const on return type disallows doing (u*k) = w ;
907 return Polynomial( f ) *= k ;
908}
909
910
911/*=============================================================================
912|
913| NAME
914|
915| Polynomial operator()
916|
917| DESCRIPTION
918|
919| Evaluate the monic polynomial f( x ) with modulo p arithmetic.
920|
921| n n-1
922| f( x ) = x + a x + ... + a 0 <= a < p
923| n-1 0 i
924|
925| EXAMPLE
926| 4
927| Let n = 4, p = 5 and f(x) = x + 3x + 3.
928|
929| By Horner's rule, f(x) = (((x + 0)x + 0)x + 3)x + 3.
930|
931| Then f(2) = (((2 + 0)2 + 0)2 + 3) = (8 + 3)2 + 3 = 1 + 2 + 3 (mod 5) = 0.
932| and f(3) = (((3 + 0)3 + 0)3 + 3)3 + 3 (mod 5) = 3
933|
934| METHOD
935|
936| By Horner's rule, f(x) = ( ... ( (x + a )x + ... )x + a .
937| n-1 0
938|
939| We evaluate recursively, f := f * x + a (mod p), starting
940| i
941| with f = 1 and i = n-1.
942|
943+============================================================================*/
944
945ppuint
946Polynomial::operator()( int x )
947{
948 ppuint val = 1 ;
949
950 for (int degree = n_- 1 ; degree >= 0 ; --degree)
951 val = mod( val * x + f_[ degree ]) ;
952
953 return( val ) ;
954}
955
956
957/*=============================================================================
958|
959| NAME
960|
961| hasLinearFactor
962|
963| DESCRIPTION
964|
965| Check if the polynomial f(x) has any linear factors.
966|
967| Polynomial f ; // A polynomial f(x) of degree n, modulo p.
968| bool hasFactor = f.hasLinearFactor() ;
969|
970| hasFactor is true if f( a ) = 0 (mod p) for a = 1, 2, ... p-1,
971| and is false otherwise.
972|
973| i.e. check if f(x) contains a linear factor (x - a). We don't need to test
974| for the root a = 0 because f(x) was chosen in main to have a non-zero
975| constant term, hence no zero root.
976|
977| EXAMPLE
978| 4 2 2
979| Let n = 4, p = 5 and f(x) = x + 3x + 3. Then f(x) = (x + 3) (x + 4x + 2)
980|
981| Then f(0) = 3 (mod 5), f(1) = 2 (mod 5), but
982| f(2) = 0 (mod 5), so we exit immediately with a true.
983|
984| 4 2
985| However, f(x) = x + 3 x + x + 1 is irreducible, so has no linear factors.
986|
987| METHOD
988|
989| Evaluate f(x) at x = 0, ..., p-1 by Horner's rule. Return instantly the
990| moment f(x) evaluates to 0.
991|
992+============================================================================*/
993
994bool
995Polynomial::hasLinearFactor()
996{
997 for (int i = 0 ; i <= p_ - 1 ; ++i)
998 if ((*this)( i ) == 0)
999 return( true ) ;
1000
1001 return( false ) ;
1002}
1003
1004
1005/*=============================================================================
1006|
1007| NAME
1008|
1009| Polynomial::isInteger
1010|
1011| DESCRIPTION
1012|
1013| Return true if a polynomial is a constant.
1014|
1015| EXAMPLE
1016|
1017| Polynomial p( "2 x ^ 2 " ) ;
1018| p.isInteger -> false
1019|
1020| Polynomial p( "2 " ) ;
1021| p.isInteger -> true
1022|
1023| METHOD
1024|
1025| A constant polynomial is zero in its first through n th degree
1026| terms. Return immediately with false if any such term is non-zero.
1027|
1028+============================================================================*/
1029
1030bool
1031Polynomial::isInteger() const
1032{
1033 // Degree 0 is constant.
1034 if (n_ == 0)
1035 return true ;
1036
1037 // Not integer if any coefficients above zero degree term are non-zero.
1038 for (int i = 1 ; i <= n_ ; ++i)
1039 if (f_[ i ] != 0)
1040 return( false ) ;
1041
1042 return( true ) ;
1043}
1044
1045
1046/*=============================================================================
1047 |
1048 | NAME
1049 |
1050 | initialTrialPoly
1051 |
1052 | DESCRIPTION
1053 |
1054 | Create an initial monic polynomial
1055 | n
1056 | f( x ) = x
1057 |
1058 | EXAMPLE
1059 | 4
1060 | Let n = 4. Set f(x) = x - 1.
1061 |
1062 |
1063 +============================================================================*/
1064
1065void TrialPolynomial::initialTrialPoly( const ppuint n, const ppuint p )
1066{
1067 (*this).setModulus(p);
1068
1069 // Allocate enough coefficients for an nth degree polynomial and
1070 // initialize all intermediate coefficients to 0.
1071
1072 if (n > numeric_limits<int>::max())
1073 {
1074 ostringstream os ;
1075 os << "Polynomial::initialTrialPoly: n = " << n << " is too large for an array index" ;
1076 throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
1077 }
1078
1079 (*this)[ static_cast<int>( n ) ] = 1 ;
1080 f_[ 0 ] = 0 ;
1081}
1082
1083
1084/*=============================================================================
1085 |
1086 | NAME
1087 |
1088 | nextTrialPoly
1089 |
1090 | DESCRIPTION
1091 |
1092 | Return the next monic polynomial in the sequence after f(x), explained
1093 | below
1094 |
1095 | EXAMPLE
1096 | 3
1097 | Let n = 3 and p = 5. Suppose f(x) = x + 4 x + 4. As a mod p number,
1098 | this is 1 0 4 4. Adding 1 gives 1 0 4 5. We reduce modulo
1099 | 5 and propagate the carry to get the number 1 0 5 0. Propagating
1100 | the carry again and reducing gives 1 1 0 0. The next polynomial
1101 | 3 2
1102 | after f(x) is x + x .
1103 |
1104 | METHOD
1105 |
1106 | Think of the polynomial coefficients as the digits of a number written
1107 | in base p. The "next" polynomial is the one you would get by adding 1
1108 | to this number in multiple precision arithmetic. Our intention is to
1109 | run through all possible monic polynomials modulo p.
1110 |
1111 | Propagate carries in digits 1 through n-2 when any digit exceeds p. No
1112 | carries take place in the n-1 st digit because our polynomial is monic.
1113 |
1114 | TODO: Find polynomials in order of Hamming weight?
1115 |
1116 +============================================================================*/
1117
1118void TrialPolynomial::nextTrialPoly()
1119{
1120 ++f_[ 0 ] ; // Add 1, i.e. increment the coefficient of the x term.
1121
1122 // Sweep through the number from right to left, propagating carries. Skip
1123 // the constant and the nth degree terms.
1124 for (int digit_num = 0 ; digit_num <= n_ - 2 ; ++digit_num)
1125 {
1126 if (f_[ digit_num ] == p_) // Propagate carry to next digit.
1127 {
1128 f_[ digit_num ] = 0 ;
1129 ++f_[ digit_num + 1 ] ;
1130 }
1131 }
1132}
1133
1134
1135
1136/*------------------------------------------------------------------------------
1137| PolyMod Implementation |
1138------------------------------------------------------------------------------*/
1139
1140/*=============================================================================
1141 |
1142 | NAME
1143 |
1144 | PolyMod default constructor
1145 |
1146 | DESCRIPTION
1147 |
1148 |
1149 | EXAMPLE
1150 |
1151 | METHOD
1152 |
1153 +============================================================================*/
1154
1155PolyMod::PolyMod()
1156 : g_()
1157 , f_()
1158 , powerTable_()
1159 , mod( f_.modulus() )
1160{
1161 constructPowerTable() ;
1162 modf() ;
1163}
1164
1165
1166/*=============================================================================
1167 |
1168 | NAME
1169 |
1170 | PolyMod destructor
1171 |
1172 | DESCRIPTION
1173 |
1174 |
1175 | EXAMPLE
1176 |
1177 | METHOD
1178 |
1179 +============================================================================*/
1180
1181PolyMod::~PolyMod()
1182{
1183// Member fields will clean up themselves.
1184}
1185
1186
1187/*=============================================================================
1188 |
1189 | NAME
1190 |
1191 | PolyMod constructor
1192 |
1193 | DESCRIPTION
1194 |
1195 | Given polynomials f( x ) and g( x ) where g is a string,
1196 | construct p( x ) = g( x ) mod f( x ).
1197 |
1198 | EXAMPLE
1199 |
1200 | METHOD
1201 |
1202 +============================================================================*/
1203
1204PolyMod::PolyMod( const string & g, const Polynomial & f )
1205 : g_( g )
1206 , f_( f )
1207 , powerTable_()
1208 , mod( f.modulus() )
1209{
1210 constructPowerTable() ;
1211 modf() ;
1212}
1213
1214
1215/*=============================================================================
1216 |
1217 | NAME
1218 |
1219 | PolyMod constructor
1220 |
1221 | DESCRIPTION
1222 |
1223 | Given polynomials f( x ) and g( x ), construct p( x ) = g( x ) mod f( x ).
1224 |
1225 | EXAMPLE
1226 |
1227 | METHOD
1228 |
1229 +============================================================================*/
1230
1231PolyMod::PolyMod( const Polynomial & g, const Polynomial & f )
1232 : g_( g )
1233 , f_( f )
1234 , powerTable_()
1235 , mod( f.modulus() )
1236{
1237 constructPowerTable() ;
1238 modf() ;
1239}
1240
1241
1242/*=============================================================================
1243 |
1244 | NAME
1245 |
1246 | PolyMod string operator
1247 |
1248 | DESCRIPTION
1249 |
1250 | Given g( x ) mod f( x ), return g( x ) as a string.
1251 |
1252 | EXAMPLE
1253 |
1254 | METHOD
1255 |
1256 +============================================================================*/
1257
1258// Operator casting to string type.
1259PolyMod::operator string() const
1260{
1261return static_cast<string>( g_ ) ;
1262}
1263
1264
1265/*=============================================================================
1266 |
1267 | NAME
1268 |
1269 | Operator << for PolyMod
1270 |
1271 | DESCRIPTION
1272 |
1273 | Given g( x ) mod f( x ), output g( x ) as a string.
1274 |
1275 | EXAMPLE
1276 |
1277 | METHOD
1278 |
1279 +============================================================================*/
1280
1281ostream & operator<<( ostream & out, const PolyMod & p )
1282{
1283 // Cast to polynomial to a string first, then output as a string.
1284 // May throw a PolynomialRangeError.
1285 out << static_cast<string>( p.g_ ) ;
1286
1287 return out ;
1288}
1289
1290
1291/*=============================================================================
1292 |
1293 | NAME
1294 |
1295 | getf
1296 |
1297 | DESCRIPTION
1298 |
1299 | Given g( x ) mod f( x ), return f( x ).
1300 |
1301 | EXAMPLE
1302 |
1303 | METHOD
1304 |
1305 +============================================================================*/
1306
1307const Polynomial PolyMod::getf() const
1308{
1309 return f_ ;
1310}
1311
1312
1313/*=============================================================================
1314 |
1315 | NAME
1316 |
1317 | getModulus
1318 |
1319 | DESCRIPTION
1320 |
1321 | Given g( x ) mod (f( x ), p) return p.
1322 |
1323 | EXAMPLE
1324 |
1325 | METHOD
1326 |
1327 +============================================================================*/
1328
1329const ppuint PolyMod::getModulus() const
1330{
1331 return f_.modulus() ;
1332}
1333
1334
1335/*=============================================================================
1336 |
1337 | NAME
1338 |
1339 | PolyMod copy constructor
1340 |
1341 | DESCRIPTION
1342 |
1343 | Copy g2 to g( x ) mod (f( x ), p)
1344 |
1345 | EXAMPLE
1346 |
1347 | METHOD
1348 |
1349 +============================================================================*/
1350
1351PolyMod::PolyMod( const PolyMod & g2 )
1352 : g_( g2.g_ )
1353 , f_( g2.f_ )
1354 , powerTable_( g2.powerTable_ )
1355 , mod( f_.modulus() )
1356{
1357}
1358
1359
1360/*=============================================================================
1361 |
1362 | NAME
1363 |
1364 | operator=
1365 |
1366 | DESCRIPTION
1367 |
1368 | PolyMod assignment operator.
1369 |
1370 +============================================================================*/
1371
1372PolyMod & PolyMod::operator=( const PolyMod & g2 )
1373{
1374 // Check for assigning to oneself: just pass back a reference to the unchanged object.
1375 if (this == &g2)
1376 return *this ;
1377
1378 g_ = g2.g_ ;
1379 g_ = g2.f_ ;
1380
1381 powerTable_ = g2.powerTable_ ;
1382 mod = g2.mod ;
1383
1384 // Return a reference to the altered object.
1385 return *this ;
1386}
1387
1388
1389/*=============================================================================
1390 |
1391 | NAME
1392 |
1393 | operator[]
1394 |
1395 | DESCRIPTION
1396 |
1397 | Bounds checked indexing operator for read only access:
1398 | coeff = p[ i ] ;
1399 |
1400 +============================================================================*/
1401
1402const ppuint PolyMod::operator[]( int i ) const
1403{
1404 // Can throw an exception.
1405 return g_[ i ] ;
1406}
1407
1408
1409/*=============================================================================
1410|
1411| NAME
1412|
1413| constructPowerTable
1414|
1415| DESCRIPTION
1416|
1417| Construct a table of powers of x:
1418|
1419| n 2n-2
1420| x (mod f(x), p) ... x (mod f(x), p)
1421|
1422|
1423| powerTable_[i][j] is the coefficient of
1424| j n+i
1425| x in x (mod f(x), p) where 0 <= i <= n-2 and 0 <= j <= n-1.
1426|
1427| EXAMPLE
1428| 4 2 4
1429| Let n = 4, p = 5 and f(x) = x + x + 2x + 3. Then x =
1430|
1431| 2 2
1432| -( x + 2x + 3) = 4 x + 3 x + 2 (mod f(x), 5), and we get
1433|
1434| 4 2
1435| x (mod f(x), 5) = 4 x + 3 x + 2 = powerTable_[0].
1436|
1437| 5 2 3 2
1438| x (mod f(x), 5) = x( 4 x + 3 x + 2) = 4 x + 3 x + 2x
1439| = powerTable_[1].
1440|
1441| 6 3 2 4 3 2
1442| x (mod f(x), 5) = x( 4 x + 3 x + 2 x) = 4x + 3 x + 2 x
1443|
1444| 2 3 2
1445| = 4 ( 4x + 3 x + 2) + 3 x + 2 x =
1446|
1447| 3 2
1448| = 3 x + 3 x + 2 x + 3 = powerTable_[2].
1449|
1450| j
1451| powerTable_[i][j]: | 0 1 2 3
1452| ---+-------------
1453| 0 | 2 3 4 0
1454| i 1 | 0 2 3 4
1455| 2 | 3 2 3 3
1456|
1457| NOTES
1458| n-1
1459| Beginning with t( x ) = x, compute the next power of x from the last
1460| n
1461| one by multiplying by x. If necessary, substitute x =
1462| n-1
1463| -( a x + ... + a ) to reduce the degree. This formula comes from
1464| n-1 0
1465| n n-1
1466| the observation f( x ) = x + a x + ... + a = 0 (mod f(x), p).
1467| n-1 0
1468|
1469+============================================================================*/
1470
1471void PolyMod::constructPowerTable()
1472{
1473 // Get hold of the degree of f(x).
1474 int n = f_.deg() ;
1475
1476 // Empty the power table.
1477 powerTable_.clear() ;
1478
1479 //
1480 // t(x) is temporary storage for x ^ k (mod f(x),p)
1481 // n <= k <= 2n-2. Its degree can go as high as
1482 // n before it is reduced again.
1483 Polynomial t ;
1484
1485
1486 // n-1
1487 // Initialize t( x ) = x mod p.
1488 t[ n-1 ] = 1 ;
1489
1490 // In Microsoft Visual Studio C++ 2008 we get garbage placed in t[ n ] in the loop
1491 // at j = n in the step
1492 // t[ j ] = t[ j-1 ] ;
1493 // Why? We first access the value of t[j-1], the compiler places it in a temporary,
1494 // we then access t[n], and this causes a resize of f_ in Polynomial.operator[],
1495 // then t[ j ] = garbage since we apparently lose the temporary. Does not happen if
1496 // we rewrite the step as
1497 // int tmp ;
1498 // tmp = t[ j-1 ] ;
1499 // t[ j ] = tmp ;
1500 // or alternatively, we prevent resizing occurring by pre-expanding:
1501 t[ n ] = 0 ; // Expand the size of t(x) now since we'll access t[n] later.
1502
1503 t.setModulus( getModulus() ) ;
1504
1505 try
1506 {
1507 // i+n
1508 // Fill the ith row of the table with x (mod f(x), p)
1509 // for i = 0 ... n-2.
1510 //
1511 for (int i = 0 ; i <= n - 2 ; ++i)
1512 {
1513 // Compute t(x) = x t(x) by shifting the coefficients
1514 // to the left and filling with zero.
1515 for (int j = n ; j >= 1 ; --j)
1516 t[ j ] = t[ j-1 ] ;
1517
1518 t[ 0 ] = 0 ;
1519
1520 // Coefficient of the x ^ n degree term of t(x).
1521 ppsint coeff = 0 ;
1522 if ( (coeff = t[ n ]) != 0)
1523 {
1524 // Zero out the x ^ n th term.
1525 t[ n ] = 0 ;
1526
1527 // n n n-1
1528 // Replace x with x (mod f(x), p) = -(a x + ... + a )
1529 // n-1 0
1530 for (int j = 0 ; j <= n-1 ; ++j)
1531
1532 t[ j ] = mod( t[ j ] +
1533 mod( -coeff * f_[ j ]) ) ;
1534 } // end if
1535
1536 // Copy t(x) into row i of power_table.
1537 powerTable_.push_back( t ) ;
1538
1539 } // end for
1540
1541 #ifdef DEBUG_PP_POLYNOMIAL
1542 cout << "PowerTable of polynomials x^n ... x^2n-2 mod f(x), p" << endl ;
1543 cout << "f(x) = " << getf() << " n = " << n << " p = " << getModulus() << endl ;
1544 for (int i = n ; i <= 2*n-2 ; ++i)
1545 cout << "powerTable[ x^" << i << " ] = " << powerTable_[ offset(i) ] << endl ;
1546 #endif
1547 }
1548 catch( bad_alloc & e )
1549 {
1550 throw PolynomialRangeError( "PolyMod::constructPowerTable had a bad_alloc memory allocation failure", __FILE__, __LINE__ ) ;
1551 }
1552
1553 // t will be automagically freed upon exit.
1554 return ;
1555}
1556
1557
1558/*=============================================================================
1559|
1560| NAME
1561|
1562| modf
1563|
1564| DESCRIPTION
1565|
1566| Reduce g(x) modulo f(x), and p.
1567|
1568| EXAMPLE
1569| 4 2
1570| Let n = 4, p = 5 and f(x) = x + x + 2x + 3.
1571|
1572| 6 3 2 4 3 2
1573| x (mod f(x), 5) = x( 4 x + 3 x + 2 x) = 4x + 3 x + 2 x
1574|
1575|
1576+============================================================================*/
1577
1578void PolyMod::modf()
1579{
1580 // Get hold of the degree of f(x).
1581 int n = f_.deg() ;
1582 int m = g_.deg() ;
1583
1584 if (m > 2 * n - 2)
1585 {
1586 ostringstream os ;
1587 os << "Error in PolyMod::modf(): degree of g(x) higher than power table can handle with deg f = " << n
1588 << " deg g = " << m << " and p = " << getModulus() ;
1589 throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
1590 }
1591
1592
1593 // i+n
1594 // Fill the ith row of the table with x (mod f(x), p)
1595 // for i = 0 ... n-2.
1596 //
1597 for (int i = n ; i <= m ; ++i)
1598 {
1599 #ifdef DEBUG_PP_POLYNOMIAL
1600 cout << "\nBefore converting, g( x ) = " << g_ << endl ;
1601 #endif
1602
1603 // Coefficient of the x ^ i degree term of g(x).
1604 ppuint coeff{ 0 } ;
1605 if ( (coeff = g_[ i ]) != 0)
1606 {
1607 // Subtract (zero out) the x ^ i th term.
1608 g_[ i ] = 0 ;
1609
1610 // i i
1611 // Replace x with x (mod f(x), p) from the power table * coeff.
1612 g_ += (powerTable_[ offset(i) ] * coeff) ;
1613 }
1614
1615 #ifdef DEBUG_PP_POLYNOMIAL
1616 cout << "\nAfter converting with coeff = " << coeff << " g( x ) = " << g_ << endl ;
1617 #endif
1618
1619 } // end for
1620
1621 return ;
1622}
1623
1624
1625/*=============================================================================
1626 |
1627 | NAME
1628 |
1629 | autoconvolve
1630 |
1631 | DESCRIPTION
1632 |
1633 | Compute a convolution-like sum for use in function coeffOfSquare,
1634 |
1635 | upper
1636 | ---
1637 | \ t t But define the sum to be 0 when lower > upper to catch
1638 | / i k-i the special cases that come up in function coeffOfSquare.
1639 | ---
1640 | i=lower
1641 |
1642 | where
1643 | n-1
1644 | Coefficients of t(x) = t x + ... + t x + t
1645 | n-1 1 0
1646 |
1647 | EXAMPLE
1648 | 3 2
1649 | Suppose t(x) = 4 x + x + 3 x + 3, lower = 1, upper = 3, n = 3,
1650 |
1651 | and p = 5. For k = 3, autoConvolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +
1652 |
1653 | t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3. For lower = 0,
1654 |
1655 | upper = -1, or for lower = 3 and upper = 2, autoConvolve = 0, no matter what
1656 |
1657 | k is.
1658 |
1659 | METHOD
1660 |
1661 | A "for" loop in the C language is not executed when its lower limit
1662 |
1663 | exceeds its upper limit, leaving sum = 0.
1664 |
1665 +============================================================================*/
1666
1667ppuint autoConvolve( const Polynomial & t, int k, int lower, int upper )
1668{
1669 ModP<ppuint,ppsint> mod( t.modulus() ) ;
1670 int deg_t = t.deg() ;
1671
1672 ppuint sum { 0 } ;
1673 for (int i = lower ; i <= upper ; ++i)
1674 {
1675 // Coeff is zero if higher or lower than degree of polynomial.
1676 ppuint coeff_ti{ 0u } ;
1677 ppuint coeff_tkmi{ 0u } ;
1678
1679 if (i <= deg_t && i >= 0)
1680 coeff_ti = t[ i ] ;
1681
1682 if (k-i <= deg_t && k-i >= 0)
1683 coeff_tkmi = t[ k - i ] ;
1684
1685 sum = mod( sum + mod( coeff_ti * coeff_tkmi )) ;
1686 }
1687
1688 return( sum ) ;
1689}
1690
1691
1692/*=============================================================================
1693 |
1694 | NAME
1695 |
1696 | convolve
1697 |
1698 | DESCRIPTION
1699 |
1700 | Compute a convolution-like sum,
1701 |
1702 | upper
1703 | ---
1704 | \ s t But define the sum to be 0 when lower > upper to catch
1705 | / i k-i the special cases
1706 | ---
1707 | i=lower
1708 |
1709 | where
1710 | n-1
1711 | Coefficients of s(x) = s x + ... + s x + s
1712 | n-1 1 0
1713 | n-1
1714 | Coefficients of t(x) = t x + ... + t x + t
1715 | n-1 1 0
1716 |
1717 | 0 <= k <= 2n - 2
1718 | 0 <= lower <= n-1
1719 | 0 <= upper <= n-1
1720 |
1721 | EXAMPLE
1722 | 3 2
1723 | Suppose s(x) = 4 x + x + 3 x + 3,
1724 |
1725 | 3 2
1726 | Suppose t(x) = 4 x + x + 3 x + 3,
1727 |
1728 |
1729 | lower = 1, upper = 3, n = 3,
1730 |
1731 | and p = 5. For k = 3, convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +
1732 |
1733 | t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3. For lower = 0,
1734 |
1735 | upper = -1, or for lower = 3 and upper = 2, convolve = 0, no matter what
1736 |
1737 | k is.
1738 |
1739 | METHOD
1740 |
1741 | A "for" loop in the C language is not executed when its lower limit
1742 |
1743 | exceeds its upper limit, leaving sum = 0.
1744 |
1745 +============================================================================*/
1746
1747ppuint convolve( const Polynomial & s, const Polynomial & t,
1748 const int k, const int lower, const int upper )
1749{
1750 ppuint sum{ 0 } ; // Convolution sum.
1751
1752 ModP<ppuint,ppsint> mod( s.modulus() ) ; // Initialize the functionoid.
1753
1754 int deg_s{ s.deg() } ;
1755 int deg_t{ t.deg() } ;
1756
1757 for (int i = lower ; i <= upper ; ++i)
1758 {
1759 // Coeff is zero if higher or lower than degree of polynomial.
1760 ppuint coeff_s{ 0u } ;
1761 ppuint coeff_t{ 0u } ;
1762
1763 if (i <= deg_s && i >= 0)
1764 coeff_s = s[ i ] ;
1765
1766 if (k-i <= deg_t && k-i >= 0)
1767 coeff_t = t[ k- i ] ;
1768
1769 sum = mod( sum + mod( coeff_s * coeff_t )) ;
1770 }
1771
1772 return( sum ) ;
1773}
1774
1775
1776/*=============================================================================
1777 |
1778 | NAME
1779 |
1780 | coeffOfSquare
1781 |
1782 | DESCRIPTION
1783 | th 2
1784 | Return the coefficient of the k power of x in g ( x ) modulo p,
1785 | given of g(x) of degree <= n-1.
1786 |
1787 | where 0 <= k <= 2n-2
1788 |
1789 | EXAMPLE
1790 | 3 2 2
1791 | Let n = 4, p = 5, and g(x) = 4 x + x + 3 x + 3. g(x) =
1792 |
1793 | 6 5
1794 | x + 3 x + 3 x + 3 x + 4, all modulo 5.
1795 |
1796 | k | 0 1 2 3 4 5 6
1797 | ----------------+---------------------
1798 | coeffOfSquare | 4 3 0 0 0 3 1
1799 |
1800 | METHOD
1801 | 2
1802 | The formulas were gotten by writing out the product g (x) explicitly.
1803 |
1804 | The sum is 0 in two cases:
1805 |
1806 | (1) when k = 0 and the limits of summation are 0 to -1
1807 |
1808 | (2) k = 2n - 2, when the limits of summation are n to n-1.
1809 |
1810 | To derive the formulas, let
1811 |
1812 | n-1
1813 | Let g(x) = g x + ... + g x + g
1814 | n-1 1 0
1815 |
1816 | Look at the formulas in coeffOfProduct for each power of x,
1817 | replacing s with t, and observe that half of the terms are
1818 | duplicates, so we can save computation time.
1819 |
1820 | Inspection yields the formulas,
1821 |
1822 | for 0 <= k <= n-1, even k,
1823 |
1824 | k/2-1
1825 | --- 2
1826 | 2 \ g g + g
1827 | / i k-i k/2
1828 | ---
1829 | i=0
1830 |
1831 | for 0 <= k <= n-1, odd k,
1832 |
1833 | (k-1)/2
1834 | ---
1835 | 2 \ g g
1836 | / i k-i
1837 | ---
1838 | i=0
1839 |
1840 | and for n <= k <= 2n-2, even k,
1841 |
1842 | n-1
1843 | --- 2
1844 | 2 \ g g + g
1845 | / i k-i k/2
1846 | ---
1847 | i=k/2+1
1848 |
1849 | and for n <= k <= 2n-2, odd k,
1850 |
1851 | n-1
1852 | ---
1853 | 2 \ g g
1854 | / i k-i
1855 | ---
1856 | i=(k+1)/2
1857 |
1858 +============================================================================*/
1859
1860ppuint coeffOfSquare( const Polynomial & g, const int k, const int n )
1861{
1862 ModP<ppuint,ppsint> mod( g.modulus() ) ; // Initialize the functionoid.
1863
1864 // 2
1865 ppuint sum { 0 } ; // kth coefficient of g( x )
1866
1867 // Coeff is zero if higher or lower than degree of polynomial.
1868 ppuint coeff_gkd2 { 0 } ;
1869 if (k/2 <= g.deg() && k/2 >= 0)
1870 coeff_gkd2 = g[ k/2 ] * g[ k/2 ] ;
1871
1872 if (0 <= k && k <= n-1)
1873 {
1874 if (k % 2 == 0) // Even k
1875 sum = mod( mod( 2 * autoConvolve( g, k, 0, k/2 - 1) ) + coeff_gkd2 ) ;
1876
1877 else // Odd k
1878 sum = mod( 2 * autoConvolve( g, k, 0, (k-1)/2)) ;
1879 }
1880 else if (n <= k && k <= 2 * n - 2)
1881 {
1882
1883 if (k % 2 == 0) // Even k
1884 sum = mod( mod( 2 * autoConvolve( g, k, k/2 + 1, n-1)) + coeff_gkd2 ) ;
1885
1886 else // Odd k
1887 sum = mod( 2 * autoConvolve( g, k, (k+1)/2, n-1)) ;
1888 }
1889
1890 return( sum ) ;
1891}
1892
1893
1894/*=============================================================================
1895 |
1896 | NAME
1897 |
1898 | coeffOfProduct
1899 |
1900 | DESCRIPTION
1901 | th
1902 | Return the coefficient of the k power of x in s( x ) t( x ) modulo p.
1903 |
1904 | EXAMPLE
1905 | 3 2 2
1906 | Let n = 4, p = 5, t(x) = 4 x + x + 4, s( x ) = 3 x + x + 2
1907 |
1908 | 5 4 3 2
1909 | then s ( x ) t( x ) = 2 x + 2 x + 4 x + 4 x + 4 x + 3
1910 |
1911 | We'll do the case k=3,
1912 |
1913 | t3 s0 + t2 s1 + t1 s2 + t0 s3 = 4 * 2 + 1 * 1 + 0 * 3 + 4 * 0 = 9 = 4 (mod 5).
1914 |
1915 | k | 0 1 2 3 4 5 6
1916 | -----------------+---------------------
1917 | coeffOfProduct | 3 4 4 4 2 2 0
1918 |
1919 | METHOD
1920 |
1921 | The formulas were gotten by writing out the product s(x) t (x) explicitly.
1922 |
1923 | The sum is 0 in two cases:
1924 |
1925 | (1) when k = 0 and the limits of summation are 0 to -1
1926 |
1927 | (2) k = 2n - 2, when the limits of summation are n to n-1.
1928 |
1929 |
1930 | To derive the formulas, let
1931 |
1932 | n-1
1933 | Let s (x) = s x + ... + s x + s and
1934 | n-1 1 0
1935 |
1936 | n-1
1937 | t (x) = t x + ... + t x + t
1938 | n-1 1 0
1939 |
1940 | and multiply out the terms, collecting like powers of x:
1941 |
1942 |
1943 | Power of x Coefficient
1944 | ==========================
1945 | 2n-2
1946 | x s t
1947 | n-1 n-1
1948 |
1949 | 2n-3
1950 | x s t + s t
1951 | n-2 n-1 n-1 n-2
1952 |
1953 | 2n-4
1954 | x s t + s t + s t
1955 | n-3 n-1 n-2 n-2 n-3 n-1
1956 |
1957 | 2n-5
1958 | x s t + s t + s t + s t
1959 | n-4 n-1 n-3 n-2 n-2 n-3 n-1 n-4
1960 |
1961 | . . .
1962 |
1963 | n
1964 | x s t + s t + ... + s t
1965 | 1 n-1 2 n-2 n-1 1
1966 |
1967 | n-1
1968 | x s t + s t + ... + s t
1969 | 0 n-1 1 n-2 n-1 0
1970 |
1971 | . . .
1972 |
1973 | 3
1974 | x s t + s t + s t + s t
1975 | 0 3 1 2 2 1 3 0
1976 |
1977 | 2
1978 | x s t + s t + s t
1979 | 0 2 1 1 2 0
1980 |
1981 |
1982 | x s t + s t
1983 | 0 1 1 0
1984 |
1985 | 1 s t
1986 | 0 0
1987 |
1988 |
1989 | Inspection yields the formulas,
1990 |
1991 |
1992 | for 0 <= k <= n-1,
1993 |
1994 | k
1995 | ---
1996 | \ s t
1997 | / i k-i
1998 | ---
1999 | i=0
2000 |
2001 |
2002 | and for n <= k <= 2n-2,
2003 |
2004 | n-1
2005 | ---
2006 | \ s t
2007 | / i k-i
2008 | ---
2009 | i=k-n+1
2010 |
2011 +============================================================================*/
2012
2013ppuint coeffOfProduct( const Polynomial & s, const Polynomial & t, const int k, const int n )
2014{
2015 // Check if p is the same for s and t, and check the degree of s and t are < n.
2016 if (s.modulus() != t.modulus() || s.deg()> n || t.deg() > n)
2017 throw PolynomialRangeError( "coeffOfProduct: degree or modulus doesn't agree for polynomials s and t",
2018 __FILE__, __LINE__ ) ;
2019
2020 ppuint sum { 0 } ; // kth coefficient of t(x) ^ 2.
2021
2022 if (0 <= k && k <= n-1)
2023 {
2024 sum = convolve( s, t, k, 0, k ) ;
2025 }
2026 else if (n <= k && k <= 2 * n - 2)
2027 {
2028 sum = convolve( s, t, k, k - n + 1, n - 1 ) ;
2029 }
2030
2031 return( sum ) ;
2032}
2033
2034
2035/*=============================================================================
2036 |
2037 | NAME
2038 |
2039 | *
2040 |
2041 | DESCRIPTION
2042 |
2043 | Compute s( x ) t( x ) (mod f(x), p)
2044 |
2045 | s(x), of degree <= n-1.
2046 | t(x), of degree <= n-1.
2047 |
2048 | Uses a precomputed table of powers of x,
2049 | powerTable contains x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.
2050 |
2051 | EXAMPLE
2052 | 3 2 2
2053 | Let n = 4 and p = 5, t( x ) = 4 x + x + 4, s( x ) = 3 x + x + 2
2054 |
2055 | 5 4 3 2
2056 | Then s( x ) t( x ) = 2 x + 2 x + 4 x + 4 x + 4 x + 3, modulo 5,
2057 |
2058 | 4 2
2059 | and after reduction modulo f( x ) = x + x + 2 x + 3, using the power
2060 |
2061 | 4 2 5 3 2
2062 | table entries for x = 4 x + 3 x + 2 and x = 4 x + 3 x + 2 x, we get
2063 |
2064 | 3 2
2065 | s( x ) t( x ) (mod f( x ), p) = 2 x + 3 x + 4 x + 2
2066 |
2067 |
2068 | METHOD
2069 |
2070 | Compute the coefficients using the function coeffOfProduct.
2071 |
2072 | The next step is to reduce s(x) t(x) modulo f(x) and p. To do so, replace
2073 |
2074 | k k
2075 | each non-zero term t x, n <= k <= 2n-2, by the term t * [ x mod f(x), p)]
2076 | k k
2077 |
2078 | which we get from the array powerTable.
2079 |
2080 +============================================================================*/
2081
2082const PolyMod operator*( const PolyMod & s,
2083 const PolyMod & t )
2084{
2085 // Do * in terms of *= to maintain consistency.
2086 // Return value optimization compiles away the copy constructor.
2087 // const on return type disallows doing (u*v) = w ;
2088 return PolyMod( s ) *= t ;
2089}
2090
2091
2092/*=============================================================================
2093 |
2094 | NAME
2095 |
2096 | *=
2097 |
2098 | DESCRIPTION
2099 |
2100 | C-like multiply by operator
2101 |
2102 +============================================================================*/
2103
2104PolyMod &
2105PolyMod::operator*=( const PolyMod & t )
2106{
2107 int i, j ; // k 2
2108 ppuint coeff; // Coefficient of x term of t(x)
2109
2110 // Temporary storage for the new t(x). Can have degree up to n.
2111 Polynomial temp ;
2112
2113 // Get hold of the degree of f(x).
2114 int n = f_.deg() ;
2115
2116 // 0 n-1
2117 // Compute the coefficients of x , ..., x. These terms do not require
2118 // reduction mod f(x) because their degree is less than n.
2119 for (i = 0 ; i <= n ; ++i)
2120 temp[ i ] = coeffOfProduct( g_, t.g_, i, n ) ;
2121
2122 // n 2n-2 k
2123 // Compute the coefficients of x , ..., x. Replace t x with
2124 // k
2125 // k
2126 // t * [ x (mod f(x), p) ] from array powerTable when t is
2127 // k k
2128 // non-zero.
2129 for (i = n ; i <= 2 * n - 2 ; ++i)
2130 if ( (coeff = coeffOfProduct( g_, t.g_, i, n)) != 0 )
2131 for (j = 0 ; j <= n - 1 ; ++j)
2132 temp[ j ] = mod( temp[ j ] +
2133 mod( coeff * powerTable_[ offset(i) ] [ j ])) ;
2134
2135 for (i = 0 ; i <= n - 1 ; ++i)
2136 g_[ i ] = temp[ i ] ;
2137
2138 // Return (reference to) the product.
2139 return *this ;
2140}
2141
2142
2143/*=============================================================================
2144 |
2145 | NAME
2146 |
2147 | timesX
2148 |
2149 | DESCRIPTION
2150 |
2151 | Compute x g(x) (mod f(x), p)
2152 |
2153 | EXAMPLE
2154 |
2155 | g.timesX( t ) ;
2156 |
2157 | EXAMPLE
2158 | 3 2
2159 | Let n = 4, p = 5, and g(x) = 2 x + 4 x + 3 x. Let f(x) =
2160 | 4 2 4 3 2
2161 | x + x + 2 x + 3. Then x t (x) = 2 x + 4 x + 3 x and
2162 | 2 3 2
2163 | x g(x) (mod f(x), p) = 2 * (4 x + 3 x + 2) + 4 x + 3 x =
2164 | 3 2
2165 | 4 x + x + x + 4.
2166 | 3 2
2167 | = 4 x + 2 x + 3 x + 2.
2168 |
2169 | METHOD
2170 |
2171 | Uses a precomputed table of powers of x.
2172 |
2173 | n-1 n-2
2174 | Multiply g(x) = g x + g x + ... + g by shifting the coefficients
2175 | n-1 n-2 0
2176 |
2177 | n
2178 | to the left. If an x term appears, eliminate it by
2179 |
2180 | substitution using powerTable.
2181 |
2182 +============================================================================*/
2183
2184void PolyMod::timesX()
2185{
2186 int n = f_.deg() ;
2187
2188 #ifdef DEBUG_PP_POLYNOMIAL
2189 cout << "timesX: g( x ) = " << g_ << endl ;
2190 #endif
2191
2192 // Multiply g(x) by x by shifting the coefficients left in the array, giving
2193 // n-1
2194 // g x + ... + g x + 0
2195 // n-2 1
2196 //
2197 // but save the coefficient g first before overwriting it.
2198 // n-1
2199 ppuint g_coeff{ g_[ n - 1 ] } ;
2200
2201 for (int i = n-1 ; i >= 1 ; --i)
2202 g_[ i ] = g_[ i-1 ] ;
2203
2204 g_[ 0 ] = 0 ;
2205
2206 // n n
2207 // Replace g x with g * [ x (mod f(x), p) ] using
2208 // n-1 n-1
2209 // n
2210 // x from powerTable
2211
2212 if (g_coeff != 0)
2213 {
2214 for (int i = 0 ; i <= n - 1 ; ++i)
2215 g_[ i ] = mod( g_[ i ] +
2216 mod( g_coeff * powerTable_[ offset(n) ] [ i ] )) ;
2217 }
2218
2219 #ifdef DEBUG_PP_POLYNOMIAL
2220 cout << "timesX: x g( x ) = " << g_ << endl ;
2221 #endif
2222}
2223
2224
2225/*=============================================================================
2226 |
2227 | NAME
2228 |
2229 | square
2230 |
2231 | DESCRIPTION
2232 |
2233 | 2
2234 | Compute g (x) (mod f(x), p)
2235 |
2236 | EXAMPLE
2237 |
2238 | g.square() ;
2239 |
2240 |
2241 | EXAMPLE
2242 | 3 2
2243 | Let n = 4, p = 5, and g(x) = 4 x + x + 4. Let f(x) =
2244 |
2245 | 4 2 2 6 5 4 3 2
2246 | x + x + 2 x + 3. Then t (x) = x + 3 x + x + 2 x + 3 x + 1
2247 |
2248 | Now subsituting powers of x modulo f(x) from the power table,
2249 |
2250 | 2 3 2
2251 | t (x) (mod f(x), p) = (3 x + 3 x + 2 x + 3) +
2252 |
2253 | 3 2 2 3 2
2254 | 3 * (4 x + 3 x + 2 x) + 4 * (4 x + 3 x + 2) + 4 x + 4 x + 3 x + 1
2255 |
2256 | 3 2
2257 | = 2 x + 4 x + x + 1.
2258 |
2259 |
2260 | METHOD
2261 |
2262 | Uses a precomputed table of powers of x.
2263 |
2264 |
2265 | 2 2n-2 n n-1
2266 | Let g (x) = g x + ... + g x + g x + ... + g .
2267 | 2n-2 n n-1 0
2268 |
2269 | Compute the coefficients g using the function coeffOfSquare.
2270 | k
2271 |
2272 | 2
2273 | The next step is to reduce g (x) modulo f(x). To do so, replace
2274 |
2275 | k k
2276 | each non-zero term g x, n <= k <= 2n-2, by the term g * [ x mod f(x), p)]
2277 | k k
2278 |
2279 | which we get from the array powerTable_.
2280 |
2281 +============================================================================*/
2282
2283void PolyMod::square()
2284{
2285 // Get hold of the degree of f(x).
2286 int n = f_.deg() ;
2287
2288 #ifdef DEBUG_PP_POLYNOMIAL
2289 cout << "square: g( x ) = " << g_ << endl ;
2290 #endif
2291
2292 // Temporary storage for the new g(x). Can have degree up to n.
2293 Polynomial t ;
2294
2295 // 0 n-1
2296 // Compute the coefficients of x , ..., x. These terms do not require
2297 //
2298 // reduction mod f(x) because their degree is less than n.
2299 for (int i = 0 ; i <= n ; ++i)
2300 t[ i ] = coeffOfSquare( g_, i, n ) ;
2301
2302 // n 2n-2 k
2303 // Compute the coefficients of x , ..., x. Replace g x with
2304 // k k
2305 // g * [ x (mod f(x), p) ] from array powerTable_ when g is
2306 // k k
2307 // non-zero.
2308 for (int i = n ; i <= 2 * n - 2 ; ++i)
2309 {
2310 ppuint coeff{ 0 } ;
2311
2312 if ( (coeff = coeffOfSquare( g_, i, n )) != 0 )
2313
2314 for (int j = 0 ; j <= n- 1 ; ++j)
2315
2316 t[ j ] = mod( t[ j ] + mod( coeff * powerTable_[ offset(i) ] [ j ])) ;
2317 }
2318
2319 for (int i = 0 ; i <= n - 1 ; ++i)
2320
2321 g_[ i ] = t[ i ] ;
2322
2323 #ifdef DEBUG_PP_POLYNOMIAL
2324 cout << "square: g( x ) ^ 2 = " << g_ << endl ;
2325 #endif
2326}
2327
2328
2329/*=============================================================================
2330 |
2331 | NAME
2332 |
2333 | power
2334 |
2335 | DESCRIPTION
2336 | m
2337 | Compute g(x) = x (mod f(x), p).
2338 |
2339 | EXAMPLE
2340 | 4 2
2341 | Let n = 4, p = 5, f(x) = x + x + 2 x + 3, and m = 156.
2342 |
2343 | 156 = 0 . . . 0 1 0 0 1 1 1 0 0 (binary representation)
2344 | |<- ignore ->| S S SX SX SX S S (exponentiation rule,
2345 | S = square, X = multiply by x)
2346 | m
2347 | x (mod f(x), p) =
2348 |
2349 | 2 2
2350 | 6 S x = x
2351 |
2352 | 4 2
2353 | 5 S x = 4 x + 3 x + 2
2354 |
2355 |
2356 | 8 3 2
2357 | 4 S x = 4 x + 4 x + 1
2358 |
2359 | 9 3 2
2360 | 4 X x = 4 x + x + 3 x + 3
2361 |
2362 |
2363 | 18 2
2364 | 3 S x = 2 x + x + 2
2365 |
2366 | 19 3 2
2367 | 3 X x = 2 x + x + 2 x
2368 |
2369 |
2370 | 38 3 2
2371 | 2 S x = 2 x + 4 x + 3 x
2372 |
2373 | 39 3 2
2374 | 2 X x = 4 x + 2 x + x + 4
2375 |
2376 |
2377 | 78 3 2
2378 | 1 S x = 4 x + 2 x + 3 x + 2
2379 |
2380 | 156
2381 | 0 S x = 3
2382 |
2383 | METHOD
2384 |
2385 | Exponentiation by repeated squaring, using precomputed table of
2386 | powers. See ART OF COMPUTER PROGRAMMING, vol. 2, 2nd Ed.,
2387 | D. E. Knuth, pgs 441-443.
2388 |
2389 | n 2n-2
2390 | x, ... , x (mod f(x), p)
2391 |
2392 | m
2393 | to find g(x) = x (mod f(x), p), expand m into binary,
2394 |
2395 | k k-1
2396 | m = a 2 + a 2 + . . . + a 2 + a
2397 | k k-1 1 0
2398 |
2399 | m
2400 | where a = 1, and split x apart into
2401 | k
2402 |
2403 | k k
2404 | 2 2 a 2 a a
2405 | m k-1 1 0
2406 | x = x x . . . x x
2407 |
2408 |
2409 | Then to raise x to the mth power, do
2410 |
2411 |
2412 | t( x ) = x
2413 |
2414 | return if m = 1
2415 |
2416 |
2417 | for i = k-1 downto 0 do
2418 |
2419 | 2
2420 | t(x) = t(x) (mod f(x), p) // Square each time.
2421 |
2422 | if a = 1 then
2423 | i
2424 |
2425 | t(x) = x t(x) (mod f(x), p) // Times x only if current bit is 1
2426 | endif
2427 |
2428 | endfor
2429 | k
2430 | 2
2431 | The initial t(x) = x gets squared k times to give x . If a = 1 for
2432 | i
2433 | 0 <= i <= k-1, we multiply by x which then gets squared i times more
2434 |
2435 | i
2436 | 2
2437 | to give x . On a binary computer, we use bit shifting and masking to
2438 |
2439 | identify the k bits { a . . . a } to the right of the leading 1
2440 | k-1 0
2441 |
2442 | bit. There are log ( m ) - 1 squarings and (number of 1 bits) - 1
2443 | 2
2444 | multiplies.
2445 |
2446 +============================================================================*/
2447
2448const PolyMod power( const PolyMod & g1, const BigInt & m )
2449{
2450 // Return if g(x) != x
2451 if (g1.f_.deg() == 1 && g1[ 0 ] == 0 && g1[ 1 ] == 1)
2452 {
2453 ostringstream os ;
2454 os << "Error in PolyMod::power(): g( x ) != x "
2455 << "with deg g = " << g1.f_.deg() << " m = " << m ;
2456 throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
2457 }
2458
2459 // Exit right away if m = 1 and return a copy of g(x).
2460 PolyMod g( g1 ) ;
2461
2462 if (m == static_cast<BigInt>( 1u ))
2463 return g ;
2464
2465 // Find the number of the leading bit.
2466 int bitNum = m.maxBitNumber() ; // Number of highest possible bit.
2467
2468 #ifdef DEBUG_PP_POLYNOMIAL
2469 cout << "initial max bitNum = " << bitNum << endl ;
2470 cout << "g( x ) = " << g << endl ;
2471 #endif
2472
2473 while (!m.testBit( bitNum ))
2474 --bitNum ;
2475
2476 #ifdef DEBUG_PP_POLYNOMIAL
2477 cout << "after skipping leading 0 bits, bitNum = " << bitNum << endl ;
2478 #endif
2479
2480 if (bitNum == -1)
2481 {
2482 ostringstream os ;
2483 os << "PolyMod::x_to_power " << "bitNum == -1 internal error in PolyMod" ;
2484 throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
2485 }
2486
2487 #ifdef DEBUG_PP_POLYNOMIAL
2488 cout << "\nAfter skipping zero bits, bitNum = " << bitNum << endl ;
2489 #endif
2490
2491 // Exponentiation by repeated squaring. Discard the leading 1 bit.
2492 // Thereafter, square for every 0 bit; square and multiply by x for
2493 // every 1 bit.
2494 while ( --bitNum >= 0 )
2495 {
2496 g.square() ;
2497
2498 if (m.testBit( bitNum ))
2499 g.timesX() ;
2500
2501 #ifdef DEBUG_PP_POLYNOMIAL
2502 cout << "S " ;
2503 if (m.testBit( bitNum ))
2504 cout << "X " ;
2505 cout << "Bit num = " << bitNum << " g( x ) = " << g << endl ;
2506 #endif
2507 }
2508
2509 #ifdef DEBUG_PP_POLYNOMIAL
2510 cout << "Out of the loop bitNum = " << bitNum << " g( x ) = " << g << endl ;
2511 #endif
2512
2513 return g ;
2514}
2515
2516
2517/*=============================================================================
2518 |
2519 | NAME
2520 |
2521 | isInteger
2522 |
2523 | DESCRIPTION
2524 |
2525 | Getter function.
2526 |
2527 +============================================================================*/
2528
2529bool PolyMod::isInteger() const
2530{
2531 return g_.isInteger() ;
2532}
2533
2534
2535
2536/*------------------------------------------------------------------------------
2537| PolyOrder Implementation |
2538------------------------------------------------------------------------------*/
2539
2540/*=============================================================================
2541 |
2542 | NAME
2543 |
2544 | PolyOrder()
2545 |
2546 | DESCRIPTION
2547 |
2548 | Set a new value of f(x) with same degree n and modulus p.
2549 |
2550 +============================================================================*/
2551
2552void PolyOrder::resetPolynomial( const Polynomial & f )
2553{
2554 f_ = f ;
2555}
2556
2557
2558/*=============================================================================
2559 |
2560 | NAME
2561 |
2562 | PolyOrder()
2563 |
2564 | DESCRIPTION
2565 |
2566 | Initialize. Mainly do the prime factoring.
2567 |
2568 +============================================================================*/
2569
2570PolyOrder::PolyOrder( const Polynomial & f )
2571 : f_( f )
2572 , n_( f.deg() )
2573 , p_( f.modulus() )
2574 , mod( f.modulus() )
2575 , p_to_n_minus_1_( BigInt( 0u ))
2576 , r_( 0u )
2577 , a_( 0u )
2578 , factors_of_p_to_n_minus_1_()
2579 , factors_of_R_()
2580 , num_prim_poly_( 0u )
2581 , max_num_poly_( 0u )
2582 , Q_( 0 )
2583 , nullity_( 0 )
2584 , statistics_()
2585{
2586 // This is the most time consuming step for large n:
2587 // n
2588 // p - 1
2589 // Compute r = -------- and factor r into the product of primes.
2590 // p - 1
2591 try
2592 {
2593 computeMaxNumPoly() ;
2594 factorR() ;
2595 computeNumPrimPoly() ;
2596 }
2597 catch( BigIntMathError & e )
2598 {
2599 ostringstream os ;
2600 os << "PolyOrder: problem computing p^n or r = (p^n - 1 )/ (p - 1), or factoring r, or finding EulerPhi( p^n - 1 )/ n "
2601 << " p = " << p_ << " n = " << n_
2602 << " [ " << e.what() << " ] " ;
2603 throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
2604 }
2605
2606 // Copy the factoring statistics, and others.
2607 statistics_ = factors_of_R_.statistics_ ;
2608 statistics_.p = p_ ;
2609 statistics_.n = n_ ;
2610 statistics_.max_num_possible_poly = max_num_poly_ ;
2611 statistics_.num_primitive_poly = num_prim_poly_ ;
2612
2613 // Prepare the Q matrix to the proper size.
2614 try
2615 {
2616 Q_.clear() ;
2617 Q_.resize( n_ ) ;
2618
2619 for (int row = 0 ; row < n_ ; ++row)
2620 {
2621 Q_[ row ].resize( n_ ) ;
2622 }
2623
2624 }
2625 // Failed to resize Q matrix.
2626 catch( length_error & e )
2627 {
2628 throw PolynomialError( "PolyOrder::PolyOrder had a length_error exception and failed to allocate the Q matrix", __FILE__, __LINE__ ) ;
2629 }
2630}
2631
2632
2633 /*=============================================================================
2634 |
2635 | NAME
2636 |
2637 | factorR
2638 |
2639 | DESCRIPTION
2640 |
2641 | This is the most time consuming step for large n due to the integer
2642 | factorization.
2643 |
2644 | n
2645 | Find the maximum number of polynomials p
2646 |
2647 | Find
2648 | n
2649 | p - 1
2650 | r = --------
2651 | p - 1
2652 |
2653 | Compute the prime factorization of r.
2654 | n
2655 | Find number of primitive polynomials = Phi( p - 1 ) / n
2656 |
2657 | EXAMPLE
2658 |
2659 | See the examples in the code below.
2660 |
2661 +============================================================================*/
2662
2663void PolyOrder::factorR()
2664{
2665 // n
2666 // p - 1
2667 p_to_n_minus_1_ = BigInt( max_num_poly_ - static_cast<BigInt>( 1u ) ) ;
2668
2669 // n
2670 // Factor p - 1 into primes.
2671 // Pass in p and n in case we can do a fast table lookup.
2672 factors_of_p_to_n_minus_1_ = Factorization<BigInt>( p_to_n_minus_1_,
2673 FactoringAlgorithm::Automatic, p_, n_ ) ;
2674
2675 // n
2676 // p - 1
2677 // r = -------
2678 // p - 1
2679 r_ = p_to_n_minus_1_ / (p_ - 1u) ;
2680
2681 #ifdef DEBUG_PP_FACTOR
2682 cout << "p = " << p_ << endl ;
2683 cout << "n = " << n_ << endl ;
2684 cout << "p^n = " << max_num_poly_ << endl ;
2685 cout << "r = (p^n-1)/(p-1) = " << r_ << endl ;
2686 cout << "p_to_n_minus_1 = " << p_to_n_minus_1_ << endl ;
2687 cout << "factorization of p^n - 1 = " << endl ;
2688 for (unsigned int i = 0 ; i < factors_of_p_to_n_minus_1_.numDistinctFactors() ; ++i)
2689 cout << factors_of_p_to_n_minus_1_.primeFactor( i ) << " ^ " << factors_of_p_to_n_minus_1_.multiplicity( i ) << endl ;
2690 #endif // DEBUG_PP_FACTOR
2691
2692 // n
2693 // Factor r by starting with the factorization of p - 1
2694
2695 // Now we have to divide out all factors of (p - 1).
2696 // e.g.
2697 //
2698 // n 8 5 2
2699 // p - 1 = 19 - 1 = 2 3 5 17 181 3833
2700 //
2701 // 2
2702 // p - 1 = 19 - 1 = 2 3
2703 //
2704 // n 8 4
2705 // p - 1 = 19 - 1 = 2 5 17 181 3833
2706 // ----- ------
2707 // p - 1 19 - 1
2708 //
2709
2710 // n
2711 // Copy over the factors of p - 1
2712 factors_of_R_ = factors_of_p_to_n_minus_1_ ;
2713
2714 // We're done if p - 1 = 1.
2715 if (p_ > 2)
2716 {
2717 // Factor p - 1 into primes.
2718 Factorization<BigInt> factors_of_p_minus_1( static_cast<BigInt>( p_ - 1u ) ) ;
2719
2720 #ifdef DEBUG_PP_FACTOR
2721 cout << "factorization of p - 1 = " << endl ;
2722 for (unsigned int i = 0 ; i < factors_of_p_minus_1.numDistinctFactors() ; ++i)
2723 cout << factors_of_p_minus_1.primeFactor( i ) << " ^ " << factors_of_p_minus_1.multiplicity( i ) << endl ;
2724 #endif // DEBUG_PP_FACTOR
2725
2726 // n n
2727 // p-1 cannot have more distinct factors than p - 1 since p - 1 | p - 1
2728 if (factors_of_p_minus_1.numDistinctFactors() > factors_of_p_to_n_minus_1_.numDistinctFactors())
2729 {
2730 ostringstream os ;
2731 os << "factorR "
2732 << " number of distinct prime factors for p-1 = " << factors_of_p_minus_1.numDistinctFactors() << " > "
2733 << " number of distinct prime factors for p^n-1 = " << factors_of_p_to_n_minus_1_.numDistinctFactors()
2734 << " which is not possible since (p-1) | (p^n - 1)" ;
2735 throw BigIntUnderflow( os.str(), __FILE__, __LINE__ ) ;
2736 }
2737
2738 // Divide out p-1, one prime factor at a time.
2739 for (int i = 0, j = 0 ; i < factors_of_p_minus_1.numDistinctFactors() ; ++i)
2740 {
2741 BigInt factor_of_p_m_1 = factors_of_p_minus_1.primeFactor( i ) ;
2742 BigInt factor_of_r = factors_of_R_.primeFactor( j ) ;
2743
2744 // Divide out the common prime factor. Advance to next prime factor in the numerator.
2745 if (factor_of_p_m_1 == factor_of_r)
2746 {
2747 factors_of_R_[ j ].count_ -= factors_of_p_minus_1[ i ].count_ ;
2748 ++j ;
2749 }
2750 // Factor in denominator < factor in numerator. Advance to next factor in denominator.
2751 else if (factor_of_p_m_1 > factor_of_r)
2752 continue ;
2753 // Factor in denominator > factor in numerator. Advance to next factor in numerator. All smaller factors in numerator should have been divided out already.
2754 else
2755 ++j ;
2756
2757 #ifdef DEBUG_PP_FACTOR
2758 cout << " i = " << i << " prime factor of p-1 = " << factor_of_p_m_1 << " j = " << j << " prime factor num = " << factor_of_r << endl ;
2759 #endif // DEBUG_PP_FACTOR
2760 }
2761
2762 #ifdef DEBUG_PP_FACTOR
2763 cout << "factorization of r = " << endl ;
2764 for (unsigned int i = 0 ; i < factors_of_R_.numDistinctFactors() ; ++i)
2765 cout << factors_of_R_.primeFactor( i ) << " ^ " << factors_of_R_.multiplicity( i ) << endl ;
2766 #endif // DEBUG_PP_FACTOR
2767 }
2768
2769 return ;
2770}
2771
2772
2773/*=============================================================================
2774 |
2775 | NAME
2776 |
2777 | computeMaxNumPoly
2778 |
2779 | DESCRIPTION
2780 | n
2781 | Maximum number of possible polynomials of degree n modulo p = p
2782 |
2783 | EXAMPLE
2784 |
2785 +============================================================================*/
2786
2787void PolyOrder::computeMaxNumPoly()
2788{
2789 max_num_poly_ = power( p_, n_ ) ;
2790
2791 return ;
2792}
2793
2794
2795/*=============================================================================
2796 |
2797 | NAME
2798 |
2799 | computeNumPrimPoly
2800 |
2801 | DESCRIPTION
2802 | n
2803 | Find number of primitive polynomials = Phi( p - 1 ) / n
2804 |
2805 | EXAMPLE
2806 |
2807 | See the examples in the code below.
2808 |
2809 +============================================================================*/
2810
2811void PolyOrder::computeNumPrimPoly()
2812{
2813 // n
2814 // Compute the number of primitive polynomials = Phi( p - 1 ) / n
2815 //
2816 // Recall Euler's totient is
2817 //
2818 // ----- ----- -----
2819 // Phi[ n ] = n | | (1 - 1/p ) = n | | (p - 1) / | | p
2820 // i i i
2821 // p = all distinct
2822 // i
2823 // prime factors of n
2824 //
2825 // For example,
2826 //
2827 // 8 5 2
2828 // 19 - 1 = 16983563040 = 2 3 5 17 181 3833
2829 //
2830 // 8 8
2831 // Phi( 19 - 1 ) = (19 - 1) (2-1) (3-1) (5-1) (17-1) (181-1) (3833-1) / (2 3 5 17 181 3833) = 4237885440
2832 //
2833 // You can check with Wolfram Alpha on the web,
2834 //
2835 // http://www.wolframalpha.com/input/?i=eulerphi%28+19^8-1%29
2836 // 8
2837 // then the number of primitive polynomials = Phi( 19 - 1) / 8 = 529735680
2838
2839 num_prim_poly_ = p_to_n_minus_1_ ;
2840 vector<BigInt> distinct_factors_of_p_to_n_minus_1_ = factors_of_p_to_n_minus_1_.getDistinctPrimeFactors() ;
2841
2842 for (auto & f : distinct_factors_of_p_to_n_minus_1_)
2843 num_prim_poly_ *= (f - static_cast<BigInt>( 1u )) ;
2844
2845 for (auto & f : distinct_factors_of_p_to_n_minus_1_)
2846 num_prim_poly_ /= f ;
2847
2848 num_prim_poly_ /= static_cast<BigInt>( n_ ) ;
2849
2850 return ;
2851}
2852
2853
2854/*=============================================================================
2855 |
2856 | NAME
2857 |
2858 | orderM
2859 |
2860 | DESCRIPTION
2861 | m
2862 | Check that x (mod f(x), p) is not an integer for m = r / p but skip
2863 | i
2864 | n
2865 | p - 1 th
2866 | this test if p | (p-1). Recall r = -------, and p = i prime in
2867 | i p - 1 i
2868 |
2869 | the factorization of r.
2870 |
2871 |
2872 | EXAMPLE
2873 | 2
2874 | Let n = 4 and p = 5. Then r = 156 = 2 * 3 * 13, and p = 2, p = 3,
2875 | 1 2
2876 |
2877 | and p = 13. m = { 156/2, 156/3, 156/13 } = { 78, 52, 12 }. We can
2878 | 3
2879 |
2880 | skip the test for m = 78 because p = 2 divides p-1 = 4. Exponentiation
2881 | 1
2882 |
2883 | 52 3 2 12
2884 | gives x = 2 x + x + 4 x, which is not an integer and x =
2885 |
2886 | 3 2
2887 | 4 x + 2 x + 4 x + 3 which is not an integer either, so we return
2888 |
2889 | true.
2890 |
2891 | METHOD
2892 |
2893 | Exponentiate x with x_to_power and test the result with is_integer.
2894 | Return right away if the result is not an integer.
2895 |
2896 +============================================================================*/
2897
2898bool PolyOrder::orderM()
2899{
2900 ppuint p{ f_.modulus() } ;
2901
2902 for (int i = 0 ; i < factors_of_R_.numDistinctFactors() ; ++i)
2903 {
2904 // Can we skip this order m test?
2905 if (!factors_of_R_.skipTest( p, i ))
2906 {
2907 BigInt m = r_ / factors_of_R_.primeFactor( i ) ;
2908
2909 Polynomial x1( "x" ) ;
2910 x1.setModulus( p ) ;
2911 PolyMod x( x1, f_ ) ;
2912
2913 PolyMod x_to_m = power( x, m ) ;
2914
2915 #ifdef DEBUG_PP_POLYNOMIAL
2916 cout << "Prime factor p[ " << i << " ] = " << factors_of_R_.primeFactor( i ) << endl ;
2917 cout << "m = " << m << endl ;
2918 cout << "x^m = " << x_to_m << endl ;
2919 #endif
2920
2921 // Early out.
2922 if (x_to_m.isInteger())
2923 return( false ) ;
2924 }
2925 }
2926
2927 return( true ) ;
2928
2929}
2930
2931
2932/*=============================================================================
2933 |
2934 | NAME
2935 |
2936 | orderR
2937 |
2938 | DESCRIPTION
2939 | n
2940 | r p - 1
2941 | Compute x (mod f(x), p) = a, where r = -------
2942 | p - 1
2943 |
2944 | If a is not an integer, return 0, else return a itself.
2945 |
2946 | EXAMPLE
2947 | 4 2
2948 | f(x) = x + x + 2 x + 3, n = 4 and p = 5. Then r = 156 and
2949 |
2950 | r 156
2951 | x = x = 3 (mod f(x), 5) = 3, so we return 3.
2952 |
2953 | 4
2954 | But for f(x) = x + x + 3, n = 4, p = 5,
2955 |
2956 | r 156 3
2957 | x = x = 3 x + 2 x + 1 (mod f(x), 5) so we return 0.
2958 |
2959 | METHOD
2960 | r
2961 | First compute g(x) = x (mod f(x), p).
2962 | Then test if g(x) is a constant polynomial.
2963 |
2964 +============================================================================*/
2965
2966ppuint PolyOrder::orderR()
2967{
2968 Polynomial x1( "x", p_ ) ;
2969 PolyMod x( x1, f_ ) ;
2970
2971 PolyMod x_to_r = power( x, r_ ) ;
2972
2973 #ifdef DEBUG_PP_POLYNOMIAL
2974 cout << "r = " << r_ << endl ;
2975 cout << "x^r = " << x_to_r << endl ;
2976 cout << "is integer = " << x_to_r.isInteger() << endl ;
2977 #endif
2978
2979 if (x_to_r.isInteger())
2980 // Return the value a = constant term of g(x).
2981 return x_to_r[ 0 ] ;
2982 else
2983 return 0 ;
2984}
2985
2986
2987/*=============================================================================
2988 |
2989 | NAME
2990 |
2991 | max_order
2992 |
2993 | DESCRIPTION
2994 | k n
2995 | Check if x = 1 (mod f(x), p) only when k = p - 1 and not for any smaller
2996 | power of k, i.e. that f(x) is a primitive polynomial.
2997 |
2998 | INPUT
2999 |
3000 | f (int *) Monic polynomial f(x).
3001 | n (int, n >= 1) Degree of f(x).
3002 | p (int) Modulo p coefficient arithmetic.
3003 |
3004 | RETURNS
3005 |
3006 | true if f( x ) is primitive.
3007 | false if it isn't.
3008 |
3009 | EXAMPLE
3010 |
3011 | 4
3012 | f( x ) = x + x + 1 is a primitive polynomial modulo p = 2,
3013 | 4
3014 | because it generates the group GF( 2 ) with the exception of
3015 | 2 15
3016 | zero. The powers {x, x , ... , x } modulo f(x), mod 2 are
3017 | 16 4
3018 | distinct and not equal to 1 until x (mod x + x + 1, 2) = 1.
3019 |
3020 | METHOD
3021 |
3022 | Confirm f(x) is primitive using the definition of primitive
3023 | polynomial as a generator of the Galois group
3024 | n n
3025 | GF( p ) by testing that p - 1 is the smallest power for which
3026 | k
3027 | x = 1 (mod f(x), p).
3028 |
3029 +============================================================================*/
3030
3031bool PolyOrder::maximal_order()
3032{
3033 // Highest possible order for x.
3034 BigInt maxOrder = power( f_.modulus(), f_.deg()) - static_cast<BigInt>( 1u ) ;
3035
3036 // Start with the polynomial x, modulo f(x) and p.
3037 Polynomial x1( "x", f_.modulus() ) ;
3038 PolyMod x( x1, f_ ) ;
3039
3040 #ifdef DEBUG_PP_SLOW_PRIMITIVITY_CHECK
3041 cout << endl << "Starting with x (mod f(x),p) = " << x << endl ;
3042 #endif
3043
3044 bool x_to_k_is_one_too_early = false ;
3045
3046 BigInt k { 1u } ;
3047 while ( k < maxOrder )
3048 {
3049 PolyMod x_to_power = power( x, k ) ; // x^k (mod f(x), p)
3050
3051 #ifdef DEBUG_PP_SLOW_PRIMITIVITY_CHECK
3052 cout << "k = " << k << " x^k (mod f(x),p) = " << x_to_power
3053 << " x^k is integer = " << (x_to_power.isInteger() ? "yes" : "no")
3054 << " x^k has constant term x[0] = " << x_to_power[0] << endl ;
3055 #endif
3056
3057 // x^k (mod f(x), p) = 1 for k < p^n-1, so x cannot be a generator of the field.
3058 if (x_to_power.isInteger() && x_to_power[0] == 1u)
3059 {
3060 #ifdef DEBUG_PP_SLOW_PRIMITIVITY_CHECK
3061 cout << "x^k = 1 for k = " << k << " break out of the loop" << endl ;
3062 #endif
3063
3064 x_to_k_is_one_too_early = true ;
3065 break ;
3066 }
3067
3068 ++k ;
3069 }
3070
3071 // Order of x < p^n-1, so definitely not primitive.
3072 if (x_to_k_is_one_too_early)
3073 return false ;
3074 // x^k != 1 for k = 1, ..., p^n-2. x could be still primitive.
3075 else
3076 {
3077 PolyMod x_to_power = power( x, maxOrder ) ;
3078
3079 #ifdef DEBUG_PP_SLOW_PRIMITIVITY_CHECK
3080 cout << "maxOrder = " << maxOrder << " x^maxOrder (mod f(x),p) = " << x_to_power
3081 << " x^maxOrder is integer = " << (x_to_power.isInteger() ? "yes" : "no")
3082 << " x^maxOrder has constant term x[0] = " << x_to_power[0] << endl ;
3083 #endif
3084
3085 if (x_to_power.isInteger() && x_to_power[0] == 1u)
3086 // x^k = 1 for k = p^n-1, so primitive by definition.
3087 return true ;
3088 else
3089 return false ;
3090 }
3091
3092 return false ;
3093}
3094
3095
3096/*=============================================================================
3097 |
3098 | NAME
3099 |
3100 | hasMultipleDistinctFactors
3101 |
3102 | DESCRIPTION
3103 |
3104 | Returns true if the monic polynomial f( x ) has two or more distinct
3105 | irreducible factors, false otherwise.
3106 |
3107 | If earlyOut is false, compute the exact nullity in findNullity() instead
3108 | of stopping when the nullity is >= 2.
3109 |
3110 | EXAMPLE
3111 |
3112 | Let n = 4, p = 5
3113 |
3114 | 4 2
3115 | f( x ) = x + x + 2 x + 3 is irreducible, so has one distinct factor.
3116 |
3117 | 4 3 2 4
3118 | f( x ) = x + 4x + x + 4x + 1 = (x + 1) has one distinct factor.
3119 |
3120 | 3 2
3121 | f( x ) = x + 3 = (x + 3x + 4)(x + 2) has two distinct irreducible factors.
3122 |
3123 | 4 3 2 2
3124 | f( x ) = x + 3x + 3x + 3x + 2 = (x + 1) (x + 2) (x + 3) has 3 distinct
3125 | irreducible factors.
3126 |
3127 | 4
3128 | f(x) = x + 4 = (x+1)(x+2)(x+3)(x+4) has 4 distinct irreducible factors.
3129 |
3130 | METHOD
3131 |
3132 | Berlekamp's method for factoring polynomials over GF( p ), modified to test
3133 | for irreducibility only.
3134 |
3135 | See my notes; I skip the polynomial GCD step which ensures polynomials
3136 | are square-free due to time constraints, but this requires a proof that
3137 | the method is still valid.
3138 |
3139 +============================================================================*/
3140
3141bool PolyOrder::hasMultipleDistinctFactors( bool earlyOut )
3142
3143{
3144 // Generate the Q-I matrix.
3145 generateQMatrix() ;
3146
3147 // Find nullity of Q-I
3148 findNullity( earlyOut ) ;
3149
3150 // If nullity_ >= 2, f( x ) is a reducible polynomial modulo p since it has
3151 // two or more distinct irreducible factors.
3152 // e
3153 // Nullity of 1 implies f( x ) = p( x ) for some power e >= 1 so we cannot
3154 // determine reducibility.
3155 if (nullity_ >= 2)
3156 return true ;
3157
3158 return false ;
3159
3160}
3161
3162
3163/*=============================================================================
3164 |
3165 | NAME
3166 |
3167 | isPrimitive
3168 |
3169 | DESCRIPTION
3170 |
3171 | Check if a given polynomial f(x) of degree n modulo p is primitive.
3172 |
3173 +============================================================================*/
3174
3175bool PolyOrder::isPrimitive()
3176{
3177 bool isPrimitive = false ;
3178 ++statistics_.num_poly_tested ;
3179
3180 try
3181 {
3182 BigInt max_num_possible_poly = power( p_, n_ ) ;
3183 ArithModP modp( p_ ) ;
3184
3185 // Constant coefficient of f(x) * (-1)^n must be a primitive root of p.
3186 if (modp.constCoeffIsPrimitiveRoot( f_[0], f_.deg() ))
3187 {
3188 ++statistics_.num_where_const_coeff_is_primitive_root ;
3189
3190 #ifdef DEBUG_PP_PRIMITIVITY
3191 cout << " (-1)^n const coeff " << f_[ 0 ] << " is primitive root of " << p_ << endl ;
3192 #endif
3193
3194 // f(x) can't have any linear factors.
3195 if (!f_.hasLinearFactor())
3196 {
3197 ++statistics_.num_free_of_linear_factors ;
3198
3199 #ifdef DEBUG_PP_PRIMITIVITY
3200 cout << " No linear factors" << endl ;
3201 #endif
3202
3203 // Do more in-depth checking.
3204
3205 // f(x) can't have two or more distinct irreducible factors.
3206 if (!hasMultipleDistinctFactors())
3207 {
3208 ++statistics_.num_irreducible_to_power ;
3209
3210 #ifdef DEBUG_PP_PRIMITIVITY
3211 cout << " One distinct irreducible factor (possibly repeated)" << endl ;
3212 #endif
3213
3214 // r
3215 // x (mod f(x), p) = a_ must be an integer.
3216 ppuint a{ orderR() } ;
3217 if (a != 0)
3218 {
3219 ++statistics_.num_order_r ;
3220
3221 #ifdef DEBUG_PP_PRIMITIVITY
3222 cout << " x^r = a (integer)" << endl ;
3223 #endif
3224
3225 // n
3226 // Check if (-1) (constant coefficient of f(x)) = a_ (mod p)
3227 //
3228 if (modp.constCoeffTest( f_[ 0 ], a, f_.deg() ))
3229 {
3230 ++statistics_.num_passing_const_coeff_test ;
3231
3232 #ifdef DEBUG_PP_PRIMITIVITY
3233 cout << " a (integer) = (-1)^n f[0]" << endl ;
3234 #endif
3235
3236 // x^m != integer for all m = r / q, q a prime divisor of r.
3237 if (orderM())
3238 {
3239 ++statistics_.num_order_m ;
3240
3241 #ifdef DEBUG_PP_PRIMITIVITY
3242 cout << " x^m != integer for m = r / prime divisor of r" << endl ;
3243 #endif
3244
3245 isPrimitive = true ;
3246 goto Exit ;
3247
3248 } // end order m
3249 } // end const coeff test
3250 } // end order r
3251 } // end can't determine if reducible
3252 } // end no linear factors
3253 } // end constant coefficient primitive.
3254 }
3255 catch( ArithModPError & e )
3256 {
3257 ostringstream os ;
3258 os << "PolyOrder::isPrimitive had a mod p arithmetic error [ " << e.what() << " ] " ;
3259 throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
3260 }
3261
3262 Exit:
3263 return isPrimitive ;
3264}
3265
3266
3267/*=============================================================================
3268 |
3269 | NAME
3270 | generateQMatrix
3271 |
3272 | DESCRIPTION
3273 |
3274 | Generate n x n matrix Q - I, where rows of Q are the powers,
3275 |
3276 | p 2p (n-1) p
3277 | 1, x (mod f(x),p), x (mod f(x), p), ... , x (mod f(x), p)
3278 |
3279 | for monic polynomial f(x).
3280 |
3281 | EXAMPLE
3282 |
3283 | 4 2
3284 | f(x) = x + x + 2 x + 3, n = 4, p = 5
3285 |
3286 | ( 1 ) ( 1 ) ( 1 0 0 0 )
3287 | ( ) ( ) ( )
3288 | ( 5 ) ( 2 3 ) ( )
3289 | Q = ( x ) ( 2x + 3x + 4 x ) ( 0 2 3 4 )
3290 | ( ) ( ) ( )
3291 | ( ) ( 2 3 ) ( )
3292 | ( 10 ) = ( 3 + 4 x + x ) = ( )
3293 | ( x ) ( ) ( 3 0 4 1 )
3294 | ( ) ( ) ( )
3295 | ( 15 ) ( 2 3 ) ( )
3296 | ( x ) ( 4x + 4x + 3x ) ( 0 4 4 3 )
3297 | ( ) ( ) ( )
3298 |
3299 | 2 3
3300 | 1 x x x
3301 | ( 0 0 0 0 )
3302 | ( 0 1 3 4 )
3303 | Q - I = ( 3 0 3 1 )
3304 | ( 0 4 4 2 )
3305 |
3306 |
3307 | The left nullspace has dimension = 1 with basis { (1 0 0 0) }.
3308 |
3309 +============================================================================*/
3310
3311void PolyOrder::generateQMatrix()
3312{
3313 // Check for invalid inputs.
3314 if (n_ < 2 || p_ < 2)
3315 throw PolynomialRangeError( "generateQMatrix has n < 2 or p < 2", __FILE__, __LINE__ ) ;
3316
3317 // Row 0 of Q = (1 0 ... 0).
3318 Q_[ 0 ][ 0 ] = 1 ;
3319 for (int i = 1 ; i < n_ ; ++i)
3320 Q_[ 0 ][ i ] = 0 ;
3321
3322 // p
3323 // Let q(x) = x (mod f(x),p)
3324 // and Q[ 1 ] = coefficients of q(x).
3325 Polynomial x1( "x", p_ ) ;
3326 PolyMod x( x1, f_ ) ;
3327 PolyMod xp = power( x, static_cast< BigInt >( p_ ) ) ;
3328
3329 #ifdef DEBUG_PP_POLYNOMIAL
3330 cout << "x ^ p (mod f(x),p) = " << xp << endl ;
3331 cout << "initial Q matrix " << printQMatrix() ;
3332 #endif
3333
3334 PolyMod q = xp ;
3335
3336 for (int i = 0 ; i < n_ ; ++i)
3337 Q_[ 1 ][ i ] = xp[ i ] ;
3338
3339 // pk
3340 // Row k of Q = x (mod f(x), p) 2 <= k <= n-1, computed by
3341 // p
3342 // multiplying each previous row by x (mod f(x),p).
3343 for (int k = 2 ; k <= n_- 1 ; ++k)
3344 {
3345 q *= xp ;
3346
3347 #ifdef DEBUG_PP_POLYNOMIAL
3348 cout << "x ^ pk (mod f(x),p) = " << q << " for k = " << k << endl ;
3349 #endif
3350
3351 for (int i = 0 ; i < n_ ; ++i)
3352 Q_[ k ][ i ] = q[ i ] ;
3353 }
3354
3355 #ifdef DEBUG_PP_POLYNOMIAL
3356 cout << "computed Q matrix " << printQMatrix() ;
3357 #endif
3358
3359 // Subtract Q - I
3360 for (int row = 0 ; row < n_ ; ++row)
3361 {
3362 Q_[ row ][ row ] = mod( Q_[ row ][ row ] - 1 ) ;
3363 }
3364
3365 #ifdef DEBUG_PP_POLYNOMIAL
3366 cout << "computed Q-I matrix " << printQMatrix() ;
3367 #endif
3368
3369 return ;
3370}
3371
3372
3373/*=============================================================================
3374 |
3375 | NAME
3376 |
3377 | printQMatrix
3378 |
3379 | DESCRIPTION
3380 |
3381 | Print the matrix to the console.
3382 |
3383 +============================================================================*/
3384
3385string PolyOrder::printQMatrix() const
3386{
3387 // Print the matrix as a string.
3388 ostringstream os ;
3389
3390 os << endl ;
3391 for (int row = 0 ; row < n_ ; ++row)
3392 {
3393 os << "( " ;
3394 for (int col = 0 ; col < n_ ; ++col)
3395 {
3396 os << setw( 4 ) << setfill( ' ' ) << Q_[ row ][ col ] ;
3397 }
3398 os << " )" << endl ;
3399 }
3400
3401 return os.str() ;
3402}
3403
3404
3405/*=============================================================================
3406 |
3407 | NAME
3408 |
3409 | findNullity
3410 |
3411 | DESCRIPTION
3412 |
3413 | Computes the nullity of Q.
3414 | If earlyOut is true, stop when the nullity is >= 2 and return 2.
3415 |
3416 | EXAMPLE
3417 |
3418 | Let p = 5 and n = 3. We use the facts that -1 = 4 (mod 5), 1/2 = 3, -1/2 = 2,
3419 | 1/3 = 2, -1/3 = 3, 1/4 = 4, -1/4 = 1.
3420 |
3421 | Consider the matrix
3422 |
3423 | ( 2 3 4 )
3424 | Q = ( 0 2 1 )
3425 | ( 3 3 3 )
3426 |
3427 | Begin with row 1. No pivotal columns have been selected yet. Scan row 1 and
3428 | pick column 1 as the pivotal column because it contains a nonzero element.
3429 |
3430 | Normalizing column 1 by multiplying by -1/pivot = -1/2 = 2 gives
3431 |
3432 | ( 4 3 4 )
3433 | ( 0 2 1 )
3434 | ( 1 3 3 )
3435 |
3436 | Now perform column reduction on column 2 by multiplying the pivotal column 1
3437 | by 3 (the column 2 element in the current row) and adding back to row 2.
3438 |
3439 | ( 4 0 4 )
3440 | ( 0 2 1 )
3441 | ( 1 1 3 )
3442 |
3443 | Column reduce column 3 by multiplying the pivotal column by 4 and adding back to row 3,
3444 |
3445 | ( 4 0 0 )
3446 | ( 0 2 1 )
3447 | ( 1 1 2 )
3448 |
3449 | For row 2, pick column 2 as the pivotal column, normalize it and reduce columns 1, then 3,
3450 |
3451 | ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 )
3452 | ( 0 2 1 ) => ( 0 4 1 ) => ( 0 4 1 ) => ( 0 4 0 )
3453 | ( 1 1 2 ) ( 1 2 2 ) ( 1 2 2 ) ( 1 2 4 )
3454 | norm. c.r. 1 c.r. 3
3455 |
3456 | For row 3, we must pick column 3 as pivotal column because we've used up columns 1 and 2,
3457 |
3458 | ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 )
3459 | ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 )
3460 | ( 1 2 4 ) ( 1 2 4 ) ( 1 2 4 ) ( 0 0 4 )
3461 | norm. c.r. 1 c.r. 2
3462 |
3463 | The nullity is zero, since we were always able to find a pivot in each row.
3464 |
3465 | METHOD
3466 |
3467 | Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth, Addison-Wesley.
3468 |
3469 | We combine operations of normalization of columns,
3470 |
3471 | ( c ) ( C )
3472 | ( c ) ( C )
3473 | ( . ) ( C )
3474 | ( . . . q . . . ) row ================> ( . . . -1 . . . ) row
3475 | ( c ) ( C )
3476 | ( c ) column times ( C )
3477 | ( c ) -1/a modulo p ( C )
3478 | pivotCol pivotCol
3479 |
3480 | and column reduction,
3481 |
3482 | ( C b ) ( C B )
3483 | ( C b ) ( C B )
3484 | ( C b ) ( C B )
3485 | ( . . . -1 . . .e . . . ) row ================> ( . . . -1 . . . 0 . . . )
3486 | ( C b ) ( C B )
3487 | ( C b ) pivotCol times ( C B )
3488 | ( C b ) e added back to col ( C B )
3489 | pivotCol col col
3490 |
3491 | to reduce the matrix to a form in which columns having pivots are zero until
3492 | the pivotal row.
3493 |
3494 | The column operations don't change the left nullspace of the
3495 | matrix.
3496 |
3497 | The matrix rank is the number of pivotal rows since they are linearly
3498 | independent. The nullity is then the number of non-pivotal rows.
3499 |
3500 +============================================================================*/
3501
3502void PolyOrder::findNullity( bool earlyOut )
3503{
3504 try
3505 {
3506 InverseModP invmod( p_ ) ;
3507
3508 #ifdef DEBUG_PP_POLYNOMIAL
3509 cout << "Q-I matrix " << printQMatrix() ;
3510 #endif
3511
3512 vector<bool>pivotInCol( n_, false ) ; // Is false if the column has no pivotal element.
3513
3514 nullity_ = 0 ;
3515 int pivotCol = -1 ; // No pivots yet.
3516
3517 // Sweep through each row.
3518 for (int row = 0 ; row < n_ ; ++row)
3519 {
3520 // Search for a pivot in this row: a non-zero element
3521 // in a column which had no previous pivot.
3522 bool found = false ;
3523 for (int col = 0 ; col < n_ ; ++col)
3524 {
3525 if (Q_[ row ][ col ] > 0 && !pivotInCol[ col ])
3526 {
3527 found = true ;
3528 pivotCol = col ;
3529 break ;
3530 }
3531 }
3532
3533 // No pivot; increase nullity by 1.
3534 if (found == false)
3535 {
3536 ++nullity_ ;
3537
3538 // Early out.
3539 if (earlyOut && nullity_ >= 2)
3540 goto EarlyOut ;
3541 }
3542 // Found a pivot, q.
3543 else
3544 {
3545 ppsint q = Q_[ row ][ pivotCol ] ;
3546
3547 // Compute -1/q (mod p)
3548 ppsint t = mod( -invmod( q )) ;
3549
3550 // Normalize the pivotal column.
3551 for (int r = 0 ; r < n_ ; ++r)
3552 {
3553 Q_[ r ][ pivotCol ] = mod( t * Q_[ r ][ pivotCol ]) ;
3554 }
3555
3556 // Do column reduction: Add C times the pivotal column to the other
3557 // columns where C = element in the other column at current row.
3558 for (int col = 0 ; col < n_ ; ++col)
3559 {
3560 if (col != pivotCol)
3561 {
3562 q = Q_[ row ][ col ] ;
3563
3564 for (int r = 0 ; r < n_ ; ++r)
3565 {
3566 t = mod( q * Q_[ r ][ pivotCol ]) ;
3567 Q_[ r ][ col ] = mod( t + Q_[ r ][ col ] ) ;
3568 }
3569 }
3570 }
3571
3572 // Record the presence of a pivot in this column.
3573 pivotInCol[ pivotCol ] = true ;
3574
3575 #ifdef DEBUG_PP_POLYNOMIAL
3576 cout << "row = " << row << " pivot = " << q << " (-1/q) = " << t << " nullity = " << nullity_
3577 << " Row reduced Q-I matrix " << printQMatrix() ;
3578 #endif
3579
3580 } // found a pivot
3581
3582 } // end for row
3583
3584 EarlyOut: ;
3585 #ifdef DEBUG_PP_POLYNOMIAL
3586 cout << "Row reduced Q-I matrix " << printQMatrix() ;
3587 #endif
3588 }
3589 catch( ArithModPError & e )
3590 {
3591 ostringstream os ;
3592 os << "PolyOrder::findNullity failed in matrix row reduction [ " << e.what() << " ] " ;
3593 throw PrimpolyError( os.str(), __FILE__, __LINE__ ) ;
3594 }
3595
3596 // Automagically free pivotInCol and mod objects.
3597
3598} // ===================== end of function findNullity =====================