My Project
Loading...
Searching...
No Matches
minpoly.cc
Go to the documentation of this file.
1/*************************************************
2 * Author: Sebastian Jambor, 2011 *
3 * GPL (e-mail from June 6, 2012, 17:00:31 MESZ) *
4 * sebastian@momo.math.rwth-aachen.de *
5 ************************************************/
6
7
8#include <cmath>
9#include <cstdlib>
10
11
12
13#include "kernel/mod2.h"
14
15#include "minpoly.h"
16
17// TODO: avoid code copying, use subclassing instead!
18
20{
21 this->n = n;
22 this->p = p;
23
24 matrix = new unsigned long *[n];
25 for(int i = 0; i < n; i++)
26 {
27 matrix[i] = new unsigned long[2 * n + 1];
28 }
29 pivots = new unsigned[n];
30 tmprow = new unsigned long[2 * n + 1];
31 rows = 0;
32}
33
35{
36 delete[]tmprow;
37 delete[]pivots;
38
39 for(int i = 0; i < n; i++)
40 {
41 delete[](matrix[i]);
42 }
43 delete[]matrix;
44}
45
47{
48 rows = 0;
49}
50
52{
53 for(int i = 0; i < n; i++)
54 if(row[i] != 0)
55 return i;
56
57 return -1;
58}
59
61{
62 for(int i = 0; i < rows; i++)
63 {
64 unsigned piv = pivots[i];
65 unsigned x = tmprow[piv];
66 // if the corresponding entry in the row is zero, there is nothing to do
67 if(x != 0)
68 {
69 // subtract tmprow[i] times the i-th row
70 for(int j = piv; j < n + rows + 1; j++)
71 {
72 if (matrix[i][j] != 0)
73 {
74 unsigned long tmp = multMod (matrix[i][j], x, p);
75 tmp = p - tmp;
76 tmprow[j] += tmp;
77 if (tmprow[j] >= p)
78 {
79 tmprow[j] -= p;
80 }
81 }
82 }
83 }
84 }
85}
86
87
89{
90 unsigned long inv = modularInverse (tmprow[i], p);
91 tmprow[i] = 1;
92 for(int j = i + 1; j < 2 * n + 1; j++)
93 tmprow[j] = multMod (tmprow[j], inv, p);
94}
95
97 unsigned long *dep)
98{
99 // Copy newRow to tmprow (we need to add RHS)
100 for(int i = 0; i < n; i++)
101 {
102 tmprow[i] = newRow[i];
103 tmprow[n + i] = 0;
104 }
105 tmprow[2 * n] = 0;
106 tmprow[n + rows] = 1;
107
108 reduceTmpRow ();
109
110 // Is tmprow reduced to zero? Then we have found a linear dependence.
111 // Otherwise, we just add tmprow to the matrix.
112 unsigned newpivot = firstNonzeroEntry (tmprow);
113 if(newpivot == -1)
114 {
115 for(int i = 0; i <= n; i++)
116 {
117 dep[i] = tmprow[n + i];
118 }
119
120 return true;
121 }
122 else
123 {
124 normalizeTmp (newpivot);
125
126 for(int i = 0; i < 2 * n + 1; i++)
127 {
128 matrix[rows][i] = tmprow[i];
129 }
130
131 pivots[rows] = newpivot;
132 rows++;
133
134 return false;
135 }
136}
137
138#if 0
139std::ostream & operator<< (std::ostream & out,
140 const LinearDependencyMatrix & mat)
141{
142 int width = ((int) log10 (mat.p)) + 1;
143
144 out << "Pivots: " << std::endl;
145 for(int j = 0; j < mat.n; j++)
146 {
147 out << std::setw (width) << mat.pivots[j] << " ";
148 }
149 out << std::endl;
150 out << "matrix:" << std::endl;
151 for(int i = 0; i < mat.rows; i++)
152 {
153 for(int j = 0; j < mat.n; j++)
154 {
155 out << std::setw (width) << mat.matrix[i][j] << " ";
156 }
157 out << "| ";
158 for(int j = mat.n; j <= 2 * mat.n; j++)
159 {
160 out << std::setw (width) << mat.matrix[i][j] << " ";
161 }
162 out << std::endl;
163 }
164 out << "tmprow: " << std::endl;
165 for(int j = 0; j < mat.n; j++)
166 {
167 out << std::setw (width) << mat.tmprow[j] << " ";
168 }
169 out << "| ";
170 for(int j = mat.n; j <= 2 * mat.n; j++)
171 {
172 out << std::setw (width) << mat.tmprow[j] << " ";
173 }
174 out << std::endl;
175
176 return out;
177}
178#endif
179
180
181NewVectorMatrix::NewVectorMatrix (unsigned n, unsigned long p)
182{
183 this->n = n;
184 this->p = p;
185
186 matrix = new unsigned long *[n];
187 for(int i = 0; i < n; i++)
188 {
189 matrix[i] = new unsigned long[n];
190 }
191
192 pivots = new unsigned[n];
193
194 nonPivots = new unsigned[n];
195
196 for (int i = 0; i < n; i++)
197 {
198 nonPivots[i] = i;
199 }
200
201 rows = 0;
202}
203
205{
206 delete nonPivots;
207 delete pivots;
208
209 for(int i = 0; i < n; i++)
210 {
211 delete[]matrix[i];
212 }
213 delete matrix;
214}
215
217{
218 for(int i = 0; i < n; i++)
219 if(row[i] != 0)
220 return i;
221
222 return -1;
223}
224
225void NewVectorMatrix::normalizeRow (unsigned long *row, unsigned i)
226{
227 unsigned long inv = modularInverse (row[i], p);
228 row[i] = 1;
229
230 for(int j = i + 1; j < n; j++)
231 {
232 row[j] = multMod (row[j], inv, p);
233 }
234}
235
236void NewVectorMatrix::insertRow (unsigned long *row)
237{
238 for(int i = 0; i < rows; i++)
239 {
240 unsigned piv = pivots[i];
241 unsigned x = row[piv];
242 // if the corresponding entry in the row is zero, there is nothing to do
243 if(x != 0)
244 {
245 // subtract row[i] times the i-th row
246 // only the non-pivot entries have to be considered
247 // (and also the first entry)
248 row[piv] = 0;
249
250 int smallestNonPivIndex = 0;
251 while (nonPivots[smallestNonPivIndex] < piv)
252 {
253 smallestNonPivIndex++;
254 }
255
256 for (int j = smallestNonPivIndex; j < n-rows; j++)
257 {
258 unsigned ind = nonPivots[j];
259 if (matrix[i][ind] != 0)
260 {
261 unsigned long tmp = multMod (matrix[i][ind], x, p);
262 tmp = p - tmp;
263 row[ind] += tmp;
264 if (row[ind] >= p)
265 {
266 row[ind] -= p;
267 }
268 }
269 }
270 }
271 }
272
273 unsigned piv = firstNonzeroEntry (row);
274
275 if(piv != -1)
276 {
277 // Normalize and insert row into the matrix.
278 // Then reduce upwards.
279 normalizeRow (row, piv);
280 for(int i = 0; i < n; i++)
281 {
282 matrix[rows][i] = row[i];
283 }
284
285
286 for (int i = 0; i < rows; i++)
287 {
288 unsigned x = matrix[i][piv];
289 // if the corresponding entry in the matrix is zero,
290 // there is nothing to do
291 if (x != 0)
292 {
293 for (int j = piv; j < n; j++)
294 {
295 if (row[j] != 0)
296 {
297 unsigned long tmp = multMod(row[j], x, p);
298 tmp = p - tmp;
299 matrix[i][j] += tmp;
300 if (matrix[i][j] >= p)
301 {
302 matrix[i][j] -= p;
303 }
304 }
305 }
306 }
307 }
308
309 pivots[rows] = piv;
310
311 // update nonPivots
312 for (int i = 0; i < n-rows; i++)
313 {
314 if (nonPivots[i] == piv)
315 {
316 // shift everything one position to the left
317 for (int j = i; j < n-rows-1; j++)
318 {
319 nonPivots[j] = nonPivots[j+1];
320 }
321
322 break;
323 }
324 }
325
326 rows++;
327 }
328}
329
330
332{
333 for(int i = 0; i < mat.rows; i++)
334 {
335 insertRow (mat.matrix[i]);
336 }
337}
338
340{
341 // This method isn't very efficient, but it is called at most a few times,
342 // so efficiency is not important.
343 if(rows == n)
344 return -1;
345
346 for(int i = 0; i < n; i++)
347 {
348 bool isPivot = false;
349 for(int j = 0; j < rows; j++)
350 {
351 if(pivots[j] == i)
352 {
353 isPivot = true;
354 break;
355 }
356 }
357
358 if(!isPivot)
359 {
360 return i;
361 }
362 }
363 abort();
364}
365
367{
368 // This method isn't very efficient, but it is called at most a few times, so efficiency is not important.
369 if(rows == n)
370 return -1;
371
372 for(int i = n-1; i >= 0; i--)
373 {
374 bool isPivot = false;
375 for(int j = 0; j < rows; j++)
376 {
377 if(pivots[j] == i)
378 {
379 isPivot = true;
380 break;
381 }
382 }
383
384 if(!isPivot)
385 {
386 return i;
387 }
388 }
389 abort();
390}
391
392
393void vectorMatrixMult (unsigned long *vec, unsigned long **mat,
394 unsigned **nonzeroIndices, unsigned *nonzeroCounts,
395 unsigned long *result, unsigned n, unsigned long p)
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}
413
414
415#if 0
416void printVec (unsigned long *vec, int n)
417{
418 for(int i = 0; i < n; i++)
419 {
420 std::cout << vec[i] << ", ";
421 }
422
423 std::cout << std::endl;
424}
425#endif
426
427
428unsigned long *computeMinimalPolynomial (unsigned long **matrix, unsigned n,
429 unsigned long p)
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}
570
571
572void rem (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
573 int degq)
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}
595
596
597void quo (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
598 int degq)
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}
645
646
647void mult (unsigned long *result, unsigned long *a, unsigned long *b,
648 unsigned long p, int dega, int degb)
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}
664
665
666int gcd (unsigned long *g, unsigned long *a, unsigned long *b,
667 unsigned long p, int dega, int degb)
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}
707
708
709int lcm (unsigned long *l, unsigned long *a, unsigned long *b,
710 unsigned long p, int dega, int degb)
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}
740
741
742// computes the gcd and the cofactors of u and v
743// gcd(u,v) = u3 = u*u1 + v*u2
744unsigned long modularInverse (long long x, long long p)
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}
775
#define swap(_i, _j)
int l
Definition: cfEzgcd.cc:100
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
unsigned long * tmprow
Definition: minpoly.h:74
bool findLinearDependency(unsigned long *newRow, unsigned long *dep)
Definition: minpoly.cc:96
void normalizeTmp(unsigned i)
Definition: minpoly.cc:88
int firstNonzeroEntry(unsigned long *row)
Definition: minpoly.cc:51
LinearDependencyMatrix(unsigned n, unsigned long p)
Definition: minpoly.cc:19
unsigned * pivots
Definition: minpoly.h:75
unsigned long n
Definition: minpoly.h:72
unsigned long ** matrix
Definition: minpoly.h:73
void normalizeRow(unsigned long *row, unsigned i)
Definition: minpoly.cc:225
unsigned * nonPivots
Definition: minpoly.h:111
unsigned p
Definition: minpoly.h:107
NewVectorMatrix(unsigned n, unsigned long p)
Definition: minpoly.cc:181
unsigned rows
Definition: minpoly.h:112
unsigned long ** matrix
Definition: minpoly.h:109
void insertRow(unsigned long *row)
Definition: minpoly.cc:236
unsigned * pivots
Definition: minpoly.h:110
int findLargestNonpivot()
Definition: minpoly.cc:366
void insertMatrix(LinearDependencyMatrix &mat)
Definition: minpoly.cc:331
int firstNonzeroEntry(unsigned long *row)
Definition: minpoly.cc:216
unsigned long n
Definition: minpoly.h:108
int findSmallestNonpivot()
Definition: minpoly.cc:339
return result
Definition: facAbsBiFact.cc:75
CanonicalForm factor
Definition: facAbsFact.cc:97
int degg
Definition: facAlgExt.cc:64
CFList tmp1
Definition: facFqBivar.cc:72
CFList tmp2
Definition: facFqBivar.cc:72
fq_nmod_poly_t * vec
Definition: facHensel.cc:108
int j
Definition: facHensel.cc:110
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
unsigned long * computeMinimalPolynomial(unsigned long **matrix, unsigned n, unsigned long p)
Definition: minpoly.cc:428
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 gcd(unsigned long *g, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition: minpoly.cc:666
int lcm(unsigned long *l, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition: minpoly.cc:709
void rem(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:572
unsigned long modularInverse(long long x, long long p)
Definition: minpoly.cc:744
static unsigned long multMod(unsigned long a, unsigned long b, unsigned long p)
Definition: minpoly.h:202
ostream & operator<<(ostream &s, const spectrum &spec)
Definition: semic.cc:249