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.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 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 
  55 using namespace std ;
  56 
  57 
  58 /*------------------------------------------------------------------------------
  59 |                                PP Include Files                              |
  60 ------------------------------------------------------------------------------*/
  61 
  62 #include "Primpoly.h"         // Global functions.
  63 #include "ppArith.h"          // Basic arithmetic functions.
  64 #include "ppBigInt.h"         // Arbitrary precision integer arithmetic.
  65 #include "ppOperationCount.h" // OperationCount collection for factoring and poly finding.
  66 #include "ppFactor.h"         // Prime factorization and Euler Phi.
  67 #include "ppPolynomial.h"     // Polynomial operations and mod polynomial operations.
  68 #include "ppParser.h"         // Parsing of polynomials and I/O services.
  69 #include "ppUnitTest.h"       // 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 
  89 Polynomial::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 
 118 Polynomial::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 
 149 Polynomial::~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 
 179 Polynomial::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 
 215 Polynomial & 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 
 283 Polynomial & 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 
 335 Polynomial::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.
 420 Polynomial::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 
 541 ostream & 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 
 576 istream & 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 
 613 bool 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 
 627 bool 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 
 666 ppuint & 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 
 714 const 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 
 755 int Polynomial::deg() const
 756 {
 757     return n_ ;
 758 }
 759 
 760 // Return the modulus p of f(x).
 761 ppuint Polynomial::modulus() const
 762 {
 763     return p_ ;
 764 }
 765 
 766 void 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 
 801 Polynomial & 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 
 857 const 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 
 878 Polynomial & 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 
 901 const 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 
 945 ppuint
 946 Polynomial::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 
 994 bool
 995 Polynomial::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 
1030 bool
1031 Polynomial::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 
1065 void Polynomial::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 
1118 void Polynomial::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 
1155 PolyMod::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 
1181 PolyMod::~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 
1204 PolyMod::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 
1231 PolyMod::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.
1259 PolyMod::operator string() const
1260 {
1261 return 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 
1281 ostream & 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 
1307 const 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 
1329 const 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 
1351 PolyMod::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 
1372 PolyMod & 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 
1402 const 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 
1471 void 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 
1578 void 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 
1667 ppuint 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 
1747 ppuint 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 
1860 ppuint 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 
2013 ppuint 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 
2082 const 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 
2104 PolyMod &
2105 PolyMod::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 
2184 void 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 
2283 void 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 
2448 const 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 
2529 bool 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 
2552 void 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 
2570 PolyOrder::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 
2663 void 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 
2787 void 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 
2811 void 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 
2898 bool 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 
2966 ppuint 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 
3031 bool PolyOrder::maximal_order()
3032 {
3033     //  Highest possible order for x.
3034     BigInt maxOrder = power( f_.modulus(), f_.deg()) - static_cast<BigInt>( 1u ) ;
3035 
3036     BigInt k { 1u } ;
3037     Polynomial x1( "x", f_.modulus() ) ;
3038     PolyMod x( x1, f_ ) ;    // g(x) = x (mod f(x), p)
3039 
3040     while ( k <= maxOrder )
3041     {
3042         PolyMod x_to_power = power( x, k ) ; // x^k
3043 
3044         if (x_to_power.isInteger() &&
3045             x_to_power[0] == 1u &&
3046             k < maxOrder)
3047         {
3048             return false ;
3049         }
3050 
3051         ++k ;
3052 
3053     }
3054 
3055     return true ;
3056 }
3057 
3058 
3059 /*=============================================================================
3060  |
3061  | NAME
3062  |
3063  |     hasMultipleDistinctFactors
3064  |
3065  | DESCRIPTION
3066  |
3067  |    Returns true if the monic polynomial f( x ) has two or more distinct
3068  |    irreducible factors, false otherwise.
3069  |
3070  |    If earlyOut is false, compute the exact nullity in findNullity() instead
3071  |    of stopping when the nullity is >= 2.
3072  |
3073  | EXAMPLE
3074  |
3075  |    Let n = 4, p = 5
3076  |
3077  |              4    2
3078  |    f( x ) = x  + x  + 2 x + 3 is irreducible, so has one distinct factor.
3079  |
3080  |              4    3   2                  4
3081  |    f( x ) = x + 4x + x + 4x + 1 = (x + 1)  has one distinct factor.
3082  |
3083  |              3         2
3084  |    f( x ) = x  + 3 = (x + 3x + 4)(x + 2) has two distinct irreducible factors.
3085  |
3086  |              4    3    2                          2
3087  |    f( x ) = x + 3x + 3x + 3x + 2 = (x + 1) (x + 2) (x + 3) has 3 distinct
3088  |    irreducible factors.
3089  |
3090  |            4
3091  |    f(x) = x + 4 = (x+1)(x+2)(x+3)(x+4) has 4 distinct irreducible factors.
3092  |
3093  | METHOD
3094  |
3095  |       Berlekamp's method for factoring polynomials over GF( p ), modified to test
3096  |       for irreducibility only.
3097  |
3098  |       See my notes;  I skip the polynomial GCD step which ensures polynomials
3099  |       are square-free due to time constraints, but this requires a proof that
3100  |       the method is still valid.
3101  |
3102  +============================================================================*/
3103 
3104 bool PolyOrder::hasMultipleDistinctFactors( bool earlyOut )
3105 
3106 {
3107     // Generate the Q-I matrix.
3108     generateQMatrix() ;
3109 
3110     // Find nullity of Q-I
3111     findNullity( earlyOut ) ;
3112 
3113     // If nullity_ >= 2, f( x ) is a reducible polynomial modulo p since it has
3114     // two or more distinct irreducible factors.
3115     //                                     e
3116     // Nullity of 1 implies f( x ) = p( x )  for some power e >= 1 so we cannot
3117     // determine reducibility.
3118     if (nullity_ >= 2)
3119         return true ;
3120 
3121     return false ;
3122 
3123 }
3124 
3125 
3126 /*=============================================================================
3127  |
3128  | NAME
3129  |
3130  |     isPrimitive
3131  |
3132  | DESCRIPTION
3133  |
3134  |     Check if a given polynomial f(x) of degree n modulo p is primitive.
3135  |
3136  +============================================================================*/
3137 
3138 bool PolyOrder::isPrimitive()
3139 {
3140     bool isPrimitive = false ;
3141     ++statistics_.num_poly_tested ;
3142 
3143     try
3144     {
3145         BigInt max_num_possible_poly = power( p_, n_ ) ;
3146         ArithModP modp( p_ ) ;
3147 
3148         // Constant coefficient of f(x) * (-1)^n must be a primitive root of p.
3149         if (modp.constCoeffIsPrimitiveRoot( f_[0], f_.deg() ))
3150         {
3151             ++statistics_.num_where_const_coeff_is_primitive_root ;
3152 
3153             #ifdef DEBUG_PP_POLYNOMIAL
3154             cout << "    (-1)^n const coeff " << f_[ 0 ] << " is primitive root of " << p_ << endl ;
3155             #endif
3156 
3157             // f(x) can't have any linear factors.
3158             if (!f_.hasLinearFactor())
3159             {
3160                 ++statistics_.num_free_of_linear_factors ;
3161 
3162                 #ifdef DEBUG_PP_POLYNOMIAL
3163                 cout << "    No linear factors" << endl ;
3164                 #endif
3165 
3166                 // Do more in-depth checking.
3167 
3168                 // f(x) can't have two or more distinct irreducible factors.
3169                 if (!hasMultipleDistinctFactors())
3170                 {
3171                     ++statistics_.num_irreducible_to_power ;
3172 
3173                     #ifdef DEBUG_PP_POLYNOMIAL
3174                     cout << "    One distinct irreducible factor (possibly repeated)" << endl ;
3175                     #endif
3176 
3177                     //  r
3178                     // x  (mod f(x), p) = a_ must be an integer.
3179                     ppuint a{ orderR() } ;
3180                     if (a != 0)
3181                     {
3182                         ++statistics_.num_order_r ;
3183 
3184                         #ifdef DEBUG_PP_POLYNOMIAL
3185                         cout << "    x^r = a (integer)" << endl ;
3186                         #endif
3187 
3188                          //              n
3189                          // Check if (-1)  (constant coefficient of f(x)) = a_ (mod p)
3190                          //
3191                          if (modp.constCoeffTest( f_[ 0 ], a, f_.deg() ))
3192                          {
3193                              ++statistics_.num_passing_const_coeff_test ;
3194 
3195                              #ifdef DEBUG_PP_POLYNOMIAL
3196                              cout << "    a (integer) = (-1)^n f[0]" << endl ;
3197                              #endif
3198 
3199                              //  x^m != integer for all m = r / q, q a prime divisor of r.
3200                              if (orderM())
3201                              {
3202                                  ++statistics_.num_order_m ;
3203 
3204                                  #ifdef DEBUG_PP_POLYNOMIAL
3205                                  cout << "    x^m != integer for m = r / prime divisor of r" << endl ;
3206                                  #endif
3207 
3208                                  isPrimitive = true ;
3209                                  goto Exit ;
3210 
3211                              } // end order m
3212                          } // end const coeff test
3213                     } // end order r
3214                 } // end can't determine if reducible
3215             } // end no linear factors
3216         } // end constant coefficient primitive.
3217     }
3218     catch( ArithModPError & e )
3219     {
3220         ostringstream os ;
3221         os << "PolyOrder::isPrimitive had a mod p arithmetic error [ " << e.what() << " ] " ;
3222         throw PolynomialRangeError( os.str(), __FILE__, __LINE__ ) ;
3223     }
3224 
3225     Exit:
3226         return isPrimitive ;
3227 }
3228 
3229 
3230 /*=============================================================================
3231  |
3232  | NAME
3233  |     generateQMatrix
3234  |
3235  | DESCRIPTION
3236  |
3237  |     Generate n x n matrix Q - I, where rows of Q are the powers,
3238  |
3239  |         p                2p                      (n-1) p
3240  |     1, x  (mod f(x),p), x  (mod f(x), p), ... , x       (mod f(x), p)
3241  |
3242  |     for monic polynomial f(x).
3243  |
3244  | EXAMPLE
3245  |
3246  |             4   2
3247  |     f(x) = x + x + 2 x + 3, n = 4, p = 5
3248  |
3249  |         (    1     )    (    1              )       ( 1    0    0    0 )
3250  |         (          )    (                   )       (                  )
3251  |         (     5    )    (         2       3 )       (                  )
3252  |     Q = (    x     )    (  2x + 3x   + 4 x  )       ( 0    2    3    4 )
3253  |         (          )    (                   )       (                  )
3254  |         (          )    (         2     3   )       (                  )
3255  |         (     10   ) =  (  3 + 4 x   + x    )   =   (                  )
3256  |         (    x     )    (                   )       ( 3    0    4    1 )
3257  |         (          )    (                   )       (                  )
3258  |         (     15   )    (          2    3   )       (                  )
3259  |         (    x     )    (   4x + 4x + 3x    )       ( 0    4    4    3 )
3260  |         (          )    (                   )       (                  )
3261  |
3262  |                                                                 2    3
3263  |                                                       1   x    x    x
3264  |              ( 0 0 0 0 )
3265  |              ( 0 1 3 4 )
3266  |     Q - I =  ( 3 0 3 1 )
3267  |              ( 0 4 4 2 )
3268  |
3269  |
3270  |     The left nullspace has dimension = 1 with basis { (1 0 0 0) }.
3271  |
3272  +============================================================================*/
3273 
3274 void PolyOrder::generateQMatrix()
3275 {
3276     // Check for invalid inputs.
3277     if (n_ < 2 || p_ < 2)
3278         throw PolynomialRangeError( "generateQMatrix has n < 2 or p < 2", __FILE__, __LINE__ ) ;
3279 
3280     // Row 0 of Q = (1 0 ... 0).
3281     Q_[ 0 ][ 0 ] = 1 ;
3282     for (int i = 1 ;  i < n_ ;  ++i)
3283         Q_[ 0 ][ i ] = 0 ;
3284 
3285     //             p
3286     // Let q(x) = x (mod f(x),p)
3287     // and Q[ 1 ] = coefficients of q(x).
3288     Polynomial x1( "x", p_ ) ;
3289     PolyMod x( x1, f_ ) ;
3290     PolyMod xp = power( x, static_cast< BigInt >( p_ ) ) ;
3291 
3292     #ifdef DEBUG_PP_POLYNOMIAL
3293         cout << "x ^ p (mod f(x),p) = " << xp << endl ;
3294     cout << "initial Q matrix " << printQMatrix() ;
3295     #endif
3296 
3297     PolyMod q = xp ;
3298 
3299     for (int i = 0 ;  i < n_ ;  ++i)
3300         Q_[ 1 ][ i ] = xp[ i ] ;
3301 
3302     //               pk
3303     // Row k of Q = x   (mod f(x), p) 2 <= k <= n-1, computed by
3304     //                                   p
3305     // multiplying each previous row by x (mod f(x),p).
3306     for (int k = 2 ;  k <= n_- 1 ;  ++k)
3307     {
3308         q *= xp ;
3309 
3310         #ifdef DEBUG_PP_POLYNOMIAL
3311         cout << "x ^ pk (mod f(x),p) = " << q << " for k = " << k << endl ;
3312         #endif
3313 
3314         for (int i = 0 ;  i < n_ ;  ++i)
3315              Q_[ k ][ i ] = q[ i ] ;
3316     }
3317 
3318     #ifdef DEBUG_PP_POLYNOMIAL
3319     cout << "computed Q matrix " << printQMatrix() ;
3320     #endif
3321 
3322     //  Subtract Q - I
3323     for (int row = 0 ;  row < n_ ;  ++row)
3324     {
3325         Q_[ row ][ row ] = mod( Q_[ row ][ row ] - 1 ) ;
3326     }
3327 
3328     #ifdef DEBUG_PP_POLYNOMIAL
3329     cout << "computed Q-I matrix " << printQMatrix() ;
3330     #endif
3331 
3332     return ;
3333 }
3334 
3335 
3336 /*=============================================================================
3337  |
3338  | NAME
3339  |
3340  |    printQMatrix
3341  |
3342  | DESCRIPTION
3343  |
3344  |    Print the matrix to the console.
3345  |
3346  +============================================================================*/
3347 
3348 string PolyOrder::printQMatrix() const
3349 {
3350     // Print the matrix as a string.
3351     ostringstream os ;
3352 
3353     os << endl ;
3354     for (int row = 0 ;  row < n_ ;  ++row)
3355     {
3356         os << "( " ;
3357         for (int col = 0 ;  col < n_ ;  ++col)
3358         {
3359             os << setw( 4 ) << setfill( ' ' ) << Q_[ row ][ col ] ;
3360         }
3361         os << " )" << endl ;
3362     }
3363 
3364     return os.str() ;
3365 }
3366 
3367 
3368 /*=============================================================================
3369  |
3370  | NAME
3371  |
3372  |    findNullity
3373  |
3374  | DESCRIPTION
3375  |
3376  |     Computes the nullity of Q.  
3377  |     If earlyOut is true, stop when the nullity is >= 2 and return 2.
3378  |
3379  | EXAMPLE
3380  |
3381  |     Let p = 5 and n = 3.  We use the facts that -1 = 4 (mod 5), 1/2 = 3, -1/2 = 2,
3382  |     1/3 = 2, -1/3 = 3, 1/4 = 4, -1/4 = 1.
3383  |
3384  |     Consider the matrix
3385  |
3386  |         ( 2 3 4 )
3387  |     Q = ( 0 2 1 )
3388  |         ( 3 3 3 )
3389  |
3390  |     Begin with row 1.  No pivotal columns have been selected yet.  Scan row 1 and
3391  |     pick column 1 as the pivotal column because it contains a nonzero element.
3392  |
3393  |     Normalizing column 1 by multiplying by -1/pivot = -1/2 = 2 gives
3394  |
3395  |         ( 4 3 4 )
3396  |         ( 0 2 1 )
3397  |         ( 1 3 3 )
3398  |
3399  |      Now perform column reduction on column 2 by multiplying the pivotal column 1
3400  |      by 3 (the column 2 element in the current row) and adding back to row 2.
3401  |
3402  |         ( 4 0 4 )
3403  |         ( 0 2 1 )
3404  |         ( 1 1 3 )
3405  |
3406  |      Column reduce column 3 by multiplying the pivotal column by 4 and adding back to row 3,
3407  |
3408  |         ( 4 0 0 )
3409  |         ( 0 2 1 )
3410  |         ( 1 1 2 )
3411  |
3412  |      For row 2, pick column 2 as the pivotal column, normalize it and reduce columns 1, then 3,
3413  |
3414  |         ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )
3415  |         ( 0 2 1 ) => ( 0 4 1 ) => ( 0 4 1 ) => ( 0 4 0 )
3416  |         ( 1 1 2 )    ( 1 2 2 )    ( 1 2 2 )    ( 1 2 4 )
3417  |                  norm.         c.r. 1      c.r. 3
3418  |
3419  |      For row 3, we must pick column 3 as pivotal column because we've used up columns 1 and 2,
3420  |
3421  |         ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )
3422  |         ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 )
3423  |         ( 1 2 4 )    ( 1 2 4 )    ( 1 2 4 )    ( 0 0 4 )
3424  |                  norm.        c.r. 1        c.r. 2
3425  |
3426  |      The nullity is zero, since we were always able to find a pivot in each row.
3427  |
3428  | METHOD
3429  |
3430  |       Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth, Addison-Wesley.
3431  |
3432  |       We combine operations of normalization of columns,
3433  |
3434  |       (       c       )                         (        C       )
3435  |       (       c       )                         (        C       )
3436  |       (       .       )                         (        C       )
3437  |       ( . . . q . . . ) row  ================>  ( . . . -1 . . . ) row
3438  |       (       c       )                         (        C       )
3439  |       (       c       )      column times       (        C       )
3440  |       (       c       )      -1/a modulo p      (        C       )
3441  |            pivotCol                                   pivotCol
3442  |
3443  |       and column reduction,
3444  |
3445  |       (        C      b       )                         (       C        B       )
3446  |       (        C      b       )                         (       C        B       )
3447  |       (        C      b       )                         (       C        B       )
3448  |       ( . . . -1 . . .e . . . ) row  ================>  ( . . . -1 . . . 0 . . . )
3449  |       (        C      b       )                         (       C        B       )
3450  |       (        C      b       )    pivotCol times       (       C        B       )
3451  |       (        C      b       )    e added back to col  (       C        B       )
3452  |            pivotCol  col                                       col
3453  |
3454  |       to reduce the matrix to a form in which columns having pivots are zero until
3455  |       the pivotal row.
3456  |
3457  |       The column operations don't change the left nullspace of the
3458  |       matrix.
3459  |
3460  |       The matrix rank is the number of pivotal rows since they are linearly
3461  |       independent.  The nullity is then the number of non-pivotal rows.
3462  |
3463  +============================================================================*/
3464 
3465 void PolyOrder::findNullity( bool earlyOut )
3466 {
3467     try
3468     {
3469         InverseModP invmod( p_ ) ;
3470 
3471         #ifdef DEBUG_PP_POLYNOMIAL
3472         cout << "Q-I matrix " << printQMatrix() ;
3473         #endif
3474 
3475         vector<bool>pivotInCol( n_, false ) ; // Is false if the column has no pivotal element.
3476 
3477         nullity_ = 0 ;
3478         int pivotCol = -1 ; // No pivots yet.
3479 
3480         // Sweep through each row.
3481         for (int row = 0 ;  row < n_ ;  ++row)
3482         {
3483             // Search for a pivot in this row:  a non-zero element
3484             // in a column which had no previous pivot.
3485             bool found = false ;
3486             for (int col = 0 ;  col < n_ ;  ++col)
3487             {
3488                 if (Q_[ row ][ col ] > 0 && !pivotInCol[ col ])
3489                 {
3490                     found = true ;
3491                     pivotCol = col ;
3492                     break ;
3493                 }
3494             }
3495 
3496             // No pivot;  increase nullity by 1.
3497             if (found == false)
3498             {
3499                 ++nullity_ ;
3500 
3501                 // Early out.
3502                 if (earlyOut && nullity_ >= 2)
3503                     goto EarlyOut ;
3504             }
3505             // Found a pivot, q.
3506             else
3507             {
3508                 ppsint q = Q_[ row ][ pivotCol ] ;
3509 
3510                 // Compute -1/q (mod p)
3511                 ppsint t = mod( -invmod( q )) ;
3512 
3513                 // Normalize the pivotal column.
3514                 for (int r = 0 ;  r < n_ ;  ++r)
3515                 {
3516                     Q_[ r ][ pivotCol ] = mod( t * Q_[ r ][ pivotCol ]) ;
3517                 }
3518 
3519                 // Do column reduction:  Add C times the pivotal column to the other
3520                 // columns where C = element in the other column at current row.
3521                 for (int col = 0 ;  col < n_ ;  ++col)
3522                 {
3523                     if (col != pivotCol)
3524                     {
3525                         q = Q_[ row ][ col ] ;
3526 
3527                         for (int r = 0 ;  r < n_ ;  ++r)
3528                         {
3529                             t = mod( q * Q_[ r ][ pivotCol ]) ;
3530                             Q_[ r ][ col ] = mod( t + Q_[ r ][ col ] ) ;
3531                         }
3532                     }
3533                 }
3534 
3535                 // Record the presence of a pivot in this column.
3536                 pivotInCol[ pivotCol ] = true ;
3537 
3538                 #ifdef DEBUG_PP_POLYNOMIAL
3539                 cout << "row = " << row << " pivot = " << q << " (-1/q) = " << t << " nullity = " << nullity_
3540                      << " Row reduced Q-I matrix " << printQMatrix() ;
3541                 #endif
3542 
3543             } // found a pivot
3544 
3545         } // end for row
3546 
3547         EarlyOut: ;
3548             #ifdef DEBUG_PP_POLYNOMIAL
3549             cout << "Row reduced Q-I matrix " << printQMatrix() ;
3550             #endif
3551     }
3552     catch( ArithModPError & e )
3553     {
3554         ostringstream os ;
3555         os << "PolyOrder::findNullity failed in matrix row reduction [ " << e.what() << " ] " ;
3556         throw PrimpolyError( os.str(), __FILE__, __LINE__ ) ;
3557     }
3558 
3559     // Automagically free pivotInCol and mod objects.
3560 
3561 } // ===================== end of function findNullity =====================