/*==============================================================================
| 
|  File Name:     
|
|     ppHelperFunc.c
|
|  Description:   
|
|     Higher level helper functions which don't belong elsewhere.
| 
|  Functions:
|
|     initial_trial_poly
|     next_trial_poly
|     const_coeff_test
|     const_coeff_is_primitive_root
|     skip_test
|     has_multi_irred_factors
|     generate_Q_matrix
|     find_nullity
|
|  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 <stdlib.h>
#include <memory.h>

#include "Primpoly.h"


/*==============================================================================
|                             initial_trial_poly                               |
================================================================================

DESCRIPTION

    Return the initial monic polynomial in the sequence for f(x).

INPUT
                   
     f (int *)                Monic polynomial f(x). 
     n      (int, n >= 1)     Degree of f(x).

RETURNS
                                             n
     f (int *)                Sets f( x ) = x  - 1

  EXAMPLE 
                             4
     Let n = 4.  Set f(x) = x  - 1.

METHOD

BUGS

    None.

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

void initial_trial_poly( int * f, int n )
{

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

int i ;

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

f[ 0 ] = -1 ;

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

    f[ i ] = 0 ;

f[ n ] = 1 ;

} /* ================= end of function initial_trial_poly ==================== */



/*==============================================================================
|                                next_trial_poly                               |
================================================================================

DESCRIPTION

    Return the next monic polynomial in the sequence after f(x).

INPUT
                   
    f (int *)           Monic polynomial f(x). 
    n (int, n >= 1)     Degree of monic polynomial f(x).
    p (int, p >= 2)     Modulo p coefficient arithmetic.

RETURNS

     f (int *)                Overwrites f(x) with the next polynomial 
                              after it in the sequence (explained below).
EXAMPLE 
                                           3
     Let n = 3 and p = 5.  Suppose f(x) = x  + 4 x + 4.  As a mod p number,
     this is 1 0 4 4.  Adding 1 gives 1 0 4 5.  We reduce modulo
     5 and propagate the carry to get the number 1 0 5 0.  Propagating
     the carry again and reducing gives 1 1 0 0.  The next polynomial
                     3    2
     after f(x) is  x  + x .

METHOD

     Think of the polynomial coefficients as the digits of a number written
     in base p.  The "next" polynomial is the one you would get by adding 1
     to this number in multiple precision arithmetic.  Our intention is to
     run through all possible monic polynomials modulo p.

     Propagate carries in digits 1 through n-2 when any digit exceeds p.  No
     carries take place in the n-1 st digit because our polynomial is monic.

BUGS

    None.

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

void
    next_trial_poly( int * f, int n, int p )
{

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

int
    digit_num ;   /*  Loop counter and digit number. */

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

++f[ 0 ] ;     /* Add 1, i.e. increment the coefficient of the x term. */

/*
     Sweep through the number from right to left, propagating carries.  Skip
     the constant and the nth degree terms.
*/

for (digit_num = 0 ;  digit_num <= n - 2 ;  ++digit_num)
{
    if (f[ digit_num ] == p)   /*  Propagate carry to next digit. */
    {
        f[ digit_num ] = 0 ;
        ++f[ digit_num + 1 ] ;
    }
}

} /* ================= end of function next_trial_poly ====================== */


/*==============================================================================
|                               const_coeff_test                               |
================================================================================

DESCRIPTION

                  n
  Test if a = (-1)  a  (mod p) where a  is the constant coefficient
                     0                0
  of polynomial f(x) and a is the number from the order_r() test.

INPUT

    f (int *)  nth degree monic mod p polynomial f(x). 
    n (int)    Its degree.
    p (int)    Modulus for coefficient arithmetic.
    a (int)

RETURNS

    YES    if we pass the test.
    NO     otherwise

EXAMPLE 
                                 11     3
    Let p = 5, n = 11, f( x ) = x  + 2 x + 1, and a = 4.
    Then return YES since

        11
    (-1)  * 1 (mod 5) = -1 (mod 5) = 4 (mod 5).

METHOD
                        n
    We test if (a - (-1) a ) mod p = 0.
                          0
BUGS

    None.

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

int const_coeff_test( int * f, int n, int p, int a )
{

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

int constantCoeff = f[ 0 ] ;

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

if (n % 2 != 0)
    constantCoeff = -constantCoeff ;    /* (-1)^n < 0 for odd n. */

return( (mod( a - constantCoeff, p ) == 0) ? YES : NO ) ;

} /* =================== end of function const_coeff_test ================ */



/*==============================================================================
|                         const_coeff_is_primitive_root                        |
================================================================================

DESCRIPTION

              n
  Test if (-1)  a  (mod p) is a primitive root of the prime p where
                 0
  a  is the constant term of the polynomial f(x).
   0

INPUT

    f (int *)  nth degree monic mod p polynomial f(x). 
    n (int)    Its degree.
    p (int)    Modulus for coefficient arithmetic.

EXAMPLE 
                                 11     3
    Let p = 7, n = 11, f( x ) = x  + 2 x + 4.
    Then return YES since

        11
    (-1)  * 4 (mod 7) = -4 (mod 7) = 3 (mod 7), and 3 is a
                                         1      2      3
    primitive root of the prime 7 since 3 = 3, 3 = 2, 3 = 6,
     4      5                                       7-1
    3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3   = 1 (mod 7).

METHOD

                   n
    We test if (-1) a  mod p is a primitive root mod p.
                     0
BUGS

    None.

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

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

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

int constantCoeff = f[ 0 ] ;

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

/* (-1)^n < 0 for odd n. */
if (n % 2 != 0)
    constantCoeff = -constantCoeff ;

return is_primitive_root( mod( constantCoeff, p ), p ) ;

} /* ============= of function const_coeff_is_primitive_root ================ */



/*==============================================================================
|                                    skip_test                                 |
================================================================================

DESCRIPTION

     Make the test p  | (p - 1).
                    i

INPUT

    i      (int)       ith prime divisor of r.

    primes (bigint *)  Array of distict prime factors of r.  primes[ i ] = p
                                                                            i 
    p      (int)        p >= 2.

RETURNS
                      
    YES    if the test succeeds for some i.
    NO     otherwise

EXAMPLE 
   
    Suppose i = 0, primes[ 0 ] = 2 and p = 5.  Return YES since 2 | 5-1.

METHOD

    Test if (p - 1) mod p = 0 for all i.
                         i
BUGS

    None.

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

int
    skip_test( int i, bigint * primes, int p )
{

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

/* 
    p  cannot divide the smaller number (p-1)
     i
*/
if ( (bigint)( p-1 ) < primes[ i ] )
    return NO ;
else
    return( ((bigint)(p-1) % primes[ i ] == 0) ? YES : NO ) ;

} /* ====================== end of function skip_test ======================= */




/*==============================================================================
|                              has_multi_irred_factors                         |
================================================================================

DESCRIPTION

   Find out if the monic polynomial f( x ) has two or more distinct 
   irreducible factors.

INPUT

    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.

RETURNS
                      
   1 if the monic polynomial f( x ) has two or more distinct rreducible factors
   0 otherwise.

EXAMPLE 
   
   Let n = 4, p = 5

             4
   f( x ) = x  + 2 is irreducible, so has one distinct factor.

             4    3   2                  4
   f( x ) = x + 4x + x + 4x + 1 = (x + 1)  has one distinct factor.

             3         2
   f( x ) = x  + 3 = (x + 3x + 4)(x + 2) has two distinct irreducible factors.

             4    3    2                          2
   f( x ) = x + 3x + 3x + 3x + 2 = (x + 1) (x + 2) (x + 3) has 3 distinct
   irreducible factors.

           4
   f(x) = x + 4 = (x+1)(x+2)(x+3)(x+4) has 4 distinct irreducible factors.

METHOD

      Berlekamp's method for factoring polynomials over GF( p ), modified to test
      for irreducibility only.

      See my notes;  I skip the polynomial GCD step which ensures polynomials
      are square-free due to time constraints, but this requires a proof that 
      the method is still valid.     

BUGS

    None.

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

int has_multi_irred_factors( int power_table[][ MAXDEGPOLY ], int n, int p )
{
    int ** Q ;
    int row ;
    int nullity = 0 ;


    /* Allocate space for the Q matrix. */
    Q = (int **) calloc( n, sizeof( int * ) ) ;

    for (row = 0 ;  row < n ;  ++row)
    {
        Q[ row ] = (int *)calloc( n, sizeof( int ) ) ;
    }


    /* Generate the Q-I matrix. */
    generate_Q_matrix( Q, power_table, n, p ) ;


    /* Find nullity of Q-I */
    nullity = find_nullity( Q, n, p ) ;

    /* Free space for the Q matrix. */
    for (row = 0 ;  row < n ;  ++row)
    {
        free( Q[ row ] ) ;
    }

    free( Q ) ;


    /* If nullity >= 2, f( x ) is a reducible polynomial modulo p since it has  */
    /* two or more distinct irreducible factors.                                */
    /*                                     e                                    */
    /* Nullity of 1 implies f( x ) = p( x )  for some power e >= 1 so we cannot */
    /* determine reducibility.                                                  */
    if (nullity >= 2)
        return 1 ;

    return 0 ;

} /* =============== end of function has_multi_irred_factors ================ */



/*==============================================================================
|                              generate_Q_matrix                               |
================================================================================

DESCRIPTION


    Generate n x n matrix Q - I, where rows of Q are the powers,

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

    for monic polynomial f(x).


INPUT

    Q (int **)             Memory is allocated for this matrix already and 
                           all entries are 0.

                            k
    power_table (int **)   x (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.

RETURNS
                      

EXAMPLE 
   
            4
    f(x) = x + 2, n = 4, p = 5

        (   1      )    (    1   )      ( 1    0    0    0 )
        (     5    )    (        )      (                  )
    Q = (    x     )    (   3x   )      ( 0    3    0    0 )
        (     10   ) =  (     2  )   =  (                  )
        (    x     )    (   4x   )      ( 0    0    4    0 )
        (     15   )    (     3  )      (                  )
        (    x     )    (   2x   )      ( 0    0    0    2 )

                                                    2    3
                                          1   x    x    x

             ( 0 0 0 0 )
             ( 0 2 0 0 )
    Q - I =  ( 0 0 3 0 )
             ( 0 0 0 1 )

METHOD

      Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth, 
      Addison-Wesley.

BUGS

    None.

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

void
    generate_Q_matrix( int ** Q, int power_table[][ MAXDEGPOLY ], int n, int p )
{

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

int xp[ MAXDEGPOLY ] ; /* x^p (mod f(x),p)         */
int q [ MAXDEGPOLY ] ; /* Current row of Q matrix. */

int row = 0 ;


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

/* Check for invalid inputs. */
if (n < 2 || p < 2 || Q == (int **)0)
{
    return ;
}



/* Row 0 of Q = 1. */
Q[ 0 ][ 0 ] = 1 ;


/*       p
   Find x  (mod f(x),p) save the value into q(x) and
   set row 1 of Q to this value.
*/
x_to_power( (bigint) p, xp, power_table, n, p ) ;

memcpy( q,     xp, n * sizeof( int ) ) ;
memcpy( Q[ 1 ], q, n * sizeof( int ) ) ;


/*               pk
   Row k of Q = x   (mod f(x), p) 2 <= k <= n-1, computed by
                                     p
   multiplying each previous row by x (mod f(x),p).
*/
for (row = 2 ;  row <= n-1 ;  ++row)
{
    product( q, xp, power_table, n, p ) ;
    memcpy( Q[ row ], q, n * sizeof( int ) ) ;
}


/*  Subtract Q - I */

for (row = 0 ;  row < n ;  ++row)
{
    Q[ row ][ row ] = mod( Q[ row ][ row ] - 1, p ) ;
}

return ;

} /* ================== end of function generate_Q_matrix =================== */




/*==============================================================================
|                                  find_nullity                                |
================================================================================

DESCRIPTION


INPUT

    Q (int **)          Matrix of integers mod p. 
    n (int, n >= 1)     Degree of monic polynomial f(x).
    p (int, p >= 2)     Modulo p coefficient arithmetic.

RETURNS

    Nullity of Q.  However if the nullity is 2 or more, we just return 2.

EXAMPLE 

    Let p = 5 and n = 3.  We use the facts that -1 = 4 (mod 5), 1/2 = 3, -1/2 = 2, 
    1/3 = 2, -1/3 = 3, 1/4 = 4, -1/4 = 1.

    Consider the matrix

        ( 2 3 4 )
    Q = ( 0 2 1 )
        ( 3 3 3 )

    Begin with row 1.  No pivotal columns have been selected yet.  Scan row 1 and
    pick column 1 as the pivotal column because it contains a nonzero element.

    Normalizing column 1 by multiplying by -1/pivot = -1/2 = 2 gives

        ( 4 3 4 )
        ( 0 2 1 )
        ( 1 3 3 )

     Now perform column reduction on column 2 by multiplying the pivotal column 1 
     by 3 (the column 2 element in the current row) and adding back to row 2.

        ( 4 0 4 )
        ( 0 2 1 )
        ( 1 1 3 )

     Column reduce column 3 by multiplying the pivotal column by 4 and adding back to row 3,

        ( 4 0 0 )
        ( 0 2 1 )
        ( 1 1 2 )

     For row 2, pick column 2 as the pivotal column, normalize it and reduce columns 1, then 3,

        ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )
        ( 0 2 1 ) => ( 0 4 1 ) => ( 0 4 1 ) => ( 0 4 0 )
        ( 1 1 2 )    ( 1 2 2 )    ( 1 2 2 )    ( 1 2 4 )
                 norm.         c.r. 1      c.r. 3

     For row 3, we must pick column 3 as pivotal column because we've used up columns 1 and 2,

        ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )    ( 4 0 0 )   
        ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 )  
        ( 1 2 4 )    ( 1 2 4 )    ( 1 2 4 )    ( 0 0 4 )  
                 norm.        c.r. 1        c.r. 2

     The nullity is zero, since we were always able to find a pivot in each row.

METHOD

      Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth, Addison-Wesley.

      We combine operations of normalization of columns,

      (       c       )                         (        C       )
      (       c       )                         (        C       )
      (       .       )                         (        C       )
      ( . . . q . . . ) row  ================>  ( . . . -1 . . . ) row
      (       c       )                         (        C       )
      (       c       )      column times       (        C       )
      (       c       )      -1/a modulo p      (        C       )
           pivotCol                                   pivotCol

      and column reduction,

      (        C      b       )                         (       C        B       )
      (        C      b       )                         (       C        B       )
      (        C      b       )                         (       C        B       )
      ( . . . -1 . . .e . . . ) row  ================>  ( . . . -1 . . . 0 . . . )
      (        C      b       )                         (       C        B       )
      (        C      b       )    pivotCol times       (       C        B       )
      (        C      b       )    e added back to col  (       C        B       )
           pivotCol  col                                       col


      to reduce the matrix to a form in which columns having pivots are zero until
      the pivotal row.

      The column operations don't change the left nullspace of the
      matrix.

      The matrix rank is the number of pivotal rows since they are linearly 
      independent.  The nullity is then the number of non-pivotal rows.
 
BUGS


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

int find_nullity( int ** Q, int n, int p )
{

int colFlag[ MAXDEGPOLY ] ; /* Is -1 if the column has no pivotal element. */
int nullity = 0 ;
int row, col ;
int r ;
int found = NO ;
int pivotCol = -1 ;
int q = 0 ;
int t = 0 ;

/*  Initialize column flags to -1. */
memset( colFlag, -1, n * sizeof( int ) ) ;

/* Sweep through each row. */
for (row = 0 ;  row < n ;  ++row)
{
    /* Search for a pivot in this row:  a non-zero element 
       in a column which had no previous pivot.
     */
    found = NO ;
    for (col = 0 ;  col < n ;  ++col)
    {
        if (Q[ row ][ col ] > 0 && colFlag[ col ] < 0)
        {
            found = YES ;
            pivotCol = col ;
            break ;
        }
    }

    /* No pivot;  increase nullity by 1. */
    if (found == NO)
    {
        /* Early out. */
        if (++nullity >= 2)
        {
            return nullity ;
        }
    }
    /* Found a pivot, q. */
    else
    {
        q = Q[ row ][ pivotCol ] ;

        /* Compute -1/q (mod p) */
        t = mod( -inverse_mod_p( q, p ), p ) ;

        /* Normalize the pivotal column. */
        for (r = 0 ;  r < n ;  ++r)
        {
            Q[ r ][ pivotCol ] = mod( t * Q[ r ][ pivotCol ], p ) ;
        }

        /* Do column reduction:  Add C times the pivotal column to the other
           columns where C = element in the other column at current row. */
        for (col = 0 ;  col < n ;  ++col)
        {
            if (col != pivotCol)
            {
                q = Q[ row ][ col ] ;

                for (r = 0 ;  r < n ;  ++r)
                {
                    t = mod( q * Q[ r ][ pivotCol ], p ) ;
                    Q[ r ][ col ] = mod( t + Q[ r ][ col ], p ) ;
                }
            }
        }

        /* Record the presence of a pivot in this column. */
        colFlag[ pivotCol ] = row ;

    } /* found a pivot */

} /* end for row */

return nullity ;

} /* ===================== end of function find_nullity ===================== */


#if 0
int test
    /*                         n         2n-2
        Precompute the powers x ,  ..., x     (mod f(x), p)
        for use in all later computations.
    */
n = 4 ; p = 5 ;
f[0] = 2 ; f[1] = 3 ; f[2] = 3 ; f[3] = 3 ; f[4] = 1 ;

construct_power_table( power_table, f, n, p ) ;
if (has_multi_irred_factors( power_table, n, p ) == 1)
    printf( "Pass\n" ) ;
else
    printf( "Fail\n" ) ;

#endif