Data Structures | Macros | Functions
minpoly.h File Reference

Go to the source code of this file.

Data Structures

class  LinearDependencyMatrix
 
class  NewVectorMatrix
 

Macros

#define ULONG64   (unsigned long)
 

Functions

unsigned long * computeMinimalPolynomial (unsigned long **matrix, unsigned n, unsigned long p)
 
unsigned long modularInverse (long long x, long long p)
 
void vectorMatrixMult (unsigned long *vec, unsigned long **mat, unsigned **nonzeroIndices, unsigned *nonzeroCounts, unsigned long *result, 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)
 
static unsigned long multMod (unsigned long a, unsigned long b, unsigned long p)
 

Macro Definition Documentation

#define ULONG64   (unsigned long)

Function Documentation

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

Definition at line 430 of file minpoly.cc.

432 {
433  LinearDependencyMatrix lindepmat (n, p);
434  NewVectorMatrix newvectormat (n, p);
435 
436  unsigned long *result = new unsigned long[n + 1];
437  unsigned long *mpvec = new unsigned long[n + 1];
438  unsigned long *tmp = new unsigned long[n + 1];
439 
440  // initialize result = 1
441  for(int i = 0; i <= n; i++)
442  {
443  result[i] = 0;
444  }
445  result[0] = 1;
446 
447  int degresult = 0;
448 
449 
450  // Store the indices where the matrix has non-zero entries.
451  // This has a huge impact on spares matrices.
452  unsigned* nonzeroCounts = new unsigned[n];
453  unsigned** nonzeroIndices = new unsigned*[n];
454  for (int i = 0; i < n; i++)
455  {
456  nonzeroIndices[i] = new unsigned[n];
457  nonzeroCounts[i] = 0;
458  for (int j = 0; j < n; j++)
459  {
460  if (matrix[j][i] != 0)
461  {
462  nonzeroIndices[i][nonzeroCounts[i]] = j;
463  nonzeroCounts[i]++;
464  }
465  }
466  }
467 
468  int i = n-1;
469 
470  unsigned long *vec = new unsigned long[n];
471  unsigned long *vecnew = new unsigned long[n];
472 
473  unsigned loopsEven = true;
474  while(i != -1)
475  {
476  for(int j = 0; j < n; j++)
477  {
478  vec[j] = 0;
479  }
480  vec[i] = 1;
481 
482  lindepmat.resetMatrix ();
483 
484  while(true)
485  {
486  bool ld = lindepmat.findLinearDependency (vec, mpvec);
487 
488  if(ld)
489  {
490  break;
491  }
492 
493  vectorMatrixMult (vec, matrix, nonzeroIndices, nonzeroCounts, vecnew, n, p);
494  unsigned long *swap = vec;
495  vec = vecnew;
496  vecnew = swap;
497  }
498 
499 
500  unsigned degmpvec = n;
501  while(mpvec[degmpvec] == 0)
502  {
503  degmpvec--;
504  }
505 
506  // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree,
507  // then we are done.
508  if(degmpvec == n)
509  {
510  unsigned long *swap = result;
511  result = mpvec;
512  mpvec = swap;
513  i = -1;
514  }
515  else
516  {
517  // otherwise, we have to compute the lcm of mpvec and prior result
518 
519  // tmp = zeropol
520  for(int j = 0; j <= n; j++)
521  {
522  tmp[j] = 0;
523  }
524  degresult = lcm (tmp, result, mpvec, p, degresult, degmpvec);
525  unsigned long *swap = result;
526  result = tmp;
527  tmp = swap;
528 
529  if(degresult == n)
530  {
531  // minimal polynomial has highest possible degree, so we are done.
532  i = -1;
533  }
534  else
535  {
536  newvectormat.insertMatrix (lindepmat);
537 
538  // choose new unit vector from the front or the end, alternating
539  // for each round. If the matrix is the companion matrix for x^n,
540  // then taking vectors from the end is best. If it is the transpose,
541  // taking vectors from the front is best.
542  // This tries to take the middle way
543  if (loopsEven)
544  {
545  i = newvectormat.findSmallestNonpivot ();
546  }
547  else
548  {
549  i = newvectormat.findLargestNonpivot ();
550  }
551  }
552  }
553 
554  loopsEven = !loopsEven;
555  }
556 
557  for (int i = 0; i < n; i++)
558  {
559  delete[] nonzeroIndices[i];
560  }
561  delete[] nonzeroIndices;
562  delete[] nonzeroCounts;
563 
564 
565  delete[]vecnew;
566  delete[]vec;
567  delete[]tmp;
568  delete[]mpvec;
569 
570  return result;
571 }
int lcm(unsigned long *l, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition: minpoly.cc:711
return P p
Definition: myNF.cc:203
fq_nmod_poly_t * vec
Definition: facHensel.cc:103
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
void vectorMatrixMult(unsigned long *vec, unsigned long **mat, unsigned **nonzeroIndices, unsigned *nonzeroCounts, unsigned long *result, unsigned n, unsigned long p)
Definition: minpoly.cc:395
#define swap(_i, _j)
return result
Definition: facAbsBiFact.cc:76
int gcd ( unsigned long *  g,
unsigned long *  a,
unsigned long *  b,
unsigned long  p,
int  dega,
int  degb 
)

Definition at line 668 of file minpoly.cc.

670 {
671  unsigned long *tmp1 = new unsigned long[dega + 1];
672  unsigned long *tmp2 = new unsigned long[degb + 1];
673  for(int i = 0; i <= dega; i++)
674  {
675  tmp1[i] = a[i];
676  }
677  for(int i = 0; i <= degb; i++)
678  {
679  tmp2[i] = b[i];
680  }
681  int degtmp1 = dega;
682  int degtmp2 = degb;
683 
684  unsigned long *swappol;
685  int swapdeg;
686 
687  while(degtmp2 >= 0)
688  {
689  rem (tmp1, tmp2, p, degtmp1, degtmp2);
690  swappol = tmp1;
691  tmp1 = tmp2;
692  tmp2 = swappol;
693 
694  swapdeg = degtmp1;
695  degtmp1 = degtmp2;
696  degtmp2 = swapdeg;
697  }
698 
699  for(int i = 0; i <= degtmp1; i++)
700  {
701  g[i] = tmp1[i];
702  }
703 
704  delete[]tmp1;
705  delete[]tmp2;
706 
707  return degtmp1;
708 }
void rem(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:574
const poly a
Definition: syzextra.cc:212
return P p
Definition: myNF.cc:203
g
Definition: cfModGcd.cc:4031
int i
Definition: cfEzgcd.cc:123
CFList tmp2
Definition: facFqBivar.cc:70
CFList tmp1
Definition: facFqBivar.cc:70
const poly b
Definition: syzextra.cc:213
int lcm ( unsigned long *  l,
unsigned long *  a,
unsigned long *  b,
unsigned long  p,
int  dega,
int  degb 
)

Definition at line 711 of file minpoly.cc.

713 {
714  unsigned long *g = new unsigned long[dega + 1];
715  // initialize entries of g to zero
716  for(int i = 0; i <= dega; i++)
717  {
718  g[i] = 0;
719  }
720 
721  int degg = gcd (g, a, b, p, dega, degb);
722 
723  if(degg > 0)
724  {
725  // non-trivial gcd, so compute a = (a/g)
726  quo (a, g, p, dega, degg);
727  }
728  mult (l, a, b, p, dega, degb);
729 
730  // normalize
731  if(l[dega + degb + 1] != 1)
732  {
733  unsigned long inv = modularInverse (l[dega + degb], p);
734  for(int i = 0; i <= dega + degb; i++)
735  {
736  l[i] = multMod (l[i], inv, p);
737  }
738  }
739 
740  return dega + degb;
741 }
const poly a
Definition: syzextra.cc:212
int gcd(unsigned long *g, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition: minpoly.cc:668
return P p
Definition: myNF.cc:203
g
Definition: cfModGcd.cc:4031
static unsigned long multMod(unsigned long a, unsigned long b, unsigned long p)
Definition: minpoly.h:204
void quo(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:599
int i
Definition: cfEzgcd.cc:123
void mult(unsigned long *result, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition: minpoly.cc:649
unsigned long modularInverse(long long x, long long p)
Definition: minpoly.cc:746
const poly b
Definition: syzextra.cc:213
int l
Definition: cfEzgcd.cc:94
unsigned long modularInverse ( long long  x,
long long  p 
)

Definition at line 746 of file minpoly.cc.

747 {
748  long long u1 = 1;
749  long long u2 = 0;
750  long long u3 = x;
751  long long v1 = 0;
752  long long v2 = 1;
753  long long v3 = p;
754 
755  long long q, t1, t2, t3;
756  while(v3 != 0)
757  {
758  q = u3 / v3;
759  t1 = u1 - q * v1;
760  t2 = u2 - q * v2;
761  t3 = u3 - q * v3;
762  u1 = v1;
763  u2 = v2;
764  u3 = v3;
765  v1 = t1;
766  v2 = t2;
767  v3 = t3;
768  }
769 
770  if(u1 < 0)
771  {
772  u1 += p;
773  }
774 
775  return (unsigned long) u1;
776 }
return P p
Definition: myNF.cc:203
Variable x
Definition: cfModGcd.cc:4023
void mult ( unsigned long *  result,
unsigned long *  a,
unsigned long *  b,
unsigned long  p,
int  dega,
int  degb 
)

Definition at line 649 of file minpoly.cc.

651 {
652  // NOTE: we assume that every entry in result is preinitialized to zero!
653 
654  for(int i = 0; i <= dega; i++)
655  {
656  for(int j = 0; j <= degb; j++)
657  {
658  result[i + j] += multMod (a[i], b[j], p);
659  if (result[i + j] >= p)
660  {
661  result[i + j] -= p;
662  }
663  }
664  }
665 }
const poly a
Definition: syzextra.cc:212
return P p
Definition: myNF.cc:203
static unsigned long multMod(unsigned long a, unsigned long b, unsigned long p)
Definition: minpoly.h:204
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
const poly b
Definition: syzextra.cc:213
return result
Definition: facAbsBiFact.cc:76
static unsigned long multMod ( unsigned long  a,
unsigned long  b,
unsigned long  p 
)
inlinestatic

Definition at line 204 of file minpoly.h.

205 {
206 #if SIZEOF_LONG == 4
207 #define ULONG64 (unsigned long long)
208 #else
209 #define ULONG64 (unsigned long)
210 #endif
211  return (unsigned long)((ULONG64 a)*(ULONG64 b) % (ULONG64 p));
212 }
const poly a
Definition: syzextra.cc:212
#define ULONG64
return P p
Definition: myNF.cc:203
const poly b
Definition: syzextra.cc:213
void quo ( unsigned long *  a,
unsigned long *  q,
unsigned long  p,
int &  dega,
int  degq 
)

Definition at line 599 of file minpoly.cc.

601 {
602  unsigned degres = dega - degq;
603  unsigned long *result = new unsigned long[degres + 1];
604 
605  // initialize to zero
606  for (int i = 0; i <= degres; i++)
607  {
608  result[i] = 0;
609  }
610 
611  while(degq <= dega)
612  {
613  unsigned d = dega - degq;
614  unsigned long inv = modularInverse (q[degq], p);
615  result[d] = multMod (a[dega], inv, p);
616  for(int i = degq; i >= 0; i--)
617  {
618  unsigned long tmp = p - multMod (result[d], q[i], p);
619  a[d + i] += tmp;
620  if (a[d + i] >= p)
621  {
622  a[d + i] -= p;
623  }
624  }
625 
626  while(dega >= 0 && a[dega] == 0)
627  {
628  dega--;
629  }
630  }
631 
632  // copy result into a
633  for(int i = 0; i <= degres; i++)
634  {
635  a[i] = result[i];
636  }
637  // set all following entries of a to zero
638  for(int i = degres + 1; i <= degq + degres; i++)
639  {
640  a[i] = 0;
641  }
642 
643  dega = degres;
644 
645  delete[]result;
646 }
const poly a
Definition: syzextra.cc:212
return P p
Definition: myNF.cc:203
static unsigned long multMod(unsigned long a, unsigned long b, unsigned long p)
Definition: minpoly.h:204
int i
Definition: cfEzgcd.cc:123
unsigned long modularInverse(long long x, long long p)
Definition: minpoly.cc:746
return result
Definition: facAbsBiFact.cc:76
void rem ( unsigned long *  a,
unsigned long *  q,
unsigned long  p,
int &  dega,
int  degq 
)

Definition at line 574 of file minpoly.cc.

576 {
577  while(degq <= dega)
578  {
579  unsigned d = dega - degq;
580  long factor = multMod (a[dega], modularInverse (q[degq], p), p);
581  for(int i = degq; i >= 0; i--)
582  {
583  long tmp = p - multMod (factor, q[i], p);
584  a[d + i] += tmp;
585  if (a[d + i] >= p)
586  {
587  a[d + i] -= p;
588  }
589  }
590 
591  while(dega >= 0 && a[dega] == 0)
592  {
593  dega--;
594  }
595  }
596 }
const poly a
Definition: syzextra.cc:212
return P p
Definition: myNF.cc:203
static unsigned long multMod(unsigned long a, unsigned long b, unsigned long p)
Definition: minpoly.h:204
int i
Definition: cfEzgcd.cc:123
CanonicalForm factor
Definition: facAbsFact.cc:101
unsigned long modularInverse(long long x, long long p)
Definition: minpoly.cc:746
void vectorMatrixMult ( unsigned long *  vec,
unsigned long **  mat,
unsigned **  nonzeroIndices,
unsigned *  nonzeroCounts,
unsigned long *  result,
unsigned  n,
unsigned long  p 
)

Definition at line 395 of file minpoly.cc.

398 {
399  unsigned long tmp;
400 
401  for(int i = 0; i < n; i++)
402  {
403  result[i] = 0;
404  for(int j = 0; j < nonzeroCounts[i]; j++)
405  {
406  tmp = multMod (vec[nonzeroIndices[i][j]], mat[nonzeroIndices[i][j]][i], p);
407  result[i] += tmp;
408  if (result[i] >= p)
409  {
410  result[i] -= p;
411  }
412  }
413  }
414 }
return P p
Definition: myNF.cc:203
static unsigned long multMod(unsigned long a, unsigned long b, unsigned long p)
Definition: minpoly.h:204
fq_nmod_poly_t * vec
Definition: facHensel.cc:103
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
return result
Definition: facAbsBiFact.cc:76