source: git/libpolys/polys/linalg_from_matpol.cc @ 2fce0e

spielwiese
Last change on this file since 2fce0e was 2fce0e, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: exposed mp_Wedge FIX: moving duplicates from linalg_from_matpol.cc to matpol.cc FIX: commented out unused(?!) stuff from linalg_from_matpol.cc in matpol.h TODO: match&clean up!
  • Property mode set to 100644
File size: 17.5 KB
Line 
1typedef int perm[100];
2static void mpReplace(int j, int n, int &sign, int *perm);
3static int mpNextperm(perm * z, int max);
4static poly mp_Leibnitz(matrix a, const ring);
5static poly minuscopy (poly p, const ring);
6static poly p_Insert(poly p1, poly p2, const ring);
7
8static void mp_PartClean(matrix, int, int, const ring);
9static void mp_FinalClean(matrix, const ring);
10static int mp_PrepareRow (matrix, int, int, const ring);
11static int mp_PreparePiv (matrix, int, int, const ring);
12static int mp_PivBar(matrix, int, int, const ring);
13static int mp_PivRow(matrix, int, int, const ring);
14static float mp_PolyWeight(poly, const ring);
15static void mp_SwapRow(matrix, int, int, int, const ring);
16static void mp_SwapCol(matrix, int, int, int, const ring);
17static void mp_ElimBar(matrix, matrix, poly, int, int, const ring);
18
19/*2
20*  prepare one step of 'Bareiss' algorithm
21*  for application in minor
22*/
23static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
24{
25  int r;
26
27  r = mp_PivBar(a,lr,lc, R);
28  if(r==0) return 0;
29  if(r<lr) mp_SwapRow(a, r, lr, lc, R);
30  return 1;
31}
32
33/*2
34*  prepare one step of 'Bareiss' algorithm
35*  for application in minor
36*/
37static int mp_PreparePiv (matrix a, int lr, int lc, const ring R)
38{
39  int c;
40
41  c = mp_PivRow(a, lr, lc, R);
42  if(c==0) return 0;
43  if(c<lc) mp_SwapCol(a, c, lr, lc, R);
44  return 1;
45}
46
47/*
48* find best row
49*/
50static int mp_PivBar(matrix a, int lr, int lc, const ring R)
51{
52  float f1, f2;
53  poly *q1;
54  int i,j,io;
55
56  io = -1;
57  f1 = 1.0e30;
58  for (i=lr-1;i>=0;i--)
59  {
60    q1 = &(a->m)[i*a->ncols];
61    f2 = 0.0;
62    for (j=lc-1;j>=0;j--)
63    {
64      if (q1[j]!=NULL)
65        f2 += mp_PolyWeight(q1[j], R);
66    }
67    if ((f2!=0.0) && (f2<f1))
68    {
69      f1 = f2;
70      io = i;
71    }
72  }
73  if (io<0) return 0;
74  else return io+1;
75}
76
77/*
78* find pivot in the last row
79*/
80static int mp_PivRow(matrix a, int lr, int lc, const ring R)
81{
82  float f1, f2;
83  poly *q1;
84  int j,jo;
85
86  jo = -1;
87  f1 = 1.0e30;
88  q1 = &(a->m)[(lr-1)*a->ncols];
89  for (j=lc-1;j>=0;j--)
90  {
91    if (q1[j]!=NULL)
92    {
93      f2 = mp_PolyWeight(q1[j], R);
94      if (f2<f1)
95      {
96        f1 = f2;
97        jo = j;
98      }
99    }
100  }
101  if (jo<0) return 0;
102  else return jo+1;
103}
104
105/*
106* weigth of a polynomial, for pivot strategy
107*/
108static float mp_PolyWeight(poly p, const ring R)
109{
110  int i;
111  float res;
112
113  if (pNext(p) == NULL)
114  {
115    res = (float)n_Size(p_GetCoeff(p, R), R);
116    for (i=rVar(R);i>0;i--)
117    {
118      if(p_GetExp(p,i, R)!=0)
119      {
120        res += 2.0;
121        break;
122      }
123    }
124  }
125  else
126  {
127    res = 0.0;
128    do
129    {
130      res += (float)n_Size(p_GetCoeff(p, R), R) + 2.0;
131      pIter(p);
132    }
133    while (p);
134  }
135  return res;
136}
137
138static void mpSwapRow(matrix a, int pos, int lr, int lc)
139{
140  poly sw;
141  int j;
142  poly* a2 = a->m;
143  poly* a1 = &a2[a->ncols*(pos-1)];
144
145  a2 = &a2[a->ncols*(lr-1)];
146  for (j=lc-1; j>=0; j--)
147  {
148    sw = a1[j];
149    a1[j] = a2[j];
150    a2[j] = sw;
151  }
152}
153
154static void mpSwapCol(matrix a, int pos, int lr, int lc)
155{
156  poly sw;
157  int j;
158  poly* a2 = a->m;
159  poly* a1 = &a2[pos-1];
160
161  a2 = &a2[lc-1];
162  for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
163  {
164    sw = a1[j];
165    a1[j] = a2[j];
166    a2[j] = sw;
167  }
168}
169
170/*
171* C++ classes for Bareiss algorithm
172*/
173class row_col_weight
174{
175  private:
176    int ym, yn;
177  public:
178    float *wrow, *wcol;
179    row_col_weight() : ym(0) {}
180    row_col_weight(int, int);
181    ~row_col_weight();
182};
183
184/*2
185*  a submatrix M of a matrix X[m,n]:
186*    0 <= i < s_m <= a_m
187*    0 <= j < s_n <= a_n
188*    M = ( Xarray[qrow[i],qcol[j]] )
189*    if a_m = a_n and s_m = s_n
190*      det(X) = sign*div^(s_m-1)*det(M)
191*    resticted pivot for elimination
192*      0 <= j < piv_s
193*/
194class mp_permmatrix
195{
196  private:
197    int       a_m, a_n, s_m, s_n, sign, piv_s;
198    int       *qrow, *qcol;
199    poly      *Xarray;
200    ring      R;
201
202    void mpInitMat();
203    poly * mpRowAdr(int);
204    poly * mpColAdr(int);
205    void mpRowWeight(float *);
206    void mpColWeight(float *);
207    void mpRowSwap(int, int);
208    void mpColSwap(int, int);
209  public:
210    mp_permmatrix() : a_m(0), R(NULL) {}
211    mp_permmatrix(matrix, const ring);
212    mp_permmatrix(mp_permmatrix *);
213    ~mp_permmatrix();
214    int mpGetRow();
215    int mpGetCol();
216    int mpGetRdim();
217    int mpGetCdim();
218    int mpGetSign();
219    void mpSetSearch(int s);
220    void mpSaveArray();
221    poly mpGetElem(int, int);
222    void mpSetElem(poly, int, int);
223    void mpDelElem(int, int);
224    void mpElimBareiss(poly);
225    int mpPivotBareiss(row_col_weight *);
226    int mpPivotRow(row_col_weight *, int);
227    void mpToIntvec(intvec *);
228    void mpRowReorder();
229    void mpColReorder();
230};
231
232
233
234/*2
235*returns the determinant of the matrix m;
236*uses Newtons formulea for symmetric functions
237*/
238poly mp_Det (matrix m, const ring R)
239{
240  int i,j,k,n;
241  poly p,q;
242  matrix a, s;
243  matrix ma[100];
244  number c=NULL, d=NULL, ONE=NULL;
245
246  n = MATROWS(m);
247  if (n != MATCOLS(m))
248  {
249    Werror("det of %d x %d matrix",n,MATCOLS(m));
250    return NULL;
251  }
252  k=rChar(R);
253  if ((k > 0) && (k <= n))
254    return mp_Leibnitz(m, R);
255  ONE = n_Init(1, R);
256  ma[1]=mp_Copy(m, R);
257  k = (n+1) / 2;
258  s = mpNew(1, n);
259  MATELEM(s,1,1) = mp_Trace(m, R);
260  for (i=2; i<=k; i++)
261  {
262    //ma[i] = mpNew(n,n);
263    ma[i]=mp_Mult(ma[i-1], ma[1], R);
264    MATELEM(s,1,i) = mp_Trace(ma[i], R);
265    p_Test(MATELEM(s,1,i), R);
266  }
267  for (i=k+1; i<=n; i++)
268  {
269    MATELEM(s,1,i) = TraceOfProd(ma[i / 2], ma[(i+1) / 2], n, R);
270    p_Test(MATELEM(s,1,i), R);
271  }
272  for (i=1; i<=k; i++)
273    id_Delete((ideal *)&(ma[i]), R);
274/* the array s contains the traces of the powers of the matrix m,
275*  these are the power sums of the eigenvalues of m */
276  a = mpNew(1,n);
277  MATELEM(a,1,1) = minuscopy(MATELEM(s,1,1), R);
278  for (i=2; i<=n; i++)
279  {
280    p = p_Copy(MATELEM(s,1,i), R);
281    for (j=i-1; j>=1; j--)
282    {
283      q = pp_Mult_qq(MATELEM(s,1,j), MATELEM(a,1,i-j), R);
284      p_Test(q, R);
285      p = p_Add_q(p,q, R);
286    }
287    // c= -1/i
288    d = n_Init(-(int)i, R);
289    c = n_Div(ONE, d, R);
290    n_Delete(&d, R);
291
292    p_Mult_nn(p, c, R);
293    p_Test(p, R);
294    MATELEM(a,1,i) = p;
295    n_Delete(&c, R);
296  }
297/* the array a contains the elementary symmetric functions of the
298*  eigenvalues of m */
299  for (i=1; i<=n-1; i++)
300  {
301    //p_Delete(&(MATELEM(a,1,i)), R);
302    p_Delete(&(MATELEM(s,1,i)), R);
303  }
304  p_Delete(&(MATELEM(s,1,n)), R);
305/* up to a sign, the determinant is the n-th elementary symmetric function */
306  if ((n/2)*2 < n)
307  {
308    d = n_Init(-1, R);
309    p_Mult_nn(MATELEM(a,1,n), d, R);
310    n_Delete(&d, R);
311  }
312  n_Delete(&ONE, R);
313  id_Delete((ideal *)&s, R);
314  poly result=MATELEM(a,1,n);
315  MATELEM(a,1,n)=NULL;
316  id_Delete((ideal *)&a, R);
317  return result;
318}
319
320
321///*2
322//*homogenize all elements of matrix (not the matrix itself)
323//*/
324//matrix mpHomogen(matrix a, int v)
325//{
326//  int i,j;
327//  poly p;
328//
329//  for (i=1;i<=MATROWS(a);i++)
330//  {
331//    for (j=1;j<=MATCOLS(a);j++)
332//    {
333//      p=pHomogen(MATELEM(a,i,j),v);
334//      p_Delete(&(MATELEM(a,i,j)), ?);
335//      MATELEM(a,i,j)=p;
336//    }
337//  }
338//  return a;
339//}
340
341
342
343/* --------------- internal stuff ------------------- */
344
345row_col_weight::row_col_weight(int i, int j)
346{
347  ym = i;
348  yn = j;
349  wrow = (float *)omAlloc(i*sizeof(float));
350  wcol = (float *)omAlloc(j*sizeof(float));
351}
352
353row_col_weight::~row_col_weight()
354{
355  if (ym!=0)
356  {
357    omFreeSize((ADDRESS)wcol, yn*sizeof(float));
358    omFreeSize((ADDRESS)wrow, ym*sizeof(float));
359  }
360}
361
362mp_permmatrix::mp_permmatrix(matrix A, const ring r) : sign(1),  R(r)
363{
364  a_m = A->nrows;
365  a_n = A->ncols;
366  this->mpInitMat();
367  Xarray = A->m;
368}
369
370mp_permmatrix::mp_permmatrix(mp_permmatrix *M)
371{
372  poly p, *athis, *aM;
373  int i, j;
374
375  a_m = M->s_m;
376  a_n = M->s_n;
377  sign = M->sign;
378  R = M->R;
379
380  this->mpInitMat();
381  Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
382  for (i=a_m-1; i>=0; i--)
383  {
384    athis = this->mpRowAdr(i);
385    aM = M->mpRowAdr(i);
386    for (j=a_n-1; j>=0; j--)
387    {
388      p = aM[M->qcol[j]];
389      if (p)
390      {
391        athis[j] = p_Copy(p, R);
392      }
393    }
394  }
395}
396
397mp_permmatrix::~mp_permmatrix()
398{
399  int k;
400
401  if (a_m != 0)
402  {
403    omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
404    omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
405    if (Xarray != NULL)
406    {
407      for (k=a_m*a_n-1; k>=0; k--)
408        p_Delete(&Xarray[k], R);
409      omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
410    }
411  }
412}
413
414int mp_permmatrix::mpGetRdim() { return s_m; }
415
416int mp_permmatrix::mpGetCdim() { return s_n; }
417
418int mp_permmatrix::mpGetSign() { return sign; }
419
420void mp_permmatrix::mpSetSearch(int s) { piv_s = s; }
421
422void mp_permmatrix::mpSaveArray() { Xarray = NULL; }
423
424poly mp_permmatrix::mpGetElem(int r, int c)
425{
426  return Xarray[a_n*qrow[r]+qcol[c]];
427}
428
429void mp_permmatrix::mpSetElem(poly p, int r, int c)
430{
431  Xarray[a_n*qrow[r]+qcol[c]] = p;
432}
433
434void mp_permmatrix::mpDelElem(int r, int c)
435{
436  p_Delete(&Xarray[a_n*qrow[r]+qcol[c]], R);
437}
438
439/*
440* the Bareiss-type elimination with division by div (div != NULL)
441*/
442void mp_permmatrix::mpElimBareiss(poly div)
443{
444  poly piv, elim, q1, q2, *ap, *a;
445  int i, j, jj;
446
447  ap = this->mpRowAdr(s_m);
448  piv = ap[qcol[s_n]];
449  for(i=s_m-1; i>=0; i--)
450  {
451    a = this->mpRowAdr(i);
452    elim = a[qcol[s_n]];
453    if (elim != NULL)
454    {
455      elim = p_Neg(elim, R);
456      for (j=s_n-1; j>=0; j--)
457      {
458        q2 = NULL;
459        jj = qcol[j];
460        if (ap[jj] != NULL)
461        {
462          q2 = SM_MULT(ap[jj], elim, div, R);
463          if (a[jj] != NULL)
464          {
465            q1 = SM_MULT(a[jj], piv, div, R);
466            p_Delete(&a[jj], R);
467            q2 = p_Add_q(q2, q1, R);
468          }
469        }
470        else if (a[jj] != NULL)
471        {
472          q2 = SM_MULT(a[jj], piv, div, R);
473        }
474        if ((q2!=NULL) && div)
475          SM_DIV(q2, div, R);
476        a[jj] = q2;
477      }
478      p_Delete(&a[qcol[s_n]], R);
479    }
480    else
481    {
482      for (j=s_n-1; j>=0; j--)
483      {
484        jj = qcol[j];
485        if (a[jj] != NULL)
486        {
487          q2 = SM_MULT(a[jj], piv, div, R);
488          p_Delete(&a[jj], R);
489          if (div)
490            SM_DIV(q2, div, R);
491          a[jj] = q2;
492        }
493      }
494    }
495  }
496}
497
498/*2
499* pivot strategy for Bareiss algorithm
500*/
501int mp_permmatrix::mpPivotBareiss(row_col_weight *C)
502{
503  poly p, *a;
504  int i, j, iopt, jopt;
505  float sum, f1, f2, fo, r, ro, lp;
506  float *dr = C->wrow, *dc = C->wcol;
507
508  fo = 1.0e20;
509  ro = 0.0;
510  iopt = jopt = -1;
511
512  s_n--;
513  s_m--;
514  if (s_m == 0)
515    return 0;
516  if (s_n == 0)
517  {
518    for(i=s_m; i>=0; i--)
519    {
520      p = this->mpRowAdr(i)[qcol[0]];
521      if (p)
522      {
523        f1 = mp_PolyWeight(p, R);
524        if (f1 < fo)
525        {
526          fo = f1;
527          if (iopt >= 0)
528            p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
529          iopt = i;
530        }
531        else
532          p_Delete(&(this->mpRowAdr(i)[qcol[0]]), R);
533      }
534    }
535    if (iopt >= 0)
536      mpReplace(iopt, s_m, sign, qrow);
537    return 0;
538  }
539  this->mpRowWeight(dr);
540  this->mpColWeight(dc);
541  sum = 0.0;
542  for(i=s_m; i>=0; i--)
543    sum += dr[i];
544  for(i=s_m; i>=0; i--)
545  {
546    r = dr[i];
547    a = this->mpRowAdr(i);
548    for(j=s_n; j>=0; j--)
549    {
550      p = a[qcol[j]];
551      if (p)
552      {
553        lp = mp_PolyWeight(p, R);
554        ro = r - lp;
555        f1 = ro * (dc[j]-lp);
556        if (f1 != 0.0)
557        {
558          f2 = lp * (sum - ro - dc[j]);
559          f2 += f1;
560        }
561        else
562          f2 = lp-r-dc[j];
563        if (f2 < fo)
564        {
565          fo = f2;
566          iopt = i;
567          jopt = j;
568        }
569      }
570    }
571  }
572  if (iopt < 0)
573    return 0;
574  mpReplace(iopt, s_m, sign, qrow);
575  mpReplace(jopt, s_n, sign, qcol);
576  return 1;
577}
578
579/*2
580* pivot strategy for Bareiss algorithm with defined row
581*/
582int mp_permmatrix::mpPivotRow(row_col_weight *C, int row)
583{
584  poly p, *a;
585  int j, iopt, jopt;
586  float sum, f1, f2, fo, r, ro, lp;
587  float *dr = C->wrow, *dc = C->wcol;
588
589  fo = 1.0e20;
590  ro = 0.0;
591  iopt = jopt = -1;
592
593  s_n--;
594  s_m--;
595  if (s_m == 0)
596    return 0;
597  if (s_n == 0)
598  {
599    p = this->mpRowAdr(row)[qcol[0]];
600    if (p)
601    {
602      f1 = mp_PolyWeight(p, R);
603      if (f1 < fo)
604      {
605        fo = f1;
606        if (iopt >= 0)
607          p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
608        iopt = row;
609      }
610      else
611        p_Delete(&(this->mpRowAdr(row)[qcol[0]]), R);
612    }
613    if (iopt >= 0)
614      mpReplace(iopt, s_m, sign, qrow);
615    return 0;
616  }
617  this->mpRowWeight(dr);
618  this->mpColWeight(dc);
619  sum = 0.0;
620  for(j=s_m; j>=0; j--)
621    sum += dr[j];
622  r = dr[row];
623  a = this->mpRowAdr(row);
624  for(j=s_n; j>=0; j--)
625  {
626    p = a[qcol[j]];
627    if (p)
628    {
629      lp = mp_PolyWeight(p, R);
630      ro = r - lp;
631      f1 = ro * (dc[j]-lp);
632      if (f1 != 0.0)
633      {
634        f2 = lp * (sum - ro - dc[j]);
635        f2 += f1;
636      }
637      else
638        f2 = lp-r-dc[j];
639      if (f2 < fo)
640      {
641        fo = f2;
642        iopt = row;
643        jopt = j;
644      }
645    }
646  }
647  if (iopt < 0)
648    return 0;
649  mpReplace(iopt, s_m, sign, qrow);
650  mpReplace(jopt, s_n, sign, qcol);
651  return 1;
652}
653
654void mp_permmatrix::mpToIntvec(intvec *v)
655{
656  int i;
657
658  for (i=v->rows()-1; i>=0; i--)
659    (*v)[i] = qcol[i]+1;
660}
661
662void mp_permmatrix::mpRowReorder()
663{
664  int k, i, i1, i2;
665
666  if (a_m > a_n)
667    k = a_m - a_n;
668  else
669    k = 0;
670  for (i=a_m-1; i>=k; i--)
671  {
672    i1 = qrow[i];
673    if (i1 != i)
674    {
675      this->mpRowSwap(i1, i);
676      i2 = 0;
677      while (qrow[i2] != i) i2++;
678      qrow[i2] = i1;
679    }
680  }
681}
682
683void mp_permmatrix::mpColReorder()
684{
685  int k, j, j1, j2;
686
687  if (a_n > a_m)
688    k = a_n - a_m;
689  else
690    k = 0;
691  for (j=a_n-1; j>=k; j--)
692  {
693    j1 = qcol[j];
694    if (j1 != j)
695    {
696      this->mpColSwap(j1, j);
697      j2 = 0;
698      while (qcol[j2] != j) j2++;
699      qcol[j2] = j1;
700    }
701  }
702}
703
704// private
705void mp_permmatrix::mpInitMat()
706{
707  int k;
708
709  s_m = a_m;
710  s_n = a_n;
711  piv_s = 0;
712  qrow = (int *)omAlloc(a_m*sizeof(int));
713  qcol = (int *)omAlloc(a_n*sizeof(int));
714  for (k=a_m-1; k>=0; k--) qrow[k] = k;
715  for (k=a_n-1; k>=0; k--) qcol[k] = k;
716}
717
718poly * mp_permmatrix::mpRowAdr(int r)
719{
720  return &(Xarray[a_n*qrow[r]]);
721}
722
723poly * mp_permmatrix::mpColAdr(int c)
724{
725  return &(Xarray[qcol[c]]);
726}
727
728void mp_permmatrix::mpRowWeight(float *wrow)
729{
730  poly p, *a;
731  int i, j;
732  float count;
733
734  for (i=s_m; i>=0; i--)
735  {
736    a = this->mpRowAdr(i);
737    count = 0.0;
738    for(j=s_n; j>=0; j--)
739    {
740      p = a[qcol[j]];
741      if (p)
742        count += mp_PolyWeight(p, R);
743    }
744    wrow[i] = count;
745  }
746}
747
748void mp_permmatrix::mpColWeight(float *wcol)
749{
750  poly p, *a;
751  int i, j;
752  float count;
753
754  for (j=s_n; j>=0; j--)
755  {
756    a = this->mpColAdr(j);
757    count = 0.0;
758    for(i=s_m; i>=0; i--)
759    {
760      p = a[a_n*qrow[i]];
761      if (p)
762        count += mp_PolyWeight(p, R);
763    }
764    wcol[j] = count;
765  }
766}
767
768void mp_permmatrix::mpRowSwap(int i1, int i2)
769{
770  poly p, *a1, *a2;
771  int j;
772
773  a1 = &(Xarray[a_n*i1]);
774  a2 = &(Xarray[a_n*i2]);
775  for (j=a_n-1; j>= 0; j--)
776  {
777    p = a1[j];
778    a1[j] = a2[j];
779    a2[j] = p;
780  }
781}
782
783void mp_permmatrix::mpColSwap(int j1, int j2)
784{
785  poly p, *a1, *a2;
786  int i, k = a_n*a_m;
787
788  a1 = &(Xarray[j1]);
789  a2 = &(Xarray[j2]);
790  for (i=0; i< k; i+=a_n)
791  {
792    p = a1[i];
793    a1[i] = a2[i];
794    a2[i] = p;
795  }
796}
797
798int mp_permmatrix::mpGetRow()
799{
800  return qrow[s_m];
801}
802
803int mp_permmatrix::mpGetCol()
804{
805  return qcol[s_n];
806}
807
808/// perform replacement for pivot strategy in Bareiss algorithm
809/// change sign of determinant
810static void mpReplace(int j, int n, int &sign, int *perm)
811{
812  int k;
813
814  if (j != n)
815  {
816    k = perm[n];
817    perm[n] = perm[j];
818    perm[j] = k;
819    sign = -sign;
820  }
821}
822
823static int mpNextperm(perm * z, int max)
824{
825  int s, i, k, t;
826  s = max;
827  do
828  {
829    s--;
830  }
831  while ((s > 0) && ((*z)[s] >= (*z)[s+1]));
832  if (s==0)
833    return 0;
834  do
835  {
836    (*z)[s]++;
837    k = 0;
838    do
839    {
840      k++;
841    }
842    while (((*z)[k] != (*z)[s]) && (k!=s));
843  }
844  while (k < s);
845  for (i=s+1; i <= max; i++)
846  {
847    (*z)[i]=0;
848    do
849    {
850      (*z)[i]++;
851      k=0;
852      do
853      {
854        k++;
855      }
856      while (((*z)[k] != (*z)[i]) && (k != i));
857    }
858    while (k < i);
859  }
860  s = max+1;
861  do
862  {
863    s--;
864  }
865  while ((s > 0) && ((*z)[s] > (*z)[s+1]));
866  t = 1;
867  for (i=1; i<max; i++)
868    for (k=i+1; k<=max; k++)
869      if ((*z)[k] < (*z)[i])
870        t = -t;
871  (*z)[0] = t;
872  return s;
873}
874
875static poly mp_Leibnitz(matrix a, const ring R)
876{
877  int i, e, n;
878  poly p, d;
879  perm z;
880
881  n = MATROWS(a);
882  memset(&z,0,(n+2)*sizeof(int));
883  p = p_One(R);
884  for (i=1; i <= n; i++)
885    p = p_Mult_q(p, p_Copy(MATELEM(a, i, i), R), R);
886  d = p;
887  for (i=1; i<= n; i++)
888    z[i] = i;
889  z[0]=1;
890  e = 1;
891  if (n!=1)
892  {
893    while (e)
894    {
895      e = mpNextperm((perm *)&z, n);
896      p = p_One(R);
897      for (i = 1; i <= n; i++)
898        p = p_Mult_q(p, p_Copy(MATELEM(a, i, z[i]), R), R);
899      if (z[0] > 0)
900        d = p_Add_q(d, p, R);
901      else
902        d = p_Sub(d, p, R);
903    }
904  }
905  return d;
906}
907
908static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
909{
910  int r=lr-1, c=lc-1;
911  poly *b = a0->m, *x = re->m;
912  poly piv, elim, q1, q2, *ap, *a, *q;
913  int i, j;
914
915  ap = &b[r*a0->ncols];
916  piv = ap[c];
917  for(j=c-1; j>=0; j--)
918    if (ap[j] != NULL) ap[j] = p_Neg(ap[j], R);
919  for(i=r-1; i>=0; i--)
920  {
921    a = &b[i*a0->ncols];
922    q = &x[i*re->ncols];
923    if (a[c] != NULL)
924    {
925      elim = a[c];
926      for (j=c-1; j>=0; j--)
927      {
928        q1 = NULL;
929        if (a[j] != NULL)
930        {
931          q1 = SM_MULT(a[j], piv, div, R);
932          if (ap[j] != NULL)
933          {
934            q2 = SM_MULT(ap[j], elim, div, R);
935            q1 = p_Add_q(q1,q2, R);
936          }
937        }
938        else if (ap[j] != NULL)
939          q1 = SM_MULT(ap[j], elim, div, R);
940        if (q1 != NULL)
941        {
942          if (div)
943            SM_DIV(q1, div, R);
944          q[j] = q1;
945        }
946      }
947    }
948    else
949    {
950      for (j=c-1; j>=0; j--)
951      {
952        if (a[j] != NULL)
953        {
954          q1 = SM_MULT(a[j], piv, div, R);
955          if (div)
956            SM_DIV(q1, div, R);
957          q[j] = q1;
958        }
959      }
960    }
961  }
962}
963
Note: See TracBrowser for help on using the repository browser.