minpoly.cc File Reference

#include <cmath>
#include <cstdlib>
#include "kernel/mod2.h"
#include "minpoly.h"

Go to the source code of this file.

Functions
void	vectorMatrixMult (unsigned long *vec, unsigned long **mat, unsigned **nonzeroIndices, unsigned *nonzeroCounts, unsigned long *result, unsigned n, unsigned long p)
unsigned long *	computeMinimalPolynomial (unsigned long **matrix, unsigned n, unsigned long p)
void	rem (unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
void	quo (unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
void	mult (unsigned long *result, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
int	gcd (unsigned long *g, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
int	lcm (unsigned long *l, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
unsigned long	modularInverse (long long x, long long p)

Function Documentation

computeMinimalPolynomial()

unsigned long * computeMinimalPolynomial	(	unsigned long **	matrix,
		unsigned	n,
		unsigned long	p
	)

Definition at line 428 of file minpoly.cc.

{
  LinearDependencyMatrix lindepmat (n, p);
  NewVectorMatrix newvectormat (n, p);
 
  unsigned long *result = new unsigned long[n + 1];
  unsigned long *mpvec = new unsigned long[n + 1];
  unsigned long *tmp = new unsigned long[n + 1];
 
  // initialize result = 1
  for(int i = 0; i <= n; i++)
  {
    result[i] = 0;
  }
  result[0] = 1;
 
  int degresult = 0;
 
 
  // Store the indices where the matrix has non-zero entries.
  // This has a huge impact on spares matrices.
  unsigned* nonzeroCounts = new unsigned[n];
  unsigned** nonzeroIndices = new unsigned*[n];
  for (int i = 0; i < n; i++)
  {
    nonzeroIndices[i] = new unsigned[n];
    nonzeroCounts[i] = 0;
    for (int j = 0; j < n; j++)
    {
      if (matrix[j][i] != 0)
      {
        nonzeroIndices[i][nonzeroCounts[i]] = j;
        nonzeroCounts[i]++;
      }
    }
  }
 
  int i = n-1;
 
  unsigned long *vec = new unsigned long[n];
  unsigned long *vecnew = new unsigned long[n];
 
  unsigned loopsEven = true;
  while(i != -1)
  {
    for(int j = 0; j < n; j++)
    {
      vec[j] = 0;
    }
    vec[i] = 1;
 
    lindepmat.resetMatrix ();
 
    while(true)
    {
      bool ld = lindepmat.findLinearDependency (vec, mpvec);
 
      if(ld)
      {
        break;
      }
 
      vectorMatrixMult (vec, matrix, nonzeroIndices, nonzeroCounts, vecnew, n, p);
      unsigned long *swap = vec;
      vec = vecnew;
      vecnew = swap;
    }
 
 
    unsigned degmpvec = n;
    while(mpvec[degmpvec] == 0)
    {
      degmpvec--;
    }
 
    // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree,
    // then we are done.
    if(degmpvec == n)
    {
      unsigned long *swap = result;
      result = mpvec;
      mpvec = swap;
      i = -1;
    }
    else
    {
      // otherwise, we have to compute the lcm of mpvec and prior result
 
      // tmp = zeropol
      for(int j = 0; j <= n; j++)
      {
        tmp[j] = 0;
      }
      degresult = lcm (tmp, result, mpvec, p, degresult, degmpvec);
      unsigned long *swap = result;
      result = tmp;
      tmp = swap;
 
      if(degresult == n)
      {
        // minimal polynomial has highest possible degree, so we are done.
        i = -1;
      }
      else
      {
        newvectormat.insertMatrix (lindepmat);
 
        // choose new unit vector from the front or the end, alternating
        // for each round. If the matrix is the companion matrix for x^n,
        // then taking vectors from the end is best. If it is the transpose,
        // taking vectors from the front is best.
        // This tries to take the middle way
        if (loopsEven)
        {
          i = newvectormat.findSmallestNonpivot ();
        }
        else
        {
          i = newvectormat.findLargestNonpivot ();
        }
      }
    }
 
    loopsEven = !loopsEven;
  }
 
  for (int i = 0; i < n; i++)
  {
    delete[] nonzeroIndices[i];
  }
  delete[] nonzeroIndices;
  delete[] nonzeroCounts;
 
 
  delete[]vecnew;
  delete[]vec;
  delete[]tmp;
  delete[]mpvec;
 
  return result;
}

gcd()

int gcd	(	unsigned long *	g,
		unsigned long *	a,
		unsigned long *	b,
		unsigned long	p,
		int	dega,
		int	degb
	)

Definition at line 666 of file minpoly.cc.

{
  unsigned long *tmp1 = new unsigned long[dega + 1];
  unsigned long *tmp2 = new unsigned long[degb + 1];
  for(int i = 0; i <= dega; i++)
  {
    tmp1[i] = a[i];
  }
  for(int i = 0; i <= degb; i++)
  {
    tmp2[i] = b[i];
  }
  int degtmp1 = dega;
  int degtmp2 = degb;
 
  unsigned long *swappol;
  int swapdeg;
 
  while(degtmp2 >= 0)
  {
    rem (tmp1, tmp2, p, degtmp1, degtmp2);
    swappol = tmp1;
    tmp1 = tmp2;
    tmp2 = swappol;
 
    swapdeg = degtmp1;
    degtmp1 = degtmp2;
    degtmp2 = swapdeg;
  }
 
  for(int i = 0; i <= degtmp1; i++)
  {
    g[i] = tmp1[i];
  }
 
  delete[]tmp1;
  delete[]tmp2;
 
  return degtmp1;
}

lcm()

int lcm	(	unsigned long *	l,
		unsigned long *	a,
		unsigned long *	b,
		unsigned long	p,
		int	dega,
		int	degb
	)

Definition at line 709 of file minpoly.cc.

{
  unsigned long *g = new unsigned long[dega + 1];
  // initialize entries of g to zero
  for(int i = 0; i <= dega; i++)
  {
    g[i] = 0;
  }
 
  int degg = gcd (g, a, b, p, dega, degb);
 
  if(degg > 0)
  {
    // non-trivial gcd, so compute a = (a/g)
    quo (a, g, p, dega, degg);
  }
  mult (l, a, b, p, dega, degb);
 
  // normalize
  if(l[dega + degb + 1] != 1)
  {
    unsigned long inv = modularInverse (l[dega + degb], p);
    for(int i = 0; i <= dega + degb; i++)
    {
      l[i] = multMod (l[i], inv, p);
    }
  }
 
  return dega + degb;
}

modularInverse()

unsigned long modularInverse	(	long long	x,
		long long	p
	)

Definition at line 744 of file minpoly.cc.

{
  long long u1 = 1;
  long long u2 = 0;
  long long u3 = x;
  long long v1 = 0;
  long long v2 = 1;
  long long v3 = p;
 
  long long q, t1, t2, t3;
  while(v3 != 0)
  {
    q = u3 / v3;
    t1 = u1 - q * v1;
    t2 = u2 - q * v2;
    t3 = u3 - q * v3;
    u1 = v1;
    u2 = v2;
    u3 = v3;
    v1 = t1;
    v2 = t2;
    v3 = t3;
  }
 
  if(u1 < 0)
  {
    u1 += p;
  }
 
  return (unsigned long) u1;
}

mult()

void mult	(	unsigned long *	result,
		unsigned long *	a,
		unsigned long *	b,
		unsigned long	p,
		int	dega,
		int	degb
	)

Definition at line 647 of file minpoly.cc.

{
  // NOTE: we assume that every entry in result is preinitialized to zero!
 
  for(int i = 0; i <= dega; i++)
  {
    for(int j = 0; j <= degb; j++)
    {
      result[i + j] += multMod (a[i], b[j], p);
      if (result[i + j] >= p)
      {
        result[i + j] -= p;
      }
    }
  }
}

quo()

void quo	(	unsigned long *	a,
		unsigned long *	q,
		unsigned long	p,
		int &	dega,
		int	degq
	)

Definition at line 597 of file minpoly.cc.

{
  unsigned degres = dega - degq;
  unsigned long *result = new unsigned long[degres + 1];
 
  // initialize to zero
  for (int i = 0; i <= degres; i++)
  {
    result[i] = 0;
  }
 
  while(degq <= dega)
  {
    unsigned d = dega - degq;
    unsigned long inv = modularInverse (q[degq], p);
    result[d] = multMod (a[dega], inv, p);
    for(int i = degq; i >= 0; i--)
    {
      unsigned long tmp = p - multMod (result[d], q[i], p);
      a[d + i] += tmp;
      if (a[d + i] >= p)
      {
        a[d + i] -= p;
      }
    }
 
    while(dega >= 0 && a[dega] == 0)
    {
      dega--;
    }
  }
 
  // copy result into a
  for(int i = 0; i <= degres; i++)
  {
    a[i] = result[i];
  }
  // set all following entries of a to zero
  for(int i = degres + 1; i <= degq + degres; i++)
  {
    a[i] = 0;
  }
 
  dega = degres;
 
  delete[]result;
}

rem()

void rem	(	unsigned long *	a,
		unsigned long *	q,
		unsigned long	p,
		int &	dega,
		int	degq
	)

Definition at line 572 of file minpoly.cc.

{
  while(degq <= dega)
  {
    unsigned d = dega - degq;
    long factor = multMod (a[dega], modularInverse (q[degq], p), p);
    for(int i = degq; i >= 0; i--)
    {
      long tmp = p - multMod (factor, q[i], p);
      a[d + i] += tmp;
      if (a[d + i] >= p)
      {
        a[d + i] -= p;
      }
    }
 
    while(dega >= 0 && a[dega] == 0)
    {
      dega--;
    }
  }
}

vectorMatrixMult()

void vectorMatrixMult	(	unsigned long *	vec,
		unsigned long **	mat,
		unsigned **	nonzeroIndices,
		unsigned *	nonzeroCounts,
		unsigned long *	result,
		unsigned	n,
		unsigned long	p
	)

Definition at line 393 of file minpoly.cc.

{
  unsigned long tmp;
 
  for(int i = 0; i < n; i++)
  {
    result[i] = 0;
    for(int j = 0; j < nonzeroCounts[i]; j++)
    {
      tmp = multMod (vec[nonzeroIndices[i][j]], mat[nonzeroIndices[i][j]][i], p);
      result[i] += tmp;
      if (result[i] >= p)
      {
        result[i] -= p;
      }
    }
  }
}