/*==============================================================================
| 
|  File Name:     
|
|     ppPolyArith.c
|
|  Description:   
|
|     Polynomial arithmetic and exponentiation routines.
| 
|  Functions:
|
|     eval_poly
|     linear_factor
|     is_integer
|     construct_power_table
|     auto_convolve
|     convolve
|     coeff_of_square
|     coeff_of_product
|     square
|     product
|     times_x
|     x_to_power
|
|  LEGAL
|
|     Primpoly Version 13.0 - A Program for Computing Primitive Polynomials.
|     Copyright (C) 1999-2017 by Sean Erik O'Connor.  All Rights Reserved.
|
|     This program is free software: you can redistribute it and/or modify
|     it under the terms of the GNU General Public License as published by
|     the Free Software Foundation, either version 3 of the License, or
|     (at your option) any later version.
|
|     This program is distributed in the hope that it will be useful,
|     but WITHOUT ANY WARRANTY; without even the implied warranty of
|     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
|     GNU General Public License for more details.
|
|     You should have received a copy of the GNU General Public License
|     along with this program.  If not, see <http://www.gnu.org/licenses/>.
|     
|     The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com
|     with !DOT! replaced by . and the !AT! replaced by @
|
==============================================================================*/

/*------------------------------------------------------------------------------
|                                Include Files                                 |
------------------------------------------------------------------------------*/

#include <stdio.h>
#include "Primpoly.h"


/*==============================================================================
|                                   eval_poly                                  |
================================================================================

DESCRIPTION

     Evaluate the polynomial f( x ) with modulo p arithmetic.

INPUT

     f (int *)      nth degree mod p polynomial
         
                              n         n-1     
                    f( x ) = x  +  a   x  + ... + a    0 <= a  < p
                                    n-1            0         i          
     x (int)        0 <= x < p

     n (int)        n >= 1

     p (int)         

RETURNS

     f( x )  (int)

EXAMPLE 
                                  4
     Let n = 4, p = 5 and f(x) = x  + 3x + 3.  Then f(x) = 
     (((x + 0)x + 0)x + 3)x + 3.  f(2) = (((2 + 0)2 + 0)2 + 3) =
     (8 + 3)2 + 3 = 1 + 2 + 3 (mod 5) = 0.

METHOD

     By Horner's rule, f(x) = ( ... ( (x + a   )x + ... )x + a .
                                            n-1               0
     We evaluate recursively, f := f * x + a (mod p), starting
                                            i 
     with f = 1 and i = n-1.

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int
    eval_poly( int * f, int x, int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    val = 1 ;   /*  The value of f(x) which is returned. */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

while ( --n >= 0 )

    val = mod( val * x + f[ n ], p ) ;

return( val ) ;


} /* ========================= end of function eval_poly ==================== */


/*==============================================================================
|                                     linear_factor                            |
================================================================================

DESCRIPTION

     Check if the polynomial f(x) has any linear factors.

INPUT

    f (int *)  nth degree monic mod p polynomial f(x). 

    n (int)    Its degree.

    p (int)    Modulus for coefficient arithmetic.

RETURNS

    YES    if f( a ) = 0 (mod p) for a = 1, 2, ... p-1.
    NO     otherwise

    i.e. check if f(x) contains a linear factor (x - a).  We don't need to test
    for the root a = 0 because f(x) was chosen in main to have a non-zero 
    constant term, hence no zero root.

EXAMPLE 
                                 4      
    Let n = 4, p = 5 and f(x) = x  + 3x + 3.  Then f(1) = 2 (mod 5), but 
    f(2) = 0 (mod 5), so we exit immediately with a yes.

METHOD

    Evaluate f(x) at x = 1, ..., p-1 by Horner's rule.  Return instantly the 
    moment f(x) evaluates to 0.

BUGS

    None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int
    linear_factor( int * f, int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int i ;  /*  Loop counter. */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

for (i = 1 ;  i <= p-1 ;  ++i)

    if (eval_poly( f, i, n, p ) == 0)

        return( YES ) ;

return( NO ) ;

} /* ====================== end of function linear_factor =================== */


/*==============================================================================
|                                 is_integer                                   |
================================================================================

DESCRIPTION

     Test if a polynomial is a constant.

INPUT

    t (int *)  mod p polynomial t(x). 
    n (int)    Its degree.

RETURNS

    YES    if t(x) is a constant polynomial.
    NO     otherwise

EXAMPLE 
                               2
    Let n = 3.  Then t(x) = 2 x  is not constant but t(x) = 2 is.

METHOD

    A constant polynomial is zero in its first through n th degree
    terms.  Return immediately with NO if any such term is non-zero.

BUGS

    None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int
    is_integer( int * t, int n )
{

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

while ( n >= 1 )

    if (t[ n-- ] != 0)

        return( NO ) ;

return( YES ) ;

} /* ====================== end of function is_integer ====================== */


/*==============================================================================
|                               construct_power_table                          |
================================================================================

DESCRIPTION

     Construct a table of powers of x:

      n                     2n-2
     x  (mod f(x), p)  ... x    (mod f(x), p)

INPUT

    f (int *)   Coefficients of f(x), a monic polynomial of degree n.
    n (int, -infinity < n < infinity)
    p (int, p > 0)

RETURNS

    power_table (int *)  power_table[i][j] is the coefficient of 
     j       n+i
    x   in  x   (mod f(x), p) where 0 <= i <= n-2 and 0 <= j <= n-1. 

EXAMPLE 
                                  4     2                     4
     Let n = 4, p = 5 and f(x) = x  +  x  +  2x  +  3.  Then x  = 

         2                  2
     -( x  +  2x  + 3) = 4 x  + 3 x + 2 (mod f(x), 5), and we get

      4                    2
     x  (mod f(x), 5) = 4 x  + 3 x + 2 = power_table[0].

      5                       2                 3      2
     x  (mod f(x), 5) = x( 4 x  + 3 x + 2) = 4 x  + 3 x  + 2x 
                      = power_table[1].

      6                       3      2           4      3      2
     x  (mod f(x), 5) = x( 4 x  + 3 x + 2 x) = 4x  + 3 x  + 2 x 

                              2                 3      2
                      = 4 ( 4x  + 3 x + 2) + 3 x  + 2 x  =

                           3     2
                      = 3 x + 3 x + 2 x + 3 = power_table[2].

                                    j
     power_table[i][j]:       | 0  1  2  3
                           ---+-------------
                            0 | 2  3  4  0 
                        i   1 | 0  2  3  4 
                            2 | 3  2  3  3 

METHOD
                              n-1
     Beginning with t( x ) = x,    compute the next power of x from the last
                                                         n
     one by multiplying by x.  If necessary, substitute x  =
             n-1
     -( a   x    + ... + a ) to reduce the degree.  This formula comes from
         n-1              0
                               n         n-1
     the observation f( x ) = x   + a   x    + ... + a    = 0 (mod f(x), p).
                                     n-1              0
BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

void
    construct_power_table( int power_table[][ MAXDEGPOLY ], int * f,
                           int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    i, j,                  /*  Loop counters.  */
    coeff,                 /*  Coefficient of x ^ n in t(x) */
    t[ MAXDEGPOLY + 1 ] ;  /*  t(x) is temporary storage for x ^ k (mod f(x),p)
                               n <= k <= 2n-2.  Its degree can go as high as
                               n before it is reduced again. */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

/*
                         n-1
    Initialize t( x ) = x 
*/

for (i = 0 ;  i <= n-2 ;  ++i)

    t[ i ] = 0 ;

t[ n-1 ] = 1 ;



/*
                                        i+n
    Fill the ith row of the table with x   (mod f(x), p).
*/

for (i = 0 ;  i <= n-2 ;  ++i)
{

    /*  Compute t(x) = x t(x)  */
    for (j = n ;  j >= 1 ;  --j)

        t[ j ] = t[ j-1 ] ;

    t[ 0 ] = 0 ;


    /*  Coefficient of the nth degree term of t(x).  */
    if ( (coeff = t[ n ]) != 0)
    {
        t[ n ] = 0 ;               /*  Zero out the x ^ n th term. */

       /*
                 n       n                        n-1
        Replace x  with x  (mod f(x), p) = -(a   x   + ... + a )
                                              n-1             0
       */

        for (j = 0 ;  j <= n-1 ;  ++j)

            t[ j ] = mod( t[ j ] +
                          mod( -coeff * f[ j ], p ),
                          p ) ;
    }  /* end if */


    /*  Copy t(x) into row i of power_table. */
    for (j = 0 ;  j <= n-1 ; ++j)

        power_table[i][j] = t[ j ] ;

} /* end for */

return ;

} /* ================== end of function construct_power_table =============== */



/*==============================================================================
|                                 autoconvolve                                  |
================================================================================

DESCRIPTION

     Compute a convolution-like sum for use in function coeff_of_square.

INPUT
                                             n-1
    t (int *)   Coefficients of t(x) = t    x    + ... + t x  + t
                                        n-1               1      0
    k (int), 0 <= k <= 2n - 2)

    lower (int, 0 <= lower <= n-1) 

    uppper (int, 0 <= upper <= n-1) 

    p (int, p > 0)

RETURNS

     upper
     ---
     \    t  t       But define the sum to be 0 when lower > upper to catch
     /     i  k-i    the special cases that come up in function coeff_of_square.
     ---
     i=lower

EXAMPLE 
                       3      2      
     Suppose t(x) = 4 x  +  x  +  3 x  +  3, lower = 1, upper = 3, n = 3,

     and p = 5.  For k = 3, auto_convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +

     t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3.  For lower = 0,

     upper = -1, or for lower = 3 and upper = 2, auto_convolve = 0, no matter what

     k is.

METHOD

    A "for" loop in the C language is not executed when its lower limit 

    exceeds its upper limit, leaving sum = 0.

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int
    auto_convolve( int * t, int k, int lower, int upper, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    sum = 0,      /* Convolution sum. */
    i ;           /* Loop counter.    */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

for (i = lower ;  i <= upper ;  ++i)

   sum = mod( sum + mod( t[ i ] * t[ k - i ], p ), p ) ;

return( sum ) ;

} /* ==================== end of function auto_convolve ===================== */


/*==============================================================================
|                                 convolve                                     |
================================================================================

DESCRIPTION

     Compute a convolution-like sum.

INPUT
                                             n-1
    s (int *)   Coefficients of s(x) = s    x    + ... + s x  + s
                                        n-1               1      0  
                                             n-1
    t (int *)   Coefficients of t(x) = t    x    + ... + t x  + t
                                        n-1               1      0

    k (int), 0 <= k <= 2n - 2)

    lower (int, 0 <= lower <= n-1) 

    uppper (int, 0 <= upper <= n-1) 

    p (int, p > 0)

RETURNS

     upper
     ---
     \    s  t       But define the sum to be 0 when lower > upper to catch
     /     i  k-i    the special cases
     ---
     i=lower

EXAMPLE 
                       3      2      
     Suppose s(x) = 4 x  +  x  +  3 x  +  3, 

                        3      2      
     Suppose t(x) = 4 x  +  x  +  3 x  +  3, 

       
     lower = 1, upper = 3, n = 3,

     and p = 5.  For k = 3, convolve = t[ 1 ] t[ 2 ] + t[ 2 ] t[ 1 ] +

     t[ 3 ] t[ 0 ] = 3 * 1 + 1 * 3 + 4 * 3 = 18 mod 5 = 3.  For lower = 0,

     upper = -1, or for lower = 3 and upper = 2, convolve = 0, no matter what

     k is.

METHOD

    A "for" loop in the C language is not executed when its lower limit 

    exceeds its upper limit, leaving sum = 0.

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int
    convolve( int * s, int * t, int k, int lower, int upper, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    sum = 0,      /* Convolution sum. */
    i ;           /* Loop counter.    */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

for (i = lower ;  i <= upper ;  ++i)

   sum = mod( sum + mod( s[ i ] * t[ k - i ], p ), p ) ;

return( sum ) ;

} /* ======================= end of function convolve ======================= */


/*==============================================================================
|                               coeff_of_square                                |
================================================================================

DESCRIPTION
                                    th                2
     Return the coefficient of the k   power of x in t( x )  modulo p.
 
INPUT

    t (int *)  Coefficients of t(x), of degree <= n-1.

    k (int, 0 <= k <= 2n-2)

    n (int, -infinity < n < infinity)

    p (int, p > 0)

RETURNS
                        th                2
    coefficient of the k   power of x in t(x) mod p.

EXAMPLE 
                                    3     2           2
    Let n = 4, p = 5, and t(x) = 4 x  +  x  + 4.  t(x) = 

             0                   1                         2
    (4 * 4) x  +  (2 * 1 * 4) x  +  (2 * 1 * 4 + 1 * 1) x  +

                             3                        4
    (2 * 1 * 1 + 2 * 4 * 4) x  + (2 * 4 * 1 + 1 * 1) x  +

                 5            6     
    (2 * 4 * 1) x  + (4 * 4) x  =  

     6      5     4     3     2
    x  + 3 x + 4 x + 4 x + 4 x  + 3 x + 1, all modulo 5.

            k        |  0  1  2  3  4  5  6
     ----------------+---------------------
     coeff_of_square |  1  3  4  4  4  3  1

METHOD
                                                         2
    The formulas were gotten by writing out the product t (x) explicitly.

    The sum is 0 in two cases:  

        (1) when k = 0 and the limits of summation are 0 to -1 

        (2) k = 2n - 2, when the limits of summation are n to n-1.

    To derive the formulas, let 

                      n-1           
    Let t (x) = t    x     +  ... + t x + t
                 n-1                 1     0

    Look at the formulas in coeff_of_product for each power of x,
    replacing s with t, and observe that half of the terms are
    duplicates, so we can save computation time.

    Inspection yields the formulas,


    for 0 <= k <= n-1, even k,

     k/2-1
      ---             2
   2  \   t  t     + t
      /    i  k-i     k/2
      ---
      i=0

    for 0 <= k <= n-1, odd k,

      (k-1)/2
      ---  
    2 \   t  t    
      /    i  k-i
      ---
      i=0


    and for n <= k <= 2n-2, even k,

      n-1
      ---            2
   2  \   t  t    + t
      /    i  k-i    k/2
      ---
      i=k/2+1

      and for n <= k <= 2n-2, odd k,

      n-1
      ---  
   2  \   t  t    
      /    i  k-i   
      ---
      i=(k+1)/2



BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int
    coeff_of_square( int * t, int k, int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    sum = 0 ;      /* kth coefficient of t(x) ^ 2. */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

if (0 <= k && k <= n-1)
{

    if (k % 2 == 0)        /* Even k */

        sum = mod( mod( 2 * auto_convolve( t, k, 0, k/2 - 1, p ), p) +
                   t[ k/2 ] * t[ k/2 ], p ) ;

     else                  /* Odd k */

         sum = mod( 2 * auto_convolve( t, k, 0, (k-1)/2, p ), p) ;
}
else if (n <= k && k <= 2 * n - 2)
{

    if (k % 2 == 0)        /* Even k */

        sum = mod( mod( 2 * auto_convolve( t, k, k/2 + 1, n-1, p ), p) +
                   t[ k/2 ] * t[ k/2 ], p ) ;

     else                  /* Odd k */

         sum = mod( 2 * auto_convolve( t, k, (k+1)/2, n-1, p ), p) ;
}

return( sum ) ;

} /* ================== end of function coeff_of_square ===================== */


/*==============================================================================
|                               coeff_of_product                               |
================================================================================

DESCRIPTION
                                    th
     Return the coefficient of the k   power of x in s( x ) t( x )  modulo p.
 
INPUT

    t (int *)  Coefficients of t(x), of degree <= n-1.

    k (int, 0 <= k <= 2n-2)

    n (int, -infinity < n < infinity)

    p (int, p > 0)

RETURNS
                        th
    coefficient of the k   power of x in s(x) t(x) mod p.

EXAMPLE  
                              3     2                  2
  Let n = 4, p = 5, t(x) = 4 x  +  x  + 4, s( x ) = 3 x  + x + 2

  We'll do the case k=3,

  t3 s0 + t2 s1 + t1 s2 + t0 s3 = 4 * 2 + 1 * 1 + 0 * 3 + 4 * 0 = 9 = 4 (mod 5).


            k         |  0  1  2  3  4  5  6
     -----------------+---------------------
     coeff_of_product |  3  4  4  4  2  2  0

METHOD
                                                         
    The formulas were gotten by writing out the product s(x) t (x) explicitly.

    The sum is 0 in two cases:  

        (1) when k = 0 and the limits of summation are 0 to -1 

        (2) k = 2n - 2, when the limits of summation are n to n-1.


    To derive the formulas, let 

                      n-1           
    Let s (x) = s    x     +  ... + s x + s  and
                 n-1                 1     0

                      n-1           
        t (x) = t    x     +  ... + t x + t 
                 n-1                 1     0

    and multiply out the terms, collecting like powers of x:


    Power of x     Coefficient
    ==========================
     2n-2
    x              s    t
                    n-1  n-1

     2n-3
    x              s    t    +  s    t
                    n-2  n-1     n-1  n-2     

     2n-4
    x              s    t    +  s    t    +  s    t
                    n-3  n-1     n-2  n-2     n-3  n-1

     2n-5
    x              s    t    +  s    t    +  s    t    +  s    t
                    n-4  n-1     n-3  n-2     n-2  n-3     n-1  n-4

     . . .

     n
    x              s  t    +  s  t    + ...  +  s    t
                    1  n-1     2  n-2            n-1  1

     n-1
    x              s  t    +  s  t    + ...  +  s    t
                    0  n-1     1  n-2            n-1  0

    . . .

     3
    x              s  t  +  s  t  +  s  t  +  s  t
                    0  3     1  2     2  1     3  0

     2
    x              s  t  +  s  t  +  s  t 
                    0  2     1  1     2  0

     
    x              s  t  +  s  t 
                    0  1     1  0 

    1              s  t
                    0  0


    Inspection yields the formulas,


    for 0 <= k <= n-1,

      k
     ---
     \   s  t
     /    i  k-i
     ---
     i=0


    and for n <= k <= 2n-2,

     n-1
     ---
     \   s  t
     /    i  k-i
     ---
     i=k-n+1

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

int
    coeff_of_product( int * s, int * t, int k, int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    sum = 0 ;      /* kth coefficient of t(x) ^ 2. */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

if (0 <= k && k <= n-1)
{
    sum = convolve( s, t, k, 0, k, p ) ;
}
else if (n <= k && k <= 2 * n - 2)
{
    sum = convolve( s, t, k, k - n + 1, n - 1, p ) ;
}

return( sum ) ;

} /* ================== end of function coeff_of_product ==================== */


/*==============================================================================
|                                  square                                      |
================================================================================

DESCRIPTION
              2
     Compute t ( x ) (mod f(x), p)

     Uses a precomputed table of powers of x.

INPUT

    t (int *)              Coefficients of t (x), of degree <= n-1.

    power_table (int **)   x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.

    n (int, -infinity < n < infinity)  Degree of f(x).

    p (int, p > 0)         Mod p coefficient arithmetic.

OUTPUT
                                             2
    t (int *)              Overwritten with t (x) (mod f(x), p)

EXAMPLE 
                                    3     2        
    Let n = 4, p = 5, and t(x) = 4 x  +  x  + 4.  Let f(x) = 

     4    2                   2       6      5      4      3     2
    x  + x  + 2 x + 3.  Then t (x) = x  + 3 x  + 4 x  + 4 x + 4 x + 3 x

           2                            3     2
    + 1.  t (x) (mod f(x), p) = 1 * (3 x + 3 x + 2 x + 3) +

            3      2                 2                3     2
    3 * (4 x  + 3 x + 2 x) + 4 * (4 x + 3 x + 2) + 4 x + 4 x + 3 x + 1

         3     2
    = 4 x + 2 x + 3 x + 2.


METHOD
         2            2n-2              n         n-1
    Let t (x) = t    x     +  ... + t  x  +  t   x   +  ... + t .
                 2n-2                n        n-1              0

    Compute the coefficients t  using the function coeff_of_square.
                              k

                                2
    The next step is to reduce t (x) modulo f(x).  To do so, replace

                           k                                      k
    each non-zero term t  x,  n <= k <= 2n-2, by the term t * [ x   mod f(x), p)]
                        k                                  k

    which we get from the array power_table.

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

void
    square( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    i, j,                     /* Loop counters. */
    coeff,                    /* Coefficient of x ^ k term of t(x) ^2 */
    temp[ MAXDEGPOLY + 1 ] ;  /* Temporary storage for the new t(x). */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

/*
                                 0        n-1
    Compute the coefficients of x , ..., x.   These terms do not require

    reduction mod f(x) because their degree is less than n.
*/

for (i = 0 ;  i <= n ;  ++i)

    temp[ i ] = coeff_of_square( t, i, n, p ) ;

/*
                                 n        2n-2             k
    Compute the coefficients of x , ..., x.    Replace t  x  with
            k                                            k
    t  * [ x  (mod f(x), p) ] from array power_table when t is 
     k                                                     k
    non-zero.

*/

for (i = n ;  i <= 2 * n - 2 ;  ++i)

    if ( (coeff = coeff_of_square( t, i, n, p )) != 0 )

        for (j = 0 ;  j <= n - 1 ;  ++j)

            temp[ j ] = mod( temp[ j ] +
                             mod( coeff * power_table[ i - n ] [ j ], p ),
                             p ) ;

for (i = 0 ;  i <= n - 1 ;  ++i)

    t[ i ] = temp[ i ] ;

} /* ======================== end of function square ======================== */



/*==============================================================================
|                                  product                                     |
================================================================================

DESCRIPTION

     Compute s( x ) t( x ) (mod f(x), p)

     Uses a precomputed table of powers of x.

INPUT

    s (int *)              Coefficients of s(x), of degree <= n-1.

    t (int *)              Coefficients of t(x), of degree <= n-1.

    power_table (int **)    x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.

    n (int, -infinity < n < infinity)  Degree of f(x).

    p (int, p > 0)         Mod p coefficient arithmetic.

OUTPUT

    s (int *)              Overwritten with s( x ) t( x ) (mod f(x), p)

EXAMPLE 
                                     3    2                 2
    Let n = 4 and p = 5, t( x ) = 4 x  + x + 4, s( x ) = 3 x + x + 2

                            5      4     3     2
    Then s( x ) t( x ) = 2 x  + 2 x + 4 x + 4 x + 4 x + 3, modulo 5, 

                                         4   2    
    and after reduction modulo f( x ) = x + x + 2 x + 3, using the power

                       4      2               5      3      2
    table entries for x  = 4 x + 3 x + 2 and x  = 4 x  + 3 x + 2 x, we get
    
                                       3      2
    s( x ) t( x ) (mod f( x ), p) = 2 x  + 3 x  + 4 x + 2     
        

METHOD

    Compute the coefficients using the function coeff_of_product.
                                 
    The next step is to reduce s(x) t(x) modulo f(x) and p.  To do so, replace

                           k                                      k
    each non-zero term t  x,  n <= k <= 2n-2, by the term t * [ x   mod f(x), p)]
                        k                                  k

    which we get from the array power_table.

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

void
product( int  * s, int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    i, j,                     /* Loop counters. */
    coeff,                    /* Coefficient of x ^ k term of t(x) ^2 */
    temp[ MAXDEGPOLY + 1 ] ;  /* Temporary storage for the new t(x). */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

/*
                                 0        n-1
    Compute the coefficients of x , ..., x.   These terms do not require

    reduction mod f(x) because their degree is less than n.
*/

for (i = 0 ;  i <= n ;  ++i)

    temp[ i ] = coeff_of_product( s, t, i, n, p ) ;

/*
                                 n        2n-2             k
    Compute the coefficients of x , ..., x.    Replace t  x  with
            k                                            k
    t  * [ x  (mod f(x), p) ] from array power_table when t is 
     k                                                     k
    non-zero.

*/

for (i = n ;  i <= 2 * n - 2 ;  ++i)

    if ( (coeff = coeff_of_product( s, t, i, n, p )) != 0 )

        for (j = 0 ;  j <= n - 1 ;  ++j)

            temp[ j ] = mod( temp[ j ] +
                             mod( coeff * power_table[ i - n ] [ j ], p ),
                             p ) ;

for (i = 0 ;  i <= n - 1 ;  ++i)

    s[ i ] = temp[ i ] ;

} /* ======================== end of function product ======================= */


/*==============================================================================
|                                  times_x                                     |
================================================================================

DESCRIPTION

     Compute x t(x) (mod f(x), p)

     Uses a precomputed table of powers of x.

INPUT
                                            2
    t (int *)              Coefficients of t (x), of degree <= n-1.

    power_table (int **)    x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.

    n (int, -infinity < n < infinity)  Degree of f(x).

    p (int, p > 0)         Mod p coefficient arithmetic.

OUTPUT

    t (int *)              Overwritten with x t(x) (mod f(x), p)

EXAMPLE 
                                    3       2        
    Let n = 4, p = 5, and t(x) = 2 x  +  4 x  + 3 x.  Let f(x) = 
     4    2                                  4      3      2
    x  + x  + 2 x + 3.  Then x t (x) = 2 x  + 4 x  + 3 x  and
                                     2                3     2
    x t(x) (mod f(x), p) = 2 * (4 x + 3 x + 2) + 4 x + 3 x  =
       3    2 
    4 x  + x + x + 4.
         3     2
    = 4 x + 2 x + 3 x + 2.

METHOD
                          n-1
    Multiply t(x) = t    x    + ... + t  by shifting the coefficients
                     n-1               0
                             n
    to the left.  If an     x   term appears, eliminate it by 

    substitution using power_table.

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

void
    times_x( int * t, int power_table[][ MAXDEGPOLY ], int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

int
    i,              /* Loop counter. */
    coeff ;         /* Coefficient of x ^ n term of x t(x). */

/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

/*
    Multiply t(x) by x.  Do it by shifting the coefficients left in the array.
*/

coeff = t[ n - 1 ] ;

for (i = n-1 ;  i >= 1 ;  --i)

    t[ i ] = t[ i-1 ] ;

t[ 0 ] = 0 ;

/*
                  n                          n
    Zero out t  x .  Replace it with t * [ x  (mod f(x), p) ] from array
              n
    power_table.
*/

if (coeff != 0)
{
    for (i = 0 ;  i <= n - 1 ;  ++i)

        t[ i ] = mod( t[ i ] +
                      mod( coeff * power_table[ 0 ] [ i ], p ),
                      p ) ;
}

} /* ======================== end of function times_x ======================= */



/*==============================================================================
|                                  x_to_power                                  |
================================================================================

DESCRIPTION
                     m
     Compute g(x) = x (mod f(x), p).

INPUT

    m (int)       1 <= m <= r.  The exponent.

    power_table (int **)    x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic.

    n (int, n >= 1)     Degree of monic polynomial f(x).

    p (int, p >= 2)     Modulo p coefficient arithmetic.

OUTPUT

    g (int *)     Polynomial of degree <= n-1.

EXAMPLE 
                              4   2 
    Let n = 4, p = 5, f(x) = x + x + 2 x + 3, and m = 156.

    156 = 0  . . . 0  1  0  0  1  1  1  0  0 (binary representation)
          |<- ignore ->| S  S SX SX SX  S  S (exponentiation rule,
                                              S = square, X = multiply by x)
     m
    x (mod f(x), p) =

         2     2
    S   x  =  x

         4       2
    S   x  =  4 x  + 3 x + 2

         8       3      2
    S   x  =  4 x  + 4 x + 1

         9       3     2
    X   x  =  4 x  +  x + 3 x + 3

         18       2  
    S   x  =  2 x  +  x + 2

         19      3     2
    X   x  =  2 x  +  x  + 2 x

         38      3       2
    S   x  =  2 x  +  4 x  +  3 x

         39      3       2
    X   x  =  4 x  +  2 x  +  x + 4

         78      3       2
    S   x  =  4 x  +  2 x  +  3 x + 2

         156 
    S   x    = 3

METHOD

    Exponentiation by repeated squaring, using precomputed table of
    powers.  See ART OF COMPUTER PROGRAMMING, vol. 2, 2nd id,
    D. E. Knuth,  pgs 441-443.

     n         2n-2
    x,  ... , x    (mod f(x), p)

                    m
    to find g(x) = x   (mod f(x), p), expand m into binary,

           k        k-1
    m = a 2  + a   2    + . . . + a 2 + a
         k      k-1                1     0

                            m
    where a = 1, and split x   apart into
           k

            k      k
           2      2  a             2 a    a
     m                k-1             1    0
    x  =  x     x           . . .  x     x


    Then to raise x to the mth power, do


    t( x ) = x

    return if m = 1


    for i = k-1 downto 0 do

                   2
        t(x) = t(x)  (mod f(x), p)       // Square each time.

        if a = 1 then
            i

            t(x) = x t(x) (mod f(x), p)  // Times x only if current bit is 1
        endif

    endfor
                                                        k
                                                       2
    The initial t(x) = x gets squared k times to give x  .  If a  = 1 for
                                                                i
    0 <= i <= k-1, we multiply by x which then gets squared i times more

              i
             2
    to give x .  On a binary computer, we use bit shifting and masking to

    identify the k bits  { a    . . .  a  } to the right of the leading 1
                            k-1         0

    bit.  There are log ( m ) - 1 squarings and (number of 1 bits) - 1
                           2
    multiplies.  

BUGS

     None.

--------------------------------------------------------------------------------
|                                Function Call                                 |
------------------------------------------------------------------------------*/

void
    x_to_power( bigint m, int * g, int power_table[][ MAXDEGPOLY ], int n, int p )
{

/*------------------------------------------------------------------------------
|                               Local Variables                                |
------------------------------------------------------------------------------*/

/*
    mask = 1000 ... 000.  That is, all bits of mask are zero except for the
    most significant bit of the computer word holding its value.  Signed or
    unsigned type for bigint shouldn't matter here, since we just mask and
    compare bits.
*/

bigint
    mask = ((bigint)1 << (NUMBITS - 1)) ;

int
    bit_count = 0,  /* Number of bits in m to the right of the leading bit. */
    i ;             /* Loop counter. */


/*------------------------------------------------------------------------------
|                                Function Body                                 |
------------------------------------------------------------------------------*/

#ifdef DEBUG_PP_PRIMPOLY
printf( "x to power = x ^ %lld\n", m ) ;
#endif

/*
    Initialize g(x) to x.  Exit right away if m = 1.
*/

for (i = 0 ;  i <= n-1 ;  ++i)

    g[ i ] = 0 ;

g[ 1 ] = 1 ;

if (m == 1)

    return ;

/*
    Advance the leading bit of m up to the word's left hand boundary.
    Count how many bits were to the right of the leading bit.
*/

while (! (m & mask))
{
    m <<= 1 ;
    ++bit_count ;
}

bit_count = NUMBITS - bit_count ;



/*
    Exponentiation by repeated squaring.  Discard the leading 1 bit.
    Thereafter, square for every 0 bit;  square and multiply by x for
    every 1 bit.
*/

while ( --bit_count > 0 )
{

    m <<= 1 ;       /*  Expose the next bit. */

    square( g, power_table, n, p ) ;

    #ifdef DEBUG_PP_PRIMPOLY
    printf( "    after squaring, poly = \n" ) ;
    write_poly( g, n-1 ) ;
    printf( "\n\n" );
    #endif

    if (m & mask)   /*  Leading bit is 1. */
    {
        times_x( g, power_table, n, p ) ;

        #ifdef DEBUG_PP_PRIMPOLY
        printf( "    after additional times x, poly = \n" ) ;
        write_poly( g, n-1 ) ;
        printf( "\n\n" );
        #endif
    }
}

} /* ===================== end of function x_to_power ======================= */