Changeset bf6a4d in git


Ignore:
Timestamp:
Jul 26, 2011, 7:37:35 PM (12 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
73ad0c1da96d24bf6cc659867bd85c7b02df2b6d
Parents:
441a2ebc3368fca13446fd60801db4c09580a110
git-author:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2011-07-26 19:37:35+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 13:01:15+01:00
Message:
add and fix: sparsmat.cc
Location:
libpolys/polys
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • libpolys/polys/Makefile.am

    r441a2e rbf6a4d  
    5050        pDebug.cc polys0.cc prCopy.cc prCopyMacros.h \
    5151        kbuckets.cc sbuckets.cc weight.cc weight0.c simpleideals.cc matpol.cc \
     52        sparsmat.cc \
    5253        ${USE_P_PROCS_STATIC_CC} ${USE_P_PROCS_DYNAMIC_CC}  mod_raw.cc \
    5354        ext_fields/algext.cc ext_fields/algext.h \
    … …  
    6869        monomials/p_polys.h monomials/maps.h polys.h prCopy.h \
    6970        kbuckets.h sbuckets.h simpleideals.h weight.h matpol.h \
    70         clapsing.h clapconv.h
     71        sparsmat.h clapsing.h clapconv.h
    7172
    7273### nobase_include_HEADERS = $(LIBPOLYSHEADERS)

  • libpolys/polys/matpol.cc

    r441a2e rbf6a4d  
    813813};
    814814
     815row_col_weight::row_col_weight(int i, int j)
     816{
     817  ym = i;
     818  yn = j;
     819  wrow = (float *)omAlloc(i*sizeof(float));
     820  wcol = (float *)omAlloc(j*sizeof(float));
     821}
     822row_col_weight::~row_col_weight()
     823{
     824  if (ym!=0)
     825  {
     826    omFreeSize((ADDRESS)wcol, yn*sizeof(float));
     827    omFreeSize((ADDRESS)wrow, ym*sizeof(float));
     828  }
     829}
     830
    815831/*2
    816832*  a submatrix M of a matrix X[m,n]:
    … …  
    829845  int       *qrow, *qcol;
    830846  poly      *Xarray;
     847  ring _R;
    831848  void mpInitMat();
    832   poly * mpRowAdr(int);
    833   poly * mpColAdr(int);
     849  poly * mpRowAdr(int r)
     850  { return &(Xarray[a_n*qrow[r]]); }
     851  poly * mpColAdr(int c)
     852  { return &(Xarray[qcol[c]]); }
    834853  void mpRowWeight(float *);
    835854  void mpColWeight(float *);
    … …  
    838857  public:
    839858  mp_permmatrix() : a_m(0) {}
    840   mp_permmatrix(matrix);
     859  mp_permmatrix(matrix, ring);
    841860  mp_permmatrix(mp_permmatrix *);
    842861  ~mp_permmatrix();
    843862  int mpGetRow();
    844863  int mpGetCol();
    845   int mpGetRdim();
    846   int mpGetCdim();
    847   int mpGetSign();
     864  int mpGetRdim() { return s_m; }
     865  int mpGetCdim() { return s_n; }
     866  int mpGetSign() { return sign; }
    848867  void mpSetSearch(int s);
    849   void mpSaveArray();
     868  void mpSaveArray() { Xarray = NULL; }
    850869  poly mpGetElem(int, int);
    851870  void mpSetElem(poly, int, int);
    … …  
    858877  void mpColReorder();
    859878};
    860 
     879mp_permmatrix::mp_permmatrix(matrix A, ring R) : sign(1)
     880{
     881  a_m = A->nrows;
     882  a_n = A->ncols;
     883  this->mpInitMat();
     884  Xarray = A->m;
     885  _R=R;
     886}
     887
     888mp_permmatrix::mp_permmatrix(mp_permmatrix *M)
     889{
     890  poly p, *athis, *aM;
     891  int i, j;
     892
     893  _R=M->_R;
     894  a_m = M->s_m;
     895  a_n = M->s_n;
     896  sign = M->sign;
     897  this->mpInitMat();
     898  Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
     899  for (i=a_m-1; i>=0; i--)
     900  {
     901    athis = this->mpRowAdr(i);
     902    aM = M->mpRowAdr(i);
     903    for (j=a_n-1; j>=0; j--)
     904    {
     905      p = aM[M->qcol[j]];
     906      if (p)
     907      {
     908        athis[j] = p_Copy(p,_R);
     909      }
     910    }
     911  }
     912}
     913
     914mp_permmatrix::~mp_permmatrix()
     915{
     916  int k;
     917
     918  if (a_m != 0)
     919  {
     920    omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
     921    omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
     922    if (Xarray != NULL)
     923    {
     924      for (k=a_m*a_n-1; k>=0; k--)
     925        p_Delete(&Xarray[k],_R);
     926      omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
     927    }
     928  }
     929}
     930
     931
     932static float mp_PolyWeight(poly p, const ring r);
     933void mp_permmatrix::mpColWeight(float *wcol)
     934{
     935  poly p, *a;
     936  int i, j;
     937  float count;
     938
     939  for (j=s_n; j>=0; j--)
     940  {
     941    a = this->mpColAdr(j);
     942    count = 0.0;
     943    for(i=s_m; i>=0; i--)
     944    {
     945      p = a[a_n*qrow[i]];
     946      if (p)
     947        count += mp_PolyWeight(p,_R);
     948    }
     949    wcol[j] = count;
     950  }
     951}
     952void mp_permmatrix::mpRowWeight(float *wrow)
     953{
     954  poly p, *a;
     955  int i, j;
     956  float count;
     957
     958  for (i=s_m; i>=0; i--)
     959  {
     960    a = this->mpRowAdr(i);
     961    count = 0.0;
     962    for(j=s_n; j>=0; j--)
     963    {
     964      p = a[qcol[j]];
     965      if (p)
     966        count += mp_PolyWeight(p,_R);
     967    }
     968    wrow[i] = count;
     969  }
     970}
     971
     972void mp_permmatrix::mpRowSwap(int i1, int i2)
     973{
     974   poly p, *a1, *a2;
     975   int j;
     976
     977   a1 = &(Xarray[a_n*i1]);
     978   a2 = &(Xarray[a_n*i2]);
     979   for (j=a_n-1; j>= 0; j--)
     980   {
     981     p = a1[j];
     982     a1[j] = a2[j];
     983     a2[j] = p;
     984   }
     985}
     986
     987void mp_permmatrix::mpColSwap(int j1, int j2)
     988{
     989   poly p, *a1, *a2;
     990   int i, k = a_n*a_m;
     991
     992   a1 = &(Xarray[j1]);
     993   a2 = &(Xarray[j2]);
     994   for (i=0; i< k; i+=a_n)
     995   {
     996     p = a1[i];
     997     a1[i] = a2[i];
     998     a2[i] = p;
     999   }
     1000}
     1001void mp_permmatrix::mpInitMat()
     1002{
     1003  int k;
     1004
     1005  s_m = a_m;
     1006  s_n = a_n;
     1007  piv_s = 0;
     1008  qrow = (int *)omAlloc(a_m*sizeof(int));
     1009  qcol = (int *)omAlloc(a_n*sizeof(int));
     1010  for (k=a_m-1; k>=0; k--) qrow[k] = k;
     1011  for (k=a_n-1; k>=0; k--) qcol[k] = k;
     1012}
     1013
     1014void mp_permmatrix::mpColReorder()
     1015{
     1016  int k, j, j1, j2;
     1017
     1018  if (a_n > a_m)
     1019    k = a_n - a_m;
     1020  else
     1021    k = 0;
     1022  for (j=a_n-1; j>=k; j--)
     1023  {
     1024    j1 = qcol[j];
     1025    if (j1 != j)
     1026    {
     1027      this->mpColSwap(j1, j);
     1028      j2 = 0;
     1029      while (qcol[j2] != j) j2++;
     1030      qcol[j2] = j1;
     1031    }
     1032  }
     1033}
     1034
     1035void mp_permmatrix::mpRowReorder()
     1036{
     1037  int k, i, i1, i2;
     1038
     1039  if (a_m > a_n)
     1040    k = a_m - a_n;
     1041  else
     1042    k = 0;
     1043  for (i=a_m-1; i>=k; i--)
     1044  {
     1045    i1 = qrow[i];
     1046    if (i1 != i)
     1047    {
     1048      this->mpRowSwap(i1, i);
     1049      i2 = 0;
     1050      while (qrow[i2] != i) i2++;
     1051      qrow[i2] = i1;
     1052    }
     1053  }
     1054}
     1055
     1056/*
     1057* perform replacement for pivot strategy in Bareiss algorithm
     1058* change sign of determinant
     1059*/
     1060static void mpReplace(int j, int n, int &sign, int *perm)
     1061{
     1062  int k;
     1063
     1064  if (j != n)
     1065  {
     1066    k = perm[n];
     1067    perm[n] = perm[j];
     1068    perm[j] = k;
     1069    sign = -sign;
     1070  }
     1071}
     1072/*2
     1073* pivot strategy for Bareiss algorithm
     1074*/
     1075int mp_permmatrix::mpPivotBareiss(row_col_weight *C)
     1076{
     1077  poly p, *a;
     1078  int i, j, iopt, jopt;
     1079  float sum, f1, f2, fo, r, ro, lp;
     1080  float *dr = C->wrow, *dc = C->wcol;
     1081
     1082  fo = 1.0e20;
     1083  ro = 0.0;
     1084  iopt = jopt = -1;
     1085
     1086  s_n--;
     1087  s_m--;
     1088  if (s_m == 0)
     1089    return 0;
     1090  if (s_n == 0)
     1091  {
     1092    for(i=s_m; i>=0; i--)
     1093    {
     1094      p = this->mpRowAdr(i)[qcol[0]];
     1095      if (p)
     1096      {
     1097        f1 = mp_PolyWeight(p,_R);
     1098        if (f1 < fo)
     1099        {
     1100          fo = f1;
     1101          if (iopt >= 0)
     1102            p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]),_R);
     1103          iopt = i;
     1104        }
     1105        else
     1106          p_Delete(&(this->mpRowAdr(i)[qcol[0]]),_R);
     1107      }
     1108    }
     1109    if (iopt >= 0)
     1110      mpReplace(iopt, s_m, sign, qrow);
     1111    return 0;
     1112  }
     1113  this->mpRowWeight(dr);
     1114  this->mpColWeight(dc);
     1115  sum = 0.0;
     1116  for(i=s_m; i>=0; i--)
     1117    sum += dr[i];
     1118  for(i=s_m; i>=0; i--)
     1119  {
     1120    r = dr[i];
     1121    a = this->mpRowAdr(i);
     1122    for(j=s_n; j>=0; j--)
     1123    {
     1124      p = a[qcol[j]];
     1125      if (p)
     1126      {
     1127        lp = mp_PolyWeight(p,_R);
     1128        ro = r - lp;
     1129        f1 = ro * (dc[j]-lp);
     1130        if (f1 != 0.0)
     1131        {
     1132          f2 = lp * (sum - ro - dc[j]);
     1133          f2 += f1;
     1134        }
     1135        else
     1136          f2 = lp-r-dc[j];
     1137        if (f2 < fo)
     1138        {
     1139          fo = f2;
     1140          iopt = i;
     1141          jopt = j;
     1142        }
     1143      }
     1144    }
     1145  }
     1146  if (iopt < 0)
     1147    return 0;
     1148  mpReplace(iopt, s_m, sign, qrow);
     1149  mpReplace(jopt, s_n, sign, qcol);
     1150  return 1;
     1151}
     1152poly mp_permmatrix::mpGetElem(int r, int c)
     1153{
     1154  return Xarray[a_n*qrow[r]+qcol[c]];
     1155}
     1156
     1157/*
     1158* the Bareiss-type elimination with division by div (div != NULL)
     1159*/
     1160void mp_permmatrix::mpElimBareiss(poly div)
     1161{
     1162  poly piv, elim, q1, q2, *ap, *a;
     1163  int i, j, jj;
     1164
     1165  ap = this->mpRowAdr(s_m);
     1166  piv = ap[qcol[s_n]];
     1167  for(i=s_m-1; i>=0; i--)
     1168  {
     1169    a = this->mpRowAdr(i);
     1170    elim = a[qcol[s_n]];
     1171    if (elim != NULL)
     1172    {
     1173      elim = p_Neg(elim,_R);
     1174      for (j=s_n-1; j>=0; j--)
     1175      {
     1176        q2 = NULL;
     1177        jj = qcol[j];
     1178        if (ap[jj] != NULL)
     1179        {
     1180          q2 = SM_MULT(ap[jj], elim, div,_R);
     1181          if (a[jj] != NULL)
     1182          {
     1183            q1 = SM_MULT(a[jj], piv, div,_R);
     1184            p_Delete(&a[jj],_R);
     1185            q2 = p_Add_q(q2, q1, _R);
     1186          }
     1187        }
     1188        else if (a[jj] != NULL)
     1189        {
     1190          q2 = SM_MULT(a[jj], piv, div, _R);
     1191        }
     1192        if ((q2!=NULL) && div)
     1193          SM_DIV(q2, div, _R);
     1194        a[jj] = q2;
     1195      }
     1196      p_Delete(&a[qcol[s_n]], _R);
     1197    }
     1198    else
     1199    {
     1200      for (j=s_n-1; j>=0; j--)
     1201      {
     1202        jj = qcol[j];
     1203        if (a[jj] != NULL)
     1204        {
     1205          q2 = SM_MULT(a[jj], piv, div, _R);
     1206          p_Delete(&a[jj], _R);
     1207          if (div)
     1208            SM_DIV(q2, div, _R);
     1209          a[jj] = q2;
     1210        }
     1211      }
     1212    }
     1213  }
     1214}
    8611215/*
    8621216* weigth of a polynomial, for pivot strategy
    … …  
    12001554  }
    12011555  matrix c = mp_Copy(a,r);
    1202   mp_permmatrix *Bareiss = new mp_permmatrix(c);
     1556  mp_permmatrix *Bareiss = new mp_permmatrix(c,r);
    12031557  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
    12041558

  • libpolys/polys/sparsmat.cc

    r441a2e rbf6a4d  
    147147  smpoly piv, oldpiv;  // pivot and previous pivot
    148148  smpoly dumm;         // allocated dummy
    149   const ring _R;
     149  ring _R;
    150150 
    151151  void smColToRow();
    … …  
    261261
    262262/* ----------------- ops with rings ------------------ */
    263 ring sm_RingChange(const ring currRing, long bound)
     263ring sm_RingChange(const ring origR, long bound)
    264264{
    265265//  *origR =currRing;
    266   ring tmpR=rCopy0(currRing,FALSE,FALSE);
     266  ring tmpR=rCopy0(origR,FALSE,FALSE);
    267267  int *ord=(int*)omAlloc0(3*sizeof(int));
    268268  int *block0=(int*)omAlloc(3*sizeof(int));
    … …  
    283283
    284284  rComplete(tmpR,1);
    285   if ((*origR)->qideal!=NULL)
    286   {
    287     tmpR->qideal= idrCopyR_NoSort((*origR)->qideal, currRing, tmpR);
    288   }
    289 //  rChangeCurrRing(tmpR);
     285  if (origR->qideal!=NULL)
     286  {
     287    tmpR->qideal= idrCopyR_NoSort(origR->qideal, origR, tmpR);
     288  }
    290289  if (TEST_OPT_PROT)
    291290    Print("[%ld:%d]", (long) tmpR->bitmask, tmpR->ExpL_Size);
    … …  
    307306*        FALSE -> same Type
    308307*/
    309 BOOLEAN smCheckDet(ideal I, int d, BOOLEAN sw)
     308BOOLEAN sm_CheckDet(ideal I, int d, BOOLEAN sw, const ring r)
    310309{
    311310  int s,t,i;
    312311  poly p;
    313312
    314   if ((d>100) || (currRing->parameter!=NULL))
     313  if ((d>100) || (rIsExtension(r)))
    315314    return sw;
    316   if (rField_is_Q(currRing)==FALSE)
     315  if (!rField_is_Q(r))
    317316    return sw;
    318317  s = t = 0;
    … …  
    324323      if (p!=NULL)
    325324      {
    326         if(!pIsConstant(p))
     325        if(!p_IsConstant(p,r))
    327326          return sw;
    328327        s++;
    329         t+=nSize(pGetCoeff(p));
     328        t+=n_Size(pGetCoeff(p),r->cf);
    330329      }
    331330    }
    … …  
    336335    {
    337336      p=I->m[i];
    338       if (!pIsConstantPoly(p))
     337      if (!p_IsConstantPoly(p,r))
    339338        return sw;
    340339      while (p!=NULL)
    341340      {
    342341        s++;
    343         t+=nSize(pGetCoeff(p));
     342        t+=n_Size(pGetCoeff(p),r->cf);
    344343        pIter(p);
    345344      }
    … …  
    408407*uses Bareiss algorithm
    409408*/
    410 poly smCallDet(ideal I)
     409poly sm_CallDet(ideal I,const ring R)
    411410{
    412411  if (I->ncols != I->rank)
    … …  
    415414    return NULL;
    416415  }
    417   int r=idRankFreeModule(I);
     416  int r=id_RankFreeModule(I,R);
    418417  if (I->ncols != r) // some 0-lines at the end
    419418  {
    420419    return NULL;
    421420  }
    422   long bound=smExpBound(I,r,r,r);
    423   number diag,h=nInit(1);
     421  long bound=sm_ExpBound(I,r,r,r,R);
     422  number diag,h=n_Init(1,R->cf);
    424423  poly res;
    425   ring origR;
    426424  ring tmpR;
    427425  sparse_mat *det;
    428426  ideal II;
    429427
    430   tmpR=smRingChange(&origR,bound);
    431   II = idrCopyR(I, origR);
    432   diag = smCleardenom(II);
    433   det = new sparse_mat(II);
    434   idDelete(&II);
     428  tmpR=sm_RingChange(R,bound);
     429  II = idrCopyR(I, R, tmpR);
     430  diag = sm_Cleardenom(II,tmpR);
     431  det = new sparse_mat(II,tmpR);
     432  id_Delete(&II,tmpR);
    435433  if (det->smGetAct() == NULL)
    436434  {
    437435    delete det;
    438     rChangeCurrRing(origR);
    439     smKillModifiedRing(tmpR);
     436    sm_KillModifiedRing(tmpR);
    440437    return NULL;
    441438  }
    442439  res=det->smDet();
    443   if(det->smGetSign()<0) res=pNeg(res);
     440  if(det->smGetSign()<0) res=p_Neg(res,tmpR);
    444441  delete det;
    445   rChangeCurrRing(origR);
    446   res = prMoveR(res, tmpR);
    447   smKillModifiedRing(tmpR);
    448   if (nEqual(diag,h) == FALSE)
    449   {
    450     pMult_nn(res,diag);
    451     pNormalize(res);
    452   }
    453   nDelete(&diag);
    454   nDelete(&h);
     442  res = prMoveR(res, tmpR, R);
     443  sm_KillModifiedRing(tmpR);
     444  if (!n_Equal(diag,h,R->cf))
     445  {
     446    p_Mult_nn(res,diag,R);
     447    p_Normalize(res,R);
     448  }
     449  n_Delete(&diag,R->cf);
     450  n_Delete(&h,R->cf);
    455451  return res;
    456452}
    457453
    458 void smCallBareiss(ideal I, int x, int y, ideal & M, intvec **iv)
    459 {
    460   int r=idRankFreeModule(I),t=r;
     454void sm_CallBareiss(ideal I, int x, int y, ideal & M, intvec **iv, const ring R)
     455{
     456  int r=id_RankFreeModule(I,R),t=r;
    461457  int c=IDELEMS(I),s=c;
    462458  long bound;
    463   ring origR;
    464459  ring tmpR;
    465460  sparse_mat *bareiss;
    … …  
    470465    s-=y;
    471466  if (t>s) t=s;
    472   bound=smExpBound(I,c,r,t);
    473   tmpR=smRingChange(&origR,bound);
    474   ideal II = idrCopyR(I, origR);
    475   bareiss = new sparse_mat(II);
     467  bound=sm_ExpBound(I,c,r,t,R);
     468  tmpR=sm_RingChange(R,bound);
     469  ideal II = idrCopyR(I, R, tmpR);
     470  bareiss = new sparse_mat(II,tmpR);
    476471  if (bareiss->smGetAct() == NULL)
    477472  {
    478473    delete bareiss;
    479     *iv=new intvec(1,pVariables);
    480     rChangeCurrRing(origR);
     474    *iv=new intvec(1,rVar(tmpR));
    481475  }
    482476  else
    483477  {
    484     idDelete(&II);
     478    id_Delete(&II,tmpR);
    485479    bareiss->smNewBareiss(x, y);
    486480    II = bareiss->smRes2Mod();
    … …  
    488482    bareiss->smToIntvec(*iv);
    489483    delete bareiss;
    490     rChangeCurrRing(origR);
    491     II = idrMoveR(II,tmpR);
    492   }
    493   smKillModifiedRing(tmpR);
     484    II = idrMoveR(II,tmpR,R);
     485  }
     486  sm_KillModifiedRing(tmpR);
    494487  M=II;
    495488}
    … …  
    535528* constructor
    536529*/
    537 sparse_mat::sparse_mat(ideal smat)
     530sparse_mat::sparse_mat(ideal smat, const ring RR)
    538531{
    539532  int i;
    540   polyset pmat;
     533  poly* pmat;
     534  _R=RR;
    541535
    542536  ncols = smat->ncols;
    543   nrows = idRankFreeModule(smat);
     537  nrows = id_RankFreeModule(smat,RR);
    544538  if (nrows <= 0)
    545539  {
    … …  
    566560  for(i=ncols; i; i--)
    567561  {
    568     m_act[i] = smPoly2Smpoly(pmat[i-1]);
     562    m_act[i] = sm_Poly2Smpoly(pmat[i-1], RR);
    569563    pmat[i-1] = NULL;
    570564  }
    … …  
    637631* transform the result to a module
    638632*/
     633
    639634ideal sparse_mat::smRes2Mod()
    640635{
    … …  
    642637  int i;
    643638
    644   for (i=crd; i; i--) res->m[i-1] = smSmpoly2Poly(m_res[i]);
    645   res->rank = idRankFreeModule(res);
     639  for (i=crd; i; i--) res->m[i-1] = sm_Smpoly2Poly(m_res[i],_R);
     640  res->rank = id_RankFreeModule(res,_R);
    646641  return res;
    647642}
    … …  
    845840      if (a->pos > tored)
    846841        break;
    847       w = a->f = smPolyWeight(a);
     842      w = a->f = sm_PolyWeight(a,_R);
    848843      wc += w;
    849844      wrw[a->pos] += w;
    … …  
    10441039  if ((c == NULL) || (r == NULL))
    10451040  {
    1046     while (r!=NULL) smElemDelete(&r);
     1041    while (r!=NULL) sm_ElemDelete(&r,_R);
    10471042    return;
    10481043  }
    … …  
    10611056        {
    10621057          res = res->n = smElemCopy(b);
    1063           res->m = ppMult_qq(b->m, w);
     1058          res->m = pp_Mult_qq(b->m, w,_R);
    10641059          res->e = 1;
    1065           res->f = smPolyWeight(res);
     1060          res->f = sm_PolyWeight(res,_R);
    10661061          b = b->n;
    10671062        } while (b != NULL);
    … …  
    10761071      {
    10771072        res = res->n = smElemCopy(b);
    1078         res->m = ppMult_qq(b->m, w);
     1073        res->m = pp_Mult_qq(b->m, w,_R);
    10791074        res->e = 1;
    1080         res->f = smPolyWeight(res);
     1075        res->f = sm_PolyWeight(res,_R);
    10811076        b = b->n;
    10821077      }
    10831078      else
    10841079      {
    1085         ha = ppMult_qq(a->m, p);
    1086         pDelete(&a->m);
    1087         hb = ppMult_qq(b->m, w);
    1088         ha = pAdd(ha, hb);
     1080        ha = pp_Mult_qq(a->m, p,_R);
     1081        p_Delete(&a->m,_R);
     1082        hb = pp_Mult_qq(b->m, w,_R);
     1083        ha = p_Add_q(ha, hb,_R);
    10891084        if (ha != NULL)
    10901085        {
    10911086          a->m = ha;
    10921087          a->e = 1;
    1093           a->f = smPolyWeight(a);
     1088          a->f = sm_PolyWeight(a,_R);
    10941089          res = res->n = a;
    10951090          a = a->n;
    … …  
    10971092        else
    10981093        {
    1099           smElemDelete(&a);
     1094          sm_ElemDelete(&a,_R);
    11001095        }
    11011096        b = b->n;
    … …  
    11081103    }
    11091104    m_act[r->pos] = dumm->n;
    1110     smElemDelete(&r);
     1105    sm_ElemDelete(&r,_R);
    11111106  } while (r != NULL);
    11121107}
    … …  
    11251120  if ((c == NULL) || (r == NULL))
    11261121  {
    1127     while(r!=NULL) smElemDelete(&r);
    1128     pDelete(&hp);
     1122    while(r!=NULL) sm_ElemDelete(&r,_R);
     1123    p_Delete(&hp,_R);
    11291124    return;
    11301125  }
    … …  
    11461141        {
    11471142          res = res->n = smElemCopy(b);
    1148           x = SM_MULT(b->m, hr, m_res[ir]->m);
     1143          x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
    11491144          b = b->n;
    1150           if(ir) SM_DIV(x, m_res[ir]->m);
     1145          if(ir) SM_DIV(x, m_res[ir]->m,_R);
    11511146          res->m = x;
    11521147          res->e = e;
    1153           res->f = smPolyWeight(res);
     1148          res->f = sm_PolyWeight(res,_R);
    11541149        } while (b != NULL);
    11551150        break;
    … …  
    11631158      {
    11641159        res = res->n = smElemCopy(b);
    1165         x = SM_MULT(b->m, hr, m_res[ir]->m);
     1160        x = SM_MULT(b->m, hr, m_res[ir]->m,_R);
    11661161        b = b->n;
    1167         if(ir) SM_DIV(x, m_res[ir]->m);
     1162        if(ir) SM_DIV(x, m_res[ir]->m,_R);
    11681163        res->m = x;
    11691164        res->e = e;
    1170         res->f = smPolyWeight(res);
     1165        res->f = sm_PolyWeight(res,_R);
    11711166      }
    11721167      else
    … …  
    11781173          if (ir > ia)
    11791174          {
    1180             x = SM_MULT(ha, m_res[ir]->m, m_res[ia]->m);
    1181             pDelete(&ha);
     1175            x = SM_MULT(ha, m_res[ir]->m, m_res[ia]->m,_R);
     1176            p_Delete(&ha,_R);
    11821177            ha = x;
    1183             if (ia) SM_DIV(ha, m_res[ia]->m);
     1178            if (ia) SM_DIV(ha, m_res[ia]->m,_R);
    11841179            ia = ir;
    11851180          }
    1186           x = SM_MULT(ha, gp, m_res[ia]->m);
    1187           pDelete(&ha);
    1188           y = SM_MULT(b->m, hr, m_res[ia]->m);
     1181          x = SM_MULT(ha, gp, m_res[ia]->m,_R);
     1182          p_Delete(&ha,_R);
     1183          y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
    11891184        }
    11901185        else if (ir >= ip)
    … …  
    11921187          if (ia < crd)
    11931188          {
    1194             x = SM_MULT(ha, m_res[crd]->m, m_res[ia]->m);
    1195             pDelete(&ha);
     1189            x = SM_MULT(ha, m_res[crd]->m, m_res[ia]->m,_R);
     1190            p_Delete(&ha,_R);
    11961191            ha = x;
    1197             SM_DIV(ha, m_res[ia]->m);
     1192            SM_DIV(ha, m_res[ia]->m,_R);
    11981193          }
    11991194          y = hp;
    12001195          if(ir > ip)
    12011196          {
    1202             y = SM_MULT(y, m_res[ir]->m, m_res[ip]->m);
    1203             if (ip) SM_DIV(y, m_res[ip]->m);
     1197            y = SM_MULT(y, m_res[ir]->m, m_res[ip]->m,_R);
     1198            if (ip) SM_DIV(y, m_res[ip]->m,_R);
    12041199          }
    12051200          ia = ir;
    1206           x = SM_MULT(ha, y, m_res[ia]->m);
    1207           if (y != hp) pDelete(&y);
    1208           pDelete(&ha);
    1209           y = SM_MULT(b->m, hr, m_res[ia]->m);
     1201          x = SM_MULT(ha, y, m_res[ia]->m,_R);
     1202          if (y != hp) p_Delete(&y,_R);
     1203          p_Delete(&ha,_R);
     1204          y = SM_MULT(b->m, hr, m_res[ia]->m,_R);
    12101205        }
    12111206        else
    12121207        {
    1213           x = SM_MULT(hr, m_res[ia]->m, m_res[ir]->m);
    1214           if (ir) SM_DIV(x, m_res[ir]->m);
    1215           y = SM_MULT(b->m, x, m_res[ia]->m);
    1216           pDelete(&x);
    1217           x = SM_MULT(ha, gp, m_res[ia]->m);
    1218           pDelete(&ha);
     1208          x = SM_MULT(hr, m_res[ia]->m, m_res[ir]->m,_R);
     1209          if (ir) SM_DIV(x, m_res[ir]->m,_R);
     1210          y = SM_MULT(b->m, x, m_res[ia]->m,_R);
     1211          p_Delete(&x,_R);
     1212          x = SM_MULT(ha, gp, m_res[ia]->m,_R);
     1213          p_Delete(&ha,_R);
    12191214        }
    1220         ha = pAdd(x, y);
     1215        ha = p_Add_q(x, y,_R);
    12211216        if (ha != NULL)
    12221217        {
    1223           if (ia) SM_DIV(ha, m_res[ia]->m);
     1218          if (ia) SM_DIV(ha, m_res[ia]->m,_R);
    12241219          a->m = ha;
    12251220          a->e = e;
    1226           a->f = smPolyWeight(a);
     1221          a->f = sm_PolyWeight(a,_R);
    12271222          res = res->n = a;
    12281223          a = a->n;
    … …  
    12311226        {
    12321227          a->m = NULL;
    1233           smElemDelete(&a);
     1228          sm_ElemDelete(&a,_R);
    12341229        }
    12351230        b = b->n;
    … …  
    12421237    }
    12431238    m_act[r->pos] = dumm->n;
    1244     smElemDelete(&r);
     1239    sm_ElemDelete(&r,_R);
    12451240  } while (r != NULL);
    1246   pDelete(&hp);
     1241  p_Delete(&hp,_R);
    12471242}
    12481243
    … …  
    12931288        {
    12941289          ap->n = a->n;
    1295           a->m = pNeg(a->m);
     1290          a->m = p_Neg(a->m,_R);
    12961291          b = b->n = a;
    12971292          b->pos = i;
    … …  
    13031298    {
    13041299      m_act[i] = a->n;
    1305       a->m = pNeg(a->m);
     1300      a->m = p_Neg(a->m,_R);
    13061301      b = b->n = a;
    13071302      b->pos = i;
    … …  
    16061601    if (f < e)
    16071602    {
    1608       ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m);
    1609       pDelete(&a->m);
    1610       if (f) SM_DIV(ha, m_res[f]->m);
     1603      ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m,_R);
     1604      p_Delete(&a->m,_R);
     1605      if (f) SM_DIV(ha, m_res[f]->m,_R);
    16111606      a->m = ha;
    1612       if (normalize) pNormalize(a->m);
     1607      if (normalize) p_Normalize(a->m,_R);
    16131608    }
    16141609    a = a->n;
    … …  
    16341629      if (f < e)
    16351630      {
    1636         ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m);
    1637         pDelete(&a->m);
    1638         if (f) SM_DIV(ha, m_res[f]->m);
     1631        ha = SM_MULT(a->m, m_res[e]->m, m_res[f]->m, _R);
     1632        p_Delete(&a->m,_R);
     1633        if (f) SM_DIV(ha, m_res[f]->m, _R);
    16391634        a->m = ha;
    16401635      }
    1641       if (normalize) pNormalize(a->m);
     1636      if (normalize) p_Normalize(a->m, _R);
    16421637      a = a->n;
    16431638    } while (a != NULL);
    … …  
    16581653    do
    16591654    {
    1660       if(smHaveDenom(a->m)) return 1;
     1655      if(sm_HaveDenom(a->m,_R)) return 1;
    16611656      a = a->n;
    16621657    } while (a != NULL);
    … …  
    16791674    do
    16801675    {
    1681       if (e == a->e) pNormalize(a->m);
     1676      if (e == a->e) p_Normalize(a->m,_R);
    16821677      a = a->n;
    16831678    } while (a != NULL);
    … …  
    16961691  {
    16971692    h = r = a->m;
    1698     h = SM_MULT(h, m_res[crd]->m, m_res[f]->m);
    1699     if (f) SM_DIV(h, m_res[f]->m);
     1693    h = SM_MULT(h, m_res[crd]->m, m_res[f]->m, _R);
     1694    if (f) SM_DIV(h, m_res[f]->m, _R);
    17001695    a->m = h;
    1701     if (normalize) pNormalize(a->m);
    1702     a->f = smPolyWeight(a);
     1696    if (normalize) p_Normalize(a->m,_R);
     1697    a->f = sm_PolyWeight(a,_R);
    17031698    return r;
    17041699  }
    … …  
    17201715    do
    17211716    {
    1722       smElemDelete(&a);
     1717      sm_ElemDelete(&a,_R);
    17231718    } while (a != NULL);
    17241719  }
    … …  
    17341729  while (a != NULL)
    17351730  {
    1736     smElemDelete(&a);
     1731    sm_ElemDelete(&a,_R);
    17371732  }
    17381733}
    … …  
    17471742  while (i != 0)
    17481743  {
    1749     smElemDelete(&m_res[i]);
     1744    sm_ElemDelete(&m_res[i],_R);
    17501745    i--;
    17511746  }
    … …  
    17941789* a destroyed, b NOT destroyed
    17951790*/
    1796 void smPolyDiv(poly a, poly b)
     1791void sm_PolyDiv(poly a, poly b, const ring R)
    17971792{
    17981793  const number x = pGetCoeff(b);
    … …  
    18051800    do
    18061801    {
    1807       if (!pLmIsConstantComp(b))
    1808       {
    1809         for (i=pVariables; i; i--)
    1810           pSubExp(a,i,pGetExp(b,i));
    1811         pSetm(a);
    1812       }
    1813       y = nDiv(pGetCoeff(a),x);
    1814       nNormalize(y);
    1815       pSetCoeff(a,y);
     1802      if (!p_LmIsConstantComp(b,R))
     1803      {
     1804        for (i=rVar(R); i; i--)
     1805          p_SubExp(a,i,p_GetExp(b,i,R),R);
     1806        p_Setm(a,R);
     1807      }
     1808      y = n_Div(pGetCoeff(a),x,R->cf);
     1809      n_Normalize(y,R->cf);
     1810      p_SetCoeff(a,y,R);
    18161811      pIter(a);
    18171812    } while (a != NULL);
    18181813    return;
    18191814  }
    1820   dummy = pInit();
     1815  dummy = p_Init(R);
    18211816  do
    18221817  {
    1823     for (i=pVariables; i; i--)
    1824       pSubExp(a,i,pGetExp(b,i));
    1825     pSetm(a);
    1826     y = nDiv(pGetCoeff(a),x);
    1827     nNormalize(y);
    1828     pSetCoeff(a,y);
    1829     yn = nNeg(nCopy(y));
     1818    for (i=rVar(R); i; i--)
     1819      p_SubExp(a,i,p_GetExp(b,i,R),R);
     1820    p_Setm(a,R);
     1821    y = n_Div(pGetCoeff(a),x,R->cf);
     1822    n_Normalize(y,R->cf);
     1823    p_SetCoeff(a,y,R);
     1824    yn = n_Neg(n_Copy(y,R->cf),R->cf);
    18301825    t = pNext(b);
    18311826    h = dummy;
    18321827    do
    18331828    {
    1834       h = pNext(h) = pInit();
     1829      h = pNext(h) = p_Init(R);
    18351830      //pSetComp(h,0);
    1836       for (i=pVariables; i; i--)
    1837         pSetExp(h,i,pGetExp(a,i)+pGetExp(t,i));
    1838       pSetm(h);
    1839       pSetCoeff0(h,nMult(yn, pGetCoeff(t)));
     1831      for (i=rVar(R); i; i--)
     1832        p_SetExp(h,i,p_GetExp(a,i,R)+p_GetExp(t,i,R),R);
     1833      p_Setm(h,R);
     1834      pSetCoeff0(h,n_Mult(yn, pGetCoeff(t),R->cf));
    18401835      pIter(t);
    18411836    } while (t != NULL);
    1842     nDelete(&yn);
     1837    n_Delete(&yn,R->cf);
    18431838    pNext(h) = NULL;
    1844     a = pNext(a) = pAdd(pNext(a), pNext(dummy));
     1839    a = pNext(a) = p_Add_q(pNext(a), pNext(dummy),R);
    18451840  } while (a!=NULL);
    1846   pLmFree(dummy);
     1841  p_LmFree(dummy, R);
    18471842}
    18481843
    … …  
    20162011  }
    20172012  res = NULL;
    2018   e = pInit();
     2013  e = p_Init();
    20192014  lead = FALSE;
    20202015  while (!lead)
    … …  
    20292024    {
    20302025      lead = TRUE;
    2031       r = ppMult_mm(a, e);
     2026      r = pp_Mult__mm(a, e);
    20322027    }
    20332028    if (lead)
    … …  
    20452040    }
    20462041    else
    2047       res = pAdd(res, r);
     2042      res = p_Add_q(res, r);
    20482043    pIter(b);
    20492044    if (b == NULL)
    … …  
    20902085      else
    20912086      {
    2092         r = ppMult_mm(a, e);
     2087        r = pp_Mult__mm(a, e);
    20932088        append = kBucket_ExtractLarger_Add_q(bucket, append, r, &la);
    20942089      }
    … …  
    21112106          smCombineChain(&pa, r);
    21122107        else
    2113           pa = pAdd(pa,r);
     2108          pa = p_Add_q(pa,r);
    21142109      }
    21152110      else
    21162111      {
    2117         r = ppMult_mm(a, e);
     2112        r = pp_Mult__mm(a, e);
    21182113        smCombineChain(&pa, r);
    21192114      }
    … …  
    21302125*  returns the part of (a*b)/exp(lead(c)) with nonegative exponents
    21312126*/
    2132 poly smMultDiv(poly a, poly b, const poly c)
     2127poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
    21332128{
    21342129  poly pa, e, res, r;
    21352130  BOOLEAN lead;
    21362131
    2137   if ((c == NULL) || pLmIsConstantComp(c))
    2138   {
    2139     return ppMult_qq(a, b);
     2132  if ((c == NULL) || p_LmIsConstantComp(c,R))
     2133  {
     2134    return pp_Mult_qq(a, b, R);
    21402135  }
    21412136  if (smSmaller(a, b))
    … …  
    21462141  }
    21472142  res = NULL;
    2148   e = pInit();
     2143  e = p_Init(R);
    21492144  lead = FALSE;
    21502145  while (!lead)
    21512146  {
    21522147    pSetCoeff0(e,pGetCoeff(b));
    2153     if (smIsNegQuot(e, b, c))
    2154     {
    2155       lead = pLmDivisibleByNoComp(e, a);
    2156       r = smSelectCopy_ExpMultDiv(a, e, b, c);
     2148    if (sm_IsNegQuot(e, b, c, R))
     2149    {
     2150      lead = p_LmDivisibleByNoComp(e, a, R);
     2151      r = sm_SelectCopy_ExpMultDiv(a, e, b, c, R);
    21572152    }
    21582153    else
    21592154    {
    21602155      lead = TRUE;
    2161       r = ppMult_mm(a, e);
     2156      r = pp_Mult_mm(a, e,R);
    21622157    }
    21632158    if (lead)
    … …  
    21652160      if (res != NULL)
    21662161      {
    2167         smFindRef(&pa, &res, r);
     2162        sm_FindRef(&pa, &res, r, R);
    21682163        if (pa == NULL)
    21692164          lead = FALSE;
    … …  
    21752170    }
    21762171    else
    2177       res = pAdd(res, r);
     2172      res = p_Add_q(res, r, R);
    21782173    pIter(b);
    21792174    if (b == NULL)
    21802175    {
    2181       pLmFree(e);
     2176      p_LmFree(e, R);
    21822177      return res;
    21832178    }
    … …  
    21862181  {
    21872182    pSetCoeff0(e,pGetCoeff(b));
    2188     if (smIsNegQuot(e, b, c))
    2189     {
    2190       r = smSelectCopy_ExpMultDiv(a, e, b, c);
    2191       if (pLmDivisibleByNoComp(e, a))
    2192         smCombineChain(&pa, r);
     2183    if (sm_IsNegQuot(e, b, c, R))
     2184    {
     2185      r = sm_SelectCopy_ExpMultDiv(a, e, b, c, R);
     2186      if (p_LmDivisibleByNoComp(e, a,R))
     2187        sm_CombineChain(&pa, r, R);
    21932188      else
    2194         pa = pAdd(pa,r);
     2189        pa = p_Add_q(pa,r,R);
    21952190    }
    21962191    else
    21972192    {
    2198       r = ppMult_mm(a, e);
    2199       smCombineChain(&pa, r);
     2193      r = pp_Mult_mm(a, e, R);
     2194      sm_CombineChain(&pa, r, R);
    22002195    }
    22012196    pIter(b);
    22022197  } while (b != NULL);
    2203   pLmFree(e);
     2198  p_LmFree(e, R);
    22042199  return res;
    22052200}
    22062201#endif
    2207 /* ------------ internals arithmetic ------------- */
    2208 static void smExactPolyDiv(poly a, poly b)
    2209 {
    2210   const number x = pGetCoeff(b);
    2211   poly tail = pNext(b), e = pInit();
    2212   poly h;
    2213   number y, yn;
    2214   int lt = pLength(tail);
    2215 
    2216   if (lt + 1 >= SM_MIN_LENGTH_BUCKET &&  !TEST_OPT_NOT_BUCKETS)
    2217   {
    2218     kBucket_pt bucket = kBucketCreate(currRing);
    2219     kBucketInit(bucket, pNext(a), 0);
    2220     int lh = 0;
    2221     do
    2222     {
    2223       y = nDiv(pGetCoeff(a), x);
    2224       nNormalize(y);
    2225       pSetCoeff(a,y);
    2226       yn = nNeg(nCopy(y));
    2227       pSetCoeff0(e,yn);
    2228       lh = lt;
    2229       if (smIsNegQuot(e, a, b))
    2230       {
    2231         h = pp_Mult_Coeff_mm_DivSelect_MultDiv(tail, lh, e, a, b);
    2232       }
    2233       else
    2234         h = ppMult_mm(tail, e);
    2235       nDelete(&yn);
    2236       kBucket_Add_q(bucket, h, &lh);
    2237 
    2238       a = pNext(a) = kBucketExtractLm(bucket);
    2239     } while (a!=NULL);
    2240     kBucketDestroy(&bucket);
    2241   }
    2242   else
    2243   {
    2244     do
    2245     {
    2246       y = nDiv(pGetCoeff(a), x);
    2247       nNormalize(y);
    2248       pSetCoeff(a,y);
    2249       yn = nNeg(nCopy(y));
    2250       pSetCoeff0(e,yn);
    2251       if (smIsNegQuot(e, a, b))
    2252         h = smSelectCopy_ExpMultDiv(tail, e, a, b);
    2253       else
    2254         h = ppMult_mm(tail, e);
    2255       nDelete(&yn);
    2256       a = pNext(a) = pAdd(pNext(a), h);
    2257     } while (a!=NULL);
    2258   }
    2259   pLmFree(e);
    2260 }
    2261 
    2262 static void smPolyDivN(poly a, const number x)
    2263 {
    2264   number y;
    2265 
    2266   do
    2267   {
    2268     y = nDiv(pGetCoeff(a),x);
    2269     nNormalize(y);
    2270     pSetCoeff(a,y);
    2271     pIter(a);
    2272   } while (a != NULL);
    2273 }
    2274 
    22752202/*n
    22762203* exact division a/b
    … …  
    22782205* a destroyed, b NOT destroyed
    22792206*/
    2280 void smSpecialPolyDiv(poly a, poly b)
     2207void sm_SpecialPolyDiv(poly a, poly b, const ring R)
    22812208{
    22822209  if (pNext(b) == NULL)
    22832210  {
    2284     smPolyDivN(a, pGetCoeff(b));
     2211    sm_PolyDivN(a, pGetCoeff(b),R);
    22852212    return;
     2213  }
     2214  sm_ExactPolyDiv(a, b, R);
     2215}
     2216
     2217
     2218/* ------------ internals arithmetic ------------- */
     2219static void sm_ExactPolyDiv(poly a, poly b, const ring R)
     2220{
     2221  const number x = pGetCoeff(b);
     2222  poly tail = pNext(b), e = p_Init(R);
     2223  poly h;
     2224  number y, yn;
     2225  int lt = pLength(tail);
     2226
     2227  if (lt + 1 >= SM_MIN_LENGTH_BUCKET &&  !TEST_OPT_NOT_BUCKETS)
     2228  {
     2229    kBucket_pt bucket = kBucketCreate(R);
     2230    kBucketInit(bucket, pNext(a), 0);
     2231    int lh = 0;
     2232    do
     2233    {
     2234      y = n_Div(pGetCoeff(a), x, R->cf);
     2235      n_Normalize(y, R->cf);
     2236      p_SetCoeff(a,y, R);
     2237      yn = n_Neg(n_Copy(y, R->cf), R->cf);
     2238      pSetCoeff0(e,yn);
     2239      lh = lt;
     2240      if (sm_IsNegQuot(e, a, b, R))
     2241      {
     2242        h = pp_Mult_Coeff_mm_DivSelect_MultDiv(tail, lh, e, a, b, R);
     2243      }
     2244      else
     2245        h = pp_Mult_mm(tail, e, R);
     2246      n_Delete(&yn, R->cf);
     2247      kBucket_Add_q(bucket, h, &lh);
     2248
     2249      a = pNext(a) = kBucketExtractLm(bucket);
     2250    } while (a!=NULL);
     2251    kBucketDestroy(&bucket);
     2252  }
     2253  else
     2254  {
     2255    do
     2256    {
     2257      y = n_Div(pGetCoeff(a), x, R->cf);
     2258      n_Normalize(y, R->cf);
     2259      p_SetCoeff(a,y, R);
     2260      yn = n_Neg(n_Copy(y, R->cf), R->cf);
     2261      pSetCoeff0(e,yn);
     2262      if (sm_IsNegQuot(e, a, b, R))
     2263        h = sm_SelectCopy_ExpMultDiv(tail, e, a, b, R);
     2264      else
     2265        h = pp_Mult_mm(tail, e, R);
     2266      n_Delete(&yn, R->cf);
     2267      a = pNext(a) = p_Add_q(pNext(a), h, R);
     2268    } while (a!=NULL);
     2269  }
     2270  p_LmFree(e, R);
     2271}
     2272
     2273// obachman --> Wilfried: check the following
     2274static BOOLEAN sm_IsNegQuot(poly a, const poly b, const poly c, const ring R)
     2275{
     2276  if (p_LmDivisibleByNoComp(c, b,R))
     2277  {
     2278    p_ExpVectorDiff(a, b, c,R);
     2279    // Hmm: here used to be a p_Setm(a): but it is unnecessary,
     2280    // if b and c are correct
     2281    return FALSE;
     2282  }
     2283  else
     2284  {
     2285    int i;
     2286    for (i=rVar(R); i>0; i--)
     2287    {
     2288      if(p_GetExp(c,i,R) > p_GetExp(b,i,R))
     2289        p_SetExp(a,i,p_GetExp(c,i,R)-p_GetExp(b,i,R),R);
     2290      else
     2291        p_SetExp(a,i,0,R);
     2292    }
     2293    // here we actually might need a p_Setm, if a is to be used in
     2294    // comparisons
     2295    return TRUE;
     2296  }
     2297}
     2298
     2299static void sm_ExpMultDiv(poly t, const poly b, const poly c, const ring R)
     2300{
     2301  int i;
     2302  p_Test(t,R);
     2303  p_LmTest(b,R);
     2304  p_LmTest(c,R);
     2305  poly bc = p_New(R);
     2306
     2307  p_ExpVectorDiff(bc, b, c, R);
     2308
     2309  while(t!=NULL)
     2310  {
     2311    p_ExpVectorAdd(t, bc, R);
     2312    pIter(t);
     2313  }
     2314  p_LmFree(bc, R);
     2315}
     2316
     2317
     2318static void sm_PolyDivN(poly a, const number x, const ring R)
     2319{
     2320  number y;
     2321
     2322  do
     2323  {
     2324    y = n_Div(pGetCoeff(a),x, R->cf);
     2325    n_Normalize(y, R->cf);
     2326    p_SetCoeff(a,y, R);
     2327    pIter(a);
     2328  } while (a != NULL);
     2329}
     2330
     2331static BOOLEAN smSmaller(poly a, poly b)
     2332{
     2333  loop
     2334  {
     2335    pIter(b);
     2336    if (b == NULL) return TRUE;
     2337    pIter(a);
     2338    if (a == NULL) return FALSE;
     2339  }
     2340}
     2341
     2342static void sm_CombineChain(poly *px, poly r, const ring R)
     2343{
     2344  poly pa = *px, pb;
     2345  number x;
     2346  int i;
     2347
     2348  loop
     2349  {
     2350    pb = pNext(pa);
     2351    if (pb == NULL)
     2352    {
     2353      pa = pNext(pa) = r;
     2354      break;
     2355    }
     2356    i = p_LmCmp(pb, r,R);
     2357    if (i > 0)
     2358      pa = pb;
     2359    else
     2360    {
     2361      if (i == 0)
     2362      {
     2363        x = n_Add(pGetCoeff(pb), pGetCoeff(r),R->cf);
     2364        p_LmDelete(&r,R);
     2365        if (n_IsZero(x,R->cf))
     2366        {
     2367          p_LmDelete(&pb,R);
     2368          pNext(pa) = p_Add_q(pb,r,R);
     2369        }
     2370        else
     2371        {
     2372          pa = pb;
     2373          p_SetCoeff(pa,x,R);
     2374          pNext(pa) = p_Add_q(pNext(pa), r, R);
     2375        }
     2376      }
     2377      else
     2378      {
     2379        pa = pNext(pa) = r;
     2380        pNext(pa) = p_Add_q(pb, pNext(pa),R);
     2381      }
     2382      break;
     2383    }
     2384  }
     2385  *px = pa;
     2386}
     2387
     2388
     2389static void sm_FindRef(poly *ref, poly *px, poly r, const ring R)
     2390{
     2391  number x;
     2392  int i;
     2393  poly pa = *px, pp = NULL;
     2394
     2395  loop
     2396  {
     2397    i = p_LmCmp(pa, r,R);
     2398    if (i > 0)
     2399    {
     2400      pp = pa;
     2401      pIter(pa);
     2402      if (pa==NULL)
     2403      {
     2404        pNext(pp) = r;
     2405        break;
     2406      }
     2407    }
     2408    else
     2409    {
     2410      if (i == 0)
     2411      {
     2412        x = n_Add(pGetCoeff(pa), pGetCoeff(r),R->cf);
     2413        p_LmDelete(&r,R);
     2414        if (n_IsZero(x,R->cf))
     2415        {
     2416          p_LmDelete(&pa,R);
     2417          if (pp!=NULL)
     2418            pNext(pp) = p_Add_q(pa,r,R);
     2419          else
     2420            *px = p_Add_q(pa,r,R);
     2421        }
     2422        else
     2423        {
     2424          pp = pa;
     2425          p_SetCoeff(pp,x,R);
     2426          pNext(pp) = p_Add_q(pNext(pp), r, R);
     2427        }
     2428      }
     2429      else
     2430      {
     2431        if (pp!=NULL)
     2432          pp = pNext(pp) = r;
     2433        else
     2434          *px = pp = r;
     2435        pNext(pp) = p_Add_q(pa, pNext(r),R);
     2436      }
     2437      break;
     2438    }
     2439  }
     2440  *ref = pp;
     2441}
     2442
     2443/* ----------------- internal 'C' stuff ------------------ */
     2444
     2445static void sm_ElemDelete(smpoly *r, const ring R)
     2446{
     2447  smpoly a = *r, b = a->n;
     2448
     2449  p_Delete(&a->m, R);
     2450  omFreeBin((void *)a,  smprec_bin);
     2451  *r = b;
     2452}
     2453
     2454static smpoly smElemCopy(smpoly a)
     2455{
     2456  smpoly r = (smpoly)omAllocBin(smprec_bin);
     2457  memcpy(r, a, sizeof(smprec));
     2458/*  r->m = pCopy(r->m); */
     2459  return r;
     2460}
     2461
     2462/*
     2463* from poly to smpoly
     2464* do not destroy p
     2465*/
     2466static smpoly sm_Poly2Smpoly(poly q, const ring R)
     2467{
     2468  poly pp;
     2469  smpoly res, a;
     2470  long x;
     2471
     2472  if (q == NULL)
     2473    return NULL;
     2474  a = res = (smpoly)omAllocBin(smprec_bin);
     2475  a->pos = x = p_GetComp(q,R);
     2476  a->m = q;
     2477  a->e = 0;
     2478  loop
     2479  {
     2480    p_SetComp(q,0,R);
     2481    pp = q;
     2482    pIter(q);
     2483    if (q == NULL)
     2484    {
     2485      a->n = NULL;
     2486      return res;
     2487    }
     2488    if (p_GetComp(q,R) != x)
     2489    {
     2490      a = a->n = (smpoly)omAllocBin(smprec_bin);
     2491      pNext(pp) = NULL;
     2492      a->pos = x = p_GetComp(q,R);
     2493      a->m = q;
     2494      a->e = 0;
     2495    }
     2496  }
     2497}
     2498
     2499/*
     2500* from smpoly to poly
     2501* destroy a
     2502*/
     2503static poly sm_Smpoly2Poly(smpoly a, const ring R)
     2504{
     2505  smpoly b;
     2506  poly res, pp, q;
     2507  long x;
     2508
     2509  if (a == NULL)
     2510    return NULL;
     2511  x = a->pos;
     2512  q = res = a->m;
     2513  loop
     2514  {
     2515    p_SetComp(q,x,R);
     2516    pp = q;
     2517    pIter(q);
     2518    if (q == NULL)
     2519      break;
     2520  }
     2521  loop
     2522  {
     2523    b = a;
     2524    a = a->n;
     2525    omFreeBin((void *)b,  smprec_bin);
     2526    if (a == NULL)
     2527      return res;
     2528    x = a->pos;
     2529    q = pNext(pp) = a->m;
     2530    loop
     2531    {
     2532      p_SetComp(q,x,R);
     2533      pp = q;
     2534      pIter(q);
     2535      if (q == NULL)
     2536        break;
     2537    }
     2538  }
     2539}
     2540
     2541/*
     2542* weigth of a polynomial, for pivot strategy
     2543*/
     2544static float sm_PolyWeight(smpoly a, const ring R)
     2545{
     2546  poly p = a->m;
     2547  int i;
     2548  float res = (float)n_Size(pGetCoeff(p),R->cf);
     2549
     2550  if (pNext(p) == NULL)
     2551  {
     2552    for(i=rVar(R); i>0; i--)
     2553    {
     2554      if (p_GetExp(p,i,R) != 0) return res+1.0;
     2555    }
     2556    return res;
     2557  }
     2558  else
     2559  {
     2560    i = 0;
     2561    res = 0.0;
     2562    do
     2563    {
     2564      i++;
     2565      res += (float)n_Size(pGetCoeff(p),R->cf);
     2566      pIter(p);
     2567    }
     2568    while (p);
     2569    return res+(float)i;
     2570  }
     2571}
     2572
     2573static BOOLEAN sm_HaveDenom(poly a, const ring R)
     2574{
     2575  BOOLEAN sw;
     2576  number x;
     2577
     2578  while (a != NULL)
     2579  {
     2580    x = n_GetDenom(pGetCoeff(a),R->cf);
     2581    sw = n_IsOne(x,R->cf);
     2582    n_Delete(&x,R->cf);
     2583    if (!sw)
     2584    {
     2585      return TRUE;
     2586    }
     2587    pIter(a);
     2588  }
     2589  return FALSE;
     2590}
     2591
     2592static number sm_Cleardenom(ideal id, const ring R)
     2593{
     2594  poly a;
     2595  number x,y,res=n_Init(1,R->cf);
     2596  BOOLEAN sw=FALSE;
     2597
     2598  for (int i=0; i<IDELEMS(id); i++)
     2599  {
     2600    a = id->m[i];
     2601    sw = sm_HaveDenom(a,R);
     2602    if (sw) break;
     2603  }
     2604  if (!sw) return res;
     2605  for (int i=0; i<IDELEMS(id); i++)
     2606  {
     2607    a = id->m[i];
     2608    if (a!=NULL)
     2609    {
     2610      x = n_Copy(pGetCoeff(a),R->cf);
     2611      p_Cleardenom(a, R);
     2612      y = n_Div(x,pGetCoeff(a),R->cf);
     2613      n_Delete(&x,R->cf);
     2614      x = n_Mult(res,y,R->cf);
     2615      n_Normalize(x,R->cf);
     2616      n_Delete(&res,R->cf);
     2617      res = x;
     2618    }
    22862619  }
    22872620  smExactPolyDiv(a, b);
    … …  
    23012634
    23022635/* declare internal 'C' stuff */
    2303 static void smNumberDelete(smnumber *);
     2636static void sm_NumberDelete(smnumber *, const ring R);
    23042637static smnumber smNumberCopy(smnumber);
    2305 static smnumber smPoly2Smnumber(poly);
    2306 static poly smSmnumber2Poly(number);
     2638static smnumber sm_Poly2Smnumber(poly, const ring);
     2639static poly sm_Smnumber2Poly(number,const ring);
    23072640static BOOLEAN smCheckSolv(ideal);
    23082641
    … …  
    23262659  smnumber piv;        // pivot
    23272660  smnumber dumm;       // allocated dummy
     2661  ring _R;
    23282662  void smColToRow();
    23292663  void smRowToCol();
    … …  
    23342668  void smAllDel();
    23352669public:
    2336   sparse_number_mat(ideal);
     2670  sparse_number_mat(ideal, const ring);
    23372671  ~sparse_number_mat();
    23382672  int smIsSing() { return sing; }
    … …  
    23492683* uses  Gauss-elimination
    23502684*/
    2351 ideal smCallSolv(ideal I)
     2685ideal sm_CallSolv(ideal I, const ring R)
    23522686{
    23532687  sparse_number_mat *linsolv;
    … …  
    23562690  ideal rr;
    23572691
    2358   if (idIsConstant(I)==FALSE)
     2692  if (id_IsConstant(I,R)==FALSE)
    23592693  {
    23602694    WerrorS("symbol in equation");
    23612695    return NULL;
    23622696  }
    2363   I->rank = idRankFreeModule(I);
     2697  I->rank = id_RankFreeModule(I,R);
    23642698  if (smCheckSolv(I)) return NULL;
    2365   tmpR=smRingChange(&origR,1);
    2366   rr=idrCopyR(I,origR);
    2367   linsolv = new sparse_number_mat(rr);
     2699  tmpR=sm_RingChange(R,1);
     2700  rr=idrCopyR(I,R, tmpR);
     2701  linsolv = new sparse_number_mat(rr,tmpR);
    23682702  rr=NULL;
    23692703  linsolv->smTriangular();
    … …  
    23762710    WerrorS("singular problem for linsolv");
    23772711  delete linsolv;
    2378   rChangeCurrRing(origR);
    23792712  if (rr!=NULL)
    2380     rr = idrMoveR(rr,tmpR);
    2381   smKillModifiedRing(tmpR);
     2713    rr = idrMoveR(rr,tmpR,R);
     2714  sm_KillModifiedRing(tmpR);
    23822715  return rr;
    23832716}
    … …  
    23862719* constructor, destroy smat
    23872720*/
    2388 sparse_number_mat::sparse_number_mat(ideal smat)
     2721sparse_number_mat::sparse_number_mat(ideal smat, const ring R)
    23892722{
    23902723  int i;
    2391   polyset pmat;
     2724  poly* pmat;
     2725  _R=R;
    23922726
    23932727  crd = sing = 0;
    … …  
    24062740  for(i=ncols; i; i--)
    24072741  {
    2408     m_act[i] = smPoly2Smnumber(pmat[i-1]);
     2742    m_act[i] = sm_Poly2Smnumber(pmat[i-1],_R);
    24092743  }
    24102744  omFreeSize((ADDRESS)pmat,smat->ncols*sizeof(poly));
    … …  
    24822816  {
    24832817    x = sol[i];
    2484     sol[i] = nDiv(x, m_res[i]->m);
    2485     nDelete(&x);
     2818    sol[i] = n_Div(x, m_res[i]->m,_R->cf);
     2819    n_Delete(&x,_R->cf);
    24862820  }
    24872821  i--;
    … …  
    24962830      if (sol[j] != NULL)
    24972831      {
    2498         z = nMult(sol[j], s->m);
     2832        z = n_Mult(sol[j], s->m,_R->cf);
    24992833        if (x != NULL)
    25002834        {
    25012835          y = x;
    2502           x = nSub(y, z);
    2503           nDelete(&y);
    2504           nDelete(&z);
     2836          x = n_Sub(y, z,_R->cf);
     2837          n_Delete(&y,_R->cf);
     2838          n_Delete(&z,_R->cf);
    25052839        }
    25062840        else
    2507           x = nNeg(z);
     2841          x = n_Neg(z,_R->cf);
    25082842      }
    25092843      s = s->n;
    … …  
    25132847      if (x != NULL)
    25142848      {
    2515         y = nAdd(x, sol[i]);
    2516         nDelete(&x);
    2517         if (nIsZero(y))
    2518         {
    2519           nDelete(&sol[i]);
     2849        y = n_Add(x, sol[i],_R->cf);
     2850        n_Delete(&x,_R->cf);
     2851        if (n_IsZero(y,_R->cf))
     2852        {
     2853          n_Delete(&sol[i],_R->cf);
    25202854          sol[i] = NULL;
    25212855        }
    … …  
    25292863    {
    25302864      x = sol[i];
    2531       sol[i] = nDiv(x, d->m);
    2532       nDelete(&x);
     2865      sol[i] = n_Div(x, d->m,_R->cf);
     2866      n_Delete(&x,_R->cf);
    25332867    }
    25342868    i--;
    … …  
    25482882  {
    25492883    j = perm[i]-1;
    2550     res->m[j] = smSmnumber2Poly(sol[i]);
     2884    res->m[j] = sm_Smnumber2Poly(sol[i],_R);
    25512885  }
    25522886  omFreeSize((ADDRESS)sol, sizeof(number)*(crd+1));
    … …  
    25652899  int i, copt, ropt;
    25662900
    2567   xo=nInit(0);
     2901  xo=n_Init(0,_R->cf);
    25682902  for (i=act; i; i--)
    25692903  {
    … …  
    25722906    {
    25732907      x = a->m;
    2574       if (nGreaterZero(x))
    2575       {
    2576         if (nGreater(x,xo))
    2577         {
    2578           nDelete(&xo);
    2579           xo = nCopy(x);
     2908      if (n_GreaterZero(x,_R->cf))
     2909      {
     2910        if (n_Greater(x,xo,_R->cf))
     2911        {
     2912          n_Delete(&xo,_R->cf);
     2913          xo = n_Copy(x,_R->cf);
    25802914          copt = i;
    25812915          ropt = a->pos;
    … …  
    25842918      else
    25852919      {
    2586         xo = nNeg(xo);
    2587         if (nGreater(xo,x))
    2588         {
    2589           nDelete(&xo);
    2590           xo = nCopy(x);
     2920        xo = n_Neg(xo,_R->cf);
     2921        if (n_Greater(xo,x,_R->cf))
     2922        {
     2923          n_Delete(&xo,_R->cf);
     2924          xo = n_Copy(x,_R->cf);
    25912925          copt = i;
    25922926          ropt = a->pos;
    25932927        }
    2594         xo = nNeg(xo);
     2928        xo = n_Neg(xo,_R->cf);
    25952929      }
    25962930      a = a->n;
    … …  
    26042938    m_act[copt] = a;
    26052939  }
    2606   nDelete(&xo);
     2940  n_Delete(&xo,_R->cf);
    26072941}
    26082942
    … …  
    26122946void sparse_number_mat::smGElim()
    26132947{
    2614   number p = nInvers(piv->m);  // pivotelement
     2948  number p = n_Invers(piv->m,_R->cf);  // pivotelement
    26152949  smnumber c = m_act[act];      // pivotcolumn
    26162950  smnumber r = red;             // row to reduce
    … …  
    26202954  if ((c == NULL) || (r == NULL))
    26212955  {
    2622     while (r!=NULL) smNumberDelete(&r);
     2956    while (r!=NULL) sm_NumberDelete(&r,_R);
    26232957    return;
    26242958  }
    … …  
    26292963    res->n = NULL;
    26302964    b = c;
    2631     w = nMult(r->m, p);
    2632     nDelete(&r->m);
     2965    w = n_Mult(r->m, p,_R->cf);
     2966    n_Delete(&r->m,_R->cf);
    26332967    r->m = w;
    26342968    loop   // combine the chains a and b: a + w*b
    … …  
    26392973        {
    26402974          res = res->n = smNumberCopy(b);
    2641           res->m = nMult(b->m, w);
     2975          res->m = n_Mult(b->m, w,_R->cf);
    26422976          b = b->n;
    26432977        } while (b != NULL);
    … …  
    26522986      {
    26532987        res = res->n = smNumberCopy(b);
    2654         res->m = nMult(b->m, w);
     2988        res->m = n_Mult(b->m, w,_R->cf);
    26552989        b = b->n;
    26562990      }
    26572991      else
    26582992      {
    2659         hb = nMult(b->m, w);
    2660         ha = nAdd(a->m, hb);
    2661         nDelete(&hb);
    2662         nDelete(&a->m);
    2663         if (nIsZero(ha))
    2664         {
    2665           smNumberDelete(&a);
     2993        hb = n_Mult(b->m, w,_R->cf);
     2994        ha = n_Add(a->m, hb,_R->cf);
     2995        n_Delete(&hb,_R->cf);
     2996        n_Delete(&a->m,_R->cf);
     2997        if (n_IsZero(ha,_R->cf))
     2998        {
     2999          sm_NumberDelete(&a,_R);
    26663000        }
    26673001        else
    … …  
    26803014    }
    26813015    m_act[r->pos] = dumm->n;
    2682     smNumberDelete(&r);
     3016    sm_NumberDelete(&r,_R);
    26833017  } while (r != NULL);
    2684   nDelete(&p);
     3018  n_Delete(&p,_R->cf);
    26853019}
    26863020
    … …  
    27313065        {
    27323066          ap->n = a->n;
    2733           a->m = nNeg(a->m);
     3067          a->m = n_Neg(a->m,_R);
    27343068          b = b->n = a;
    27353069          b->pos = i;
    … …  
    27413075    {
    27423076      m_act[i] = a->n;
    2743       a->m = nNeg(a->m);
     3077      a->m = n_Neg(a->m,_R);
    27443078      b = b->n = a;
    27453079      b->pos = i;
    … …  
    28383172    a = m_act[i];
    28393173    while (a != NULL)
    2840       smNumberDelete(&a);
     3174      sm_NumberDelete(&a,_R);
    28413175  }
    28423176  for (i=crd; i; i--)
    … …  
    28443178    a = m_res[i];
    28453179    while (a != NULL)
    2846       smNumberDelete(&a);
     3180      sm_NumberDelete(&a,_R);
    28473181  }
    28483182  if (act)
    … …  
    28523186      a = m_row[i];
    28533187      while (a != NULL)
    2854         smNumberDelete(&a);
     3188        sm_NumberDelete(&a,_R);
    28553189    }
    28563190  }
    … …  
    28593193/* ----------------- internal 'C' stuff ------------------ */
    28603194
    2861 static void smNumberDelete(smnumber *r)
     3195static void sm_NumberDelete(smnumber *r, const ring R)
    28623196{
    28633197  smnumber a = *r, b = a->n;
    28643198
    2865   nDelete(&a->m);
     3199  n_Delete(&a->m,R->cf);
    28663200  omFreeBin((ADDRESS)a,  smnrec_bin);
    28673201  *r = b;
    … …  
    28793213* do not destroy p
    28803214*/
    2881 static smnumber smPoly2Smnumber(poly q)
     3215static smnumber sm_Poly2Smnumber(poly q, const ring R)
    28823216{
    28833217  smnumber a, res;
    … …  
    28873221    return NULL;
    28883222  a = res = (smnumber)omAllocBin(smnrec_bin);
    2889   a->pos = pGetComp(p);
     3223  a->pos = p_GetComp(p,R);
    28903224  a->m = pGetCoeff(p);
    2891   nNew(&pGetCoeff(p));
     3225  n_New(&pGetCoeff(p),R->cf);
    28923226  loop
    28933227  {
    … …  
    28953229    if (p == NULL)
    28963230    {
    2897       pDelete(&q);
     3231      p_Delete(&q,R);
    28983232      a->n = NULL;
    28993233      return res;
    29003234    }
    29013235    a = a->n = (smnumber)omAllocBin(smnrec_bin);
    2902     a->pos = pGetComp(p);
     3236    a->pos = p_GetComp(p,R);
    29033237    a->m = pGetCoeff(p);
    2904     nNew(&pGetCoeff(p));
     3238    n_New(&pGetCoeff(p),R->cf);
    29053239  }
    29063240}
    … …  
    29103244* destroy a
    29113245*/
    2912 static poly smSmnumber2Poly(number a)
     3246static poly sm_Smnumber2Poly(number a, const ring R)
    29133247{
    29143248  poly res;
    29153249
    29163250  if (a == NULL) return NULL;
    2917   res = pInit();
     3251  res = p_Init(R);
    29183252  pSetCoeff0(res, a);
    29193253  return res;
Note: See TracChangeset for help on using the changeset viewer.