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 =====================