Changeset 845729b in git for libpolys/polys/matpol.cc


Ignore:
Timestamp:
Apr 18, 2011, 8:44:58 PM (13 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
Children:
4f684a7f137cff04239a5a3fa6eb3dafd2cc5084
Parents:
1820e70c46d10bc89d249ced6bc28edfccc0922e
git-author:
Oleksandr Motsak <motsak@mathematik.uni-kl.de>2011-04-18 20:44:58+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:31:17+01:00
Message:
Removed linear algebra stuff (Bareiss/Det etc) from matpol.cc
TODO: matpol -> simplematrices?
File:
1 edited

Legend:

Unmodified
Added
Removed

  • libpolys/polys/matpol.cc

    r1820e7 r845729b  
    3636#include "prCopy.h"
    3737
    38 #include "sparsmat.h"
     38// #include "sparsmat.h"
    3939   
    4040//omBin ip_smatrix_bin = omGetSpecBin(sizeof(ip_smatrix));
    … …  
    4242/*0 implementation*/
    4343
    44 
    45 typedef int perm[100];
    46 static void mpReplace(int j, int n, int &sign, int *perm);
    47 static int mpNextperm(perm * z, int max);
    48 static poly mp_Leibnitz(matrix a, const ring);
    49 static poly minuscopy (poly p, const ring);
    50 static poly p_Insert(poly p1, poly p2, const ring);
    5144static poly mp_Exdiv ( poly m, poly d, poly vars, const ring);
    5245static poly mp_Select (poly fro, poly what, const ring);
    53 
    54 static void mp_PartClean(matrix, int, int, const ring);
    55 static void mp_FinalClean(matrix, const ring);
    56 static int mp_PrepareRow (matrix, int, int, const ring);
    57 static int mp_PreparePiv (matrix, int, int, const ring);
    58 static int mp_PivBar(matrix, int, int, const ring);
    59 static int mp_PivRow(matrix, int, int, const ring);
    60 static float mp_PolyWeight(poly, const ring);
    61 static void mp_SwapRow(matrix, int, int, int, const ring);
    62 static void mp_SwapCol(matrix, int, int, int, const ring);
    63 static void mp_ElimBar(matrix, matrix, poly, int, int, const ring);
    6446
    6547/// create a r x c zero-matrix
    … …  
    334316}
    335317
    336 /*
    337 * C++ classes for Bareiss algorithm
    338 */
    339 class row_col_weight
    340 {
    341   private:
    342   int ym, yn;
    343   public:
    344   float *wrow, *wcol;
    345   row_col_weight() : ym(0) {}
    346   row_col_weight(int, int);
    347   ~row_col_weight();
    348 };
    349 
    350 /*2
    351 *  a submatrix M of a matrix X[m,n]:
    352 *    0 <= i < s_m <= a_m
    353 *    0 <= j < s_n <= a_n
    354 *    M = ( Xarray[qrow[i],qcol[j]] )
    355 *    if a_m = a_n and s_m = s_n
    356 *      det(X) = sign*div^(s_m-1)*det(M)
    357 *    resticted pivot for elimination
    358 *      0 <= j < piv_s
    359 */
    360 class mp_permmatrix
    361 {
    362   private:
    363   int       a_m, a_n, s_m, s_n, sign, piv_s;
    364   int       *qrow, *qcol;
    365   poly      *Xarray;
    366   ring      R;
    367  
    368   void mpInitMat();
    369   poly * mpRowAdr(int);
    370   poly * mpColAdr(int);
    371   void mpRowWeight(float *);
    372   void mpColWeight(float *);
    373   void mpRowSwap(int, int);
    374   void mpColSwap(int, int);
    375   public:
    376   mp_permmatrix() : a_m(0), R(NULL) {}
    377   mp_permmatrix(matrix, const ring);
    378   mp_permmatrix(mp_permmatrix *);
    379   ~mp_permmatrix();
    380   int mpGetRow();
    381   int mpGetCol();
    382   int mpGetRdim();
    383   int mpGetCdim();
    384   int mpGetSign();
    385   void mpSetSearch(int s);
    386   void mpSaveArray();
    387   poly mpGetElem(int, int);
    388   void mpSetElem(poly, int, int);
    389   void mpDelElem(int, int);
    390   void mpElimBareiss(poly);
    391   int mpPivotBareiss(row_col_weight *);
    392   int mpPivotRow(row_col_weight *, int);
    393   void mpToIntvec(intvec *);
    394   void mpRowReorder();
    395   void mpColReorder();
    396 };
    397 
    398 #ifndef SIZE_OF_SYSTEM_PAGE
    399 #define SIZE_OF_SYSTEM_PAGE 4096
    400 #endif
    401 
    402 /*
    403 /// entries of a are minors and go to result (only if not in R)
    404 void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
    405                      ideal R, const ring R)
    406 {
    407   poly *q1;
    408   int e=IDELEMS(result);
    409   int i,j;
    410 
    411   if (R != NULL)
    412   {
    413     for (i=r-1;i>=0;i--)
    414     {
    415       q1 = &(a->m)[i*a->ncols];
    416       for (j=c-1;j>=0;j--)
    417       {
    418         if (q1[j]!=NULL) q1[j] = kNF(R,currQuotient,q1[j]);
    419       }
    420     }
    421   }
    422   for (i=r-1;i>=0;i--)
    423   {
    424     q1 = &(a->m)[i*a->ncols];
    425     for (j=c-1;j>=0;j--)
    426     {
    427       if (q1[j]!=NULL)
    428       {
    429         if (elems>=e)
    430         {
    431           if(e<SIZE_OF_SYSTEM_PAGE)
    432           {
    433             pEnlargeSet(&(result->m),e,e);
    434             e += e;
    435           }
    436           else
    437           {
    438             pEnlargeSet(&(result->m),e,SIZE_OF_SYSTEM_PAGE);
    439             e += SIZE_OF_SYSTEM_PAGE;
    440           }
    441           IDELEMS(result) =e;
    442         }
    443         result->m[elems] = q1[j];
    444         q1[j] = NULL;
    445         elems++;
    446       }
    447     }
    448   }
    449 }
    450 
    451 /// produces recursively the ideal of all arxar-minors of a
    452 void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
    453               poly barDiv, ideal R, const ring R)
    454 {
    455   int k;
    456   int kr=lr-1,kc=lc-1;
    457   matrix nextLevel=mpNew(kr,kc);
    458 
    459   loop
    460   {
    461 // --- look for an optimal row and bring it to last position ------------
    462     if(mpPrepareRow(a,lr,lc)==0) break;
    463 // --- now take all pivots from the last row ------------
    464     k = lc;
    465     loop
    466     {
    467       if(mpPreparePiv(a,lr,k)==0) break;
    468       mpElimBar(a,nextLevel,barDiv,lr,k);
    469       k--;
    470       if (ar>1)
    471       {
    472         mpRecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R);
    473         mpPartClean(nextLevel,kr,k);
    474       }
    475       else mpMinorToResult(result,elems,nextLevel,kr,k,R);
    476       if (ar>k-1) break;
    477     }
    478     if (ar>=kr) break;
    479 // --- now we have to take out the last row...------------
    480     lr = kr;
    481     kr--;
    482   }
    483   mpFinalClean(nextLevel);
    484 }
    485 */
    486 
    487 
    488 /*2
    489 *returns the determinant of the matrix m;
    490 *uses Bareiss algorithm
    491 */
    492 poly mp_DetBareiss (matrix a, const ring R)
    493 {
    494   int s;
    495   poly div, res;
    496   if (MATROWS(a) != MATCOLS(a))
    497   {
    498     Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
    499     return NULL;
    500   }
    501   matrix c = mp_Copy(a, R);
    502   mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
    503   row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
    504 
    505   /* Bareiss */
    506   div = NULL;
    507   while(Bareiss->mpPivotBareiss(&w))
    508   {
    509     Bareiss->mpElimBareiss(div);
    510     div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
    511   }
    512   Bareiss->mpRowReorder();
    513   Bareiss->mpColReorder();
    514   Bareiss->mpSaveArray();
    515   s = Bareiss->mpGetSign();
    516   delete Bareiss;
    517 
    518   /* result */
    519   res = MATELEM(c,1,1);
    520   MATELEM(c,1,1) = NULL;
    521   id_Delete((ideal *)&c, R);
    522   if (s < 0)
    523     res = p_Neg(res, R);
    524   return res;
    525 }
    526 
    527 /*2
    528 *returns the determinant of the matrix m;
    529 *uses Newtons formulea for symmetric functions
    530 */
    531 poly mp_Det (matrix m, const ring R)
    532 {
    533   int i,j,k,n;
    534   poly p,q;
    535   matrix a, s;
    536   matrix ma[100];
    537   number c=NULL, d=NULL, ONE=NULL;
    538 
    539   n = MATROWS(m);
    540   if (n != MATCOLS(m))
    541   {
    542     Werror("det of %d x %d matrix",n,MATCOLS(m));
    543     return NULL;
    544   }
    545   k=rChar(R);
    546   if ((k > 0) && (k <= n))
    547     return mp_Leibnitz(m, R);
    548   ONE = n_Init(1, R);
    549   ma[1]=mp_Copy(m, R);
    550   k = (n+1) / 2;
    551   s = mpNew(1, n);
    552   MATELEM(s,1,1) = mp_Trace(m, R);
    553   for (i=2; i<=k; i++)
    554   {
    555     //ma[i] = mpNew(n,n);
    556     ma[i]=mp_Mult(ma[i-1], ma[1], R);
    557     MATELEM(s,1,i) = mp_Trace(ma[i], R);
    558     p_Test(MATELEM(s,1,i), R);
    559   }
    560   for (i=k+1; i<=n; i++)
    561   {
    562     MATELEM(s,1,i) = TraceOfProd(ma[i / 2], ma[(i+1) / 2], n, R);
    563     p_Test(MATELEM(s,1,i), R);
    564   }
    565   for (i=1; i<=k; i++)
    566     id_Delete((ideal *)&(ma[i]), R);
    567 /* the array s contains the traces of the powers of the matrix m,
    568 *  these are the power sums of the eigenvalues of m */
    569   a = mpNew(1,n);
    570   MATELEM(a,1,1) = minuscopy(MATELEM(s,1,1), R);
    571   for (i=2; i<=n; i++)
    572   {
    573     p = p_Copy(MATELEM(s,1,i), R);
    574     for (j=i-1; j>=1; j--)
    575     {
    576       q = pp_Mult_qq(MATELEM(s,1,j), MATELEM(a,1,i-j), R);
    577       p_Test(q, R);
    578       p = p_Add_q(p,q, R);
    579     }
    580     // c= -1/i
    581     d = n_Init(-(int)i, R);
    582     c = n_Div(ONE, d, R);
    583     n_Delete(&d, R);
    584 
    585     p_Mult_nn(p, c, R);
    586     p_Test(p, R);
    587     MATELEM(a,1,i) = p;
    588     n_Delete(&c, R);
    589   }
    590 /* the array a contains the elementary symmetric functions of the
    591 *  eigenvalues of m */
    592   for (i=1; i<=n-1; i++)
    593   {
    594     //p_Delete(&(MATELEM(a,1,i)), R);
    595     p_Delete(&(MATELEM(s,1,i)), R);
    596   }
    597   p_Delete(&(MATELEM(s,1,n)), R);
    598 /* up to a sign, the determinant is the n-th elementary symmetric function */
    599   if ((n/2)*2 < n)
    600   {
    601     d = n_Init(-1, R);
    602     p_Mult_nn(MATELEM(a,1,n), d, R);
    603     n_Delete(&d, R);
    604   }
    605   n_Delete(&ONE, R);
    606   id_Delete((ideal *)&s, R);
    607   poly result=MATELEM(a,1,n);
    608   MATELEM(a,1,n)=NULL;
    609   id_Delete((ideal *)&a, R);
    610   return result;
    611 }
    612 
    613 /*2
    614 * compute all ar-minors of the matrix a
    615 */
    616 matrix mp_Wedge(matrix a, int ar, const ring R)
    617 {
    618   int     i,j,k,l;
    619   int *rowchoise,*colchoise;
    620   BOOLEAN rowch,colch;
    621   matrix result;
    622   matrix tmp;
    623   poly p;
    624 
    625   i = binom(a->nrows,ar);
    626   j = binom(a->ncols,ar);
    627 
    628   rowchoise=(int *)omAlloc(ar*sizeof(int));
    629   colchoise=(int *)omAlloc(ar*sizeof(int));
    630   result =mpNew(i,j);
    631   tmp=mpNew(ar,ar);
    632   l = 1; /* k,l:the index in result*/
    633   idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
    634   while (!rowch)
    635   {
    636     k=1;
    637     idInitChoise(ar,1,a->ncols,&colch,colchoise);
    638     while (!colch)
    639     {
    640       for (i=1; i<=ar; i++)
    641       {
    642         for (j=1; j<=ar; j++)
    643         {
    644           MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
    645         }
    646       }
    647       p = mp_DetBareiss(tmp, R);
    648       if ((k+l) & 1) p=p_Neg(p, R);
    649       MATELEM(result,l,k) = p;
    650       k++;
    651       idGetNextChoise(ar,a->ncols,&colch,colchoise);
    652     }
    653     idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
    654     l++;
    655   }
    656 
    657   /*delete the matrix tmp*/
    658   for (i=1; i<=ar; i++)
    659   {
    660     for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
    661   }
    662   id_Delete((ideal *) &tmp, R);
    663   return (result);
    664 }
    665 
    666 ///*2
    667 //*homogenize all elements of matrix (not the matrix itself)
    668 //*/
    669 //matrix mpHomogen(matrix a, int v)
    670 //{
    671 //  int i,j;
    672 //  poly p;
    673 //
    674 //  for (i=1;i<=MATROWS(a);i++)
    675 //  {
    676 //    for (j=1;j<=MATCOLS(a);j++)
    677 //    {
    678 //      p=pHomogen(MATELEM(a,i,j),v);
    679 //      p_Delete(&(MATELEM(a,i,j)), ?);
    680 //      MATELEM(a,i,j)=p;
    681 //    }
    682 //  }
    683 //  return a;
    684 //}
     318// #ifndef SIZE_OF_SYSTEM_PAGE
     319// #define SIZE_OF_SYSTEM_PAGE 4096
     320// #endif
    685321
    686322/*2
    … …  
    974610}
    975611
    976 /* --------------- internal stuff ------------------- */
    977 
    978 row_col_weight::row_col_weight(int i, int j)
    979 {
    980   ym = i;
    981   yn = j;
    982   wrow = (float *)omAlloc(i*sizeof(float));
    983   wcol = (float *)omAlloc(j*sizeof(float));
    984 }
    985 
    986 row_col_weight::~row_col_weight()
    987 {
    988   if (ym!=0)
    989   {
    990     omFreeSize((ADDRESS)wcol, yn*sizeof(float));
    991     omFreeSize((ADDRESS)wrow, ym*sizeof(float));
    992   }
    993 }
    994 
    995 mp_permmatrix::mp_permmatrix(matrix A, const ring r) : sign(1),  R(r)
    996 {
    997   a_m = A->nrows;
    998   a_n = A->ncols;
    999   this->mpInitMat();
    1000   Xarray = A->m;
    1001 }
    1002 
    1003 mp_permmatrix::mp_permmatrix(mp_permmatrix *M)
    1004 {
    1005   poly p, *athis, *aM;
    1006   int i, j;
    1007 
    1008   a_m = M->s_m;
    1009   a_n = M->s_n;
    1010   sign = M->sign;
    1011   R = M->R;
    1012  
    1013   this->mpInitMat();
    1014   Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
    1015   for (i=a_m-1; i>=0; i--)
    1016   {
    1017     athis = this->mpRowAdr(i);
    1018     aM = M->mpRowAdr(i);
    1019     for (j=a_n-1; j>=0; j--)
    1020     {
    1021       p = aM[M->qcol[j]];
    1022       if (p)
    1023       {
    1024         athis[j] = p_Copy(p, R);
    1025       }
    1026     }
    1027   }
    1028 }
    1029 
    1030 mp_permmatrix::~mp_permmatrix()
    1031 {
    1032   int k;
    1033 
    1034   if (a_m != 0)
    1035   {
    1036     omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
    1037     omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
    1038     if (Xarray != NULL)
    1039     {
    1040       for (k=a_m*a_n-1; k>=0; k--)
    1041         p_Delete(&Xarray[k], R);
    1042       omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
    1043     }
    1044   }
    1045 }
    1046 
    1047 int mp_permmatrix::mpGetRdim() { return s_m; }
    1048 
    1049 int mp_permmatrix::mpGetCdim() { return s_n; }
    1050 
    1051 int mp_permmatrix::mpGetSign() { return sign; }
    1052 
    1053 void mp_permmatrix::mpSetSearch(int s) { piv_s = s; }
    1054 
    1055 void mp_permmatrix::mpSaveArray() { Xarray = NULL; }
    1056 
    1057 poly mp_permmatrix::mpGetElem(int r, int c)
    1058 {
    1059   return Xarray[a_n*qrow[r]+qcol[c]];
    1060 }
    1061 
    1062 void mp_permmatrix::mpSetElem(poly p, int r, int c)
    1063 {
    1064   Xarray[a_n*qrow[r]+qcol[c]] = p;
    1065 }
    1066 
    1067 void mp_permmatrix::mpDelElem(int r, int c)
    1068 {
    1069   p_Delete(&Xarray[a_n*qrow[r]+qcol[c]], R);
    1070 }
    1071 
    1072 /*
    1073 * the Bareiss-type elimination with division by div (div != NULL)
    1074 */
    1075 void mp_permmatrix::mpElimBareiss(poly div)
    1076 {
    1077   poly piv, elim, q1, q2, *ap, *a;
    1078   int i, j, jj;
    1079 
    1080   ap = this->mpRowAdr(s_m);
    1081   piv = ap[qcol[s_n]];
    1082   for(i=s_m-1; i>=0; i--)
    1083   {
    1084     a = this->mpRowAdr(i);
    1085     elim = a[qcol[s_n]];
    1086     if (elim != NULL)
    1087     {
    1088       elim = p_Neg(elim, R);
    1089       for (j=s_n-1; j>=0; j--)
    1090       {
    1091         q2 = NULL;
    1092         jj = qcol[j];
    1093         if (ap[jj] != NULL)
    1094         {
    1095           q2 = SM_MULT(ap[jj], elim, div, R);
    1096           if (a[jj] != NULL)
    1097           {
    1098             q1 = SM_MULT(a[jj], piv, div, R);
    1099             p_Delete(&a[jj], R);
    1100             q2 = p_Add_q(q2, q1, R);
    1101           }
    1102         }
    1103         else if (a[jj] != NULL)
    1104         {
    1105           q2 = SM_MULT(a[jj], piv, div, R);
    1106         }
    1107         if ((q2!=NULL) && div)
    1108           SM_DIV(q2, div, R);
    1109         a[jj] = q2;
    1110       }
    1111       p_Delete(&a[qcol[s_n]], R);
    1112     }
    1113     else
    1114     {
    1115       for (j=s_n-1; j>=0; j--)
    1116       {
    1117         jj = qcol[j];
    1118         if (a[jj] != NULL)
    1119         {
    1120           q2 = SM_MULT(a[jj], piv, div, R);
    1121           p_Delete(&a[jj], R);
    1122           if (div)
    1123             SM_DIV(q2, div, R);
    1124           a[jj] = q2;
    1125         }
    1126       }
    1127     }
    1128   }
    1129 }
    1130 
    1131 /*2
    1132 * pivot strategy for Bareiss algorithm
    1133 */
    1134 int mp_permmatrix::mpPivotBareiss(row_col_weight *C)
    1135 {
    1136   poly p, *a;
    1137   int i, j, iopt, jopt;
    1138   float sum, f1, f2, fo, r, ro, lp;
    1139   float *dr = C->wrow, *dc = C->wcol;
    1140 
    1141   fo = 1.0e20;
    1142   ro = 0.0;
    1143   iopt = jopt = -1;
    1144 
    1145   s_n--;
    1146   s_m--;
    1147   if (s_m == 0)
    1148     return 0;
    1149   if (s_n == 0)
    1150   {
    1151     for(i=s_m; i>=0; i--)
    1152     {
    1153       p = this->mpRowAdr(i)[qcol[0]];
    1154       if (p)
    1155       {
    1156         f1 = mp_PolyWeight(p, R);
    1157         if (f1 < fo)
    1158         {
    1159           fo = f1;
    1160           if (iopt >= 0)
    1161             p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
    1162           iopt = i;
    1163         }
    1164         else
    1165           p_Delete(&(this->mpRowAdr(i)[qcol[0]]), R);
    1166       }
    1167     }
    1168     if (iopt >= 0)
    1169       mpReplace(iopt, s_m, sign, qrow);
    1170     return 0;
    1171   }
    1172   this->mpRowWeight(dr);
    1173   this->mpColWeight(dc);
    1174   sum = 0.0;
    1175   for(i=s_m; i>=0; i--)
    1176     sum += dr[i];
    1177   for(i=s_m; i>=0; i--)
    1178   {
    1179     r = dr[i];
    1180     a = this->mpRowAdr(i);
    1181     for(j=s_n; j>=0; j--)
    1182     {
    1183       p = a[qcol[j]];
    1184       if (p)
    1185       {
    1186         lp = mp_PolyWeight(p, R);
    1187         ro = r - lp;
    1188         f1 = ro * (dc[j]-lp);
    1189         if (f1 != 0.0)
    1190         {
    1191           f2 = lp * (sum - ro - dc[j]);
    1192           f2 += f1;
    1193         }
    1194         else
    1195           f2 = lp-r-dc[j];
    1196         if (f2 < fo)
    1197         {
    1198           fo = f2;
    1199           iopt = i;
    1200           jopt = j;
    1201         }
    1202       }
    1203     }
    1204   }
    1205   if (iopt < 0)
    1206     return 0;
    1207   mpReplace(iopt, s_m, sign, qrow);
    1208   mpReplace(jopt, s_n, sign, qcol);
    1209   return 1;
    1210 }
    1211 
    1212 /*2
    1213 * pivot strategy for Bareiss algorithm with defined row
    1214 */
    1215 int mp_permmatrix::mpPivotRow(row_col_weight *C, int row)
    1216 {
    1217   poly p, *a;
    1218   int j, iopt, jopt;
    1219   float sum, f1, f2, fo, r, ro, lp;
    1220   float *dr = C->wrow, *dc = C->wcol;
    1221 
    1222   fo = 1.0e20;
    1223   ro = 0.0;
    1224   iopt = jopt = -1;
    1225 
    1226   s_n--;
    1227   s_m--;
    1228   if (s_m == 0)
    1229     return 0;
    1230   if (s_n == 0)
    1231   {
    1232     p = this->mpRowAdr(row)[qcol[0]];
    1233     if (p)
    1234     {
    1235       f1 = mp_PolyWeight(p, R);
    1236       if (f1 < fo)
    1237       {
    1238         fo = f1;
    1239         if (iopt >= 0)
    1240         p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
    1241         iopt = row;
    1242       }
    1243       else
    1244         p_Delete(&(this->mpRowAdr(row)[qcol[0]]), R);
    1245     }
    1246     if (iopt >= 0)
    1247       mpReplace(iopt, s_m, sign, qrow);
    1248     return 0;
    1249   }
    1250   this->mpRowWeight(dr);
    1251   this->mpColWeight(dc);
    1252   sum = 0.0;
    1253   for(j=s_m; j>=0; j--)
    1254     sum += dr[j];
    1255   r = dr[row];
    1256   a = this->mpRowAdr(row);
    1257   for(j=s_n; j>=0; j--)
    1258   {
    1259     p = a[qcol[j]];
    1260     if (p)
    1261     {
    1262       lp = mp_PolyWeight(p, R);
    1263       ro = r - lp;
    1264       f1 = ro * (dc[j]-lp);
    1265       if (f1 != 0.0)
    1266       {
    1267         f2 = lp * (sum - ro - dc[j]);
    1268         f2 += f1;
    1269       }
    1270       else
    1271         f2 = lp-r-dc[j];
    1272       if (f2 < fo)
    1273       {
    1274         fo = f2;
    1275         iopt = row;
    1276         jopt = j;
    1277       }
    1278     }
    1279   }
    1280   if (iopt < 0)
    1281     return 0;
    1282   mpReplace(iopt, s_m, sign, qrow);
    1283   mpReplace(jopt, s_n, sign, qcol);
    1284   return 1;
    1285 }
    1286 
    1287 void mp_permmatrix::mpToIntvec(intvec *v)
    1288 {
    1289   int i;
    1290 
    1291   for (i=v->rows()-1; i>=0; i--)
    1292     (*v)[i] = qcol[i]+1;
    1293 }
    1294 
    1295 void mp_permmatrix::mpRowReorder()
    1296 {
    1297   int k, i, i1, i2;
    1298 
    1299   if (a_m > a_n)
    1300     k = a_m - a_n;
    1301   else
    1302     k = 0;
    1303   for (i=a_m-1; i>=k; i--)
    1304   {
    1305     i1 = qrow[i];
    1306     if (i1 != i)
    1307     {
    1308       this->mpRowSwap(i1, i);
    1309       i2 = 0;
    1310       while (qrow[i2] != i) i2++;
    1311       qrow[i2] = i1;
    1312     }
    1313   }
    1314 }
    1315 
    1316 void mp_permmatrix::mpColReorder()
    1317 {
    1318   int k, j, j1, j2;
    1319 
    1320   if (a_n > a_m)
    1321     k = a_n - a_m;
    1322   else
    1323     k = 0;
    1324   for (j=a_n-1; j>=k; j--)
    1325   {
    1326     j1 = qcol[j];
    1327     if (j1 != j)
    1328     {
    1329       this->mpColSwap(j1, j);
    1330       j2 = 0;
    1331       while (qcol[j2] != j) j2++;
    1332       qcol[j2] = j1;
    1333     }
    1334   }
    1335 }
    1336 
    1337 // private
    1338 void mp_permmatrix::mpInitMat()
    1339 {
    1340   int k;
    1341 
    1342   s_m = a_m;
    1343   s_n = a_n;
    1344   piv_s = 0;
    1345   qrow = (int *)omAlloc(a_m*sizeof(int));
    1346   qcol = (int *)omAlloc(a_n*sizeof(int));
    1347   for (k=a_m-1; k>=0; k--) qrow[k] = k;
    1348   for (k=a_n-1; k>=0; k--) qcol[k] = k;
    1349 }
    1350 
    1351 poly * mp_permmatrix::mpRowAdr(int r)
    1352 {
    1353   return &(Xarray[a_n*qrow[r]]);
    1354 }
    1355 
    1356 poly * mp_permmatrix::mpColAdr(int c)
    1357 {
    1358   return &(Xarray[qcol[c]]);
    1359 }
    1360 
    1361 void mp_permmatrix::mpRowWeight(float *wrow)
    1362 {
    1363   poly p, *a;
    1364   int i, j;
    1365   float count;
    1366 
    1367   for (i=s_m; i>=0; i--)
    1368   {
    1369     a = this->mpRowAdr(i);
    1370     count = 0.0;
    1371     for(j=s_n; j>=0; j--)
    1372     {
    1373       p = a[qcol[j]];
    1374       if (p)
    1375         count += mp_PolyWeight(p, R);
    1376     }
    1377     wrow[i] = count;
    1378   }
    1379 }
    1380 
    1381 void mp_permmatrix::mpColWeight(float *wcol)
    1382 {
    1383   poly p, *a;
    1384   int i, j;
    1385   float count;
    1386 
    1387   for (j=s_n; j>=0; j--)
    1388   {
    1389     a = this->mpColAdr(j);
    1390     count = 0.0;
    1391     for(i=s_m; i>=0; i--)
    1392     {
    1393       p = a[a_n*qrow[i]];
    1394       if (p)
    1395         count += mp_PolyWeight(p, R);
    1396     }
    1397     wcol[j] = count;
    1398   }
    1399 }
    1400 
    1401 void mp_permmatrix::mpRowSwap(int i1, int i2)
    1402 {
    1403    poly p, *a1, *a2;
    1404    int j;
    1405 
    1406    a1 = &(Xarray[a_n*i1]);
    1407    a2 = &(Xarray[a_n*i2]);
    1408    for (j=a_n-1; j>= 0; j--)
    1409    {
    1410      p = a1[j];
    1411      a1[j] = a2[j];
    1412      a2[j] = p;
    1413    }
    1414 }
    1415 
    1416 void mp_permmatrix::mpColSwap(int j1, int j2)
    1417 {
    1418    poly p, *a1, *a2;
    1419    int i, k = a_n*a_m;
    1420 
    1421    a1 = &(Xarray[j1]);
    1422    a2 = &(Xarray[j2]);
    1423    for (i=0; i< k; i+=a_n)
    1424    {
    1425      p = a1[i];
    1426      a1[i] = a2[i];
    1427      a2[i] = p;
    1428    }
    1429 }
    1430 
    1431 int mp_permmatrix::mpGetRow()
    1432 {
    1433   return qrow[s_m];
    1434 }
    1435 
    1436 int mp_permmatrix::mpGetCol()
    1437 {
    1438   return qcol[s_n];
    1439 }
    1440 
    1441 /// perform replacement for pivot strategy in Bareiss algorithm
    1442 /// change sign of determinant
    1443 static void mpReplace(int j, int n, int &sign, int *perm)
    1444 {
    1445   int k;
    1446 
    1447   if (j != n)
    1448   {
    1449     k = perm[n];
    1450     perm[n] = perm[j];
    1451     perm[j] = k;
    1452     sign = -sign;
    1453   }
    1454 }
    1455 
    1456 static int mpNextperm(perm * z, int max)
    1457 {
    1458   int s, i, k, t;
    1459   s = max;
    1460   do
    1461   {
    1462     s--;
    1463   }
    1464   while ((s > 0) && ((*z)[s] >= (*z)[s+1]));
    1465   if (s==0)
    1466     return 0;
    1467   do
    1468   {
    1469     (*z)[s]++;
    1470     k = 0;
    1471     do
    1472     {
    1473       k++;
    1474     }
    1475     while (((*z)[k] != (*z)[s]) && (k!=s));
    1476   }
    1477   while (k < s);
    1478   for (i=s+1; i <= max; i++)
    1479   {
    1480     (*z)[i]=0;
    1481     do
    1482     {
    1483       (*z)[i]++;
    1484       k=0;
    1485       do
    1486       {
    1487         k++;
    1488       }
    1489       while (((*z)[k] != (*z)[i]) && (k != i));
    1490     }
    1491     while (k < i);
    1492   }
    1493   s = max+1;
    1494   do
    1495   {
    1496     s--;
    1497   }
    1498   while ((s > 0) && ((*z)[s] > (*z)[s+1]));
    1499   t = 1;
    1500   for (i=1; i<max; i++)
    1501     for (k=i+1; k<=max; k++)
    1502       if ((*z)[k] < (*z)[i])
    1503         t = -t;
    1504   (*z)[0] = t;
    1505   return s;
    1506 }
    1507 
    1508 static poly mp_Leibnitz(matrix a, const ring R)
    1509 {
    1510   int i, e, n;
    1511   poly p, d;
    1512   perm z;
    1513 
    1514   n = MATROWS(a);
    1515   memset(&z,0,(n+2)*sizeof(int));
    1516   p = p_One(R);
    1517   for (i=1; i <= n; i++)
    1518     p = p_Mult_q(p, p_Copy(MATELEM(a, i, i), R), R);
    1519   d = p;
    1520   for (i=1; i<= n; i++)
    1521     z[i] = i;
    1522   z[0]=1;
    1523   e = 1;
    1524   if (n!=1)
    1525   {
    1526     while (e)
    1527     {
    1528       e = mpNextperm((perm *)&z, n);
    1529       p = p_One(R);
    1530       for (i = 1; i <= n; i++)
    1531         p = p_Mult_q(p, p_Copy(MATELEM(a, i, z[i]), R), R);
    1532       if (z[0] > 0)
    1533         d = p_Add_q(d, p, R);
    1534       else
    1535         d = p_Sub(d, p, R);
    1536     }
    1537   }
    1538   return d;
    1539 }
    1540 
    1541612static poly minuscopy (poly p, const ring R)
    1542613{
    … …  
    1664735}
    1665736
    1666 /*2
    1667 *  prepare one step of 'Bareiss' algorithm
    1668 *  for application in minor
    1669 */
    1670 static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
    1671 {
    1672   int r;
    1673 
    1674   r = mp_PivBar(a,lr,lc, R);
    1675   if(r==0) return 0;
    1676   if(r<lr) mp_SwapRow(a, r, lr, lc, R);
    1677   return 1;
    1678 }
    1679 
    1680 /*2
    1681 *  prepare one step of 'Bareiss' algorithm
    1682 *  for application in minor
    1683 */
    1684 static int mp_PreparePiv (matrix a, int lr, int lc, const ring R)
    1685 {
    1686   int c;
    1687 
    1688   c = mp_PivRow(a, lr, lc, R);
    1689   if(c==0) return 0;
    1690   if(c<lc) mp_SwapCol(a, c, lr, lc, R);
    1691   return 1;
    1692 }
    1693 
    1694 /*
    1695 * find best row
    1696 */
    1697 static int mp_PivBar(matrix a, int lr, int lc, const ring R)
    1698 {
    1699   float f1, f2;
    1700   poly *q1;
    1701   int i,j,io;
    1702 
    1703   io = -1;
    1704   f1 = 1.0e30;
    1705   for (i=lr-1;i>=0;i--)
    1706   {
    1707     q1 = &(a->m)[i*a->ncols];
    1708     f2 = 0.0;
    1709     for (j=lc-1;j>=0;j--)
    1710     {
    1711       if (q1[j]!=NULL)
    1712         f2 += mp_PolyWeight(q1[j], R);
    1713     }
    1714     if ((f2!=0.0) && (f2<f1))
    1715     {
    1716       f1 = f2;
    1717       io = i;
    1718     }
    1719   }
    1720   if (io<0) return 0;
    1721   else return io+1;
    1722 }
    1723 
    1724 /*
    1725 * find pivot in the last row
    1726 */
    1727 static int mp_PivRow(matrix a, int lr, int lc, const ring R)
    1728 {
    1729   float f1, f2;
    1730   poly *q1;
    1731   int j,jo;
    1732 
    1733   jo = -1;
    1734   f1 = 1.0e30;
    1735   q1 = &(a->m)[(lr-1)*a->ncols];
    1736   for (j=lc-1;j>=0;j--)
    1737   {
    1738     if (q1[j]!=NULL)
    1739     {
    1740       f2 = mp_PolyWeight(q1[j], R);
    1741       if (f2<f1)
    1742       {
    1743         f1 = f2;
    1744         jo = j;
    1745       }
    1746     }
    1747   }
    1748   if (jo<0) return 0;
    1749   else return jo+1;
    1750 }
    1751 
    1752 /*
    1753 * weigth of a polynomial, for pivot strategy
    1754 */
    1755 static float mp_PolyWeight(poly p, const ring R)
    1756 {
    1757   int i;
    1758   float res;
    1759 
    1760   if (pNext(p) == NULL)
    1761   {
    1762     res = (float)n_Size(p_GetCoeff(p, R), R);
    1763     for (i=rVar(R);i>0;i--)
    1764     {
    1765       if(p_GetExp(p,i, R)!=0)
    1766       {
    1767         res += 2.0;
    1768         break;
    1769       }
    1770     }
    1771   }
    1772   else
    1773   {
    1774     res = 0.0;
    1775     do
    1776     {
    1777       res += (float)n_Size(p_GetCoeff(p, R), R) + 2.0;
    1778       pIter(p);
    1779     }
    1780     while (p);
    1781   }
    1782   return res;
    1783 }
    1784 
    1785 static void mpSwapRow(matrix a, int pos, int lr, int lc)
    1786 {
    1787   poly sw;
    1788   int j;
    1789   poly* a2 = a->m;
    1790   poly* a1 = &a2[a->ncols*(pos-1)];
    1791 
    1792   a2 = &a2[a->ncols*(lr-1)];
    1793   for (j=lc-1; j>=0; j--)
    1794   {
    1795     sw = a1[j];
    1796     a1[j] = a2[j];
    1797     a2[j] = sw;
    1798   }
    1799 }
    1800 
    1801 static void mpSwapCol(matrix a, int pos, int lr, int lc)
    1802 {
    1803   poly sw;
    1804   int j;
    1805   poly* a2 = a->m;
    1806   poly* a1 = &a2[pos-1];
    1807 
    1808   a2 = &a2[lc-1];
    1809   for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
    1810   {
    1811     sw = a1[j];
    1812     a1[j] = a2[j];
    1813     a2[j] = sw;
    1814   }
    1815 }
    1816 
    1817 static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
    1818 {
    1819   int r=lr-1, c=lc-1;
    1820   poly *b = a0->m, *x = re->m;
    1821   poly piv, elim, q1, q2, *ap, *a, *q;
    1822   int i, j;
    1823 
    1824   ap = &b[r*a0->ncols];
    1825   piv = ap[c];
    1826   for(j=c-1; j>=0; j--)
    1827     if (ap[j] != NULL) ap[j] = p_Neg(ap[j], R);
    1828   for(i=r-1; i>=0; i--)
    1829   {
    1830     a = &b[i*a0->ncols];
    1831     q = &x[i*re->ncols];
    1832     if (a[c] != NULL)
    1833     {
    1834       elim = a[c];
    1835       for (j=c-1; j>=0; j--)
    1836       {
    1837         q1 = NULL;
    1838         if (a[j] != NULL)
    1839         {
    1840           q1 = SM_MULT(a[j], piv, div, R);
    1841           if (ap[j] != NULL)
    1842           {
    1843             q2 = SM_MULT(ap[j], elim, div, R);
    1844             q1 = p_Add_q(q1,q2, R);
    1845           }
    1846         }
    1847         else if (ap[j] != NULL)
    1848           q1 = SM_MULT(ap[j], elim, div, R);
    1849         if (q1 != NULL)
    1850         {
    1851           if (div)
    1852             SM_DIV(q1, div, R);
    1853           q[j] = q1;
    1854         }
    1855       }
    1856     }
    1857     else
    1858     {
    1859       for (j=c-1; j>=0; j--)
    1860       {
    1861         if (a[j] != NULL)
    1862         {
    1863           q1 = SM_MULT(a[j], piv, div, R);
    1864           if (div)
    1865             SM_DIV(q1, div, R);
    1866           q[j] = q1;
    1867         }
    1868       }
    1869     }
    1870   }
    1871 }
    1872 
    1873737BOOLEAN mp_IsDiagUnit(matrix U, const ring R)
    1874738{
Note: See TracChangeset for help on using the changeset viewer.