My Project
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

◆ ULONG64

#define ULONG64   (unsigned long)

Function Documentation

◆ computeMinimalPolynomial()

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

Definition at line 428 of file minpoly.cc.

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

◆ 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.

668 {
669  unsigned long *tmp1 = new unsigned long[dega + 1];
670  unsigned long *tmp2 = new unsigned long[degb + 1];
671  for(int i = 0; i <= dega; i++)
672  {
673  tmp1[i] = a[i];
674  }
675  for(int i = 0; i <= degb; i++)
676  {
677  tmp2[i] = b[i];
678  }
679  int degtmp1 = dega;
680  int degtmp2 = degb;
681 
682  unsigned long *swappol;
683  int swapdeg;
684 
685  while(degtmp2 >= 0)
686  {
687  rem (tmp1, tmp2, p, degtmp1, degtmp2);
688  swappol = tmp1;
689  tmp1 = tmp2;
690  tmp2 = swappol;
691 
692  swapdeg = degtmp1;
693  degtmp1 = degtmp2;
694  degtmp2 = swapdeg;
695  }
696 
697  for(int i = 0; i <= degtmp1; i++)
698  {
699  g[i] = tmp1[i];
700  }
701 
702  delete[]tmp1;
703  delete[]tmp2;
704 
705  return degtmp1;
706 }
g
Definition: cfModGcd.cc:4092
CanonicalForm b
Definition: cfModGcd.cc:4105
CFList tmp1
Definition: facFqBivar.cc:72
CFList tmp2
Definition: facFqBivar.cc:72
void rem(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:572

◆ 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.

711 {
712  unsigned long *g = new unsigned long[dega + 1];
713  // initialize entries of g to zero
714  for(int i = 0; i <= dega; i++)
715  {
716  g[i] = 0;
717  }
718 
719  int degg = gcd (g, a, b, p, dega, degb);
720 
721  if(degg > 0)
722  {
723  // non-trivial gcd, so compute a = (a/g)
724  quo (a, g, p, dega, degg);
725  }
726  mult (l, a, b, p, dega, degb);
727 
728  // normalize
729  if(l[dega + degb + 1] != 1)
730  {
731  unsigned long inv = modularInverse (l[dega + degb], p);
732  for(int i = 0; i <= dega + degb; i++)
733  {
734  l[i] = multMod (l[i], inv, p);
735  }
736  }
737 
738  return dega + degb;
739 }
int l
Definition: cfEzgcd.cc:100
int degg
Definition: facAlgExt.cc:64
void mult(unsigned long *result, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition: minpoly.cc:647
void quo(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:597
unsigned long modularInverse(long long x, long long p)
Definition: minpoly.cc:744
int gcd(unsigned long *g, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition: minpoly.cc:666
static unsigned long multMod(unsigned long a, unsigned long b, unsigned long p)
Definition: minpoly.h:202

◆ modularInverse()

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

Definition at line 744 of file minpoly.cc.

745 {
746  long long u1 = 1;
747  long long u2 = 0;
748  long long u3 = x;
749  long long v1 = 0;
750  long long v2 = 1;
751  long long v3 = p;
752 
753  long long q, t1, t2, t3;
754  while(v3 != 0)
755  {
756  q = u3 / v3;
757  t1 = u1 - q * v1;
758  t2 = u2 - q * v2;
759  t3 = u3 - q * v3;
760  u1 = v1;
761  u2 = v2;
762  u3 = v3;
763  v1 = t1;
764  v2 = t2;
765  v3 = t3;
766  }
767 
768  if(u1 < 0)
769  {
770  u1 += p;
771  }
772 
773  return (unsigned long) u1;
774 }
Variable x
Definition: cfModGcd.cc:4084

◆ 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.

649 {
650  // NOTE: we assume that every entry in result is preinitialized to zero!
651 
652  for(int i = 0; i <= dega; i++)
653  {
654  for(int j = 0; j <= degb; j++)
655  {
656  result[i + j] += multMod (a[i], b[j], p);
657  if (result[i + j] >= p)
658  {
659  result[i + j] -= p;
660  }
661  }
662  }
663 }

◆ multMod()

static unsigned long multMod ( unsigned long  a,
unsigned long  b,
unsigned long  p 
)
inlinestatic

Definition at line 202 of file minpoly.h.

203 {
204 #if SIZEOF_LONG == 4
205 #define ULONG64 (unsigned long long)
206 #else
207 #define ULONG64 (unsigned long)
208 #endif
209  return (unsigned long)((ULONG64 a)*(ULONG64 b) % (ULONG64 p));
210 }
#define ULONG64

◆ quo()

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

Definition at line 597 of file minpoly.cc.

599 {
600  unsigned degres = dega - degq;
601  unsigned long *result = new unsigned long[degres + 1];
602 
603  // initialize to zero
604  for (int i = 0; i <= degres; i++)
605  {
606  result[i] = 0;
607  }
608 
609  while(degq <= dega)
610  {
611  unsigned d = dega - degq;
612  unsigned long inv = modularInverse (q[degq], p);
613  result[d] = multMod (a[dega], inv, p);
614  for(int i = degq; i >= 0; i--)
615  {
616  unsigned long tmp = p - multMod (result[d], q[i], p);
617  a[d + i] += tmp;
618  if (a[d + i] >= p)
619  {
620  a[d + i] -= p;
621  }
622  }
623 
624  while(dega >= 0 && a[dega] == 0)
625  {
626  dega--;
627  }
628  }
629 
630  // copy result into a
631  for(int i = 0; i <= degres; i++)
632  {
633  a[i] = result[i];
634  }
635  // set all following entries of a to zero
636  for(int i = degres + 1; i <= degq + degres; i++)
637  {
638  a[i] = 0;
639  }
640 
641  dega = degres;
642 
643  delete[]result;
644 }

◆ rem()

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

Definition at line 572 of file minpoly.cc.

574 {
575  while(degq <= dega)
576  {
577  unsigned d = dega - degq;
578  long factor = multMod (a[dega], modularInverse (q[degq], p), p);
579  for(int i = degq; i >= 0; i--)
580  {
581  long tmp = p - multMod (factor, q[i], p);
582  a[d + i] += tmp;
583  if (a[d + i] >= p)
584  {
585  a[d + i] -= p;
586  }
587  }
588 
589  while(dega >= 0 && a[dega] == 0)
590  {
591  dega--;
592  }
593  }
594 }
CanonicalForm factor
Definition: facAbsFact.cc:97

◆ 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.

396 {
397  unsigned long tmp;
398 
399  for(int i = 0; i < n; i++)
400  {
401  result[i] = 0;
402  for(int j = 0; j < nonzeroCounts[i]; j++)
403  {
404  tmp = multMod (vec[nonzeroIndices[i][j]], mat[nonzeroIndices[i][j]][i], p);
405  result[i] += tmp;
406  if (result[i] >= p)
407  {
408  result[i] -= p;
409  }
410  }
411  }
412 }