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.4 - A Program for Computing Primitive Polynomials.
  29|     Copyright (C) 1999-2025 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 seanerikoconnor!AT!gmail!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
  61DESCRIPTION
  62
  63     Evaluate the polynomial f( x ) with modulo p arithmetic.
  64
  65INPUT
  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
  78RETURNS
  79
  80     f( x )  (int)
  81
  82EXAMPLE 
  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
  88METHOD
  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
  96BUGS
  97
  98     None.
  99
 100--------------------------------------------------------------------------------
 101|                                Function Call                                 |
 102------------------------------------------------------------------------------*/
 103
 104int
 105    eval_poly( int * f, int x, int n, int p )
 106{
 107
 108/*------------------------------------------------------------------------------
 109|                               Local Variables                                |
 110------------------------------------------------------------------------------*/
 111
 112int
 113    val = 1 ;   /*  The value of f(x) which is returned. */
 114
 115/*------------------------------------------------------------------------------
 116|                                Function Body                                 |
 117------------------------------------------------------------------------------*/
 118
 119while ( --n >= 0 )
 120
 121    val = mod( val * x + f[ n ], p ) ;
 122
 123return( val ) ;
 124
 125
 126} /* ========================= end of function eval_poly ==================== */
 127
 128
 129/*==============================================================================
 130|                                     linear_factor                            |
 131================================================================================
 132
 133DESCRIPTION
 134
 135     Check if the polynomial f(x) has any linear factors.
 136
 137INPUT
 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
 145RETURNS
 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
 154EXAMPLE 
 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
 159METHOD
 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
 164BUGS
 165
 166    None.
 167
 168--------------------------------------------------------------------------------
 169|                                Function Call                                 |
 170------------------------------------------------------------------------------*/
 171
 172int 
 173    linear_factor( int * f, int n, int p )
 174{
 175
 176/*------------------------------------------------------------------------------
 177|                               Local Variables                                |
 178------------------------------------------------------------------------------*/
 179
 180int i ;  /*  Loop counter. */
 181
 182/*------------------------------------------------------------------------------
 183|                                Function Body                                 |
 184------------------------------------------------------------------------------*/
 185
 186for (i = 1 ;  i <= p-1 ;  ++i)
 187
 188    if (eval_poly( f, i, n, p ) == 0)
 189
 190        return( YES ) ;
 191
 192return( NO ) ;
 193
 194} /* ====================== end of function linear_factor =================== */
 195
 196
 197/*==============================================================================
 198|                                 is_integer                                   |
 199================================================================================
 200
 201DESCRIPTION
 202
 203     Test if a polynomial is a constant.
 204
 205INPUT
 206
 207    t (int *)  mod p polynomial t(x). 
 208    n (int)    Its degree.
 209
 210RETURNS
 211
 212    YES    if t(x) is a constant polynomial.
 213    NO     otherwise
 214
 215EXAMPLE 
 216                               2
 217    Let n = 3.  Then t(x) = 2 x  is not constant but t(x) = 2 is.
 218
 219METHOD
 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
 224BUGS
 225
 226    None.
 227
 228--------------------------------------------------------------------------------
 229|                                Function Call                                 |
 230------------------------------------------------------------------------------*/
 231
 232int 
 233    is_integer( int * t, int n )
 234{
 235
 236/*------------------------------------------------------------------------------
 237|                                Function Body                                 |
 238------------------------------------------------------------------------------*/
 239
 240while ( n >= 1 )
 241
 242    if (t[ n-- ] != 0) 
 243
 244        return( NO ) ;
 245
 246return( YES ) ;
 247
 248} /* ====================== end of function is_integer ====================== */
 249
 250
 251/*==============================================================================
 252|                               construct_power_table                          |
 253================================================================================
 254
 255DESCRIPTION
 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
 262INPUT
 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
 268RETURNS
 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
 274EXAMPLE 
 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
 304METHOD
 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
 315BUGS
 316
 317     None.
 318
 319--------------------------------------------------------------------------------
 320|                                Function Call                                 |
 321------------------------------------------------------------------------------*/
 322
 323void 
 324    construct_power_table( int power_table[][ MAXDEGPOLY ], int * f, 
 325                           int n, int p )
 326{
 327
 328/*------------------------------------------------------------------------------
 329|                               Local Variables                                |
 330------------------------------------------------------------------------------*/
 331
 332int 
 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
 348for (i = 0 ;  i <= n-2 ;  ++i)
 349
 350    t[ i ] = 0 ;
 351
 352t[ 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
 361for (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
 398return ;
 399
 400} /* ================== end of function construct_power_table =============== */
 401
 402
 403
 404/*==============================================================================
 405|                                 autoconvolve                                  |
 406================================================================================
 407
 408DESCRIPTION
 409
 410     Compute a convolution-like sum for use in function coeff_of_square.
 411
 412INPUT
 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
 424RETURNS
 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
 433EXAMPLE 
 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
 445METHOD
 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
 451BUGS
 452
 453     None.
 454
 455--------------------------------------------------------------------------------
 456|                                Function Call                                 |
 457------------------------------------------------------------------------------*/
 458
 459int
 460    auto_convolve( int * t, int k, int lower, int upper, int p )
 461{
 462
 463/*------------------------------------------------------------------------------
 464|                               Local Variables                                |
 465------------------------------------------------------------------------------*/
 466
 467int 
 468    sum = 0,      /* Convolution sum. */
 469    i ;           /* Loop counter.    */
 470
 471/*------------------------------------------------------------------------------
 472|                                Function Body                                 |
 473------------------------------------------------------------------------------*/
 474
 475for (i = lower ;  i <= upper ;  ++i)
 476
 477   sum = mod( sum + mod( t[ i ] * t[ k - i ], p ), p ) ;
 478
 479return( sum ) ;
 480
 481} /* ==================== end of function auto_convolve ===================== */
 482
 483
 484/*==============================================================================
 485|                                 convolve                                     |
 486================================================================================
 487
 488DESCRIPTION
 489
 490     Compute a convolution-like sum.
 491
 492INPUT
 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
 508RETURNS
 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
 517EXAMPLE 
 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
 535METHOD
 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
 541BUGS
 542
 543     None.
 544
 545--------------------------------------------------------------------------------
 546|                                Function Call                                 |
 547------------------------------------------------------------------------------*/
 548
 549int
 550    convolve( int * s, int * t, int k, int lower, int upper, int p )
 551{
 552
 553/*------------------------------------------------------------------------------
 554|                               Local Variables                                |
 555------------------------------------------------------------------------------*/
 556
 557int 
 558    sum = 0,      /* Convolution sum. */
 559    i ;           /* Loop counter.    */
 560
 561/*------------------------------------------------------------------------------
 562|                                Function Body                                 |
 563------------------------------------------------------------------------------*/
 564
 565for (i = lower ;  i <= upper ;  ++i)
 566
 567   sum = mod( sum + mod( s[ i ] * t[ k - i ], p ), p ) ;
 568
 569return( sum ) ;
 570
 571} /* ======================= end of function convolve ======================= */
 572
 573
 574/*==============================================================================
 575|                               coeff_of_square                                |
 576================================================================================
 577
 578DESCRIPTION
 579                                    th                2
 580     Return the coefficient of the k   power of x in t( x )  modulo p.
 581 
 582INPUT
 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
 592RETURNS
 593                        th                2
 594    coefficient of the k   power of x in t(x) mod p.
 595
 596EXAMPLE 
 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
 616METHOD
 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
 678BUGS
 679
 680     None.
 681
 682--------------------------------------------------------------------------------
 683|                                Function Call                                 |
 684------------------------------------------------------------------------------*/
 685
 686int 
 687    coeff_of_square( int * t, int k, int n, int p )
 688{
 689
 690/*------------------------------------------------------------------------------
 691|                               Local Variables                                |
 692------------------------------------------------------------------------------*/
 693
 694int 
 695    sum = 0 ;      /* kth coefficient of t(x) ^ 2. */
 696
 697/*------------------------------------------------------------------------------
 698|                                Function Body                                 |
 699------------------------------------------------------------------------------*/
 700
 701if (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}
 713else 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
 726return( sum ) ;
 727
 728} /* ================== end of function coeff_of_square ===================== */
 729
 730
 731/*==============================================================================
 732|                               coeff_of_product                               |
 733================================================================================
 734
 735DESCRIPTION
 736                                    th
 737     Return the coefficient of the k   power of x in s( x ) t( x )  modulo p.
 738 
 739INPUT
 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
 749RETURNS
 750                        th
 751    coefficient of the k   power of x in s(x) t(x) mod p.
 752
 753EXAMPLE  
 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
 766METHOD
 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
 858BUGS
 859
 860     None.
 861
 862--------------------------------------------------------------------------------
 863|                                Function Call                                 |
 864------------------------------------------------------------------------------*/
 865
 866int 
 867    coeff_of_product( int * s, int * t, int k, int n, int p )
 868{
 869
 870/*------------------------------------------------------------------------------
 871|                               Local Variables                                |
 872------------------------------------------------------------------------------*/
 873
 874int 
 875    sum = 0 ;      /* kth coefficient of t(x) ^ 2. */
 876
 877/*------------------------------------------------------------------------------
 878|                                Function Body                                 |
 879------------------------------------------------------------------------------*/
 880
 881if (0 <= k && k <= n-1)
 882{
 883    sum = convolve( s, t, k, 0, k, p ) ;
 884}
 885else if (n <= k && k <= 2 * n - 2)
 886{
 887    sum = convolve( s, t, k, k - n + 1, n - 1, p ) ;
 888}
 889
 890return( sum ) ;
 891
 892} /* ================== end of function coeff_of_product ==================== */
 893
 894
 895/*==============================================================================
 896|                                  square                                      |
 897================================================================================
 898
 899DESCRIPTION
 900              2
 901     Compute t ( x ) (mod f(x), p)
 902
 903     Uses a precomputed table of powers of x.
 904
 905INPUT
 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
 915OUTPUT
 916                                             2
 917    t (int *)              Overwritten with t (x) (mod f(x), p)
 918
 919EXAMPLE 
 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
 936METHOD
 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
 953BUGS
 954
 955     None.
 956
 957--------------------------------------------------------------------------------
 958|                                Function Call                                 |
 959------------------------------------------------------------------------------*/
 960
 961void 
 962    square( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
 963{
 964
 965/*------------------------------------------------------------------------------
 966|                               Local Variables                                |
 967------------------------------------------------------------------------------*/
 968
 969int 
 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
 985for (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
 999for (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
1009for (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
1021DESCRIPTION
1022
1023     Compute s( x ) t( x ) (mod f(x), p)
1024
1025     Uses a precomputed table of powers of x.
1026
1027INPUT
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
1039OUTPUT
1040
1041    s (int *)              Overwritten with s( x ) t( x ) (mod f(x), p)
1042
1043EXAMPLE 
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
1060METHOD
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
1072BUGS
1073
1074     None.
1075
1076--------------------------------------------------------------------------------
1077|                                Function Call                                 |
1078------------------------------------------------------------------------------*/
1079
1080void 
1081product( int  * s, int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
1082{
1083
1084/*------------------------------------------------------------------------------
1085|                               Local Variables                                |
1086------------------------------------------------------------------------------*/
1087
1088int 
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
1104for (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
1118for (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
1128for (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
1139DESCRIPTION
1140
1141     Compute x t(x) (mod f(x), p)
1142
1143     Uses a precomputed table of powers of x.
1144
1145INPUT
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
1155OUTPUT
1156
1157    t (int *)              Overwritten with x t(x) (mod f(x), p)
1158
1159EXAMPLE 
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
1171METHOD
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
1180BUGS
1181
1182     None.
1183
1184--------------------------------------------------------------------------------
1185|                                Function Call                                 |
1186------------------------------------------------------------------------------*/
1187
1188void 
1189    times_x( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
1190{
1191
1192/*------------------------------------------------------------------------------
1193|                               Local Variables                                |
1194------------------------------------------------------------------------------*/
1195
1196int 
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
1208coeff = t[ n - 1 ] ;
1209
1210for (i = n-1 ;  i >= 1 ;  --i)
1211
1212    t[ i ] = t[ i-1 ] ;
1213
1214t[ 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
1223if (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
1240DESCRIPTION
1241                     m
1242     Compute g(x) = x (mod f(x), p).
1243
1244INPUT
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
1254OUTPUT
1255
1256    g (int *)     Polynomial of degree <= n-1.
1257
1258EXAMPLE 
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
1298METHOD
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
1361BUGS
1362
1363     None.
1364
1365--------------------------------------------------------------------------------
1366|                                Function Call                                 |
1367------------------------------------------------------------------------------*/
1368
1369void 
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
1384bigint 
1385    mask = ((bigint)1 << (NUMBITS - 1)) ;
1386
1387int 
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
1397printf( "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
1404for (i = 0 ;  i <= n-1 ;  ++i)
1405
1406    g[ i ] = 0 ;
1407
1408g[ 1 ] = 1 ;
1409
1410if (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
1419while (! (m & mask))
1420{
1421    m <<= 1 ;
1422    ++bit_count ;
1423}
1424
1425bit_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
1435while ( --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 ======================= */