source: git/Singular/minpoly.cc @ 665ca8

spielwiese
Last change on this file since 665ca8 was f613a24, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix: Singular/minpoly.{cc|h} bugfixes and optimizations on behalf of Sebastian Jambor version sent in via e-mail on June 6, 2012, 17:00:31 MESZ for master
  • Property mode set to 100644
File size: 14.9 KB
Line 
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 "config.h"
10#include<kernel/mod2.h>
11
12//#include<iomanip>
13
14#include "minpoly.h"
15
16// TODO: avoid code copying, use subclassing instead!
17
18LinearDependencyMatrix::LinearDependencyMatrix (unsigned n, unsigned long p)
19{
20  this->n = n;
21  this->p = p;
22
23  matrix = new unsigned long *[n];
24  for(int i = 0; i < n; i++)
25  {
26    matrix[i] = new unsigned long[2 * n + 1];
27  }
28  pivots = new unsigned[n];
29  tmprow = new unsigned long[2 * n + 1];
30  rows = 0;
31}
32
33LinearDependencyMatrix::~LinearDependencyMatrix ()
34{
35  delete[]tmprow;
36  delete[]pivots;
37
38  for(int i = 0; i < n; i++)
39  {
40    delete[](matrix[i]);
41  }
42  delete[]matrix;
43}
44
45void LinearDependencyMatrix::resetMatrix ()
46{
47  rows = 0;
48}
49
50int LinearDependencyMatrix::firstNonzeroEntry (unsigned long *row)
51{
52  for(int i = 0; i < n; i++)
53    if(row[i] != 0)
54      return i;
55
56  return -1;
57}
58
59void LinearDependencyMatrix::reduceTmpRow ()
60{
61  for(int i = 0; i < rows; i++)
62  {
63    unsigned piv = pivots[i];
64    unsigned x = tmprow[piv];
65    // if the corresponding entry in the row is zero, there is nothing to do
66    if(x != 0)
67    {
68      // subtract tmprow[i] times the i-th row
69      for(int j = piv; j < n + rows + 1; j++)
70      {
71        if (matrix[i][j] != 0)
72        {
73          unsigned long tmp = multMod (matrix[i][j], x, p);
74          tmp = p - tmp;
75          tmprow[j] += tmp;
76          if (tmprow[j] >= p)
77          {
78            tmprow[j] -= p;
79          }
80        }
81      }
82    }
83  }
84}
85
86
87void LinearDependencyMatrix::normalizeTmp (unsigned i)
88{
89  unsigned long inv = modularInverse (tmprow[i], p);
90  tmprow[i] = 1;
91  for(int j = i + 1; j < 2 * n + 1; j++)
92    tmprow[j] = multMod (tmprow[j], inv, p);
93}
94
95bool LinearDependencyMatrix::findLinearDependency (unsigned long *newRow,
96                                                   unsigned long *dep)
97{
98  // Copy newRow to tmprow (we need to add RHS)
99  for(int i = 0; i < n; i++)
100  {
101    tmprow[i] = newRow[i];
102    tmprow[n + i] = 0;
103  }
104  tmprow[2 * n] = 0;
105  tmprow[n + rows] = 1;
106
107  reduceTmpRow ();
108
109  // Is tmprow reduced to zero? Then we have found a linear dependence.
110  // Otherwise, we just add tmprow to the matrix.
111  unsigned newpivot = firstNonzeroEntry (tmprow);
112  if(newpivot == -1)
113  {
114    for(int i = 0; i <= n; i++)
115    {
116      dep[i] = tmprow[n + i];
117    }
118
119    return true;
120  }
121  else
122  {
123    normalizeTmp (newpivot);
124
125    for(int i = 0; i < 2 * n + 1; i++)
126    {
127      matrix[rows][i] = tmprow[i];
128    }
129
130    pivots[rows] = newpivot;
131    rows++;
132
133    return false;
134  }
135}
136
137#if 0
138std::ostream & operator<< (std::ostream & out,
139                           const LinearDependencyMatrix & mat)
140{
141  int width = ((int) log10 (mat.p)) + 1;
142
143  out << "Pivots: " << std::endl;
144  for(int j = 0; j < mat.n; j++)
145  {
146    out << std::setw (width) << mat.pivots[j] << " ";
147  }
148  out << std::endl;
149  out << "matrix:" << std::endl;
150  for(int i = 0; i < mat.rows; i++)
151  {
152    for(int j = 0; j < mat.n; j++)
153    {
154      out << std::setw (width) << mat.matrix[i][j] << " ";
155    }
156    out << "| ";
157    for(int j = mat.n; j <= 2 * mat.n; j++)
158    {
159      out << std::setw (width) << mat.matrix[i][j] << " ";
160    }
161    out << std::endl;
162  }
163  out << "tmprow: " << std::endl;
164  for(int j = 0; j < mat.n; j++)
165  {
166    out << std::setw (width) << mat.tmprow[j] << " ";
167  }
168  out << "| ";
169  for(int j = mat.n; j <= 2 * mat.n; j++)
170  {
171    out << std::setw (width) << mat.tmprow[j] << " ";
172  }
173  out << std::endl;
174
175  return out;
176}
177#endif
178
179
180NewVectorMatrix::NewVectorMatrix (unsigned n, unsigned long p)
181{
182  this->n = n;
183  this->p = p;
184
185  matrix = new unsigned long *[n];
186  for(int i = 0; i < n; i++)
187  {
188    matrix[i] = new unsigned long[n];
189  }
190
191  pivots = new unsigned[n];
192
193  nonPivots = new unsigned[n];
194
195  for (int i = 0; i < n; i++)
196  {
197     nonPivots[i] = i;
198  }
199
200  rows = 0;
201}
202
203NewVectorMatrix::~NewVectorMatrix ()
204{
205  delete nonPivots;
206  delete pivots;
207
208  for(int i = 0; i < n; i++)
209  {
210    delete[]matrix[i];
211  }
212  delete matrix;
213}
214
215int NewVectorMatrix::firstNonzeroEntry (unsigned long *row)
216{
217  for(int i = 0; i < n; i++)
218    if(row[i] != 0)
219      return i;
220
221  return -1;
222}
223
224void NewVectorMatrix::normalizeRow (unsigned long *row, unsigned i)
225{
226  unsigned long inv = modularInverse (row[i], p);
227  row[i] = 1;
228
229  for(int j = i + 1; j < n; j++)
230  {
231    row[j] = multMod (row[j], inv, p);
232  }
233}
234
235void NewVectorMatrix::insertRow (unsigned long *row)
236{
237  for(int i = 0; i < rows; i++)
238  {
239    unsigned piv = pivots[i];
240    unsigned x = row[piv];
241    // if the corresponding entry in the row is zero, there is nothing to do
242    if(x != 0)
243    {
244      // subtract row[i] times the i-th row
245      // only the non-pivot entries have to be considered
246      // (and also the first entry)
247      row[piv] = 0;
248
249      int smallestNonPivIndex = 0;
250      while (nonPivots[smallestNonPivIndex] < piv)
251      {
252        smallestNonPivIndex++;
253      }
254
255      for (int j = smallestNonPivIndex; j < n-rows; j++)
256      {
257        unsigned ind = nonPivots[j];
258        if (matrix[i][ind] != 0)
259        {
260          unsigned long tmp = multMod (matrix[i][ind], x, p);
261          tmp = p - tmp;
262          row[ind] += tmp;
263          if (row[ind] >= p)
264          {
265            row[ind] -= p;
266          }
267        }
268      }
269    }
270  }
271
272  unsigned piv = firstNonzeroEntry (row);
273
274  if(piv != -1)
275  {
276    // Normalize and insert row into the matrix.
277    // Then reduce upwards.
278    normalizeRow (row, piv);
279    for(int i = 0; i < n; i++)
280    {
281      matrix[rows][i] = row[i];
282    }
283
284
285    for (int i = 0; i < rows; i++)
286    {
287      unsigned x = matrix[i][piv];
288      // if the corresponding entry in the matrix is zero,
289      // there is nothing to do
290      if (x != 0)
291      {
292        for (int j = piv; j < n; j++)
293        {
294          if (row[j] != 0)
295          {
296            unsigned long tmp = multMod(row[j], x, p);
297            tmp = p - tmp;
298            matrix[i][j] += tmp;
299            if (matrix[i][j] >= p)
300            {
301              matrix[i][j] -= p;
302            }
303          }
304        }
305      }
306    }
307
308    pivots[rows] = piv;
309
310    // update nonPivots
311    for (int i = 0; i < n-rows; i++)
312    {
313      if (nonPivots[i] == piv)
314      {
315        // shift everything one position to the left
316        for (int j = i; j < n-rows-1; j++)
317        {
318          nonPivots[j] = nonPivots[j+1];
319        }
320
321        break;
322      }
323    }
324
325    rows++;
326  }
327}
328
329
330void NewVectorMatrix::insertMatrix (LinearDependencyMatrix & mat)
331{
332  for(int i = 0; i < mat.rows; i++)
333  {
334    insertRow (mat.matrix[i]);
335  }
336}
337
338int NewVectorMatrix::findSmallestNonpivot ()
339{
340  // This method isn't very efficient, but it is called at most a few times,
341  // so efficiency is not important.
342  if(rows == n)
343    return -1;
344
345  for(int i = 0; i < n; i++)
346  {
347    bool isPivot = false;
348    for(int j = 0; j < rows; j++)
349    {
350      if(pivots[j] == i)
351      {
352        isPivot = true;
353        break;
354      }
355    }
356
357    if(!isPivot)
358    {
359      return i;
360    }
361  }
362}
363
364int NewVectorMatrix::findLargestNonpivot ()
365{
366  // This method isn't very efficient, but it is called at most a few times, so efficiency is not important.
367  if(rows == n)
368    return -1;
369
370  for(int i = n-1; i >= 0; i--)
371  {
372    bool isPivot = false;
373    for(int j = 0; j < rows; j++)
374    {
375      if(pivots[j] == i)
376      {
377        isPivot = true;
378        break;
379      }
380    }
381
382    if(!isPivot)
383    {
384      return i;
385    }
386  }
387}
388
389
390void vectorMatrixMult (unsigned long *vec, unsigned long **mat,
391                       unsigned **nonzeroIndices, unsigned *nonzeroCounts,
392                       unsigned long *result, unsigned n, unsigned long p)
393{
394  unsigned long tmp;
395
396  for(int i = 0; i < n; i++)
397  {
398    result[i] = 0;
399    for(int j = 0; j < nonzeroCounts[i]; j++)
400    {
401      tmp = multMod (vec[nonzeroIndices[i][j]], mat[nonzeroIndices[i][j]][i], p);
402      result[i] += tmp;
403      if (result[i] >= p)
404      {
405        result[i] -= p;
406      }
407    }
408  }
409}
410
411
412#if 0
413void printVec (unsigned long *vec, int n)
414{
415  for(int i = 0; i < n; i++)
416  {
417    std::cout << vec[i] << ", ";
418  }
419
420  std::cout << std::endl;
421}
422#endif
423
424
425unsigned long *computeMinimalPolynomial (unsigned long **matrix, unsigned n,
426                                         unsigned long p)
427{
428  LinearDependencyMatrix lindepmat (n, p);
429  NewVectorMatrix newvectormat (n, p);
430
431  unsigned long *result = new unsigned long[n + 1];
432  unsigned long *mpvec = new unsigned long[n + 1];
433  unsigned long *tmp = new unsigned long[n + 1];
434
435  // initialize result = 1
436  for(int i = 0; i <= n; i++)
437  {
438    result[i] = 0;
439  }
440  result[0] = 1;
441
442  int degresult = 0;
443
444
445  // Store the indices where the matrix has non-zero entries.
446  // This has a huge impact on spares matrices.
447  unsigned* nonzeroCounts = new unsigned[n];
448  unsigned** nonzeroIndices = new unsigned*[n];
449  for (int i = 0; i < n; i++)
450  {
451    nonzeroIndices[i] = new unsigned[n];
452    nonzeroCounts[i] = 0;
453    for (int j = 0; j < n; j++)
454    {
455      if (matrix[j][i] != 0)
456      {
457        nonzeroIndices[i][nonzeroCounts[i]] = j;
458        nonzeroCounts[i]++;
459      }
460    }
461  }
462
463  int i = n-1;
464
465  unsigned long *vec = new unsigned long[n];
466  unsigned long *vecnew = new unsigned long[n];
467
468  unsigned loopsEven = true;
469  while(i != -1)
470  {
471    for(int j = 0; j < n; j++)
472    {
473      vec[j] = 0;
474    }
475    vec[i] = 1;
476
477    lindepmat.resetMatrix ();
478
479    while(true)
480    {
481      bool ld = lindepmat.findLinearDependency (vec, mpvec);
482
483      if(ld)
484      {
485        break;
486      }
487
488      vectorMatrixMult (vec, matrix, nonzeroIndices, nonzeroCounts, vecnew, n, p);
489      unsigned long *swap = vec;
490      vec = vecnew;
491      vecnew = swap;
492    }
493
494
495    unsigned degmpvec = n;
496    while(mpvec[degmpvec] == 0)
497    {
498      degmpvec--;
499    }
500
501    // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree,
502    // then we are done.
503    if(degmpvec == n)
504    {
505      unsigned long *swap = result;
506      result = mpvec;
507      mpvec = swap;
508      i = -1;
509    }
510    else
511    {
512      // otherwise, we have to compute the lcm of mpvec and prior result
513
514      // tmp = zeropol
515      for(int j = 0; j <= n; j++)
516      {
517        tmp[j] = 0;
518      }
519      degresult = lcm (tmp, result, mpvec, p, degresult, degmpvec);
520      unsigned long *swap = result;
521      result = tmp;
522      tmp = swap;
523
524      if(degresult == n)
525      {
526        // minimal polynomial has highest possible degree, so we are done.
527        i = -1;
528      }
529      else
530      {
531        newvectormat.insertMatrix (lindepmat);
532
533        // choose new unit vector from the front or the end, alternating
534        // for each round. If the matrix is the companion matrix for x^n,
535        // then taking vectors from the end is best. If it is the transpose,
536        // taking vectors from the front is best.
537        // This tries to take the middle way
538        if (loopsEven)
539        {
540          i = newvectormat.findSmallestNonpivot ();
541        }
542        else
543        {
544          i = newvectormat.findLargestNonpivot ();
545        }
546      }
547    }
548
549    loopsEven = !loopsEven;
550  }
551
552  for (int i = 0; i < n; i++)
553  {
554    delete[] nonzeroIndices[i];
555  }
556  delete[] nonzeroIndices;
557  delete[] nonzeroCounts;
558
559
560  delete[]vecnew;
561  delete[]vec;
562  delete[]tmp;
563  delete[]mpvec;
564
565  return result;
566}
567
568
569void rem (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
570          int degq)
571{
572  while(degq <= dega)
573  {
574    unsigned d = dega - degq;
575    long factor = multMod (a[dega], modularInverse (q[degq], p), p);
576    for(int i = degq; i >= 0; i--)
577    {
578      long tmp = p - multMod (factor, q[i], p);
579      a[d + i] += tmp;
580      if (a[d + i] >= p)
581      {
582        a[d + i] -= p;
583      }
584    }
585
586    while(dega >= 0 && a[dega] == 0)
587    {
588      dega--;
589    }
590  }
591}
592
593
594void quo (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
595          int degq)
596{
597  unsigned degres = dega - degq;
598  unsigned long *result = new unsigned long[degres + 1];
599
600  // initialize to zero
601  for (int i = 0; i <= degres; i++)
602  {
603    result[i] = 0;
604  }
605
606  while(degq <= dega)
607  {
608    unsigned d = dega - degq;
609    unsigned long inv = modularInverse (q[degq], p);
610    result[d] = multMod (a[dega], inv, p);
611    for(int i = degq; i >= 0; i--)
612    {
613      unsigned long tmp = p - multMod (result[d], q[i], p);
614      a[d + i] += tmp;
615      if (a[d + i] >= p)
616      {
617        a[d + i] -= p;
618      }
619    }
620
621    while(dega >= 0 && a[dega] == 0)
622    {
623      dega--;
624    }
625  }
626
627  // copy result into a
628  for(int i = 0; i <= degres; i++)
629  {
630    a[i] = result[i];
631  }
632  // set all following entries of a to zero
633  for(int i = degres + 1; i <= degq + degres; i++)
634  {
635    a[i] = 0;
636  }
637
638  dega = degres;
639
640  delete[]result;
641}
642
643
644void mult (unsigned long *result, unsigned long *a, unsigned long *b,
645           unsigned long p, int dega, int degb)
646{
647  // NOTE: we assume that every entry in result is preinitialized to zero!
648
649  for(int i = 0; i <= dega; i++)
650  {
651    for(int j = 0; j <= degb; j++)
652    {
653      result[i + j] += multMod (a[i], b[j], p);
654      if (result[i + j] >= p)
655      {
656        result[i + j] -= p;
657      }
658    }
659  }
660}
661
662
663int gcd (unsigned long *g, unsigned long *a, unsigned long *b,
664         unsigned long p, int dega, int degb)
665{
666  unsigned long *tmp1 = new unsigned long[dega + 1];
667  unsigned long *tmp2 = new unsigned long[degb + 1];
668  for(int i = 0; i <= dega; i++)
669  {
670    tmp1[i] = a[i];
671  }
672  for(int i = 0; i <= degb; i++)
673  {
674    tmp2[i] = b[i];
675  }
676  int degtmp1 = dega;
677  int degtmp2 = degb;
678
679  unsigned long *swappol;
680  int swapdeg;
681
682  while(degtmp2 >= 0)
683  {
684    rem (tmp1, tmp2, p, degtmp1, degtmp2);
685    swappol = tmp1;
686    tmp1 = tmp2;
687    tmp2 = swappol;
688
689    swapdeg = degtmp1;
690    degtmp1 = degtmp2;
691    degtmp2 = swapdeg;
692  }
693
694  for(int i = 0; i <= degtmp1; i++)
695  {
696    g[i] = tmp1[i];
697  }
698
699  delete[]tmp1;
700  delete[]tmp2;
701
702  return degtmp1;
703}
704
705
706int lcm (unsigned long *l, unsigned long *a, unsigned long *b,
707         unsigned long p, int dega, int degb)
708{
709  unsigned long *g = new unsigned long[dega + 1];
710  // initialize entries of g to zero
711  for(int i = 0; i <= dega; i++)
712  {
713    g[i] = 0;
714  }
715
716  int degg = gcd (g, a, b, p, dega, degb);
717
718  if(degg > 0)
719  {
720    // non-trivial gcd, so compute a = (a/g)
721    quo (a, g, p, dega, degg);
722  }
723  mult (l, a, b, p, dega, degb);
724
725  // normalize
726  if(l[dega + degb + 1] != 1)
727  {
728    unsigned long inv = modularInverse (l[dega + degb], p);
729    for(int i = 0; i <= dega + degb; i++)
730    {
731      l[i] = multMod (l[i], inv, p);
732    }
733  }
734
735  return dega + degb;
736}
737
738
739// computes the gcd and the cofactors of u and v
740// gcd(u,v) = u3 = u*u1 + v*u2
741unsigned long modularInverse (long long x, long long p)
742{
743  long long u1 = 1;
744  long long u2 = 0;
745  long long u3 = x;
746  long long v1 = 0;
747  long long v2 = 1;
748  long long v3 = p;
749
750  long long q, t1, t2, t3;
751  while(v3 != 0)
752  {
753    q = u3 / v3;
754    t1 = u1 - q * v1;
755    t2 = u2 - q * v2;
756    t3 = u3 - q * v3;
757    u1 = v1;
758    u2 = v2;
759    u3 = v3;
760    v1 = t1;
761    v2 = t2;
762    v3 = t3;
763  }
764
765  if(u1 < 0)
766  {
767    u1 += p;
768  }
769
770  return (unsigned long) u1;
771}
772
Note: See TracBrowser for help on using the repository browser.