1 /*==============================================================================
  2 |
  3 |  NAME
  4 |
  5 |     ppPolynomial.h
  6 |
  7 |  DESCRIPTION
  8 |
  9 |     Header file for all the polynomial 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.1 - A Program for Computing Primitive Polynomials.
 17 |     Copyright (C) 1999-2021 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 !DOT! replaced by . and the !AT! replaced by @
 34 |
 35 ==============================================================================*/
 36 
 37 // Wrap this header file to prevent duplication if it is included
 38 // accidentally more than once.
 39 #ifndef __POLYNOMIAL_H__
 40 #define __POLYNOMIAL_H__
 41 
 42 /*=============================================================================
 43 |
 44 | NAME
 45 |
 46 |     PolynomialError
 47 |     PolynomialRangeError
 48 |
 49 | DESCRIPTION
 50 |
 51 |     Exception classes for classes Polynomial and PolyMod
 52 |     derived from the STL exception class runtime_error.
 53 |
 54 +============================================================================*/
 55 
 56 // General purpose exception for a polynomial.
 57 class PolynomialError : public runtime_error
 58 {
 59     public:
 60         // Throw with error message, file name and line number.
 61         PolynomialError( const string & description, const string & file, const int & line )
 62         : runtime_error( description + " in file " + file + " at line " + to_string(line) )
 63         {
 64         } ;
 65 
 66         // Throw with an error message.
 67         PolynomialError( const string & description )
 68             : runtime_error( description )
 69         {
 70         } ;
 71 
 72         // Throw with default error message.
 73         PolynomialError()
 74             : runtime_error( "Polynomial error:  " )
 75         {
 76         } ;
 77 
 78 } ; // end class PolynomialError
 79 
 80 // Typically thrown when we parse a bad string which is not a proper polynomial.
 81 class PolynomialRangeError : public PolynomialError
 82 {
 83     public:
 84         // Throw with error message, file name and line number.
 85         PolynomialRangeError( const string & description, const string & file, const int & line )
 86         : PolynomialError( description + " in file " + file + " at line " + to_string(line) )
 87         {
 88         } ;
 89 
 90         // Throw with an error message.
 91         PolynomialRangeError( const string & description )
 92             : PolynomialError( description )
 93         {
 94         } ;
 95 
 96         // Throw with default error message.
 97         PolynomialRangeError()
 98             : PolynomialError( "Polynomial range error:  " )
 99         {
100         } ;
101 
102 } ; // end class PolynomialRangeError
103 
104 
105 /*=============================================================================
106 |
107 | NAME
108 |
109 |     Polynomial
110 |
111 | DESCRIPTION
112 |
113 |     Represents the monic polynomial f( x ) of degree n, with
114 |     coefficients in GF( p ).  Precisely,
115 |
116 |                   n         n-1
117 |         f( x ) = x  +  a   x    + ... a x + a     0 <= a  < p
118 |                         n-1            1     0          i
119 |
120 |     The constant polynomial is of degree 0.
121 |
122 |     We support these basic operations:
123 |
124 |         Polynomial f ;                Construct f(x) = 0 mod 2.
125 |         <destructor>
126 |         Polynomial f( f2 ) ;          Copy, f(x) = f2(x)
127 |         Polynomial f = f2 ;           Overwrite, f(x) = f2(x)
128 |         Polynomial f( "x^2+1, 3" ) ;  Polynomial from string.
129 |         String s = f ;                Poly to string.
130 |         stream << p                   Read the poly.
131 |         p >> stream                   Print the poly.
132 |         p1 == p2, p1 != p2            Test for equality.
133 |         f[ i ] = coeff                Set poly coefficient.
134 |         f.deg()                       Return the degree of f(x).
135 |         f.modulus()                   Return the modulus p of f(x).
136 |         int coeff = f[ i ]            Get poly coefficient.
137 |         f(x) + g(x)                   Add poly mod p.
138 |         f(x)                          Evaluate f(x) for integer x
139 |         f.hasLinearFactor()           Check if f(x) has any linear factors.
140 |         f.isInteger()                 Is f(x) = constant?
141 |         f.firstTrialPoly()            Set f(x) = x^n - 1
142 |         f.nextTrialPoly()             Set f(x) = next poly in sequence.
143 |
144 |     Exceptions:  throw PolynomialError or one of its subclasses such as 
145 |     PolynomialRangeError.
146 |
147 | NOTES
148 |
149 |     The member functions and friends are documented in ppPolynomial.cpp
150 |
151 +============================================================================*/
152 
153 class Polynomial
154 {
155     public:
156         // Default constructor which sets f(x) = 0 modulo 2.
157         Polynomial() ;
158 
159         // Constructor for a polynomial from a vector of integers.
160         explicit Polynomial( const vector<ppuint> ) ;
161 
162         // Destructor.
163         virtual ~Polynomial() ;
164 
165         // Copy constructor.
166         Polynomial( const Polynomial & g ) ;
167 
168         // Assignment.
169         virtual Polynomial & operator=( const Polynomial & g ) ;
170 
171         // String to polynomial assignment.
172         virtual Polynomial & operator=( string s ) ;
173 
174         // Construct from string:
175         // Polynomial p( "x^2 + 2 x + 1, 3" ) ;
176         // If modulus isn't specified, use the one in specified in the
177         // polynomial string.
178         Polynomial( string s, ppuint p = 0 ) ;
179 
180         // Operator casting to string type.
181         operator string() const ;
182 
183         // Equality tests.
184         friend bool operator==( const Polynomial & p1, const Polynomial & p2 ) ;
185         friend bool operator!=( const Polynomial & p1, const Polynomial & p2 ) ;
186 
187         // cout << p and cin >> p
188         friend ostream & operator<<( ostream & out, const Polynomial & p ) ;
189         friend istream & operator>>( istream & in,        Polynomial & p ) ;
190 
191         // Bounds checked indexing operator which allows an lvalue:
192         //  p[ i ] = coeff ;
193         ppuint & operator[]( int i ) ;
194 
195         // Bounds checked indexing operator for read only access:
196         // coeff = p[ i ] ;
197         const ppuint operator[]( int i ) const ;
198 
199         // Return the degree n of f(x).
200         int deg() const ;
201 
202         // Return the modulus p of f(x).
203         ppuint modulus() const ;
204 
205         // Set the modulus p of f(x).
206         void setModulus( const ppuint p ) ;
207 
208         // Addition modulo p:  f(x) + g(x) mod p
209         friend const Polynomial operator+( const Polynomial & f, const Polynomial & g ) ;
210 
211         // Addition.
212         Polynomial & operator+=( const Polynomial & g ) ;
213 
214         // Scalar multiple.
215         friend const Polynomial
216            operator*( const Polynomial & f, const ppuint k ) ;
217 
218         Polynomial & operator*=( const ppuint k ) ;
219 
220         // Evaluation:  f( x ) mod p
221         ppuint operator()( int x ) ;
222 
223         // Does f(x) have any linear factor?
224         bool hasLinearFactor() ;
225 
226         // Is f(x) an integer?
227         bool isInteger() const ;
228 
229         // First trial polynomial.  Set
230         //                 n
231         //       f( x ) = x  - 1
232         void initialTrialPoly( const ppuint n, const ppuint p ) ;
233 
234         // Update f( x ) := next polynomial in sequence.
235         void nextTrialPoly() ;
236 
237     // Private data accessible by member functions only, and
238     // derived classes for convenience.
239     protected:
240 
241         vector<ppuint>     f_ ;   // Polynomial coefficients:
242                                   // f_[0] = const coeff ... f_[n] = nth degree coeff.
243                                   // then f_.size() = n+1
244         int                 n_ ;  // Degree of the polynomial.
245         ppuint              p_ ;  // Coefficients are in GF( p ).
246         ModP<ppuint,ppsint> mod ; // modulo p functionoid.
247 } ;
248 
249 /*=============================================================================
250 |
251 | NAME
252 |
253 |     PolyMod
254 |
255 | DESCRIPTION
256 |
257 |     Represents the monic polynomial g( x ) of degree m with coefficients in GF( p ),
258 |
259 |                   m         m-1
260 |         g( x ) = x  +  a   x  + ... + a    0 <= a  < p
261 |                         m-1            0         i
262 |
263 |     modulo f( x ) for monic polynomial f( x ) of degree n over GF( p ),
264 |
265 |                   n         n-1
266 |         f( x ) = x  +  a   x  + ... + a    0 <= a  < p
267 |                         n-1            0         i
268 |
269 |     We support these basic operations:
270 |     
271 |         PolyMod p() ;             Set g(x)=0 and f(x) = 0 mod 2
272 |                                   Destructor.   
273 |         PolyMod p( g, f )         Constructor from polynomials g(x) and f(x)
274 |         PolyMod p( p2 )           Copy p(x) = p2(x).
275 |         PolyMod p = p2            Assign p(x) = p2(x).
276 |         p.timesX() ;              p(x) := x p( x ) (mod f( x ), p)
277 |                                            2
278 |         p.square() ;              p(x) := p( x ) (mod f( x ), p)
279 |         p *= p2                   Do p(x) = p(x) * p2(x) (mod f( x ), p)
280 |         p1 * p2                   Do p1(x) * p2(x) (mod f( x ), p)
281 |                                           m 
282 |         p.power( m )              p(x) = x  (mod f(x), p)
283 |     
284 |     Exceptions:  throw PolynomialError or one of its subclasses such as 
285 |     PolynomialRangeError.
286 |
287 | NOTES
288 |
289 |     The member functions and friends are documented in ppPolynomial.h
290 |
291 +============================================================================*/
292 
293 // Friends of PolyMod
294 ppuint autoConvolve( const Polynomial & t, const int k, const int lower, const int upper ) ;
295 
296 ppuint coeffOfSquare( const Polynomial & g, const int k, const int n )  ;
297 
298 ppuint coeffOfProduct( const Polynomial & s, const Polynomial & t, const int k, const int n )   ;
299 
300 ppuint convolve( const Polynomial & s, const Polynomial & t, const int k, const int lower, const int upper ) ;
301 
302 class PolyMod
303 {
304     public:
305         // Set g( x ) = 0 modulo f(x) = 0 and p = 2
306         PolyMod() ;
307 
308         // Destructor.
309         virtual ~PolyMod() ;
310 
311         // Construct from polynomials g(x) and f(x).
312         PolyMod( const Polynomial & g, const Polynomial & f ) ;
313 
314         // Construct from string g and polynomial f(x).
315         PolyMod( const string & g, const Polynomial & f ) ;
316 
317         // Operator casting g(x) to string type.
318         operator string() const ;
319 
320          // cout << p prints g(x) to output stream
321         friend ostream & operator<<( ostream & out, const PolyMod & p ) ;
322 
323         // Copy g( x ) = g2( x ).
324         PolyMod( const PolyMod & g2 ) ;
325 
326         // Assign g( x ) = g2( x ) 
327         virtual PolyMod & operator=( const PolyMod & g2 ) ;
328 
329         // Bounds checked indexing operator for read only access:
330         // coeff = p[ i ] ;
331         const ppuint operator[]( int i ) const ;
332 
333         // Multiply by x:  g(x) := g(x) x (mod f( x ), p)
334         void timesX() ;
335 
336         // Squaring:             2
337         //           g(x) := g(x) (mod f( x ), p)
338         void square() ;
339 
340         // Multiplication:  g(x) := g(x) g2(x) (mod f( x ), p)
341         PolyMod & operator*=( const PolyMod & g2 ) ;
342 
343         // Multiplication:  g(x) := s(x) t(x) (mod f( x ), p)
344         friend const PolyMod operator*( const PolyMod & s, const PolyMod & t ) ;
345 
346         // Special exponentiation:  g(x) ^ m (mod f(x), p)
347         // for g(x) = x only for now!
348         friend const PolyMod power( const PolyMod & g, const BigInt & m ) ;
349 
350         bool isInteger() const ;
351 
352         //-----------------< Helper functions >-------------------------------
353         const Polynomial getf() const ;
354 
355         const ppuint getModulus() const ;
356 
357 
358     private:
359         // Polynomial g(x).
360         Polynomial g_ ;
361 
362         // Modulus polynomial f(x) and modulus p.
363         Polynomial f_ ;
364 
365         // Two dimensional precomputed power table.
366         //
367         //  Precompute the n-1 polynomials
368         //
369         //       n      2n-2
370         //      x  ... x     (mod f(x), p)
371         //
372         //                          n+i
373         //      powerTable_[ i ] = x   (mod f(x), p)
374         //
375         vector<Polynomial> powerTable_ ;
376 
377         ModP<ppuint,ppsint>  mod ; //  modulo p functionoid.
378 
379     // Helper functions.  Note the friend functions are really public due to C++ rules.
380     protected:
381         // Offset into powerTable.
382         int inline offset( const int i ) const { return i - f_.deg() ; }
383 
384         // Power table of the polynomial.
385         void constructPowerTable() ;
386 
387         // Reduce g( x ) mod f( x ) and p
388         void modf() ;
389 
390         // Autoconvolution product.
391         friend ppuint autoConvolve( const Polynomial & t, const int k, const int lower, const int upper )   ;
392 
393         // Convolution product.
394         friend ppuint convolve( const Polynomial & s, const Polynomial & t,
395                              const int k, const int lower, const int upper )  ;
396 
397         //               2
398         // kth coeff of g (x) for degree n.
399         friend ppuint coeffOfSquare( const Polynomial & g, const int k, const int n ) ;
400 
401         // Coeff of s(x) t(x) for degree n.
402         friend ppuint coeffOfProduct( const Polynomial & s, const Polynomial & t, const int k, const int n ) ;
403 } ;
404 
405 
406 
407 /*=============================================================================
408 |
409 | NAME
410 |
411 |     PolyOrder
412 |
413 | DESCRIPTION
414 |
415 |     Order tests on the monic polynomial f( x ) of degree n,
416 |     with coefficients in GF( p ).
417 |
418 |                   n         n-1
419 |         f( x ) = x  +  a   x  + ... + a    0 <= a  < p
420 |                         n-1            0         i
421 |
422 |     Exceptions:  throw PolynomialError or one of its subclasses such as
423 |     PolynomialRangeError.
424 |
425 | NOTES
426 |
427 |     The member functions and friends are documented in ppPolynomial.h
428 |
429 +============================================================================*/
430 
431 class PolyOrder
432 {
433     public:
434         // Do tests on the nth degree polynomial f(x) modulo p.
435         explicit PolyOrder( const Polynomial & f ) ;
436 
437         //               n
438         //              p  - 1
439         //  Compute r = -------- and factor r into the product of primes
440         //              p - 1
441         void factorR() ;
442 
443         void computeMaxNumPoly() ;
444 
445         void computeNumPrimPoly() ;
446 
447         void resetPolynomial( const Polynomial &f ) ;
448 
449         //             m                                          r
450         // Check that x  (mod f(x), p) is not an integer for m = ---
451         //                                                        p
452         //                                                         i
453         // but skip this test if p  | (p-1).
454         //                        i
455         bool orderM() ;
456 
457         //           r
458         // Check if x  (mod f(x), p) = a, for integer a.
459         // If a is not an integer, return 0.
460         ppuint orderR() ;
461 
462         //           k                                 n
463         // Check if x  = 1 (mod f(x), p) only for k = p - 1
464         // Note this is O( p^n ); very expensive.
465         bool maximal_order() ;
466 
467         // Check if the monic polynomial f( x ) has 2 or more distinct factors.
468         // Uses x_to_power().
469         bool hasMultipleDistinctFactors( bool earlyOut = true ) ;
470 
471         inline BigInt getR() const { return r_ ; } ;
472 
473         inline BigInt getNumPrimPoly() const { return num_prim_poly_ ; } ;
474 
475         inline BigInt getMaxNumPoly() const { return max_num_poly_ ; } ;
476 
477         // Test if a given polynomial f(x) is primitive.
478         bool isPrimitive() ;
479 
480         // Test function.
481         string printQMatrix() const ;
482 
483         inline int getNullity() const { return nullity_ ; } ;
484 
485         inline Factorization<BigInt> getFactorsOfR() const { return factors_of_R_ ; }
486 
487     // Allow direct access to this simple data type for convenience.
488     public:
489         OperationCount statistics_ ;
490 
491     protected:
492         // Polynomial f(x) modulo p of degree n which is to be tested.
493         Polynomial f_ ;
494         ppuint     n_ ;
495         ppuint     p_ ;
496         ModP<ppuint,ppsint> mod ;
497 
498         //  n
499         // p  - 1
500         BigInt p_to_n_minus_1_ ;
501 
502         //           n
503         //          p  - 1 
504         //   r =   -------
505         //          p  - 1
506         BigInt r_ ;
507 
508         // Constant factor a (see notes).
509         ppuint a_ ;
510 
511         //                   n
512         // Factorization of p  - 1
513         Factorization<BigInt> factors_of_p_to_n_minus_1_ ;
514 
515         // Factorization of r.
516         Factorization<BigInt> factors_of_R_ ;
517 
518         // Number of possible primitive polynomials.
519         BigInt num_prim_poly_ ;
520 
521         // Total number of possible polynomials.
522         BigInt max_num_poly_ ;
523 
524         // Two dimensional Q matrix for irreducibility testing.
525         vector< vector<ppsint> > Q_ ;
526 
527         int nullity_ ;
528 
529     // Helper functions. 
530     protected:
531         void generateQMatrix() ;
532 
533         void findNullity( bool earlyOut = true ) ;
534 
535 } ;
536 
537 #endif // __POLYNOMIAL_H__ --- End of wrapper for header file.