My Project
Loading...
Searching...
No Matches
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:4078
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:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
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:4082

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