1 /*==============================================================================
   2 | 
   3 |  File Name:     
   4 |
   5 |     ppPolyArith.c
   6 |
   7 |  Description:   
   8 |
   9 |     Polynomial arithmetic and exponentiation routines.
  10 | 
  11 |  Functions:
  12 |
  13 |     eval_poly
  14 |     linear_factor
  15 |     is_integer
  16 |     construct_power_table
  17 |     auto_convolve
  18 |     convolve
  19 |     coeff_of_square
  20 |     coeff_of_product
  21 |     square
  22 |     product
  23 |     times_x
  24 |     x_to_power
  25 |
  26 |  LEGAL
  27 |
  28 |     Primpoly Version 16.1 - A Program for Computing Primitive Polynomials.
  29 |     Copyright (C) 1999-2021 by Sean Erik O'Connor.  All Rights Reserved.
  30 |
  31 |     This program is free software: you can redistribute it and/or modify
  32 |     it under the terms of the GNU General Public License as published by
  33 |     the Free Software Foundation, either version 3 of the License, or
  34 |     (at your option) any later version.
  35 |
  36 |     This program is distributed in the hope that it will be useful,
  37 |     but WITHOUT ANY WARRANTY; without even the implied warranty of
  38 |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  39 |     GNU General Public License for more details.
  40 |
  41 |     You should have received a copy of the GNU General Public License
  42 |     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  43 |     
  44 |     The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com
  45 |     with !DOT! replaced by . and the !AT! replaced by @
  46 |
  47 ==============================================================================*/
  48 
  49 /*------------------------------------------------------------------------------
  50 |                                Include Files                                 |
  51 ------------------------------------------------------------------------------*/
  52 
  53 #include <stdio.h>
  54 #include "Primpoly.h"
  55 
  56 
  57 /*==============================================================================
  58 |                                   eval_poly                                  |
  59 ================================================================================
  60 
  61 DESCRIPTION
  62 
  63      Evaluate the polynomial f( x ) with modulo p arithmetic.
  64 
  65 INPUT
  66 
  67      f (int *)      nth degree mod p polynomial
  68          
  69                               n         n-1     
  70                     f( x ) = x  +  a   x  + ... + a    0 <= a  < p
  71                                     n-1            0         i          
  72      x (int)        0 <= x < p
  73 
  74      n (int)        n >= 1
  75 
  76      p (int)         
  77 
  78 RETURNS
  79 
  80      f( x )  (int)
  81 
  82 EXAMPLE 
  83                                   4
  84      Let n = 4, p = 5 and f(x) = x  + 3x + 3.  Then f(x) = 
  85      (((x + 0)x + 0)x + 3)x + 3.  f(2) = (((2 + 0)2 + 0)2 + 3) =
  86      (8 + 3)2 + 3 = 1 + 2 + 3 (mod 5) = 0.
  87 
  88 METHOD
  89 
  90      By Horner's rule, f(x) = ( ... ( (x + a   )x + ... )x + a .
  91                                             n-1               0
  92      We evaluate recursively, f := f * x + a (mod p), starting
  93                                             i 
  94      with f = 1 and i = n-1.
  95 
  96 BUGS
  97 
  98      None.
  99 
 100 --------------------------------------------------------------------------------
 101 |                                Function Call                                 |
 102 ------------------------------------------------------------------------------*/
 103 
 104 int
 105     eval_poly( int * f, int x, int n, int p )
 106 {
 107 
 108 /*------------------------------------------------------------------------------
 109 |                               Local Variables                                |
 110 ------------------------------------------------------------------------------*/
 111 
 112 int
 113     val = 1 ;   /*  The value of f(x) which is returned. */
 114 
 115 /*------------------------------------------------------------------------------
 116 |                                Function Body                                 |
 117 ------------------------------------------------------------------------------*/
 118 
 119 while ( --n >= 0 )
 120 
 121     val = mod( val * x + f[ n ], p ) ;
 122 
 123 return( val ) ;
 124 
 125 
 126 } /* ========================= end of function eval_poly ==================== */
 127 
 128 
 129 /*==============================================================================
 130 |                                     linear_factor                            |
 131 ================================================================================
 132 
 133 DESCRIPTION
 134 
 135      Check if the polynomial f(x) has any linear factors.
 136 
 137 INPUT
 138 
 139     f (int *)  nth degree monic mod p polynomial f(x). 
 140 
 141     n (int)    Its degree.
 142 
 143     p (int)    Modulus for coefficient arithmetic.
 144 
 145 RETURNS
 146 
 147     YES    if f( a ) = 0 (mod p) for a = 1, 2, ... p-1.
 148     NO     otherwise
 149 
 150     i.e. check if f(x) contains a linear factor (x - a).  We don't need to test
 151     for the root a = 0 because f(x) was chosen in main to have a non-zero 
 152     constant term, hence no zero root.
 153 
 154 EXAMPLE 
 155                                  4      
 156     Let n = 4, p = 5 and f(x) = x  + 3x + 3.  Then f(1) = 2 (mod 5), but 
 157     f(2) = 0 (mod 5), so we exit immediately with a yes.
 158 
 159 METHOD
 160 
 161     Evaluate f(x) at x = 1, ..., p-1 by Horner's rule.  Return instantly the 
 162     moment f(x) evaluates to 0.
 163 
 164 BUGS
 165 
 166     None.
 167 
 168 --------------------------------------------------------------------------------
 169 |                                Function Call                                 |
 170 ------------------------------------------------------------------------------*/
 171 
 172 int
 173     linear_factor( int * f, int n, int p )
 174 {
 175 
 176 /*------------------------------------------------------------------------------
 177 |                               Local Variables                                |
 178 ------------------------------------------------------------------------------*/
 179 
 180 int i ;  /*  Loop counter. */
 181 
 182 /*------------------------------------------------------------------------------
 183 |                                Function Body                                 |
 184 ------------------------------------------------------------------------------*/
 185 
 186 for (i = 1 ;  i <= p-1 ;  ++i)
 187 
 188     if (eval_poly( f, i, n, p ) == 0)
 189 
 190         return( YES ) ;
 191 
 192 return( NO ) ;
 193 
 194 } /* ====================== end of function linear_factor =================== */
 195 
 196 
 197 /*==============================================================================
 198 |                                 is_integer                                   |
 199 ================================================================================
 200 
 201 DESCRIPTION
 202 
 203      Test if a polynomial is a constant.
 204 
 205 INPUT
 206 
 207     t (int *)  mod p polynomial t(x). 
 208     n (int)    Its degree.
 209 
 210 RETURNS
 211 
 212     YES    if t(x) is a constant polynomial.
 213     NO     otherwise
 214 
 215 EXAMPLE 
 216                                2
 217     Let n = 3.  Then t(x) = 2 x  is not constant but t(x) = 2 is.
 218 
 219 METHOD
 220 
 221     A constant polynomial is zero in its first through n th degree
 222     terms.  Return immediately with NO if any such term is non-zero.
 223 
 224 BUGS
 225 
 226     None.
 227 
 228 --------------------------------------------------------------------------------
 229 |                                Function Call                                 |
 230 ------------------------------------------------------------------------------*/
 231 
 232 int
 233     is_integer( int * t, int n )
 234 {
 235 
 236 /*------------------------------------------------------------------------------
 237 |                                Function Body                                 |
 238 ------------------------------------------------------------------------------*/
 239 
 240 while ( n >= 1 )
 241 
 242     if (t[ n-- ] != 0)
 243 
 244         return( NO ) ;
 245 
 246 return( YES ) ;
 247 
 248 } /* ====================== end of function is_integer ====================== */
 249 
 250 
 251 /*==============================================================================
 252 |                               construct_power_table                          |
 253 ================================================================================
 254 
 255 DESCRIPTION
 256 
 257      Construct a table of powers of x:
 258 
 259       n                     2n-2
 260      x  (mod f(x), p)  ... x    (mod f(x), p)
 261 
 262 INPUT
 263 
 264     f (int *)   Coefficients of f(x), a monic polynomial of degree n.
 265     n (int, -infinity < n < infinity)
 266     p (int, p > 0)
 267 
 268 RETURNS
 269 
 270     power_table (int *)  power_table[i][j] is the coefficient of 
 271      j       n+i
 272     x   in  x   (mod f(x), p) where 0 <= i <= n-2 and 0 <= j <= n-1. 
 273 
 274 EXAMPLE 
 275                                   4     2                     4
 276      Let n = 4, p = 5 and f(x) = x  +  x  +  2x  +  3.  Then x  = 
 277 
 278          2                  2
 279      -( x  +  2x  + 3) = 4 x  + 3 x + 2 (mod f(x), 5), and we get
 280 
 281       4                    2
 282      x  (mod f(x), 5) = 4 x  + 3 x + 2 = power_table[0].
 283 
 284       5                       2                 3      2
 285      x  (mod f(x), 5) = x( 4 x  + 3 x + 2) = 4 x  + 3 x  + 2x 
 286                       = power_table[1].
 287 
 288       6                       3      2           4      3      2
 289      x  (mod f(x), 5) = x( 4 x  + 3 x + 2 x) = 4x  + 3 x  + 2 x 
 290 
 291                               2                 3      2
 292                       = 4 ( 4x  + 3 x + 2) + 3 x  + 2 x  =
 293 
 294                            3     2
 295                       = 3 x + 3 x + 2 x + 3 = power_table[2].
 296 
 297                                     j
 298      power_table[i][j]:       | 0  1  2  3
 299                            ---+-------------
 300                             0 | 2  3  4  0 
 301                         i   1 | 0  2  3  4 
 302                             2 | 3  2  3  3 
 303 
 304 METHOD
 305                               n-1
 306      Beginning with t( x ) = x,    compute the next power of x from the last
 307                                                          n
 308      one by multiplying by x.  If necessary, substitute x  =
 309              n-1
 310      -( a   x    + ... + a ) to reduce the degree.  This formula comes from
 311          n-1              0
 312                                n         n-1
 313      the observation f( x ) = x   + a   x    + ... + a    = 0 (mod f(x), p).
 314                                      n-1              0
 315 BUGS
 316 
 317      None.
 318 
 319 --------------------------------------------------------------------------------
 320 |                                Function Call                                 |
 321 ------------------------------------------------------------------------------*/
 322 
 323 void
 324     construct_power_table( int power_table[][ MAXDEGPOLY ], int * f,
 325                            int n, int p )
 326 {
 327 
 328 /*------------------------------------------------------------------------------
 329 |                               Local Variables                                |
 330 ------------------------------------------------------------------------------*/
 331 
 332 int
 333     i, j,                  /*  Loop counters.  */
 334     coeff,                 /*  Coefficient of x ^ n in t(x) */
 335     t[ MAXDEGPOLY + 1 ] ;  /*  t(x) is temporary storage for x ^ k (mod f(x),p)
 336                                n <= k <= 2n-2.  Its degree can go as high as
 337                                n before it is reduced again. */
 338 
 339 /*------------------------------------------------------------------------------
 340 |                                Function Body                                 |
 341 ------------------------------------------------------------------------------*/
 342 
 343 /*
 344                          n-1
 345     Initialize t( x ) = x 
 346 */
 347 
 348 for (i = 0 ;  i <= n-2 ;  ++i)
 349 
 350     t[ i ] = 0 ;
 351 
 352 t[ n-1 ] = 1 ;
 353 
 354 
 355 
 356 /*
 357                                         i+n
 358     Fill the ith row of the table with x   (mod f(x), p).
 359 */
 360 
 361 for (i = 0 ;  i <= n-2 ;  ++i)
 362 {
 363 
 364     /*  Compute t(x) = x t(x)  */
 365     for (j = n ;  j >= 1 ;  --j)
 366 
 367         t[ j ] = t[ j-1 ] ;
 368 
 369     t[ 0 ] = 0 ;
 370 
 371 
 372     /*  Coefficient of the nth degree term of t(x).  */
 373     if ( (coeff = t[ n ]) != 0)
 374     {
 375         t[ n ] = 0 ;               /*  Zero out the x ^ n th term. */
 376 
 377        /*
 378                  n       n                        n-1
 379         Replace x  with x  (mod f(x), p) = -(a   x   + ... + a )
 380                                               n-1             0
 381        */
 382 
 383         for (j = 0 ;  j <= n-1 ;  ++j)
 384 
 385             t[ j ] = mod( t[ j ] +
 386                           mod( -coeff * f[ j ], p ),
 387                           p ) ;
 388     }  /* end if */
 389 
 390 
 391     /*  Copy t(x) into row i of power_table. */
 392     for (j = 0 ;  j <= n-1 ; ++j)
 393 
 394         power_table[i][j] = t[ j ] ;
 395 
 396 } /* end for */
 397 
 398 return ;
 399 
 400 } /* ================== end of function construct_power_table =============== */
 401 
 402 
 403 
 404 /*==============================================================================
 405 |                                 autoconvolve                                  |
 406 ================================================================================
 407 
 408 DESCRIPTION
 409 
 410      Compute a convolution-like sum for use in function coeff_of_square.
 411 
 412 INPUT
 413                                              n-1
 414     t (int *)   Coefficients of t(x) = t    x    + ... + t x  + t
 415                                         n-1               1      0
 416     k (int), 0 <= k <= 2n - 2)
 417 
 418     lower (int, 0 <= lower <= n-1) 
 419 
 420     uppper (int, 0 <= upper <= n-1) 
 421 
 422     p (int, p > 0)
 423 
 424 RETURNS
 425 
 426      upper
 427      ---
 428      \    t  t       But define the sum to be 0 when lower > upper to catch
 429      /     i  k-i    the special cases that come up in function coeff_of_square.
 430      ---
 431      i=lower
 432 
 433 EXAMPLE 
 434                        3      2      
 435      Suppose t(x) = 4 x  +  x  +  3 x  +  3, lower = 1, upper = 3, n = 3,
 436 
 437      and p = 5.  For k = 3, auto_convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +
 438 
 439      t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3.  For lower = 0,
 440 
 441      upper = -1, or for lower = 3 and upper = 2, auto_convolve = 0, no matter what
 442 
 443      k is.
 444 
 445 METHOD
 446 
 447     A "for" loop in the C language is not executed when its lower limit 
 448 
 449     exceeds its upper limit, leaving sum = 0.
 450 
 451 BUGS
 452 
 453      None.
 454 
 455 --------------------------------------------------------------------------------
 456 |                                Function Call                                 |
 457 ------------------------------------------------------------------------------*/
 458 
 459 int
 460     auto_convolve( int * t, int k, int lower, int upper, int p )
 461 {
 462 
 463 /*------------------------------------------------------------------------------
 464 |                               Local Variables                                |
 465 ------------------------------------------------------------------------------*/
 466 
 467 int
 468     sum = 0,      /* Convolution sum. */
 469     i ;           /* Loop counter.    */
 470 
 471 /*------------------------------------------------------------------------------
 472 |                                Function Body                                 |
 473 ------------------------------------------------------------------------------*/
 474 
 475 for (i = lower ;  i <= upper ;  ++i)
 476 
 477    sum = mod( sum + mod( t[ i ] * t[ k - i ], p ), p ) ;
 478 
 479 return( sum ) ;
 480 
 481 } /* ==================== end of function auto_convolve ===================== */
 482 
 483 
 484 /*==============================================================================
 485 |                                 convolve                                     |
 486 ================================================================================
 487 
 488 DESCRIPTION
 489 
 490      Compute a convolution-like sum.
 491 
 492 INPUT
 493                                              n-1
 494     s (int *)   Coefficients of s(x) = s    x    + ... + s x  + s
 495                                         n-1               1      0  
 496                                              n-1
 497     t (int *)   Coefficients of t(x) = t    x    + ... + t x  + t
 498                                         n-1               1      0
 499 
 500     k (int), 0 <= k <= 2n - 2)
 501 
 502     lower (int, 0 <= lower <= n-1) 
 503 
 504     uppper (int, 0 <= upper <= n-1) 
 505 
 506     p (int, p > 0)
 507 
 508 RETURNS
 509 
 510      upper
 511      ---
 512      \    s  t       But define the sum to be 0 when lower > upper to catch
 513      /     i  k-i    the special cases
 514      ---
 515      i=lower
 516 
 517 EXAMPLE 
 518                        3      2      
 519      Suppose s(x) = 4 x  +  x  +  3 x  +  3, 
 520 
 521                         3      2      
 522      Suppose t(x) = 4 x  +  x  +  3 x  +  3, 
 523 
 524        
 525      lower = 1, upper = 3, n = 3,
 526 
 527      and p = 5.  For k = 3, convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +
 528 
 529      t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3.  For lower = 0,
 530 
 531      upper = -1, or for lower = 3 and upper = 2, convolve = 0, no matter what
 532 
 533      k is.
 534 
 535 METHOD
 536 
 537     A "for" loop in the C language is not executed when its lower limit 
 538 
 539     exceeds its upper limit, leaving sum = 0.
 540 
 541 BUGS
 542 
 543      None.
 544 
 545 --------------------------------------------------------------------------------
 546 |                                Function Call                                 |
 547 ------------------------------------------------------------------------------*/
 548 
 549 int
 550     convolve( int * s, int * t, int k, int lower, int upper, int p )
 551 {
 552 
 553 /*------------------------------------------------------------------------------
 554 |                               Local Variables                                |
 555 ------------------------------------------------------------------------------*/
 556 
 557 int
 558     sum = 0,      /* Convolution sum. */
 559     i ;           /* Loop counter.    */
 560 
 561 /*------------------------------------------------------------------------------
 562 |                                Function Body                                 |
 563 ------------------------------------------------------------------------------*/
 564 
 565 for (i = lower ;  i <= upper ;  ++i)
 566 
 567    sum = mod( sum + mod( s[ i ] * t[ k - i ], p ), p ) ;
 568 
 569 return( sum ) ;
 570 
 571 } /* ======================= end of function convolve ======================= */
 572 
 573 
 574 /*==============================================================================
 575 |                               coeff_of_square                                |
 576 ================================================================================
 577 
 578 DESCRIPTION
 579                                     th                2
 580      Return the coefficient of the k   power of x in t( x )  modulo p.
 581  
 582 INPUT
 583 
 584     t (int *)  Coefficients of t(x), of degree <= n-1.
 585 
 586     k (int, 0 <= k <= 2n-2)
 587 
 588     n (int, -infinity < n < infinity)
 589 
 590     p (int, p > 0)
 591 
 592 RETURNS
 593                         th                2
 594     coefficient of the k   power of x in t(x) mod p.
 595 
 596 EXAMPLE 
 597                                     3     2           2
 598     Let n = 4, p = 5, and t(x) = 4 x  +  x  + 4.  t(x) = 
 599 
 600              0                   1                         2
 601     (4 * 4) x  +  (2 * 1 * 4) x  +  (2 * 1 * 4 + 1 * 1) x  +
 602 
 603                              3                        4
 604     (2 * 1 * 1 + 2 * 4 * 4) x  + (2 * 4 * 1 + 1 * 1) x  +
 605 
 606                  5            6     
 607     (2 * 4 * 1) x  + (4 * 4) x  =  
 608 
 609      6      5     4     3     2
 610     x  + 3 x + 4 x + 4 x + 4 x  + 3 x + 1, all modulo 5.
 611 
 612             k        |  0  1  2  3  4  5  6
 613      ----------------+---------------------
 614      coeff_of_square |  1  3  4  4  4  3  1
 615 
 616 METHOD
 617                                                          2
 618     The formulas were gotten by writing out the product t (x) explicitly.
 619 
 620     The sum is 0 in two cases:  
 621 
 622         (1) when k = 0 and the limits of summation are 0 to -1 
 623 
 624         (2) k = 2n - 2, when the limits of summation are n to n-1.
 625 
 626     To derive the formulas, let 
 627 
 628                       n-1           
 629     Let t (x) = t    x     +  ... + t x + t
 630                  n-1                 1     0
 631 
 632     Look at the formulas in coeff_of_product for each power of x,
 633     replacing s with t, and observe that half of the terms are
 634     duplicates, so we can save computation time.
 635 
 636     Inspection yields the formulas,
 637 
 638 
 639     for 0 <= k <= n-1, even k,
 640 
 641      k/2-1
 642       ---             2
 643    2  \   t  t     + t
 644       /    i  k-i     k/2
 645       ---
 646       i=0
 647 
 648     for 0 <= k <= n-1, odd k,
 649 
 650       (k-1)/2
 651       ---  
 652     2 \   t  t    
 653       /    i  k-i
 654       ---
 655       i=0
 656 
 657 
 658     and for n <= k <= 2n-2, even k,
 659 
 660       n-1
 661       ---            2
 662    2  \   t  t    + t
 663       /    i  k-i    k/2
 664       ---
 665       i=k/2+1
 666 
 667       and for n <= k <= 2n-2, odd k,
 668 
 669       n-1
 670       ---  
 671    2  \   t  t    
 672       /    i  k-i   
 673       ---
 674       i=(k+1)/2
 675 
 676 
 677 
 678 BUGS
 679 
 680      None.
 681 
 682 --------------------------------------------------------------------------------
 683 |                                Function Call                                 |
 684 ------------------------------------------------------------------------------*/
 685 
 686 int
 687     coeff_of_square( int * t, int k, int n, int p )
 688 {
 689 
 690 /*------------------------------------------------------------------------------
 691 |                               Local Variables                                |
 692 ------------------------------------------------------------------------------*/
 693 
 694 int
 695     sum = 0 ;      /* kth coefficient of t(x) ^ 2. */
 696 
 697 /*------------------------------------------------------------------------------
 698 |                                Function Body                                 |
 699 ------------------------------------------------------------------------------*/
 700 
 701 if (0 <= k && k <= n-1)
 702 {
 703 
 704     if (k % 2 == 0)        /* Even k */
 705 
 706         sum = mod( mod( 2 * auto_convolve( t, k, 0, k/2 - 1, p ), p) +
 707                    t[ k/2 ] * t[ k/2 ], p ) ;
 708 
 709      else                  /* Odd k */
 710 
 711          sum = mod( 2 * auto_convolve( t, k, 0, (k-1)/2, p ), p) ;
 712 }
 713 else if (n <= k && k <= 2 * n - 2)
 714 {
 715 
 716     if (k % 2 == 0)        /* Even k */
 717 
 718         sum = mod( mod( 2 * auto_convolve( t, k, k/2 + 1, n-1, p ), p) +
 719                    t[ k/2 ] * t[ k/2 ], p ) ;
 720 
 721      else                  /* Odd k */
 722 
 723          sum = mod( 2 * auto_convolve( t, k, (k+1)/2, n-1, p ), p) ;
 724 }
 725 
 726 return( sum ) ;
 727 
 728 } /* ================== end of function coeff_of_square ===================== */
 729 
 730 
 731 /*==============================================================================
 732 |                               coeff_of_product                               |
 733 ================================================================================
 734 
 735 DESCRIPTION
 736                                     th
 737      Return the coefficient of the k   power of x in s( x ) t( x )  modulo p.
 738  
 739 INPUT
 740 
 741     t (int *)  Coefficients of t(x), of degree <= n-1.
 742 
 743     k (int, 0 <= k <= 2n-2)
 744 
 745     n (int, -infinity < n < infinity)
 746 
 747     p (int, p > 0)
 748 
 749 RETURNS
 750                         th
 751     coefficient of the k   power of x in s(x) t(x) mod p.
 752 
 753 EXAMPLE  
 754                               3     2                  2
 755   Let n = 4, p = 5, t(x) = 4 x  +  x  + 4, s( x ) = 3 x  + x + 2
 756 
 757   We'll do the case k=3,
 758 
 759   t3 s0 + t2 s1 + t1 s2 + t0 s3 = 4 * 2 + 1 * 1 + 0 * 3 + 4 * 0 = 9 = 4 (mod 5).
 760 
 761 
 762             k         |  0  1  2  3  4  5  6
 763      -----------------+---------------------
 764      coeff_of_product |  3  4  4  4  2  2  0
 765 
 766 METHOD
 767                                                          
 768     The formulas were gotten by writing out the product s(x) t (x) explicitly.
 769 
 770     The sum is 0 in two cases:  
 771 
 772         (1) when k = 0 and the limits of summation are 0 to -1 
 773 
 774         (2) k = 2n - 2, when the limits of summation are n to n-1.
 775 
 776 
 777     To derive the formulas, let 
 778 
 779                       n-1           
 780     Let s (x) = s    x     +  ... + s x + s  and
 781                  n-1                 1     0
 782 
 783                       n-1           
 784         t (x) = t    x     +  ... + t x + t 
 785                  n-1                 1     0
 786 
 787     and multiply out the terms, collecting like powers of x:
 788 
 789 
 790     Power of x     Coefficient
 791     ==========================
 792      2n-2
 793     x              s    t
 794                     n-1  n-1
 795 
 796      2n-3
 797     x              s    t    +  s    t
 798                     n-2  n-1     n-1  n-2     
 799 
 800      2n-4
 801     x              s    t    +  s    t    +  s    t
 802                     n-3  n-1     n-2  n-2     n-3  n-1
 803 
 804      2n-5
 805     x              s    t    +  s    t    +  s    t    +  s    t
 806                     n-4  n-1     n-3  n-2     n-2  n-3     n-1  n-4
 807 
 808      . . .
 809 
 810      n
 811     x              s  t    +  s  t    + ...  +  s    t
 812                     1  n-1     2  n-2            n-1  1
 813 
 814      n-1
 815     x              s  t    +  s  t    + ...  +  s    t
 816                     0  n-1     1  n-2            n-1  0
 817 
 818     . . .
 819 
 820      3
 821     x              s  t  +  s  t  +  s  t  +  s  t
 822                     0  3     1  2     2  1     3  0
 823 
 824      2
 825     x              s  t  +  s  t  +  s  t 
 826                     0  2     1  1     2  0
 827 
 828      
 829     x              s  t  +  s  t 
 830                     0  1     1  0 
 831 
 832     1              s  t
 833                     0  0
 834 
 835 
 836     Inspection yields the formulas,
 837 
 838 
 839     for 0 <= k <= n-1,
 840 
 841       k
 842      ---
 843      \   s  t
 844      /    i  k-i
 845      ---
 846      i=0
 847 
 848 
 849     and for n <= k <= 2n-2,
 850 
 851      n-1
 852      ---
 853      \   s  t
 854      /    i  k-i
 855      ---
 856      i=k-n+1
 857 
 858 BUGS
 859 
 860      None.
 861 
 862 --------------------------------------------------------------------------------
 863 |                                Function Call                                 |
 864 ------------------------------------------------------------------------------*/
 865 
 866 int
 867     coeff_of_product( int * s, int * t, int k, int n, int p )
 868 {
 869 
 870 /*------------------------------------------------------------------------------
 871 |                               Local Variables                                |
 872 ------------------------------------------------------------------------------*/
 873 
 874 int
 875     sum = 0 ;      /* kth coefficient of t(x) ^ 2. */
 876 
 877 /*------------------------------------------------------------------------------
 878 |                                Function Body                                 |
 879 ------------------------------------------------------------------------------*/
 880 
 881 if (0 <= k && k <= n-1)
 882 {
 883     sum = convolve( s, t, k, 0, k, p ) ;
 884 }
 885 else if (n <= k && k <= 2 * n - 2)
 886 {
 887     sum = convolve( s, t, k, k - n + 1, n - 1, p ) ;
 888 }
 889 
 890 return( sum ) ;
 891 
 892 } /* ================== end of function coeff_of_product ==================== */
 893 
 894 
 895 /*==============================================================================
 896 |                                  square                                      |
 897 ================================================================================
 898 
 899 DESCRIPTION
 900               2
 901      Compute t ( x ) (mod f(x), p)
 902 
 903      Uses a precomputed table of powers of x.
 904 
 905 INPUT
 906 
 907     t (int *)              Coefficients of t (x), of degree <= n-1.
 908 
 909     power_table (int **)   x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.
 910 
 911     n (int, -infinity < n < infinity)  Degree of f(x).
 912 
 913     p (int, p > 0)         Mod p coefficient arithmetic.
 914 
 915 OUTPUT
 916                                              2
 917     t (int *)              Overwritten with t (x) (mod f(x), p)
 918 
 919 EXAMPLE 
 920                                     3     2        
 921     Let n = 4, p = 5, and t(x) = 4 x  +  x  + 4.  Let f(x) = 
 922 
 923      4    2                   2       6      5      4      3     2
 924     x  + x  + 2 x + 3.  Then t (x) = x  + 3 x  + 4 x  + 4 x + 4 x + 3 x
 925 
 926            2                            3     2
 927     + 1.  t (x) (mod f(x), p) = 1 * (3 x + 3 x + 2 x + 3) +
 928 
 929             3      2                 2                3     2
 930     3 * (4 x  + 3 x + 2 x) + 4 * (4 x + 3 x + 2) + 4 x + 4 x + 3 x + 1
 931 
 932          3     2
 933     = 4 x + 2 x + 3 x + 2.
 934 
 935 
 936 METHOD
 937          2            2n-2              n         n-1
 938     Let t (x) = t    x     +  ... + t  x  +  t   x   +  ... + t .
 939                  2n-2                n        n-1              0
 940 
 941     Compute the coefficients t  using the function coeff_of_square.
 942                               k
 943 
 944                                 2
 945     The next step is to reduce t (x) modulo f(x).  To do so, replace
 946 
 947                            k                                      k
 948     each non-zero term t  x,  n <= k <= 2n-2, by the term t * [ x   mod f(x), p)]
 949                         k                                  k
 950 
 951     which we get from the array power_table.
 952 
 953 BUGS
 954 
 955      None.
 956 
 957 --------------------------------------------------------------------------------
 958 |                                Function Call                                 |
 959 ------------------------------------------------------------------------------*/
 960 
 961 void
 962     square( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
 963 {
 964 
 965 /*------------------------------------------------------------------------------
 966 |                               Local Variables                                |
 967 ------------------------------------------------------------------------------*/
 968 
 969 int
 970     i, j,                     /* Loop counters. */
 971     coeff,                    /* Coefficient of x ^ k term of t(x) ^2 */
 972     temp[ MAXDEGPOLY + 1 ] ;  /* Temporary storage for the new t(x). */
 973 
 974 /*------------------------------------------------------------------------------
 975 |                                Function Body                                 |
 976 ------------------------------------------------------------------------------*/
 977 
 978 /*
 979                                  0        n-1
 980     Compute the coefficients of x , ..., x.   These terms do not require
 981 
 982     reduction mod f(x) because their degree is less than n.
 983 */
 984 
 985 for (i = 0 ;  i <= n ;  ++i)
 986 
 987     temp[ i ] = coeff_of_square( t, i, n, p ) ;
 988 
 989 /*
 990                                  n        2n-2             k
 991     Compute the coefficients of x , ..., x.    Replace t  x  with
 992             k                                            k
 993     t  * [ x  (mod f(x), p) ] from array power_table when t is 
 994      k                                                     k
 995     non-zero.
 996 
 997 */
 998 
 999 for (i = n ;  i <= 2 * n - 2 ;  ++i)
1000 
1001     if ( (coeff = coeff_of_square( t, i, n, p )) != 0 )
1002 
1003         for (j = 0 ;  j <= n - 1 ;  ++j)
1004 
1005             temp[ j ] = mod( temp[ j ] +
1006                              mod( coeff * power_table[ i - n ] [ j ], p ),
1007                              p ) ;
1008 
1009 for (i = 0 ;  i <= n - 1 ;  ++i)
1010 
1011     t[ i ] = temp[ i ] ;
1012 
1013 } /* ======================== end of function square ======================== */
1014 
1015 
1016 
1017 /*==============================================================================
1018 |                                  product                                     |
1019 ================================================================================
1020 
1021 DESCRIPTION
1022 
1023      Compute s( x ) t( x ) (mod f(x), p)
1024 
1025      Uses a precomputed table of powers of x.
1026 
1027 INPUT
1028 
1029     s (int *)              Coefficients of s(x), of degree <= n-1.
1030 
1031     t (int *)              Coefficients of t(x), of degree <= n-1.
1032 
1033     power_table (int **)    x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.
1034 
1035     n (int, -infinity < n < infinity)  Degree of f(x).
1036 
1037     p (int, p > 0)         Mod p coefficient arithmetic.
1038 
1039 OUTPUT
1040 
1041     s (int *)              Overwritten with s( x ) t( x ) (mod f(x), p)
1042 
1043 EXAMPLE 
1044                                      3    2                 2
1045     Let n = 4 and p = 5, t( x ) = 4 x  + x + 4, s( x ) = 3 x + x + 2
1046 
1047                             5      4     3     2
1048     Then s( x ) t( x ) = 2 x  + 2 x + 4 x + 4 x + 4 x + 3, modulo 5, 
1049 
1050                                          4   2    
1051     and after reduction modulo f( x ) = x + x + 2 x + 3, using the power
1052 
1053                        4      2               5      3      2
1054     table entries for x  = 4 x + 3 x + 2 and x  = 4 x  + 3 x + 2 x, we get
1055     
1056                                        3      2
1057     s( x ) t( x ) (mod f( x ), p) = 2 x  + 3 x  + 4 x + 2     
1058         
1059 
1060 METHOD
1061 
1062     Compute the coefficients using the function coeff_of_product.
1063                                  
1064     The next step is to reduce s(x) t(x) modulo f(x) and p.  To do so, replace
1065 
1066                            k                                      k
1067     each non-zero term t  x,  n <= k <= 2n-2, by the term t * [ x   mod f(x), p)]
1068                         k                                  k
1069 
1070     which we get from the array power_table.
1071 
1072 BUGS
1073 
1074      None.
1075 
1076 --------------------------------------------------------------------------------
1077 |                                Function Call                                 |
1078 ------------------------------------------------------------------------------*/
1079 
1080 void
1081 product( int  * s, int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
1082 {
1083 
1084 /*------------------------------------------------------------------------------
1085 |                               Local Variables                                |
1086 ------------------------------------------------------------------------------*/
1087 
1088 int
1089     i, j,                     /* Loop counters. */
1090     coeff,                    /* Coefficient of x ^ k term of t(x) ^2 */
1091     temp[ MAXDEGPOLY + 1 ] ;  /* Temporary storage for the new t(x). */
1092 
1093 /*------------------------------------------------------------------------------
1094 |                                Function Body                                 |
1095 ------------------------------------------------------------------------------*/
1096 
1097 /*
1098                                  0        n-1
1099     Compute the coefficients of x , ..., x.   These terms do not require
1100 
1101     reduction mod f(x) because their degree is less than n.
1102 */
1103 
1104 for (i = 0 ;  i <= n ;  ++i)
1105 
1106     temp[ i ] = coeff_of_product( s, t, i, n, p ) ;
1107 
1108 /*
1109                                  n        2n-2             k
1110     Compute the coefficients of x , ..., x.    Replace t  x  with
1111             k                                            k
1112     t  * [ x  (mod f(x), p) ] from array power_table when t is 
1113      k                                                     k
1114     non-zero.
1115 
1116 */
1117 
1118 for (i = n ;  i <= 2 * n - 2 ;  ++i)
1119 
1120     if ( (coeff = coeff_of_product( s, t, i, n, p )) != 0 )
1121 
1122         for (j = 0 ;  j <= n - 1 ;  ++j)
1123 
1124             temp[ j ] = mod( temp[ j ] +
1125                              mod( coeff * power_table[ i - n ] [ j ], p ),
1126                              p ) ;
1127 
1128 for (i = 0 ;  i <= n - 1 ;  ++i)
1129 
1130     s[ i ] = temp[ i ] ;
1131 
1132 } /* ======================== end of function product ======================= */
1133 
1134 
1135 /*==============================================================================
1136 |                                  times_x                                     |
1137 ================================================================================
1138 
1139 DESCRIPTION
1140 
1141      Compute x t(x) (mod f(x), p)
1142 
1143      Uses a precomputed table of powers of x.
1144 
1145 INPUT
1146                                             2
1147     t (int *)              Coefficients of t (x), of degree <= n-1.
1148 
1149     power_table (int **)    x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.
1150 
1151     n (int, -infinity < n < infinity)  Degree of f(x).
1152 
1153     p (int, p > 0)         Mod p coefficient arithmetic.
1154 
1155 OUTPUT
1156 
1157     t (int *)              Overwritten with x t(x) (mod f(x), p)
1158 
1159 EXAMPLE 
1160                                     3       2        
1161     Let n = 4, p = 5, and t(x) = 2 x  +  4 x  + 3 x.  Let f(x) = 
1162      4    2                                  4      3      2
1163     x  + x  + 2 x + 3.  Then x t (x) = 2 x  + 4 x  + 3 x  and
1164                                      2                3     2
1165     x t(x) (mod f(x), p) = 2 * (4 x + 3 x + 2) + 4 x + 3 x  =
1166        3    2 
1167     4 x  + x + x + 4.
1168          3     2
1169     = 4 x + 2 x + 3 x + 2.
1170 
1171 METHOD
1172                           n-1
1173     Multiply t(x) = t    x    + ... + t  by shifting the coefficients
1174                      n-1               0
1175                              n
1176     to the left.  If an     x   term appears, eliminate it by 
1177 
1178     substitution using power_table.
1179 
1180 BUGS
1181 
1182      None.
1183 
1184 --------------------------------------------------------------------------------
1185 |                                Function Call                                 |
1186 ------------------------------------------------------------------------------*/
1187 
1188 void
1189     times_x( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
1190 {
1191 
1192 /*------------------------------------------------------------------------------
1193 |                               Local Variables                                |
1194 ------------------------------------------------------------------------------*/
1195 
1196 int
1197     i,              /* Loop counter. */
1198     coeff ;         /* Coefficient of x ^ n term of x t(x). */
1199 
1200 /*------------------------------------------------------------------------------
1201 |                                Function Body                                 |
1202 ------------------------------------------------------------------------------*/
1203 
1204 /*
1205     Multiply t(x) by x.  Do it by shifting the coefficients left in the array.
1206 */
1207 
1208 coeff = t[ n - 1 ] ;
1209 
1210 for (i = n-1 ;  i >= 1 ;  --i)
1211 
1212     t[ i ] = t[ i-1 ] ;
1213 
1214 t[ 0 ] = 0 ;
1215 
1216 /*
1217                   n                          n
1218     Zero out t  x .  Replace it with t * [ x  (mod f(x), p) ] from array
1219               n
1220     power_table.
1221 */
1222 
1223 if (coeff != 0)
1224 {
1225     for (i = 0 ;  i <= n - 1 ;  ++i)
1226 
1227         t[ i ] = mod( t[ i ] +
1228                       mod( coeff * power_table[ 0 ] [ i ], p ),
1229                       p ) ;
1230 }
1231 
1232 } /* ======================== end of function times_x ======================= */
1233 
1234 
1235 
1236 /*==============================================================================
1237 |                                  x_to_power                                  |
1238 ================================================================================
1239 
1240 DESCRIPTION
1241                      m
1242      Compute g(x) = x (mod f(x), p).
1243 
1244 INPUT
1245 
1246     m (int)       1 <= m <= r.  The exponent.
1247 
1248     power_table (int **)    x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.
1249 
1250     n (int, n >= 1)     Degree of monic polynomial f(x).
1251 
1252     p (int, p >= 2)     Modulo p coefficient arithmetic.
1253 
1254 OUTPUT
1255 
1256     g (int *)     Polynomial of degree <= n-1.
1257 
1258 EXAMPLE 
1259                               4   2 
1260     Let n = 4, p = 5, f(x) = x + x + 2 x + 3, and m = 156.
1261 
1262     156 = 0  . . . 0  1  0  0  1  1  1  0  0 (binary representation)
1263           |<- ignore ->| S  S SX SX SX  S  S (exponentiation rule,
1264                                               S = square, X = multiply by x)
1265      m
1266     x (mod f(x), p) =
1267 
1268          2     2
1269     S   x  =  x
1270 
1271          4       2
1272     S   x  =  4 x  + 3 x + 2
1273 
1274          8       3      2
1275     S   x  =  4 x  + 4 x + 1
1276 
1277          9       3     2
1278     X   x  =  4 x  +  x + 3 x + 3
1279 
1280          18       2  
1281     S   x  =  2 x  +  x + 2
1282 
1283          19      3     2
1284     X   x  =  2 x  +  x  + 2 x
1285 
1286          38      3       2
1287     S   x  =  2 x  +  4 x  +  3 x
1288 
1289          39      3       2
1290     X   x  =  4 x  +  2 x  +  x + 4
1291 
1292          78      3       2
1293     S   x  =  4 x  +  2 x  +  3 x + 2
1294 
1295          156 
1296     S   x    = 3
1297 
1298 METHOD
1299 
1300     Exponentiation by repeated squaring, using precomputed table of
1301     powers.  See ART OF COMPUTER PROGRAMMING, vol. 2, 2nd id,
1302     D. E. Knuth,  pgs 441-443.
1303 
1304      n         2n-2
1305     x,  ... , x    (mod f(x), p)
1306 
1307                     m
1308     to find g(x) = x   (mod f(x), p), expand m into binary,
1309 
1310            k        k-1
1311     m = a 2  + a   2    + . . . + a 2 + a
1312          k      k-1                1     0
1313 
1314                             m
1315     where a = 1, and split x   apart into
1316            k
1317 
1318             k      k
1319            2      2  a             2 a    a
1320      m                k-1             1    0
1321     x  =  x     x           . . .  x     x
1322 
1323 
1324     Then to raise x to the mth power, do
1325 
1326 
1327     t( x ) = x
1328 
1329     return if m = 1
1330 
1331 
1332     for i = k-1 downto 0 do
1333 
1334                    2
1335         t(x) = t(x)  (mod f(x), p)       // Square each time.
1336 
1337         if a = 1 then
1338             i
1339 
1340             t(x) = x t(x) (mod f(x), p)  // Times x only if current bit is 1
1341         endif
1342 
1343     endfor
1344                                                         k
1345                                                        2
1346     The initial t(x) = x gets squared k times to give x  .  If a  = 1 for
1347                                                                 i
1348     0 <= i <= k-1, we multiply by x which then gets squared i times more
1349 
1350               i
1351              2
1352     to give x .  On a binary computer, we use bit shifting and masking to
1353 
1354     identify the k bits  { a    . . .  a  } to the right of the leading 1
1355                             k-1         0
1356 
1357     bit.  There are log ( m ) - 1 squarings and (number of 1 bits) - 1
1358                            2
1359     multiplies.  
1360 
1361 BUGS
1362 
1363      None.
1364 
1365 --------------------------------------------------------------------------------
1366 |                                Function Call                                 |
1367 ------------------------------------------------------------------------------*/
1368 
1369 void
1370     x_to_power( bigint m, int * g, int power_table[][ MAXDEGPOLY ], int n, int p )
1371 {
1372 
1373 /*------------------------------------------------------------------------------
1374 |                               Local Variables                                |
1375 ------------------------------------------------------------------------------*/
1376 
1377 /*
1378     mask = 1000 ... 000.  That is, all bits of mask are zero except for the
1379     most significant bit of the computer word holding its value.  Signed or
1380     unsigned type for bigint shouldn't matter here, since we just mask and
1381     compare bits.
1382 */
1383 
1384 bigint
1385     mask = ((bigint)1 << (NUMBITS - 1)) ;
1386 
1387 int
1388     bit_count = 0,  /* Number of bits in m to the right of the leading bit. */
1389     i ;             /* Loop counter. */
1390 
1391 
1392 /*------------------------------------------------------------------------------
1393 |                                Function Body                                 |
1394 ------------------------------------------------------------------------------*/
1395 
1396 #ifdef DEBUG_PP_PRIMPOLY
1397 printf( "x to power = x ^ %lld\n", m ) ;
1398 #endif
1399 
1400 /*
1401     Initialize g(x) to x.  Exit right away if m = 1.
1402 */
1403 
1404 for (i = 0 ;  i <= n-1 ;  ++i)
1405 
1406     g[ i ] = 0 ;
1407 
1408 g[ 1 ] = 1 ;
1409 
1410 if (m == 1)
1411 
1412     return ;
1413 
1414 /*
1415     Advance the leading bit of m up to the word's left hand boundary.
1416     Count how many bits were to the right of the leading bit.
1417 */
1418 
1419 while (! (m & mask))
1420 {
1421     m <<= 1 ;
1422     ++bit_count ;
1423 }
1424 
1425 bit_count = NUMBITS - bit_count ;
1426 
1427 
1428 
1429 /*
1430     Exponentiation by repeated squaring.  Discard the leading 1 bit.
1431     Thereafter, square for every 0 bit;  square and multiply by x for
1432     every 1 bit.
1433 */
1434 
1435 while ( --bit_count > 0 )
1436 {
1437 
1438     m <<= 1 ;       /*  Expose the next bit. */
1439 
1440     square( g, power_table, n, p ) ;
1441 
1442     #ifdef DEBUG_PP_PRIMPOLY
1443     printf( "    after squaring, poly = \n" ) ;
1444     write_poly( g, n-1 ) ;
1445     printf( "\n\n" );
1446     #endif
1447 
1448     if (m & mask)   /*  Leading bit is 1. */
1449     {
1450         times_x( g, power_table, n, p ) ;
1451 
1452         #ifdef DEBUG_PP_PRIMPOLY
1453         printf( "    after additional times x, poly = \n" ) ;
1454         write_poly( g, n-1 ) ;
1455         printf( "\n\n" );
1456         #endif
1457     }
1458 }
1459 
1460 } /* ===================== end of function x_to_power ======================= */