Changeset 441a2e in git


Ignore:
Timestamp:
Jul 26, 2011, 4:40:43 PM (12 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')
Children:
bf6a4d5920ef44ffb48d0fc83b4d53589111d60c
Parents:
b28bafedafd8d93bd2ff41d9466d46678bec1341
git-author:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2011-07-26 16:40:43+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:53:35+01:00
Message:
restore: mp_RecMin, mp_detBareiss -> matpol.cc
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • Singular/eigenval_ip.cc

    rb28bafe r441a2e  
    180180
    181181      intvec *m0;
    182       ideal e0=singclap_factorize(mpDetBareiss(M0),&m0,2);
     182      ideal e0=singclap_factorize(mp_DetBareiss(M,currRing)),&m0,2);
    183183      if (e0==NULL)
    184184      {

  • Singular/iparith.cc

    rb28bafe r441a2e  
    50235023static BOOLEAN jjmpDetBareiss(leftv res, leftv v)
    50245024{
    5025   res->data = (char *)mpDetBareiss((matrix)v->Data());
     5025  res->data = (char *)mp_DetBareiss((matrix)v->Data(),currRing);
    50265026  return FALSE;
    50275027}

  • kernel/ideals.cc

    rb28bafe r441a2e  
    543543        p = pCopy(tempstd->m[j]);
    544544      else
    545         p = prCopyR(tempstd->m[j], syz_ring);
     545        p = prCopyR(tempstd->m[j], syz_ring,currRing);
    546546      p_Shift(&p,-syzComp-isIdeal,currRing);
    547547      result->m[k] = p;
    … …  
    694694  assume(currRing != NULL);
    695695  ring orig_ring=currRing;
    696   ring syz_ring=rAssure_SyzComp(orig_ring.TRUE);
     696  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
    697697  rChangeCurrRing(syz_ring);
    698698
    … …  
    702702  if (orig_ring != syz_ring)
    703703  {
    704     s_h1=idrCopyR_NoSort(h1,orig_ring);
     704    s_h1=idrCopyR_NoSort(h1,orig_ring,syz_ring);
    705705  }
    706706  else
    … …  
    774774  )
    775775  {
    776     ring dp_C_ring = rCurrRingAssure_dp_C();
     776    ring dp_C_ring = rAssure_dp_C(syz_ring);
    777777    if (dp_C_ring != syz_ring)
     778    {
     779      rChangeCurrRing(dp_C_ring);
    778780      e = idrMoveR_NoSort(e, syz_ring, dp_C_ring);
     781    }
    779782    resolvente res = sySchreyerResolvente(e,-1,&length,TRUE, TRUE);
    780783    intvec * dummy = syBetti(res,length,&reg, *w);
    … …  
    824827  if (orig_ring != syz_ring)
    825828  {
    826     s_h1=idrCopyR_NoSort(h1,orig_ring);
     829    s_h1=idrCopyR_NoSort(h1,orig_ring, syz_ring);
    827830  }
    828831  else
    … …  
    898901
    899902  if (orig_ring != syz_ring)
    900     s_h1 = idrCopyR_NoSort(h1,orig_ring);
     903    s_h1 = idrCopyR_NoSort(h1,orig_ring,syz_ring);
    901904  else
    902905    s_h1 = h1;
    … …  
    11061109  if (orig_ring != syz_ring)
    11071110  {
    1108     s_mod = idrCopyR_NoSort(mod,orig_ring);
    1109     s_temp = idrCopyR_NoSort(submod,orig_ring);
     1111    s_mod = idrCopyR_NoSort(mod,orig_ring,syz_ring);
     1112    s_temp = idrCopyR_NoSort(submod,orig_ring,syz_ring);
    11101113  }
    11111114  else
    … …  
    12881291      if(pDivisibleBy(Q->m[j],p))
    12891292      {
    1290         poly p0=pDivideM(pHead(p),pHead(Q->m[j]));
     1293        poly p0=p_DivideM(pHead(p),pHead(Q->m[j]),currRing);
    12911294        if(w==NULL)
    12921295          p=pJet(pSub(p,ppMult_mm(Q->m[j],p0)),N);
    … …  
    14601463
    14611464  ring orig_ring=currRing;
    1462   ring syz_ring=rAssure_SyzComp(orig_ring.TRUE);
     1465  ring syz_ring=rAssure_SyzComp(orig_ring,TRUE);
    14631466  rChangeCurrRing(syz_ring);
    14641467  rSetSyzComp(kmax-1,syz_ring);
    … …  
    18061809          }
    18071810        }
    1808         p = mpDetBareiss(tmp);
     1811        p = mp_DetBareiss(tmp,currRing);
    18091812        if (p!=NULL)
    18101813        {
    … …  
    18701873        }
    18711874      }
    1872       p = mpDetBareiss(tmp);
     1875      p = mp_DetBareiss(tmp,vcurrRing);
    18731876      if (p!=NULL)
    18741877      {
    … …  
    19331936    return NULL;
    19341937  }
    1935   h = idMatrix2Module(mpCopy(a));
    1936   bound = smExpBound(h,c,r,ar);
     1938  h = idMatrix2Module(mp_Copy(a,currRing));
     1939  bound = sm_ExpBound(h,c,r,ar,currRing);
    19371940  idDelete(&h);
    1938   tmpR=smRingChange(&origR,bound);
     1941  tmpR=sm_RingChange(origR,bound);
    19391942  b = mpNew(r,c);
    19401943  for (i=r*c-1;i>=0;i--)
    19411944  {
    19421945    if (a->m[i])
    1943       b->m[i] = prCopyR(a->m[i],origR);
     1946      b->m[i] = prCopyR(a->m[i],origR,currRing);
    19441947  }
    19451948  if (R!=NULL)
    19461949  {
    1947     R = idrCopyR(R,origR);
     1950    R = idrCopyR(R,origR,currRing);
    19481951    //if (ar>1) // otherwise done in mpMinorToResult
    19491952    //{
    … …  
    19601963  idSkipZeroes(result);
    19611964  rChangeCurrRing(origR);
    1962   result = idrMoveR(result,tmpR);
     1965  result = idrMoveR(result,tmpR,origR);
    19631966  smKillModifiedRing(tmpR);
    19641967  idTest(result);

  • libpolys/polys/matpol.cc

    rb28bafe r441a2e  
    3636#include "prCopy.h"
    3737
    38 // #include "sparsmat.h"
     38#include "sparsmat.h"
    3939   
    4040//omBin ip_smatrix_bin = omGetSpecBin(sizeof(ip_smatrix));
    … …  
    798798  id_Delete((ideal *) a, r);
    799799}
     800
     801/*
     802* C++ classes for Bareiss algorithm
     803*/
     804class row_col_weight
     805{
     806  private:
     807  int ym, yn;
     808  public:
     809  float *wrow, *wcol;
     810  row_col_weight() : ym(0) {}
     811  row_col_weight(int, int);
     812  ~row_col_weight();
     813};
     814
     815/*2
     816*  a submatrix M of a matrix X[m,n]:
     817*    0 <= i < s_m <= a_m
     818*    0 <= j < s_n <= a_n
     819*    M = ( Xarray[qrow[i],qcol[j]] )
     820*    if a_m = a_n and s_m = s_n
     821*      det(X) = sign*div^(s_m-1)*det(M)
     822*    resticted pivot for elimination
     823*      0 <= j < piv_s
     824*/
     825class mp_permmatrix
     826{
     827  private:
     828  int       a_m, a_n, s_m, s_n, sign, piv_s;
     829  int       *qrow, *qcol;
     830  poly      *Xarray;
     831  void mpInitMat();
     832  poly * mpRowAdr(int);
     833  poly * mpColAdr(int);
     834  void mpRowWeight(float *);
     835  void mpColWeight(float *);
     836  void mpRowSwap(int, int);
     837  void mpColSwap(int, int);
     838  public:
     839  mp_permmatrix() : a_m(0) {}
     840  mp_permmatrix(matrix);
     841  mp_permmatrix(mp_permmatrix *);
     842  ~mp_permmatrix();
     843  int mpGetRow();
     844  int mpGetCol();
     845  int mpGetRdim();
     846  int mpGetCdim();
     847  int mpGetSign();
     848  void mpSetSearch(int s);
     849  void mpSaveArray();
     850  poly mpGetElem(int, int);
     851  void mpSetElem(poly, int, int);
     852  void mpDelElem(int, int);
     853  void mpElimBareiss(poly);
     854  int mpPivotBareiss(row_col_weight *);
     855  int mpPivotRow(row_col_weight *, int);
     856  void mpToIntvec(intvec *);
     857  void mpRowReorder();
     858  void mpColReorder();
     859};
     860
     861/*
     862* weigth of a polynomial, for pivot strategy
     863*/
     864static float mp_PolyWeight(poly p, const ring r)
     865{
     866  int i;
     867  float res;
     868
     869  if (pNext(p) == NULL)
     870  {
     871    res = (float)n_Size(pGetCoeff(p),r->cf);
     872    for (i=rVar(r);i>0;i--)
     873    {
     874      if(p_GetExp(p,i,r)!=0)
     875      {
     876        res += 2.0;
     877        break;
     878      }
     879    }
     880  }
     881  else
     882  {
     883    res = 0.0;
     884    do
     885    {
     886      res += (float)n_Size(pGetCoeff(p),r->cf)+2.0;
     887      pIter(p);
     888    }
     889    while (p);
     890  }
     891  return res;
     892}
     893/*
     894* find best row
     895*/
     896static int mp_PivBar(matrix a, int lr, int lc, const ring r)
     897{
     898  float f1, f2;
     899  poly *q1;
     900  int i,j,io;
     901
     902  io = -1;
     903  f1 = 1.0e30;
     904  for (i=lr-1;i>=0;i--)
     905  {
     906    q1 = &(a->m)[i*a->ncols];
     907    f2 = 0.0;
     908    for (j=lc-1;j>=0;j--)
     909    {
     910      if (q1[j]!=NULL)
     911        f2 += mp_PolyWeight(q1[j],r);
     912    }
     913    if ((f2!=0.0) && (f2<f1))
     914    {
     915      f1 = f2;
     916      io = i;
     917    }
     918  }
     919  if (io<0) return 0;
     920  else return io+1;
     921}
     922
     923static void mpSwapRow(matrix a, int pos, int lr, int lc)
     924{
     925  poly sw;
     926  int j;
     927  poly* a2 = a->m;
     928  poly* a1 = &a2[a->ncols*(pos-1)];
     929
     930  a2 = &a2[a->ncols*(lr-1)];
     931  for (j=lc-1; j>=0; j--)
     932  {
     933    sw = a1[j];
     934    a1[j] = a2[j];
     935    a2[j] = sw;
     936  }
     937}
     938
     939/*2
     940*  prepare one step of 'Bareiss' algorithm
     941*  for application in minor
     942*/
     943static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
     944{
     945  int r;
     946
     947  r = mp_PivBar(a,lr,lc,R);
     948  if(r==0) return 0;
     949  if(r<lr) mpSwapRow(a, r, lr, lc);
     950  return 1;
     951}
     952
     953/*
     954* find pivot in the last row
     955*/
     956static int mp_PivRow(matrix a, int lr, int lc, const ring r)
     957{
     958  float f1, f2;
     959  poly *q1;
     960  int j,jo;
     961
     962  jo = -1;
     963  f1 = 1.0e30;
     964  q1 = &(a->m)[(lr-1)*a->ncols];
     965  for (j=lc-1;j>=0;j--)
     966  {
     967    if (q1[j]!=NULL)
     968    {
     969      f2 = mp_PolyWeight(q1[j],r);
     970      if (f2<f1)
     971      {
     972        f1 = f2;
     973        jo = j;
     974      }
     975    }
     976  }
     977  if (jo<0) return 0;
     978  else return jo+1;
     979}
     980
     981static void mpSwapCol(matrix a, int pos, int lr, int lc)
     982{
     983  poly sw;
     984  int j;
     985  poly* a2 = a->m;
     986  poly* a1 = &a2[pos-1];
     987
     988  a2 = &a2[lc-1];
     989  for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
     990  {
     991    sw = a1[j];
     992    a1[j] = a2[j];
     993    a2[j] = sw;
     994  }
     995}
     996
     997/*2
     998*  prepare one step of 'Bareiss' algorithm
     999*  for application in minor
     1000*/
     1001static int mp_PreparePiv (matrix a, int lr, int lc,const ring r)
     1002{
     1003  int c;
     1004
     1005  c = mp_PivRow(a, lr, lc,r);
     1006  if(c==0) return 0;
     1007  if(c<lc) mpSwapCol(a, c, lr, lc);
     1008  return 1;
     1009}
     1010
     1011static inline BOOLEAN smSmaller(poly a, poly b)
     1012{
     1013  loop
     1014  {
     1015    pIter(b);
     1016    if (b == NULL) return TRUE;
     1017    pIter(a);
     1018    if (a == NULL) return FALSE;
     1019  }
     1020}
     1021
     1022static BOOLEAN sm_IsNegQuot(poly a, const poly b, const poly c, const ring R)
     1023{
     1024  if (p_LmDivisibleByNoComp(c, b, R))
     1025  {
     1026    p_ExpVectorDiff(a, b, c, R);
     1027    // Hmm: here used to be a pSetm(a): but it is unnecessary,
     1028    // if b and c are correct
     1029    return FALSE;
     1030  }
     1031  else
     1032  {
     1033    int i;
     1034    for (i=rVar(R); i>0; i--)
     1035    {
     1036      if(p_GetExp(c,i,R) > p_GetExp(b,i,R))
     1037        p_SetExp(a,i,p_GetExp(c,i,R)-p_GetExp(b,i,R),R);
     1038      else
     1039        p_SetExp(a,i,0,R);
     1040    }
     1041    // here we actually might need a pSetm, if a is to be used in
     1042    // comparisons
     1043    return TRUE;
     1044  }
     1045}
     1046
     1047static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
     1048{
     1049  int r=lr-1, c=lc-1;
     1050  poly *b = a0->m, *x = re->m;
     1051  poly piv, elim, q1, q2, *ap, *a, *q;
     1052  int i, j;
     1053
     1054  ap = &b[r*a0->ncols];
     1055  piv = ap[c];
     1056  for(j=c-1; j>=0; j--)
     1057    if (ap[j] != NULL) ap[j] = p_Neg(ap[j],R);
     1058  for(i=r-1; i>=0; i--)
     1059  {
     1060    a = &b[i*a0->ncols];
     1061    q = &x[i*re->ncols];
     1062    if (a[c] != NULL)
     1063    {
     1064      elim = a[c];
     1065      for (j=c-1; j>=0; j--)
     1066      {
     1067        q1 = NULL;
     1068        if (a[j] != NULL)
     1069        {
     1070          q1 = sm_MultDiv(a[j], piv, div,R);
     1071          if (ap[j] != NULL)
     1072          {
     1073            q2 = sm_MultDiv(ap[j], elim, div, R);
     1074            q1 = p_Add_q(q1,q2,R);
     1075          }
     1076        }
     1077        else if (ap[j] != NULL)
     1078          q1 = sm_MultDiv(ap[j], elim, div, R);
     1079        if (q1 != NULL)
     1080        {
     1081          if (div)
     1082            sm_SpecialPolyDiv(q1, div,R);
     1083          q[j] = q1;
     1084        }
     1085      }
     1086    }
     1087    else
     1088    {
     1089      for (j=c-1; j>=0; j--)
     1090      {
     1091        if (a[j] != NULL)
     1092        {
     1093          q1 = sm_MultDiv(a[j], piv, div, R);
     1094          if (div)
     1095            sm_SpecialPolyDiv(q1, div, R);
     1096          q[j] = q1;
     1097        }
     1098      }
     1099    }
     1100  }
     1101}
     1102
     1103/*2
     1104* entries of a are minors and go to result (only if not in R)
     1105*/
     1106void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
     1107                     ideal R, const ring ri)
     1108{
     1109  poly *q1;
     1110  int e=IDELEMS(result);
     1111  int i,j;
     1112
     1113  if (R != NULL)
     1114  {
     1115    for (i=r-1;i>=0;i--)
     1116    {
     1117      q1 = &(a->m)[i*a->ncols];
     1118      //for (j=c-1;j>=0;j--)
     1119      //{
     1120      //  if (q1[j]!=NULL) q1[j] = kNF(R,currQuotient,q1[j]);
     1121      //}
     1122    }
     1123  }
     1124  for (i=r-1;i>=0;i--)
     1125  {
     1126    q1 = &(a->m)[i*a->ncols];
     1127    for (j=c-1;j>=0;j--)
     1128    {
     1129      if (q1[j]!=NULL)
     1130      {
     1131        if (elems>=e)
     1132        {
     1133          pEnlargeSet(&(result->m),e,e);
     1134          e += e;
     1135          IDELEMS(result) =e;
     1136        }
     1137        result->m[elems] = q1[j];
     1138        q1[j] = NULL;
     1139        elems++;
     1140      }
     1141    }
     1142  }
     1143}
     1144
     1145static void mpFinalClean(matrix a)
     1146{
     1147  omFreeSize((ADDRESS)a->m,a->nrows*a->ncols*sizeof(poly));
     1148  omFreeBin((ADDRESS)a, ip_smatrix_bin);
     1149}
     1150
     1151/*2
     1152* produces recursively the ideal of all arxar-minors of a
     1153*/
     1154void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
     1155              poly barDiv, ideal R, const ring r)
     1156{
     1157  int k;
     1158  int kr=lr-1,kc=lc-1;
     1159  matrix nextLevel=mpNew(kr,kc);
     1160
     1161  loop
     1162  {
     1163/*--- look for an optimal row and bring it to last position ------------*/
     1164    if(mp_PrepareRow(a,lr,lc,r)==0) break;
     1165/*--- now take all pivots from the last row ------------*/
     1166    k = lc;
     1167    loop
     1168    {
     1169      if(mp_PreparePiv(a,lr,k,r)==0) break;
     1170      mp_ElimBar(a,nextLevel,barDiv,lr,k,r);
     1171      k--;
     1172      if (ar>1)
     1173      {
     1174        mp_RecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R,r);
     1175        mp_PartClean(nextLevel,kr,k, r);
     1176      }
     1177      else mp_MinorToResult(result,elems,nextLevel,kr,k,R,r);
     1178      if (ar>k-1) break;
     1179    }
     1180    if (ar>=kr) break;
     1181/*--- now we have to take out the last row...------------*/
     1182    lr = kr;
     1183    kr--;
     1184  }
     1185  mpFinalClean(nextLevel);
     1186}
     1187
     1188/*2
     1189*returns the determinant of the matrix m;
     1190*uses Bareiss algorithm
     1191*/
     1192poly mp_DetBareiss (matrix a, const ring r)
     1193{
     1194  int s;
     1195  poly div, res;
     1196  if (MATROWS(a) != MATCOLS(a))
     1197  {
     1198    Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
     1199    return NULL;
     1200  }
     1201  matrix c = mp_Copy(a,r);
     1202  mp_permmatrix *Bareiss = new mp_permmatrix(c);
     1203  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
     1204
     1205  /* Bareiss */
     1206  div = NULL;
     1207  while(Bareiss->mpPivotBareiss(&w))
     1208  {
     1209    Bareiss->mpElimBareiss(div);
     1210    div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
     1211  }
     1212  Bareiss->mpRowReorder();
     1213  Bareiss->mpColReorder();
     1214  Bareiss->mpSaveArray();
     1215  s = Bareiss->mpGetSign();
     1216  delete Bareiss;
     1217
     1218  /* result */
     1219  res = MATELEM(c,1,1);
     1220  MATELEM(c,1,1) = NULL;
     1221  id_Delete((ideal *)&c,r);
     1222  if (s < 0)
     1223    res = p_Neg(res,r);
     1224  return res;
     1225}

  • libpolys/polys/matpol.h

    rb28bafe r441a2e  
    6161// BOOLEAN mpJacobi(leftv res,leftv a);
    6262// BOOLEAN mpKoszul(leftv res,leftv b/*in*/, leftv c/*ip*/, leftv id=NULL);
    63 // poly mp_DetBareiss (matrix a, const ring r);
     63poly mp_DetBareiss (matrix a, const ring r);
    6464
    6565//matrix mp_Homogen(matrix a, int v, const ring r);
    … …  
    8181
    8282/// for minors with Bareiss
    83 // void mp_RecMin(int, ideal, int &, matrix, int, int, poly, ideal, const ring r);
     83void mp_RecMin(int, ideal, int &, matrix, int, int, poly, ideal, const ring r);
    8484// void mp_MinorToResult(ideal, int &, matrix, int, int, ideal, const ring r);
    8585
Note: See TracChangeset for help on using the changeset viewer.