source: git/Singular/minpoly.cc @ 3ec6bba

spielwiese
Last change on this file since 3ec6bba was f5b40a, checked in by Hans Schoenemann <hannes@…>, 13 years ago
use kernel/mod2.h everythere
  • Property mode set to 100644
File size: 11.9 KB
Line 
1/***********************************************
2 * Copyright (C) 2011 Sebastian Jambor         *
3 * sebastian@momo.math.rwth-aachen.de          *
4 ***********************************************/
5
6
7#include<cmath>
8#include<kernel/mod2.h>
9//#include<iomanip>
10
11#include "minpoly.h"
12
13// TODO: avoid code copying, use subclassing instead!
14
15LinearDependencyMatrix::LinearDependencyMatrix (unsigned n, unsigned long p)
16{
17  this->n = n;
18  this->p = p;
19
20  matrix = new unsigned long *[n];
21  for(int i = 0; i < n; i++)
22  {
23    matrix[i] = new unsigned long[2 * n + 1];
24  }
25  pivots = new unsigned[n];
26  tmprow = new unsigned long[2 * n + 1];
27  rows = 0;
28}
29
30LinearDependencyMatrix::~LinearDependencyMatrix ()
31{
32  delete[]tmprow;
33  delete[]pivots;
34
35  for(int i = 0; i < n; i++)
36  {
37    delete[](matrix[i]);
38  }
39  delete[]matrix;
40}
41
42void LinearDependencyMatrix::resetMatrix ()
43{
44  rows = 0;
45}
46
47int LinearDependencyMatrix::firstNonzeroEntry (unsigned long *row)
48{
49  for(int i = 0; i < n; i++)
50    if(row[i] != 0)
51      return i;
52
53  return -1;
54}
55
56void LinearDependencyMatrix::reduceTmpRow ()
57{
58  for(int i = 0; i < rows; i++)
59  {
60    unsigned piv = pivots[i];
61    unsigned x = tmprow[piv];
62    // if the corresponding entry in the row is zero, there is nothing to do
63    if(x != 0)
64    {
65      // subtract tmprow[i] times the i-th row
66      for(int j = piv; j < n + rows + 1; j++)
67      {
68        unsigned long tmp = multMod (matrix[i][j], x, p);
69        tmp = p - tmp;
70        tmprow[j] += tmp;
71        tmprow[j] %= p;
72      }
73    }
74  }
75}
76
77
78void LinearDependencyMatrix::normalizeTmp (unsigned i)
79{
80  unsigned long inv = modularInverse (tmprow[i], p);
81  tmprow[i] = 1;
82  for(int j = i + 1; j < 2 * n + 1; j++)
83    tmprow[j] = multMod (tmprow[j], inv, p);
84}
85
86bool LinearDependencyMatrix::findLinearDependency (unsigned long *newRow,
87                                                   unsigned long *dep)
88{
89  // Copy newRow to tmprow (we need to add RHS)
90  for(int i = 0; i < n; i++)
91  {
92    tmprow[i] = newRow[i];
93    tmprow[n + i] = 0;
94  }
95  tmprow[2 * n] = 0;
96  tmprow[n + rows] = 1;
97
98  reduceTmpRow ();
99
100  // Is tmprow reduced to zero? Then we have found a linear dependence.
101  // Otherwise, we just add tmprow to the matrix.
102  unsigned newpivot = firstNonzeroEntry (tmprow);
103  if(newpivot == -1)
104  {
105    for(int i = 0; i <= n; i++)
106    {
107      dep[i] = tmprow[n + i];
108    }
109
110    return true;
111  }
112  else
113  {
114    normalizeTmp (newpivot);
115
116    for(int i = 0; i < 2 * n + 1; i++)
117    {
118      matrix[rows][i] = tmprow[i];
119    }
120
121    pivots[rows] = newpivot;
122    rows++;
123
124    return false;
125  }
126}
127
128#if 0
129std::ostream & operator<< (std::ostream & out,
130                           const LinearDependencyMatrix & mat)
131{
132  int width = ((int) log10 (mat.p)) + 1;
133
134  out << "Pivots: " << std::endl;
135  for(int j = 0; j < mat.n; j++)
136  {
137    out << std::setw (width) << mat.pivots[j] << " ";
138  }
139  out << std::endl;
140  out << "matrix:" << std::endl;
141  for(int i = 0; i < mat.rows; i++)
142  {
143    for(int j = 0; j < mat.n; j++)
144    {
145      out << std::setw (width) << mat.matrix[i][j] << " ";
146    }
147    out << "| ";
148    for(int j = mat.n; j <= 2 * mat.n; j++)
149    {
150      out << std::setw (width) << mat.matrix[i][j] << " ";
151    }
152    out << std::endl;
153  }
154  out << "tmprow: " << std::endl;
155  for(int j = 0; j < mat.n; j++)
156  {
157    out << std::setw (width) << mat.tmprow[j] << " ";
158  }
159  out << "| ";
160  for(int j = mat.n; j <= 2 * mat.n; j++)
161  {
162    out << std::setw (width) << mat.tmprow[j] << " ";
163  }
164  out << std::endl;
165
166  return out;
167}
168#endif
169
170
171NewVectorMatrix::NewVectorMatrix (unsigned n, unsigned long p)
172{
173  this->n = n;
174  this->p = p;
175
176  matrix = new unsigned long *[n];
177  for(int i = 0; i < n; i++)
178  {
179    matrix[i] = new unsigned long[n];
180  }
181
182  pivots = new unsigned[n];
183  rows = 0;
184}
185
186NewVectorMatrix::~NewVectorMatrix ()
187{
188  delete pivots;
189
190  for(int i = 0; i < n; i++)
191  {
192    delete[]matrix[i];
193  }
194  delete matrix;
195}
196
197int NewVectorMatrix::firstNonzeroEntry (unsigned long *row)
198{
199  for(int i = 0; i < n; i++)
200    if(row[i] != 0)
201      return i;
202
203  return -1;
204}
205
206void NewVectorMatrix::normalizeRow (unsigned long *row, unsigned i)
207{
208  unsigned long inv = modularInverse (row[i], p);
209  row[i] = 1;
210
211  for(int j = i + 1; j < n; j++)
212  {
213    row[j] = multMod (row[j], inv, p);
214  }
215}
216
217void NewVectorMatrix::insertRow (unsigned long *row)
218{
219  for(int i = 0; i < rows; i++)
220  {
221    unsigned piv = pivots[i];
222    unsigned x = row[piv];
223    // if the corresponding entry in the row is zero, there is nothing to do
224    if(x != 0)
225    {
226      // subtract row[i] times the i-th row
227      for(int j = piv; j < n; j++)
228      {
229        unsigned long tmp = multMod (matrix[i][j], x, p);
230        tmp = p - tmp;
231        row[j] += tmp;
232        row[j] %= p;
233      }
234    }
235  }
236
237  unsigned piv = firstNonzeroEntry (row);
238
239  if(piv != -1)
240  {
241    // normalize and insert row into the matrix
242    normalizeRow (row, piv);
243    for(int i = 0; i < n; i++)
244    {
245      matrix[rows][i] = row[i];
246    }
247
248    pivots[rows] = piv;
249    rows++;
250  }
251}
252
253
254void NewVectorMatrix::insertMatrix (LinearDependencyMatrix & mat)
255{
256  // The matrix in LinearDependencyMatrix is already in reduced form.
257  // Thus, if the current matrix is empty, we can simply copy the other matrix.
258  if(rows == 0)
259  {
260    for(int i = 0; i < mat.rows; i++)
261    {
262      for(int j = 0; j < n; j++)
263      {
264        matrix[i][j] = mat.matrix[i][j];
265
266        rows = mat.rows;
267        for(int i = 0; i < rows; i++)
268        {
269          pivots[i] = mat.pivots[i];
270        }
271      }
272    }
273  }
274  else
275  {
276    for(int i = 0; i < mat.rows; i++)
277    {
278      insertRow (mat.matrix[i]);
279    }
280  }
281}
282
283int NewVectorMatrix::findSmallestNonpivot ()
284{
285  // This method isn't very efficient, but it is called at most a few times, so efficiency is not important.
286  if(rows == n)
287    return -1;
288
289  for(int i = 0; i < n; i++)
290  {
291    bool isPivot = false;
292    for(int j = 0; j < rows; j++)
293    {
294      if(pivots[j] == i)
295      {
296        isPivot = true;
297        break;
298      }
299    }
300
301    if(!isPivot)
302    {
303      return i;
304    }
305  }
306}
307
308
309void vectorMatrixMult (unsigned long *vec, unsigned long **mat,
310                       unsigned long *result, unsigned n, unsigned long p)
311{
312  unsigned long tmp;
313  for(int i = 0; i < n; i++)
314  {
315    result[i] = 0;
316    for(int j = 0; j < n; j++)
317    {
318      tmp = multMod (vec[j], mat[j][i], p);
319      result[i] += tmp;
320      result[i] %= p;
321    }
322  }
323}
324
325
326#if 0
327void printVec (unsigned long *vec, int n)
328{
329  for(int i = 0; i < n; i++)
330  {
331    std::cout << vec[i] << ", ";
332  }
333
334  std::cout << std::endl;
335}
336#endif
337
338
339unsigned long *computeMinimalPolynomial (unsigned long **matrix, unsigned n,
340                                         unsigned long p)
341{
342  LinearDependencyMatrix lindepmat (n, p);
343  NewVectorMatrix newvectormat (n, p);
344
345  unsigned long *result = new unsigned long[n + 1];
346  unsigned long *mpvec = new unsigned long[n + 1];
347  unsigned long *tmp = new unsigned long[n + 1];
348
349  // initialize result = 1
350  for(int i = 0; i <= n; i++)
351  {
352    result[i] = 0;
353  }
354  result[0] = 1;
355
356  int degresult = 0;
357
358
359  int i = 0;
360
361  unsigned long *vec = new unsigned long[n];
362  unsigned long *vecnew = new unsigned long[n];
363
364  while(i != -1)
365  {
366    for(int j = 0; j < n; j++)
367    {
368      vec[j] = 0;
369    }
370    vec[i] = 1;
371
372    lindepmat.resetMatrix ();
373
374    while(true)
375    {
376      bool ld = lindepmat.findLinearDependency (vec, mpvec);
377
378      if(ld)
379      {
380        break;
381      }
382
383      vectorMatrixMult (vec, matrix, vecnew, n, p);
384      unsigned long *swap = vec;
385      vec = vecnew;
386      vecnew = swap;
387    }
388
389
390    unsigned degmpvec = n;
391    while(mpvec[degmpvec] == 0)
392    {
393      degmpvec--;
394    }
395
396    // if the dimension is already maximal, i.e., the minimal polynomial of vec has the highest possible degree,
397    // then we are done.
398    if(degmpvec == n)
399    {
400      unsigned long *swap = result;
401      result = mpvec;
402      mpvec = swap;
403      i = -1;
404    }
405    else
406    {
407      // otherwise, we have to compute the lcm of mpvec and prior result
408
409      // tmp = zeropol
410      for(int j = 0; j <= n; j++)
411      {
412        tmp[j] = 0;
413      }
414      degresult = lcm (tmp, result, mpvec, p, degresult, degmpvec);
415      unsigned long *swap = result;
416      result = tmp;
417      tmp = swap;
418
419      if(degresult == n)
420      {
421        // minimal polynomial has highest possible degree, so we are done.
422        i = -1;
423      }
424      else
425      {
426        newvectormat.insertMatrix (lindepmat);
427        i = newvectormat.findSmallestNonpivot ();
428      }
429    }
430  }
431
432  // TODO: take lcms of the different monomials!
433
434  delete[]vecnew;
435  delete[]vec;
436  delete[]tmp;
437  delete[]mpvec;
438
439  return result;
440}
441
442
443void rem (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
444          int degq)
445{
446  while(degq <= dega)
447  {
448    unsigned d = dega - degq;
449    long factor = multMod (a[dega], modularInverse (q[degq], p), p);
450    for(int i = degq; i >= 0; i--)
451    {
452      long tmp = p - multMod (factor, q[i], p);
453      a[d + i] += tmp;
454      a[d + i] %= p;
455    }
456
457    while(a[dega] == 0 && dega >= 0)
458    {
459      dega--;
460    }
461  }
462}
463
464
465void quo (unsigned long *a, unsigned long *q, unsigned long p, int &dega,
466          int degq)
467{
468  unsigned degres = dega - degq;
469  unsigned long *result = new unsigned long[degres + 1];
470
471  while(degq <= dega)
472  {
473    unsigned d = dega - degq;
474    long inv = modularInverse (q[degq], p);
475    result[d] = multMod (a[dega], inv, p);
476    for(int i = degq; i >= 0; i--)
477    {
478      long tmp = p - multMod (result[d], q[i], p);
479      a[d + i] += tmp;
480      a[d + i] %= p;
481    }
482
483    while(a[dega] == 0 && dega >= 0)
484    {
485      dega--;
486    }
487  }
488
489  // copy result into a
490  for(int i = 0; i <= degres; i++)
491  {
492    a[i] = result[i];
493  }
494  // set all following entries of a to zero
495  for(int i = degres + 1; i <= degq + degres; i++)
496  {
497    a[i] = 0;
498  }
499
500  dega = degres;
501
502  delete[]result;
503}
504
505
506void mult (unsigned long *result, unsigned long *a, unsigned long *b,
507           unsigned long p, int dega, int degb)
508{
509  // NOTE: we assume that every entry in result is preinitialized to zero!
510
511  for(int i = 0; i <= dega; i++)
512  {
513    for(int j = 0; j <= degb; j++)
514    {
515      result[i + j] += multMod (a[i], b[j], p);
516      result[i + j] %= p;
517    }
518  }
519}
520
521
522int gcd (unsigned long *g, unsigned long *a, unsigned long *b,
523         unsigned long p, int dega, int degb)
524{
525  unsigned long *tmp1 = new unsigned long[dega + 1];
526  unsigned long *tmp2 = new unsigned long[degb + 1];
527  for(int i = 0; i <= dega; i++)
528  {
529    tmp1[i] = a[i];
530  }
531  for(int i = 0; i <= degb; i++)
532  {
533    tmp2[i] = b[i];
534  }
535  int degtmp1 = dega;
536  int degtmp2 = degb;
537
538  unsigned long *swappol;
539  int swapdeg;
540
541  while(degtmp2 >= 0)
542  {
543    rem (tmp1, tmp2, p, degtmp1, degtmp2);
544    swappol = tmp1;
545    tmp1 = tmp2;
546    tmp2 = swappol;
547
548    swapdeg = degtmp1;
549    degtmp1 = degtmp2;
550    degtmp2 = swapdeg;
551  }
552
553  for(int i = 0; i <= degtmp1; i++)
554  {
555    g[i] = tmp1[i];
556  }
557
558  delete[]tmp1;
559  delete[]tmp2;
560
561  return degtmp1;
562}
563
564
565int lcm (unsigned long *l, unsigned long *a, unsigned long *b,
566         unsigned long p, int dega, int degb)
567{
568  unsigned long *g = new unsigned long[dega + 1];
569  // initialize entries of g to zero
570  for(int i = 0; i <= dega; i++)
571  {
572    g[i] = 0;
573  }
574
575  int degg = gcd (g, a, b, p, dega, degb);
576
577  if(degg > 0)
578  {
579    // non-trivial gcd, so compute a = (a/g)
580    quo (a, g, p, dega, degg);
581  }
582  mult (l, a, b, p, dega, degb);
583
584  // normalize
585  if(l[dega + degb + 1] != 1)
586  {
587    unsigned long inv = modularInverse (l[dega + degb], p);
588    for(int i = 0; i <= dega + degb; i++)
589    {
590      l[i] = multMod (l[i], inv, p);
591    }
592  }
593
594  return dega + degb;
595}
596
597
598// computes the gcd and the cofactors of u and v
599// gcd(u,v) = u3 = u*u1 + v*u2
600unsigned long modularInverse (long long x, long long p)
601{
602  long long u1 = 1;
603  long long u2 = 0;
604  long long u3 = x;
605  long long v1 = 0;
606  long long v2 = 1;
607  long long v3 = p;
608
609  long long q, t1, t2, t3;
610  while(v3 != 0)
611  {
612    q = u3 / v3;
613    t1 = u1 - q * v1;
614    t2 = u2 - q * v2;
615    t3 = u3 - q * v3;
616    u1 = v1;
617    u2 = v2;
618    u3 = v3;
619    v1 = t1;
620    v2 = t2;
621    v3 = t3;
622  }
623
624  if(u1 < 0)
625  {
626    u1 += p;
627  }
628
629  return (unsigned long) u1;
630}
631
Note: See TracBrowser for help on using the repository browser.