source: git/libpolys/polys/linalg_from_matpol.cc @ 18dab28

spielwiese
Last change on this file since 18dab28 was 845729b, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
Removed linear algebra stuff (Bareiss/Det etc) from matpol.cc TODO: matpol -> simplematrices?
  • Property mode set to 100644
File size: 21.3 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/*
172/// entries of a are minors and go to result (only if not in R)
173void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
174                     ideal R, const ring R)
175{
176  poly *q1;
177  int e=IDELEMS(result);
178  int i,j;
179
180  if (R != NULL)
181  {
182    for (i=r-1;i>=0;i--)
183    {
184      q1 = &(a->m)[i*a->ncols];
185      for (j=c-1;j>=0;j--)
186      {
187        if (q1[j]!=NULL) q1[j] = kNF(R,currQuotient,q1[j]);
188      }
189    }
190  }
191  for (i=r-1;i>=0;i--)
192  {
193    q1 = &(a->m)[i*a->ncols];
194    for (j=c-1;j>=0;j--)
195    {
196      if (q1[j]!=NULL)
197      {
198        if (elems>=e)
199        {
200          if(e<SIZE_OF_SYSTEM_PAGE)
201          {
202            pEnlargeSet(&(result->m),e,e);
203            e += e;
204          }
205          else
206          {
207            pEnlargeSet(&(result->m),e,SIZE_OF_SYSTEM_PAGE);
208            e += SIZE_OF_SYSTEM_PAGE;
209          }
210          IDELEMS(result) =e;
211        }
212        result->m[elems] = q1[j];
213        q1[j] = NULL;
214        elems++;
215      }
216    }
217  }
218}
219
220/// produces recursively the ideal of all arxar-minors of a
221void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
222              poly barDiv, ideal R, const ring R)
223{
224  int k;
225  int kr=lr-1,kc=lc-1;
226  matrix nextLevel=mpNew(kr,kc);
227
228  loop
229  {
230// --- look for an optimal row and bring it to last position ------------
231    if(mpPrepareRow(a,lr,lc)==0) break;
232// --- now take all pivots from the last row ------------
233    k = lc;
234    loop
235    {
236      if(mpPreparePiv(a,lr,k)==0) break;
237      mpElimBar(a,nextLevel,barDiv,lr,k);
238      k--;
239      if (ar>1)
240      {
241        mpRecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R);
242        mpPartClean(nextLevel,kr,k);
243      }
244      else mpMinorToResult(result,elems,nextLevel,kr,k,R);
245      if (ar>k-1) break;
246    }
247    if (ar>=kr) break;
248// --- now we have to take out the last row...------------
249    lr = kr;
250    kr--;
251  }
252  mpFinalClean(nextLevel);
253}
254*/
255
256
257
258/*
259* C++ classes for Bareiss algorithm
260*/
261class row_col_weight
262{
263  private:
264    int ym, yn;
265  public:
266    float *wrow, *wcol;
267    row_col_weight() : ym(0) {}
268    row_col_weight(int, int);
269    ~row_col_weight();
270};
271
272/*2
273*  a submatrix M of a matrix X[m,n]:
274*    0 <= i < s_m <= a_m
275*    0 <= j < s_n <= a_n
276*    M = ( Xarray[qrow[i],qcol[j]] )
277*    if a_m = a_n and s_m = s_n
278*      det(X) = sign*div^(s_m-1)*det(M)
279*    resticted pivot for elimination
280*      0 <= j < piv_s
281*/
282class mp_permmatrix
283{
284  private:
285    int       a_m, a_n, s_m, s_n, sign, piv_s;
286    int       *qrow, *qcol;
287    poly      *Xarray;
288    ring      R;
289
290    void mpInitMat();
291    poly * mpRowAdr(int);
292    poly * mpColAdr(int);
293    void mpRowWeight(float *);
294    void mpColWeight(float *);
295    void mpRowSwap(int, int);
296    void mpColSwap(int, int);
297  public:
298    mp_permmatrix() : a_m(0), R(NULL) {}
299    mp_permmatrix(matrix, const ring);
300    mp_permmatrix(mp_permmatrix *);
301    ~mp_permmatrix();
302    int mpGetRow();
303    int mpGetCol();
304    int mpGetRdim();
305    int mpGetCdim();
306    int mpGetSign();
307    void mpSetSearch(int s);
308    void mpSaveArray();
309    poly mpGetElem(int, int);
310    void mpSetElem(poly, int, int);
311    void mpDelElem(int, int);
312    void mpElimBareiss(poly);
313    int mpPivotBareiss(row_col_weight *);
314    int mpPivotRow(row_col_weight *, int);
315    void mpToIntvec(intvec *);
316    void mpRowReorder();
317    void mpColReorder();
318};
319
320
321
322/*2
323*returns the determinant of the matrix m;
324*uses Bareiss algorithm
325*/
326poly mp_DetBareiss (matrix a, const ring R)
327{
328  int s;
329  poly div, res;
330  if (MATROWS(a) != MATCOLS(a))
331  {
332    Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
333    return NULL;
334  }
335  matrix c = mp_Copy(a, R);
336  mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
337  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
338
339  /* Bareiss */
340  div = NULL;
341  while(Bareiss->mpPivotBareiss(&w))
342  {
343    Bareiss->mpElimBareiss(div);
344    div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
345  }
346  Bareiss->mpRowReorder();
347  Bareiss->mpColReorder();
348  Bareiss->mpSaveArray();
349  s = Bareiss->mpGetSign();
350  delete Bareiss;
351
352  /* result */
353  res = MATELEM(c,1,1);
354  MATELEM(c,1,1) = NULL;
355  id_Delete((ideal *)&c, R);
356  if (s < 0)
357    res = p_Neg(res, R);
358  return res;
359}
360
361
362
363/*2
364*returns the determinant of the matrix m;
365*uses Newtons formulea for symmetric functions
366*/
367poly mp_Det (matrix m, const ring R)
368{
369  int i,j,k,n;
370  poly p,q;
371  matrix a, s;
372  matrix ma[100];
373  number c=NULL, d=NULL, ONE=NULL;
374
375  n = MATROWS(m);
376  if (n != MATCOLS(m))
377  {
378    Werror("det of %d x %d matrix",n,MATCOLS(m));
379    return NULL;
380  }
381  k=rChar(R);
382  if ((k > 0) && (k <= n))
383    return mp_Leibnitz(m, R);
384  ONE = n_Init(1, R);
385  ma[1]=mp_Copy(m, R);
386  k = (n+1) / 2;
387  s = mpNew(1, n);
388  MATELEM(s,1,1) = mp_Trace(m, R);
389  for (i=2; i<=k; i++)
390  {
391    //ma[i] = mpNew(n,n);
392    ma[i]=mp_Mult(ma[i-1], ma[1], R);
393    MATELEM(s,1,i) = mp_Trace(ma[i], R);
394    p_Test(MATELEM(s,1,i), R);
395  }
396  for (i=k+1; i<=n; i++)
397  {
398    MATELEM(s,1,i) = TraceOfProd(ma[i / 2], ma[(i+1) / 2], n, R);
399    p_Test(MATELEM(s,1,i), R);
400  }
401  for (i=1; i<=k; i++)
402    id_Delete((ideal *)&(ma[i]), R);
403/* the array s contains the traces of the powers of the matrix m,
404*  these are the power sums of the eigenvalues of m */
405  a = mpNew(1,n);
406  MATELEM(a,1,1) = minuscopy(MATELEM(s,1,1), R);
407  for (i=2; i<=n; i++)
408  {
409    p = p_Copy(MATELEM(s,1,i), R);
410    for (j=i-1; j>=1; j--)
411    {
412      q = pp_Mult_qq(MATELEM(s,1,j), MATELEM(a,1,i-j), R);
413      p_Test(q, R);
414      p = p_Add_q(p,q, R);
415    }
416    // c= -1/i
417    d = n_Init(-(int)i, R);
418    c = n_Div(ONE, d, R);
419    n_Delete(&d, R);
420
421    p_Mult_nn(p, c, R);
422    p_Test(p, R);
423    MATELEM(a,1,i) = p;
424    n_Delete(&c, R);
425  }
426/* the array a contains the elementary symmetric functions of the
427*  eigenvalues of m */
428  for (i=1; i<=n-1; i++)
429  {
430    //p_Delete(&(MATELEM(a,1,i)), R);
431    p_Delete(&(MATELEM(s,1,i)), R);
432  }
433  p_Delete(&(MATELEM(s,1,n)), R);
434/* up to a sign, the determinant is the n-th elementary symmetric function */
435  if ((n/2)*2 < n)
436  {
437    d = n_Init(-1, R);
438    p_Mult_nn(MATELEM(a,1,n), d, R);
439    n_Delete(&d, R);
440  }
441  n_Delete(&ONE, R);
442  id_Delete((ideal *)&s, R);
443  poly result=MATELEM(a,1,n);
444  MATELEM(a,1,n)=NULL;
445  id_Delete((ideal *)&a, R);
446  return result;
447}
448
449/*2
450* compute all ar-minors of the matrix a
451*/
452matrix mp_Wedge(matrix a, int ar, const ring R)
453{
454  int     i,j,k,l;
455  int *rowchoise,*colchoise;
456  BOOLEAN rowch,colch;
457  matrix result;
458  matrix tmp;
459  poly p;
460
461  i = binom(a->nrows,ar);
462  j = binom(a->ncols,ar);
463
464  rowchoise=(int *)omAlloc(ar*sizeof(int));
465  colchoise=(int *)omAlloc(ar*sizeof(int));
466  result =mpNew(i,j);
467  tmp=mpNew(ar,ar);
468  l = 1; /* k,l:the index in result*/
469  idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
470  while (!rowch)
471  {
472    k=1;
473    idInitChoise(ar,1,a->ncols,&colch,colchoise);
474    while (!colch)
475    {
476      for (i=1; i<=ar; i++)
477      {
478        for (j=1; j<=ar; j++)
479        {
480          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
481        }
482      }
483      p = mp_DetBareiss(tmp, R);
484      if ((k+l) & 1) p=p_Neg(p, R);
485      MATELEM(result,l,k) = p;
486      k++;
487      idGetNextChoise(ar,a->ncols,&colch,colchoise);
488    }
489    idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
490    l++;
491  }
492
493  /*delete the matrix tmp*/
494  for (i=1; i<=ar; i++)
495  {
496    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
497  }
498  id_Delete((ideal *) &tmp, R);
499  return (result);
500}
501
502///*2
503//*homogenize all elements of matrix (not the matrix itself)
504//*/
505//matrix mpHomogen(matrix a, int v)
506//{
507//  int i,j;
508//  poly p;
509//
510//  for (i=1;i<=MATROWS(a);i++)
511//  {
512//    for (j=1;j<=MATCOLS(a);j++)
513//    {
514//      p=pHomogen(MATELEM(a,i,j),v);
515//      p_Delete(&(MATELEM(a,i,j)), ?);
516//      MATELEM(a,i,j)=p;
517//    }
518//  }
519//  return a;
520//}
521
522
523
524/* --------------- internal stuff ------------------- */
525
526row_col_weight::row_col_weight(int i, int j)
527{
528  ym = i;
529  yn = j;
530  wrow = (float *)omAlloc(i*sizeof(float));
531  wcol = (float *)omAlloc(j*sizeof(float));
532}
533
534row_col_weight::~row_col_weight()
535{
536  if (ym!=0)
537  {
538    omFreeSize((ADDRESS)wcol, yn*sizeof(float));
539    omFreeSize((ADDRESS)wrow, ym*sizeof(float));
540  }
541}
542
543mp_permmatrix::mp_permmatrix(matrix A, const ring r) : sign(1),  R(r)
544{
545  a_m = A->nrows;
546  a_n = A->ncols;
547  this->mpInitMat();
548  Xarray = A->m;
549}
550
551mp_permmatrix::mp_permmatrix(mp_permmatrix *M)
552{
553  poly p, *athis, *aM;
554  int i, j;
555
556  a_m = M->s_m;
557  a_n = M->s_n;
558  sign = M->sign;
559  R = M->R;
560
561  this->mpInitMat();
562  Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
563  for (i=a_m-1; i>=0; i--)
564  {
565    athis = this->mpRowAdr(i);
566    aM = M->mpRowAdr(i);
567    for (j=a_n-1; j>=0; j--)
568    {
569      p = aM[M->qcol[j]];
570      if (p)
571      {
572        athis[j] = p_Copy(p, R);
573      }
574    }
575  }
576}
577
578mp_permmatrix::~mp_permmatrix()
579{
580  int k;
581
582  if (a_m != 0)
583  {
584    omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
585    omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
586    if (Xarray != NULL)
587    {
588      for (k=a_m*a_n-1; k>=0; k--)
589        p_Delete(&Xarray[k], R);
590      omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
591    }
592  }
593}
594
595int mp_permmatrix::mpGetRdim() { return s_m; }
596
597int mp_permmatrix::mpGetCdim() { return s_n; }
598
599int mp_permmatrix::mpGetSign() { return sign; }
600
601void mp_permmatrix::mpSetSearch(int s) { piv_s = s; }
602
603void mp_permmatrix::mpSaveArray() { Xarray = NULL; }
604
605poly mp_permmatrix::mpGetElem(int r, int c)
606{
607  return Xarray[a_n*qrow[r]+qcol[c]];
608}
609
610void mp_permmatrix::mpSetElem(poly p, int r, int c)
611{
612  Xarray[a_n*qrow[r]+qcol[c]] = p;
613}
614
615void mp_permmatrix::mpDelElem(int r, int c)
616{
617  p_Delete(&Xarray[a_n*qrow[r]+qcol[c]], R);
618}
619
620/*
621* the Bareiss-type elimination with division by div (div != NULL)
622*/
623void mp_permmatrix::mpElimBareiss(poly div)
624{
625  poly piv, elim, q1, q2, *ap, *a;
626  int i, j, jj;
627
628  ap = this->mpRowAdr(s_m);
629  piv = ap[qcol[s_n]];
630  for(i=s_m-1; i>=0; i--)
631  {
632    a = this->mpRowAdr(i);
633    elim = a[qcol[s_n]];
634    if (elim != NULL)
635    {
636      elim = p_Neg(elim, R);
637      for (j=s_n-1; j>=0; j--)
638      {
639        q2 = NULL;
640        jj = qcol[j];
641        if (ap[jj] != NULL)
642        {
643          q2 = SM_MULT(ap[jj], elim, div, R);
644          if (a[jj] != NULL)
645          {
646            q1 = SM_MULT(a[jj], piv, div, R);
647            p_Delete(&a[jj], R);
648            q2 = p_Add_q(q2, q1, R);
649          }
650        }
651        else if (a[jj] != NULL)
652        {
653          q2 = SM_MULT(a[jj], piv, div, R);
654        }
655        if ((q2!=NULL) && div)
656          SM_DIV(q2, div, R);
657        a[jj] = q2;
658      }
659      p_Delete(&a[qcol[s_n]], R);
660    }
661    else
662    {
663      for (j=s_n-1; j>=0; j--)
664      {
665        jj = qcol[j];
666        if (a[jj] != NULL)
667        {
668          q2 = SM_MULT(a[jj], piv, div, R);
669          p_Delete(&a[jj], R);
670          if (div)
671            SM_DIV(q2, div, R);
672          a[jj] = q2;
673        }
674      }
675    }
676  }
677}
678
679/*2
680* pivot strategy for Bareiss algorithm
681*/
682int mp_permmatrix::mpPivotBareiss(row_col_weight *C)
683{
684  poly p, *a;
685  int i, j, iopt, jopt;
686  float sum, f1, f2, fo, r, ro, lp;
687  float *dr = C->wrow, *dc = C->wcol;
688
689  fo = 1.0e20;
690  ro = 0.0;
691  iopt = jopt = -1;
692
693  s_n--;
694  s_m--;
695  if (s_m == 0)
696    return 0;
697  if (s_n == 0)
698  {
699    for(i=s_m; i>=0; i--)
700    {
701      p = this->mpRowAdr(i)[qcol[0]];
702      if (p)
703      {
704        f1 = mp_PolyWeight(p, R);
705        if (f1 < fo)
706        {
707          fo = f1;
708          if (iopt >= 0)
709            p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
710          iopt = i;
711        }
712        else
713          p_Delete(&(this->mpRowAdr(i)[qcol[0]]), R);
714      }
715    }
716    if (iopt >= 0)
717      mpReplace(iopt, s_m, sign, qrow);
718    return 0;
719  }
720  this->mpRowWeight(dr);
721  this->mpColWeight(dc);
722  sum = 0.0;
723  for(i=s_m; i>=0; i--)
724    sum += dr[i];
725  for(i=s_m; i>=0; i--)
726  {
727    r = dr[i];
728    a = this->mpRowAdr(i);
729    for(j=s_n; j>=0; j--)
730    {
731      p = a[qcol[j]];
732      if (p)
733      {
734        lp = mp_PolyWeight(p, R);
735        ro = r - lp;
736        f1 = ro * (dc[j]-lp);
737        if (f1 != 0.0)
738        {
739          f2 = lp * (sum - ro - dc[j]);
740          f2 += f1;
741        }
742        else
743          f2 = lp-r-dc[j];
744        if (f2 < fo)
745        {
746          fo = f2;
747          iopt = i;
748          jopt = j;
749        }
750      }
751    }
752  }
753  if (iopt < 0)
754    return 0;
755  mpReplace(iopt, s_m, sign, qrow);
756  mpReplace(jopt, s_n, sign, qcol);
757  return 1;
758}
759
760/*2
761* pivot strategy for Bareiss algorithm with defined row
762*/
763int mp_permmatrix::mpPivotRow(row_col_weight *C, int row)
764{
765  poly p, *a;
766  int j, iopt, jopt;
767  float sum, f1, f2, fo, r, ro, lp;
768  float *dr = C->wrow, *dc = C->wcol;
769
770  fo = 1.0e20;
771  ro = 0.0;
772  iopt = jopt = -1;
773
774  s_n--;
775  s_m--;
776  if (s_m == 0)
777    return 0;
778  if (s_n == 0)
779  {
780    p = this->mpRowAdr(row)[qcol[0]];
781    if (p)
782    {
783      f1 = mp_PolyWeight(p, R);
784      if (f1 < fo)
785      {
786        fo = f1;
787        if (iopt >= 0)
788          p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
789        iopt = row;
790      }
791      else
792        p_Delete(&(this->mpRowAdr(row)[qcol[0]]), R);
793    }
794    if (iopt >= 0)
795      mpReplace(iopt, s_m, sign, qrow);
796    return 0;
797  }
798  this->mpRowWeight(dr);
799  this->mpColWeight(dc);
800  sum = 0.0;
801  for(j=s_m; j>=0; j--)
802    sum += dr[j];
803  r = dr[row];
804  a = this->mpRowAdr(row);
805  for(j=s_n; j>=0; j--)
806  {
807    p = a[qcol[j]];
808    if (p)
809    {
810      lp = mp_PolyWeight(p, R);
811      ro = r - lp;
812      f1 = ro * (dc[j]-lp);
813      if (f1 != 0.0)
814      {
815        f2 = lp * (sum - ro - dc[j]);
816        f2 += f1;
817      }
818      else
819        f2 = lp-r-dc[j];
820      if (f2 < fo)
821      {
822        fo = f2;
823        iopt = row;
824        jopt = j;
825      }
826    }
827  }
828  if (iopt < 0)
829    return 0;
830  mpReplace(iopt, s_m, sign, qrow);
831  mpReplace(jopt, s_n, sign, qcol);
832  return 1;
833}
834
835void mp_permmatrix::mpToIntvec(intvec *v)
836{
837  int i;
838
839  for (i=v->rows()-1; i>=0; i--)
840    (*v)[i] = qcol[i]+1;
841}
842
843void mp_permmatrix::mpRowReorder()
844{
845  int k, i, i1, i2;
846
847  if (a_m > a_n)
848    k = a_m - a_n;
849  else
850    k = 0;
851  for (i=a_m-1; i>=k; i--)
852  {
853    i1 = qrow[i];
854    if (i1 != i)
855    {
856      this->mpRowSwap(i1, i);
857      i2 = 0;
858      while (qrow[i2] != i) i2++;
859      qrow[i2] = i1;
860    }
861  }
862}
863
864void mp_permmatrix::mpColReorder()
865{
866  int k, j, j1, j2;
867
868  if (a_n > a_m)
869    k = a_n - a_m;
870  else
871    k = 0;
872  for (j=a_n-1; j>=k; j--)
873  {
874    j1 = qcol[j];
875    if (j1 != j)
876    {
877      this->mpColSwap(j1, j);
878      j2 = 0;
879      while (qcol[j2] != j) j2++;
880      qcol[j2] = j1;
881    }
882  }
883}
884
885// private
886void mp_permmatrix::mpInitMat()
887{
888  int k;
889
890  s_m = a_m;
891  s_n = a_n;
892  piv_s = 0;
893  qrow = (int *)omAlloc(a_m*sizeof(int));
894  qcol = (int *)omAlloc(a_n*sizeof(int));
895  for (k=a_m-1; k>=0; k--) qrow[k] = k;
896  for (k=a_n-1; k>=0; k--) qcol[k] = k;
897}
898
899poly * mp_permmatrix::mpRowAdr(int r)
900{
901  return &(Xarray[a_n*qrow[r]]);
902}
903
904poly * mp_permmatrix::mpColAdr(int c)
905{
906  return &(Xarray[qcol[c]]);
907}
908
909void mp_permmatrix::mpRowWeight(float *wrow)
910{
911  poly p, *a;
912  int i, j;
913  float count;
914
915  for (i=s_m; i>=0; i--)
916  {
917    a = this->mpRowAdr(i);
918    count = 0.0;
919    for(j=s_n; j>=0; j--)
920    {
921      p = a[qcol[j]];
922      if (p)
923        count += mp_PolyWeight(p, R);
924    }
925    wrow[i] = count;
926  }
927}
928
929void mp_permmatrix::mpColWeight(float *wcol)
930{
931  poly p, *a;
932  int i, j;
933  float count;
934
935  for (j=s_n; j>=0; j--)
936  {
937    a = this->mpColAdr(j);
938    count = 0.0;
939    for(i=s_m; i>=0; i--)
940    {
941      p = a[a_n*qrow[i]];
942      if (p)
943        count += mp_PolyWeight(p, R);
944    }
945    wcol[j] = count;
946  }
947}
948
949void mp_permmatrix::mpRowSwap(int i1, int i2)
950{
951  poly p, *a1, *a2;
952  int j;
953
954  a1 = &(Xarray[a_n*i1]);
955  a2 = &(Xarray[a_n*i2]);
956  for (j=a_n-1; j>= 0; j--)
957  {
958    p = a1[j];
959    a1[j] = a2[j];
960    a2[j] = p;
961  }
962}
963
964void mp_permmatrix::mpColSwap(int j1, int j2)
965{
966  poly p, *a1, *a2;
967  int i, k = a_n*a_m;
968
969  a1 = &(Xarray[j1]);
970  a2 = &(Xarray[j2]);
971  for (i=0; i< k; i+=a_n)
972  {
973    p = a1[i];
974    a1[i] = a2[i];
975    a2[i] = p;
976  }
977}
978
979int mp_permmatrix::mpGetRow()
980{
981  return qrow[s_m];
982}
983
984int mp_permmatrix::mpGetCol()
985{
986  return qcol[s_n];
987}
988
989/// perform replacement for pivot strategy in Bareiss algorithm
990/// change sign of determinant
991static void mpReplace(int j, int n, int &sign, int *perm)
992{
993  int k;
994
995  if (j != n)
996  {
997    k = perm[n];
998    perm[n] = perm[j];
999    perm[j] = k;
1000    sign = -sign;
1001  }
1002}
1003
1004static int mpNextperm(perm * z, int max)
1005{
1006  int s, i, k, t;
1007  s = max;
1008  do
1009  {
1010    s--;
1011  }
1012  while ((s > 0) && ((*z)[s] >= (*z)[s+1]));
1013  if (s==0)
1014    return 0;
1015  do
1016  {
1017    (*z)[s]++;
1018    k = 0;
1019    do
1020    {
1021      k++;
1022    }
1023    while (((*z)[k] != (*z)[s]) && (k!=s));
1024  }
1025  while (k < s);
1026  for (i=s+1; i <= max; i++)
1027  {
1028    (*z)[i]=0;
1029    do
1030    {
1031      (*z)[i]++;
1032      k=0;
1033      do
1034      {
1035        k++;
1036      }
1037      while (((*z)[k] != (*z)[i]) && (k != i));
1038    }
1039    while (k < i);
1040  }
1041  s = max+1;
1042  do
1043  {
1044    s--;
1045  }
1046  while ((s > 0) && ((*z)[s] > (*z)[s+1]));
1047  t = 1;
1048  for (i=1; i<max; i++)
1049    for (k=i+1; k<=max; k++)
1050      if ((*z)[k] < (*z)[i])
1051        t = -t;
1052  (*z)[0] = t;
1053  return s;
1054}
1055
1056static poly mp_Leibnitz(matrix a, const ring R)
1057{
1058  int i, e, n;
1059  poly p, d;
1060  perm z;
1061
1062  n = MATROWS(a);
1063  memset(&z,0,(n+2)*sizeof(int));
1064  p = p_One(R);
1065  for (i=1; i <= n; i++)
1066    p = p_Mult_q(p, p_Copy(MATELEM(a, i, i), R), R);
1067  d = p;
1068  for (i=1; i<= n; i++)
1069    z[i] = i;
1070  z[0]=1;
1071  e = 1;
1072  if (n!=1)
1073  {
1074    while (e)
1075    {
1076      e = mpNextperm((perm *)&z, n);
1077      p = p_One(R);
1078      for (i = 1; i <= n; i++)
1079        p = p_Mult_q(p, p_Copy(MATELEM(a, i, z[i]), R), R);
1080      if (z[0] > 0)
1081        d = p_Add_q(d, p, R);
1082      else
1083        d = p_Sub(d, p, R);
1084    }
1085  }
1086  return d;
1087}
1088
1089static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
1090{
1091  int r=lr-1, c=lc-1;
1092  poly *b = a0->m, *x = re->m;
1093  poly piv, elim, q1, q2, *ap, *a, *q;
1094  int i, j;
1095
1096  ap = &b[r*a0->ncols];
1097  piv = ap[c];
1098  for(j=c-1; j>=0; j--)
1099    if (ap[j] != NULL) ap[j] = p_Neg(ap[j], R);
1100  for(i=r-1; i>=0; i--)
1101  {
1102    a = &b[i*a0->ncols];
1103    q = &x[i*re->ncols];
1104    if (a[c] != NULL)
1105    {
1106      elim = a[c];
1107      for (j=c-1; j>=0; j--)
1108      {
1109        q1 = NULL;
1110        if (a[j] != NULL)
1111        {
1112          q1 = SM_MULT(a[j], piv, div, R);
1113          if (ap[j] != NULL)
1114          {
1115            q2 = SM_MULT(ap[j], elim, div, R);
1116            q1 = p_Add_q(q1,q2, R);
1117          }
1118        }
1119        else if (ap[j] != NULL)
1120          q1 = SM_MULT(ap[j], elim, div, R);
1121        if (q1 != NULL)
1122        {
1123          if (div)
1124            SM_DIV(q1, div, R);
1125          q[j] = q1;
1126        }
1127      }
1128    }
1129    else
1130    {
1131      for (j=c-1; j>=0; j--)
1132      {
1133        if (a[j] != NULL)
1134        {
1135          q1 = SM_MULT(a[j], piv, div, R);
1136          if (div)
1137            SM_DIV(q1, div, R);
1138          q[j] = q1;
1139        }
1140      }
1141    }
1142  }
1143}
1144
Note: See TracBrowser for help on using the repository browser.