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