Ignore:
Timestamp:
Apr 12, 2011, 3:00:18 PM (13 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
4aa8610313f7dca2da0233c6313149b9a6164584
Parents:
1377c9221c63abafb14abbd71a5934af8237629a
git-author:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2011-04-12 15:00:18+02:00
git-committer:
Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:12:28+01:00
Message:
same missing include for simpleideals.*
removed std stuff from simpleideals.cc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • libpolys/polys/simpleideals.cc

    r1377c9 rf71e8c5  
    1111#include "config.h"
    1212#include <misc/auxiliary.h>
     13#include <misc/options.h>
     14#include <omalloc/omalloc.h>
     15#include <polys/monomials/p_polys.h>
    1316
    1417#include "simpleideals.h"
     
    4447  {
    4548    Print("Module of rank %ld,real rank %ld and %d generators.\n",
    46           id->rank,idRankFreeModule(id, lmRing, tailRing),IDELEMS(id));
     49          id->rank,id_RankFreeModule(id, lmRing, tailRing),IDELEMS(id));
    4750
    4851    int j = (id->ncols*id->nrows) - 1;
     
    5861/* index of generator with leading term in ground ring (if any);
    5962   otherwise -1 */
    60 int idPosConstant(ideal id)
     63int id_PosConstant(ideal id, const ring r)
    6164{
    6265  int k;
    6366  for (k = IDELEMS(id)-1; k>=0; k--)
    6467  {
    65     if (p_LmIsConstantComp(id->m[k], currRing) == TRUE)
     68    if (p_LmIsConstantComp(id->m[k], r) == TRUE)
    6669      return k;
    6770  }
     
    7275* initialise the maximal ideal (at 0)
    7376*/
    74 ideal idMaxIdeal (void)
     77ideal id_MaxIdeal (const ring r)
    7578{
    7679  int l;
    7780  ideal hh=NULL;
    7881
    79   hh=idInit(pVariables,1);
    80   for (l=0; l<pVariables; l++)
    81   {
    82     hh->m[l] = pOne();
    83     pSetExp(hh->m[l],l+1,1);
    84     pSetm(hh->m[l]);
     82  hh=idInit(rVar(r),1);
     83  for (l=0; l<rVar(r); l++)
     84  {
     85    hh->m[l] = p_One(r);
     86    p_SetExp(hh->m[l],l+1,1,r);
     87    p_Setm(hh->m[l],r);
    8588  }
    8689  return hh;
     
    174177* (Note that the copied polynomials may be zero.)
    175178*/
    176 ideal idCopyFirstK (const ideal ide, const int k)
     179ideal id_CopyFirstK (const ideal ide, const int k,const ring r)
    177180{
    178181  ideal newI = idInit(k, 0);
    179182  for (int i = 0; i < k; i++)
    180     newI->m[i] = pCopy(ide->m[i]);
     183    newI->m[i] = p_Copy(ide->m[i],r);
    181184  return newI;
    182185}
     
    229232#else
    230233          if (pComparePolys(id->m[i], id->m[j])) pDelete(&id->m[j]);
    231 #endif
     234#endif         
    232235        }
    233236      }
     
    331334#ifdef HAVE_RINGS
    332335          }
    333 #endif
     336#endif         
    334337        }
    335338      }
     
    404407* for idSort: compare a and b revlex inclusive module comp.
    405408*/
    406 static int pComp_RevLex(poly a, poly b,BOOLEAN nolex)
     409static int p_Comp_RevLex(poly a, poly b,BOOLEAN nolex, const ring r)
    407410{
    408411  if (b==NULL) return 1;
    409412  if (a==NULL) return -1;
    410413
    411   if (nolex)
     414  if (nolex) 
    412415  {
    413416    int r=pLmCmp(a,b);
     
    418421    return r;
    419422  }
    420   int l=pVariables;
     423  int l=rVar(r);
    421424  while ((l>0) && (pGetExp(a,l)==pGetExp(b,l))) l--;
    422425  if (l==0)
     
    758761
    759762/*2
    760 *returns a minimized set of generators of h1
    761 */
    762 ideal idMinBase (ideal h1)
    763 {
    764   ideal h2, h3,h4,e;
    765   int j,k;
    766   int i,l,ll;
    767   intvec * wth;
    768   BOOLEAN homog;
    769 
    770   homog = idHomModule(h1,currQuotient,&wth);
    771   if (rHasGlobalOrdering_currRing())
    772   {
    773     if(!homog)
    774     {
    775       WarnS("minbase applies only to the local or homogeneous case");
    776       e=idCopy(h1);
    777       return e;
    778     }
    779     else
    780     {
    781       ideal re=kMin_std(h1,currQuotient,(tHomog)homog,&wth,h2,NULL,0,3);
    782       idDelete(&re);
    783       return h2;
    784     }
    785   }
    786   e=idInit(1,h1->rank);
    787   if (idIs0(h1))
    788   {
    789     return e;
    790   }
    791   pEnlargeSet(&(e->m),IDELEMS(e),15);
    792   IDELEMS(e) = 16;
    793   h2 = kStd(h1,currQuotient,isNotHomog,NULL);
    794   h3 = idMaxIdeal();
    795   h4=idMult(h2,h3);
    796   idDelete(&h3);
    797   h3=kStd(h4,currQuotient,isNotHomog,NULL);
    798   k = IDELEMS(h3);
    799   while ((k > 0) && (h3->m[k-1] == NULL)) k--;
    800   j = -1;
    801   l = IDELEMS(h2);
    802   while ((l > 0) && (h2->m[l-1] == NULL)) l--;
    803   for (i=l-1; i>=0; i--)
    804   {
    805     if (h2->m[i] != NULL)
    806     {
    807       ll = 0;
    808       while ((ll < k) && ((h3->m[ll] == NULL)
    809       || !pDivisibleBy(h3->m[ll],h2->m[i])))
    810         ll++;
    811       if (ll >= k)
    812       {
    813         j++;
    814         if (j > IDELEMS(e)-1)
    815         {
    816           pEnlargeSet(&(e->m),IDELEMS(e),16);
    817           IDELEMS(e) += 16;
    818         }
    819         e->m[j] = pCopy(h2->m[i]);
    820       }
    821     }
    822   }
    823   idDelete(&h2);
    824   idDelete(&h3);
    825   idDelete(&h4);
    826   if (currQuotient!=NULL)
    827   {
    828     h3=idInit(1,e->rank);
    829     h2=kNF(h3,currQuotient,e);
    830     idDelete(&h3);
    831     idDelete(&e);
    832     e=h2;
    833   }
    834   idSkipZeroes(e);
    835   return e;
    836 }
    837 
    838 /*2
    839763*the minimal index of used variables - 1
    840764*/
     
    11561080
    11571081  ideal res=idElimination(h,t);
    1158   // cleanup
     1082  // cleanup 
    11591083  idDelete(&h);
    1160   if (res!=NULL) res=idrMoveR(res,r,origRing);
     1084  res=idrMoveR(res,r,origRing);
    11611085  rChangeCurrRing(origRing);
    11621086  rKill(r);
    11631087  return res;
    1164 }
    1165 /*2
    1166 * h3 := h1 intersect h2
    1167 */
    1168 ideal idSect (ideal h1,ideal h2)
    1169 {
    1170   int i,j,k,length;
    1171   int flength = idRankFreeModule(h1);
    1172   int slength = idRankFreeModule(h2);
    1173   int rank=si_min(flength,slength);
    1174   if ((idIs0(h1)) || (idIs0(h2)))  return idInit(1,rank);
    1175 
    1176   ideal first,second,temp,temp1,result;
    1177   poly p,q;
    1178 
    1179   if (IDELEMS(h1)<IDELEMS(h2))
    1180   {
    1181     first = h1;
    1182     second = h2;
    1183   }
    1184   else
    1185   {
    1186     first = h2;
    1187     second = h1;
    1188     int t=flength; flength=slength; slength=t;
    1189   }
    1190   length  = si_max(flength,slength);
    1191   if (length==0)
    1192   {
    1193     if ((currQuotient==NULL)
    1194     && (currRing->OrdSgn==1)
    1195     && (!rIsPluralRing(currRing))
    1196     && ((TEST_V_INTERSECT_ELIM) || (!TEST_V_INTERSECT_SYZ)))
    1197       return idSectWithElim(first,second);
    1198     else length = 1;
    1199   }
    1200   if (TEST_OPT_PROT) PrintS("intersect by syzygy methods\n");
    1201   j = IDELEMS(first);
    1202 
    1203   ring orig_ring=currRing;
    1204   ring syz_ring=rCurrRingAssure_SyzComp();
    1205   rSetSyzComp(length);
    1206 
    1207   while ((j>0) && (first->m[j-1]==NULL)) j--;
    1208   temp = idInit(j /*IDELEMS(first)*/+IDELEMS(second),length+j);
    1209   k = 0;
    1210   for (i=0;i<j;i++)
    1211   {
    1212     if (first->m[i]!=NULL)
    1213     {
    1214       if (syz_ring==orig_ring)
    1215         temp->m[k] = pCopy(first->m[i]);
    1216       else
    1217         temp->m[k] = prCopyR(first->m[i], orig_ring);
    1218       q = pOne();
    1219       pSetComp(q,i+1+length);
    1220       pSetmComp(q);
    1221       if (flength==0) pShift(&(temp->m[k]),1);
    1222       p = temp->m[k];
    1223       while (pNext(p)!=NULL) pIter(p);
    1224       pNext(p) = q;
    1225       k++;
    1226     }
    1227   }
    1228   for (i=0;i<IDELEMS(second);i++)
    1229   {
    1230     if (second->m[i]!=NULL)
    1231     {
    1232       if (syz_ring==orig_ring)
    1233         temp->m[k] = pCopy(second->m[i]);
    1234       else
    1235         temp->m[k] = prCopyR(second->m[i], orig_ring);
    1236       if (slength==0) pShift(&(temp->m[k]),1);
    1237       k++;
    1238     }
    1239   }
    1240   intvec *w=NULL;
    1241   temp1 = kStd(temp,currQuotient,testHomog,&w,NULL,length);
    1242   if (w!=NULL) delete w;
    1243   idDelete(&temp);
    1244   if(syz_ring!=orig_ring)
    1245     rChangeCurrRing(orig_ring);
    1246 
    1247   result = idInit(IDELEMS(temp1),rank);
    1248   j = 0;
    1249   for (i=0;i<IDELEMS(temp1);i++)
    1250   {
    1251     if ((temp1->m[i]!=NULL)
    1252     && (p_GetComp(temp1->m[i],syz_ring)>length))
    1253     {
    1254       if(syz_ring==orig_ring)
    1255       {
    1256         p = temp1->m[i];
    1257       }
    1258       else
    1259       {
    1260         p = prMoveR(temp1->m[i], syz_ring);
    1261       }
    1262       temp1->m[i]=NULL;
    1263       while (p!=NULL)
    1264       {
    1265         q = pNext(p);
    1266         pNext(p) = NULL;
    1267         k = pGetComp(p)-1-length;
    1268         pSetComp(p,0);
    1269         pSetmComp(p);
    1270         /* Warning! multiply only from the left! it's very important for Plural */
    1271         result->m[j] = pAdd(result->m[j],pMult(p,pCopy(first->m[k])));
    1272         p = q;
    1273       }
    1274       j++;
    1275     }
    1276   }
    1277   if(syz_ring!=orig_ring)
    1278   {
    1279     rChangeCurrRing(syz_ring);
    1280     idDelete(&temp1);
    1281     rChangeCurrRing(orig_ring);
    1282     rKill(syz_ring);
    1283   }
    1284   else
    1285   {
    1286     idDelete(&temp1);
    1287   }
    1288 
    1289   idSkipZeroes(result);
    1290   if (TEST_OPT_RETURN_SB)
    1291   {
    1292      w=NULL;
    1293      temp1=kStd(result,currQuotient,testHomog,&w);
    1294      if (w!=NULL) delete w;
    1295      idDelete(&result);
    1296      idSkipZeroes(temp1);
    1297      return temp1;
    1298   }
    1299   else //temp1=kInterRed(result,currQuotient);
    1300     return result;
    1301 }
    1302 
    1303 /*2
    1304 * ideal/module intersection for a list of objects
    1305 * given as 'resolvente'
    1306 */
    1307 ideal idMultSect(resolvente arg, int length)
    1308 {
    1309   int i,j=0,k=0,syzComp,l,maxrk=-1,realrki;
    1310   ideal bigmat,tempstd,result;
    1311   poly p;
    1312   int isIdeal=0;
    1313   intvec * w=NULL;
    1314 
    1315   /* find 0-ideals and max rank -----------------------------------*/
    1316   for (i=0;i<length;i++)
    1317   {
    1318     if (!idIs0(arg[i]))
    1319     {
    1320       realrki=idRankFreeModule(arg[i]);
    1321       k++;
    1322       j += IDELEMS(arg[i]);
    1323       if (realrki>maxrk) maxrk = realrki;
    1324     }
    1325     else
    1326     {
    1327       if (arg[i]!=NULL)
    1328       {
    1329         return idInit(1,arg[i]->rank);
    1330       }
    1331     }
    1332   }
    1333   if (maxrk == 0)
    1334   {
    1335     isIdeal = 1;
    1336     maxrk = 1;
    1337   }
    1338   /* init -----------------------------------------------------------*/
    1339   j += maxrk;
    1340   syzComp = k*maxrk;
    1341 
    1342   ring orig_ring=currRing;
    1343   ring syz_ring=rCurrRingAssure_SyzComp();
    1344   rSetSyzComp(syzComp);
    1345 
    1346   bigmat = idInit(j,(k+1)*maxrk);
    1347   /* create unit matrices ------------------------------------------*/
    1348   for (i=0;i<maxrk;i++)
    1349   {
    1350     for (j=0;j<=k;j++)
    1351     {
    1352       p = pOne();
    1353       pSetComp(p,i+1+j*maxrk);
    1354       pSetmComp(p);
    1355       bigmat->m[i] = pAdd(bigmat->m[i],p);
    1356     }
    1357   }
    1358   /* enter given ideals ------------------------------------------*/
    1359   i = maxrk;
    1360   k = 0;
    1361   for (j=0;j<length;j++)
    1362   {
    1363     if (arg[j]!=NULL)
    1364     {
    1365       for (l=0;l<IDELEMS(arg[j]);l++)
    1366       {
    1367         if (arg[j]->m[l]!=NULL)
    1368         {
    1369           if (syz_ring==orig_ring)
    1370             bigmat->m[i] = pCopy(arg[j]->m[l]);
    1371           else
    1372             bigmat->m[i] = prCopyR(arg[j]->m[l], orig_ring);
    1373           pShift(&(bigmat->m[i]),k*maxrk+isIdeal);
    1374           i++;
    1375         }
    1376       }
    1377       k++;
    1378     }
    1379   }
    1380   /* std computation --------------------------------------------*/
    1381   tempstd = kStd(bigmat,currQuotient,testHomog,&w,NULL,syzComp);
    1382   if (w!=NULL) delete w;
    1383   idDelete(&bigmat);
    1384 
    1385   if(syz_ring!=orig_ring)
    1386     rChangeCurrRing(orig_ring);
    1387 
    1388   /* interprete result ----------------------------------------*/
    1389   result = idInit(IDELEMS(tempstd),maxrk);
    1390   k = 0;
    1391   for (j=0;j<IDELEMS(tempstd);j++)
    1392   {
    1393     if ((tempstd->m[j]!=NULL) && (p_GetComp(tempstd->m[j],syz_ring)>syzComp))
    1394     {
    1395       if (syz_ring==orig_ring)
    1396         p = pCopy(tempstd->m[j]);
    1397       else
    1398         p = prCopyR(tempstd->m[j], syz_ring);
    1399       pShift(&p,-syzComp-isIdeal);
    1400       result->m[k] = p;
    1401       k++;
    1402     }
    1403   }
    1404   /* clean up ----------------------------------------------------*/
    1405   if(syz_ring!=orig_ring)
    1406     rChangeCurrRing(syz_ring);
    1407   idDelete(&tempstd);
    1408   if(syz_ring!=orig_ring)
    1409   {
    1410     rChangeCurrRing(orig_ring);
    1411     rKill(syz_ring);
    1412   }
    1413   idSkipZeroes(result);
    1414   return result;
    1415 }
    1416 
    1417 /*2
    1418 *computes syzygies of h1,
    1419 *if quot != NULL it computes in the quotient ring modulo "quot"
    1420 *works always in a ring with ringorder_s
    1421 */
    1422 static ideal idPrepare (ideal  h1, tHomog hom, int syzcomp, intvec **w)
    1423 {
    1424   ideal   h2, h3;
    1425   int     i;
    1426   int     j,jj=0,k;
    1427   poly    p,q;
    1428 
    1429   if (idIs0(h1)) return NULL;
    1430   k = idRankFreeModule(h1);
    1431   h2=idCopy(h1);
    1432   i = IDELEMS(h2)-1;
    1433   if (k == 0)
    1434   {
    1435     for (j=0; j<=i; j++) pShift(&(h2->m[j]),1);
    1436     k = 1;
    1437   }
    1438   if (syzcomp<k)
    1439   {
    1440     Warn("syzcomp too low, should be %d instead of %d",k,syzcomp);
    1441     syzcomp = k;
    1442     rSetSyzComp(k);
    1443   }
    1444   h2->rank = syzcomp+i+1;
    1445 
    1446   //if (hom==testHomog)
    1447   //{
    1448   //  if(idHomIdeal(h1,currQuotient))
    1449   //  {
    1450   //    hom=TRUE;
    1451   //  }
    1452   //}
    1453 
    1454 #if MYTEST
    1455 #ifdef RDEBUG
    1456   Print("Prepare::h2: ");
    1457   idPrint(h2);
    1458 
    1459   for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);
    1460 
    1461 #endif
    1462 #endif
    1463 
    1464   for (j=0; j<=i; j++)
    1465   {
    1466     p = h2->m[j];
    1467     q = pOne();
    1468     pSetComp(q,syzcomp+1+j);
    1469     pSetmComp(q);
    1470     if (p!=NULL)
    1471     {
    1472       while (pNext(p)) pIter(p);
    1473       p->next = q;
    1474     }
    1475     else
    1476       h2->m[j]=q;
    1477   }
    1478 
    1479 #ifdef PDEBUG
    1480   for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);
    1481 
    1482 #if MYTEST
    1483 #ifdef RDEBUG
    1484   Print("Prepare::Input: ");
    1485   idPrint(h2);
    1486 
    1487   Print("Prepare::currQuotient: ");
    1488   idPrint(currQuotient);
    1489 #endif
    1490 #endif
    1491 
    1492 #endif
    1493 
    1494   idTest(h2);
    1495 
    1496   h3 = kStd(h2,currQuotient,hom,w,NULL,syzcomp);
    1497 
    1498 #if MYTEST
    1499 #ifdef RDEBUG
    1500   Print("Prepare::Output: ");
    1501   idPrint(h3);
    1502   for(j=0;j<IDELEMS(h2);j++) pTest(h3->m[j]);
    1503 #endif
    1504 #endif
    1505 
    1506 
    1507   idDelete(&h2);
    1508   return h3;
    1509 }
    1510 
    1511 /*2
    1512 * compute the syzygies of h1 in R/quot,
    1513 * weights of components are in w
    1514 * if setRegularity, return the regularity in deg
    1515 * do not change h1,  w
    1516 */
    1517 ideal idSyzygies (ideal  h1, tHomog h,intvec **w, BOOLEAN setSyzComp,
    1518                   BOOLEAN setRegularity, int *deg)
    1519 {
    1520   ideal s_h1;
    1521   poly  p;
    1522   int   j, k, length=0,reg;
    1523   BOOLEAN isMonomial=TRUE;
    1524   int ii, idElemens_h1;
    1525 
    1526   assume(h1 != NULL);
    1527 
    1528   idElemens_h1=IDELEMS(h1);
    1529 #ifdef PDEBUG
    1530   for(ii=0;ii<idElemens_h1 /*IDELEMS(h1)*/;ii++) pTest(h1->m[ii]);
    1531 #endif
    1532   if (idIs0(h1))
    1533   {
    1534     ideal result=idFreeModule(idElemens_h1/*IDELEMS(h1)*/);
    1535     int curr_syz_limit=rGetCurrSyzLimit();
    1536     if (curr_syz_limit>0)
    1537     for (ii=0;ii<idElemens_h1/*IDELEMS(h1)*/;ii++)
    1538     {
    1539       if (h1->m[ii]!=NULL)
    1540         pShift(&h1->m[ii],curr_syz_limit);
    1541     }
    1542     return result;
    1543   }
    1544   int slength=(int)idRankFreeModule(h1);
    1545   k=si_max(1,slength /*idRankFreeModule(h1)*/);
    1546 
    1547   assume(currRing != NULL);
    1548   ring orig_ring=currRing;
    1549   ring syz_ring=rCurrRingAssure_SyzComp();
    1550 
    1551   if (setSyzComp)
    1552     rSetSyzComp(k);
    1553 
    1554   if (orig_ring != syz_ring)
    1555   {
    1556     s_h1=idrCopyR_NoSort(h1,orig_ring);
    1557   }
    1558   else
    1559   {
    1560     s_h1 = h1;
    1561   }
    1562 
    1563   idTest(s_h1);
    1564 
    1565   ideal s_h3=idPrepare(s_h1,h,k,w); // main (syz) GB computation
    1566 
    1567   if (s_h3==NULL)
    1568   {
    1569     return idFreeModule( idElemens_h1 /*IDELEMS(h1)*/);
    1570   }
    1571 
    1572   if (orig_ring != syz_ring)
    1573   {
    1574     idDelete(&s_h1);
    1575     for (j=0; j<IDELEMS(s_h3); j++)
    1576     {
    1577       if (s_h3->m[j] != NULL)
    1578       {
    1579         if (p_MinComp(s_h3->m[j],syz_ring) > k)
    1580           pShift(&s_h3->m[j], -k);
    1581         else
    1582           pDelete(&s_h3->m[j]);
    1583       }
    1584     }
    1585     idSkipZeroes(s_h3);
    1586     s_h3->rank -= k;
    1587     rChangeCurrRing(orig_ring);
    1588     s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
    1589     rKill(syz_ring);
    1590     #ifdef HAVE_PLURAL
    1591     if (rIsPluralRing(currRing))
    1592     {
    1593       idDelMultiples(s_h3);
    1594       idSkipZeroes(s_h3);
    1595     }
    1596     #endif
    1597     idTest(s_h3);
    1598     return s_h3;
    1599   }
    1600 
    1601   ideal e = idInit(IDELEMS(s_h3), s_h3->rank);
    1602 
    1603   for (j=IDELEMS(s_h3)-1; j>=0; j--)
    1604   {
    1605     if (s_h3->m[j] != NULL)
    1606     {
    1607       if (p_MinComp(s_h3->m[j],syz_ring) <= k)
    1608       {
    1609         e->m[j] = s_h3->m[j];
    1610         isMonomial=isMonomial && (pNext(s_h3->m[j])==NULL);
    1611         pDelete(&pNext(s_h3->m[j]));
    1612         s_h3->m[j] = NULL;
    1613       }
    1614     }
    1615   }
    1616 
    1617   idSkipZeroes(s_h3);
    1618   idSkipZeroes(e);
    1619 
    1620   if ((deg != NULL)
    1621   && (!isMonomial)
    1622   && (!TEST_OPT_NOTREGULARITY)
    1623   && (setRegularity)
    1624   && (h==isHomog)
    1625   && (!rIsPluralRing(currRing))
    1626   )
    1627   {
    1628     ring dp_C_ring = rCurrRingAssure_dp_C();
    1629     if (dp_C_ring != syz_ring)
    1630       e = idrMoveR_NoSort(e, syz_ring);
    1631     resolvente res = sySchreyerResolvente(e,-1,&length,TRUE, TRUE);
    1632     intvec * dummy = syBetti(res,length,&reg, *w);
    1633     *deg = reg+2;
    1634     delete dummy;
    1635     for (j=0;j<length;j++)
    1636     {
    1637       if (res[j]!=NULL) idDelete(&(res[j]));
    1638     }
    1639     omFreeSize((ADDRESS)res,length*sizeof(ideal));
    1640     idDelete(&e);
    1641     if (dp_C_ring != syz_ring)
    1642     {
    1643       rChangeCurrRing(syz_ring);
    1644       rKill(dp_C_ring);
    1645     }
    1646   }
    1647   else
    1648   {
    1649     idDelete(&e);
    1650   }
    1651   idTest(s_h3);
    1652   if (currQuotient != NULL)
    1653   {
    1654     ideal ts_h3=kStd(s_h3,currQuotient,h,w);
    1655     idDelete(&s_h3);
    1656     s_h3 = ts_h3;
    1657   }
    1658   return s_h3;
    1659 }
    1660 
    1661 /*2
    1662 */
    1663 ideal idXXX (ideal  h1, int k)
    1664 {
    1665   ideal s_h1;
    1666   int j;
    1667   intvec *w=NULL;
    1668 
    1669   assume(currRing != NULL);
    1670   ring orig_ring=currRing;
    1671   ring syz_ring=rCurrRingAssure_SyzComp();
    1672 
    1673   rSetSyzComp(k);
    1674 
    1675   if (orig_ring != syz_ring)
    1676   {
    1677     s_h1=idrCopyR_NoSort(h1,orig_ring);
    1678   }
    1679   else
    1680   {
    1681     s_h1 = h1;
    1682   }
    1683 
    1684   ideal s_h3=kStd(s_h1,NULL,testHomog,&w,NULL,k);
    1685 
    1686   if (s_h3==NULL)
    1687   {
    1688     return idFreeModule(IDELEMS(h1));
    1689   }
    1690 
    1691   if (orig_ring != syz_ring)
    1692   {
    1693     idDelete(&s_h1);
    1694     idSkipZeroes(s_h3);
    1695     rChangeCurrRing(orig_ring);
    1696     s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
    1697     rKill(syz_ring);
    1698     idTest(s_h3);
    1699     return s_h3;
    1700   }
    1701 
    1702   idSkipZeroes(s_h3);
    1703   idTest(s_h3);
    1704   return s_h3;
    1705 }
    1706 
    1707 /*
    1708 *computes a standard basis for h1 and stores the transformation matrix
    1709 * in ma
    1710 */
    1711 ideal idLiftStd (ideal  h1, matrix* ma, tHomog hi, ideal * syz)
    1712 {
    1713   int   i, j, k, t, inputIsIdeal=idRankFreeModule(h1);
    1714   poly  p=NULL, q, qq;
    1715   intvec *w=NULL;
    1716 
    1717   idDelete((ideal*)ma);
    1718   BOOLEAN lift3=FALSE;
    1719   if (syz!=NULL) { lift3=TRUE; idDelete(syz); }
    1720   if (idIs0(h1))
    1721   {
    1722     *ma=mpNew(1,0);
    1723     if (lift3)
    1724     {
    1725       *syz=idFreeModule(IDELEMS(h1));
    1726       int curr_syz_limit=rGetCurrSyzLimit();
    1727       if (curr_syz_limit>0)
    1728       for (int ii=0;ii<IDELEMS(h1);ii++)
    1729       {
    1730         if (h1->m[ii]!=NULL)
    1731           pShift(&h1->m[ii],curr_syz_limit);
    1732       }
    1733     }
    1734     return idInit(1,h1->rank);
    1735   }
    1736 
    1737   BITSET save_verbose=verbose;
    1738 
    1739   k=si_max(1,(int)idRankFreeModule(h1));
    1740 
    1741   if ((k==1) && (!lift3)) verbose |=Sy_bit(V_IDLIFT);
    1742 
    1743   ring orig_ring = currRing;
    1744   ring syz_ring = rCurrRingAssure_SyzComp();
    1745   rSetSyzComp(k);
    1746 
    1747   ideal s_h1=h1;
    1748 
    1749   if (orig_ring != syz_ring)
    1750     s_h1 = idrCopyR_NoSort(h1,orig_ring);
    1751   else
    1752     s_h1 = h1;
    1753 
    1754   ideal s_h3=idPrepare(s_h1,hi,k,&w); // main (syz) GB computation
    1755 
    1756   ideal s_h2 = idInit(IDELEMS(s_h3), s_h3->rank);
    1757 
    1758   if (lift3) (*syz)=idInit(IDELEMS(s_h3),IDELEMS(h1));
    1759 
    1760   if (w!=NULL) delete w;
    1761   i = 0;
    1762 
    1763   // now sort the result, SB : leave in s_h3
    1764   //                      T:  put in s_h2
    1765   //                      syz: put in *syz
    1766   for (j=0; j<IDELEMS(s_h3); j++)
    1767   {
    1768     if (s_h3->m[j] != NULL)
    1769     {
    1770       //if (p_MinComp(s_h3->m[j],syz_ring) <= k)
    1771       if (pGetComp(s_h3->m[j]) <= k) // syz_ring == currRing
    1772       {
    1773         i++;
    1774         q = s_h3->m[j];
    1775         while (pNext(q) != NULL)
    1776         {
    1777           if (pGetComp(pNext(q)) > k)
    1778           {
    1779             s_h2->m[j] = pNext(q);
    1780             pNext(q) = NULL;
    1781           }
    1782           else
    1783           {
    1784             pIter(q);
    1785           }
    1786         }
    1787         if (!inputIsIdeal) pShift(&(s_h3->m[j]), -1);
    1788       }
    1789       else
    1790       {
    1791         // we a syzygy here:
    1792         if (lift3)
    1793         {
    1794           pShift(&s_h3->m[j], -k);
    1795           (*syz)->m[j]=s_h3->m[j];
    1796           s_h3->m[j]=NULL;
    1797         }
    1798         else
    1799           pDelete(&(s_h3->m[j]));
    1800       }
    1801     }
    1802   }
    1803   idSkipZeroes(s_h3);
    1804   //extern char * iiStringMatrix(matrix im, int dim,char ch);
    1805   //PrintS("SB: ----------------------------------------\n");
    1806   //PrintS(iiStringMatrix((matrix)s_h3,k,'\n'));
    1807   //PrintLn();
    1808   //PrintS("T: ----------------------------------------\n");
    1809   //PrintS(iiStringMatrix((matrix)s_h2,h1->rank,'\n'));
    1810   //PrintLn();
    1811 
    1812   if (lift3) idSkipZeroes(*syz);
    1813 
    1814   j = IDELEMS(s_h1);
    1815 
    1816 
    1817   if (syz_ring!=orig_ring)
    1818   {
    1819     idDelete(&s_h1);
    1820     rChangeCurrRing(orig_ring);
    1821   }
    1822 
    1823   *ma = mpNew(j,i);
    1824 
    1825   i = 1;
    1826   for (j=0; j<IDELEMS(s_h2); j++)
    1827   {
    1828     if (s_h2->m[j] != NULL)
    1829     {
    1830       q = prMoveR( s_h2->m[j], syz_ring);
    1831       s_h2->m[j] = NULL;
    1832 
    1833       while (q != NULL)
    1834       {
    1835         p = q;
    1836         pIter(q);
    1837         pNext(p) = NULL;
    1838         t=pGetComp(p);
    1839         pSetComp(p,0);
    1840         pSetmComp(p);
    1841         MATELEM(*ma,t-k,i) = pAdd(MATELEM(*ma,t-k,i),p);
    1842       }
    1843       i++;
    1844     }
    1845   }
    1846   idDelete(&s_h2);
    1847 
    1848   for (i=0; i<IDELEMS(s_h3); i++)
    1849   {
    1850     s_h3->m[i] = prMoveR_NoSort(s_h3->m[i], syz_ring);
    1851   }
    1852   if (lift3)
    1853   {
    1854     for (i=0; i<IDELEMS(*syz); i++)
    1855     {
    1856       (*syz)->m[i] = prMoveR_NoSort((*syz)->m[i], syz_ring);
    1857     }
    1858   }
    1859 
    1860   if (syz_ring!=orig_ring) rKill(syz_ring);
    1861   verbose = save_verbose;
    1862   return s_h3;
    1863 }
    1864 
    1865 static void idPrepareStd(ideal s_temp, int k)
    1866 {
    1867   int j,rk=idRankFreeModule(s_temp);
    1868   poly p,q;
    1869 
    1870   if (rk == 0)
    1871   {
    1872     for (j=0; j<IDELEMS(s_temp); j++)
    1873     {
    1874       if (s_temp->m[j]!=NULL) pSetCompP(s_temp->m[j],1);
    1875     }
    1876     k = si_max(k,1);
    1877   }
    1878   for (j=0; j<IDELEMS(s_temp); j++)
    1879   {
    1880     if (s_temp->m[j]!=NULL)
    1881     {
    1882       p = s_temp->m[j];
    1883       q = pOne();
    1884       //pGetCoeff(q)=nNeg(pGetCoeff(q));   //set q to -1
    1885       pSetComp(q,k+1+j);
    1886       pSetmComp(q);
    1887       while (pNext(p)) pIter(p);
    1888       pNext(p) = q;
    1889     }
    1890   }
    1891 }
    1892 
    1893 /*2
    1894 *computes a representation of the generators of submod with respect to those
    1895 * of mod
    1896 */
    1897 
    1898 ideal idLift(ideal mod, ideal submod,ideal *rest, BOOLEAN goodShape,
    1899              BOOLEAN isSB, BOOLEAN divide, matrix *unit)
    1900 {
    1901   int lsmod =idRankFreeModule(submod), i, j, k;
    1902   int comps_to_add=0;
    1903   poly p;
    1904 
    1905   if (idIs0(submod))
    1906   {
    1907     if (unit!=NULL)
    1908     {
    1909       *unit=mpNew(1,1);
    1910       MATELEM(*unit,1,1)=pOne();
    1911     }
    1912     if (rest!=NULL)
    1913     {
    1914       *rest=idInit(1,mod->rank);
    1915     }
    1916     return idInit(1,mod->rank);
    1917   }
    1918   if (idIs0(mod)) /* and not idIs0(submod) */
    1919   {
    1920     WerrorS("2nd module does not lie in the first");
    1921     return NULL;
    1922   }
    1923   if (unit!=NULL)
    1924   {
    1925     comps_to_add = IDELEMS(submod);
    1926     while ((comps_to_add>0) && (submod->m[comps_to_add-1]==NULL))
    1927       comps_to_add--;
    1928   }
    1929   k=si_max(idRankFreeModule(mod),idRankFreeModule(submod));
    1930   if  ((k!=0) && (lsmod==0)) lsmod=1;
    1931   k=si_max(k,(int)mod->rank);
    1932   if (k<submod->rank) { WarnS("rk(submod) > rk(mod) ?");k=submod->rank; }
    1933 
    1934   ring orig_ring=currRing;
    1935   ring syz_ring=rCurrRingAssure_SyzComp();
    1936   rSetSyzComp(k);
    1937 
    1938   ideal s_mod, s_temp;
    1939   if (orig_ring != syz_ring)
    1940   {
    1941     s_mod = idrCopyR_NoSort(mod,orig_ring);
    1942     s_temp = idrCopyR_NoSort(submod,orig_ring);
    1943   }
    1944   else
    1945   {
    1946     s_mod = mod;
    1947     s_temp = idCopy(submod);
    1948   }
    1949   ideal s_h3;
    1950   if (isSB)
    1951   {
    1952     s_h3 = idCopy(s_mod);
    1953     idPrepareStd(s_h3, k+comps_to_add);
    1954   }
    1955   else
    1956   {
    1957     s_h3 = idPrepare(s_mod,(tHomog)FALSE,k+comps_to_add,NULL);
    1958   }
    1959   if (!goodShape)
    1960   {
    1961     for (j=0;j<IDELEMS(s_h3);j++)
    1962     {
    1963       if ((s_h3->m[j] != NULL) && (pMinComp(s_h3->m[j]) > k))
    1964         pDelete(&(s_h3->m[j]));
    1965     }
    1966   }
    1967   idSkipZeroes(s_h3);
    1968   if (lsmod==0)
    1969   {
    1970     for (j=IDELEMS(s_temp);j>0;j--)
    1971     {
    1972       if (s_temp->m[j-1]!=NULL)
    1973         pShift(&(s_temp->m[j-1]),1);
    1974     }
    1975   }
    1976   if (unit!=NULL)
    1977   {
    1978     for(j = 0;j<comps_to_add;j++)
    1979     {
    1980       p = s_temp->m[j];
    1981       if (p!=NULL)
    1982       {
    1983         while (pNext(p)!=NULL) pIter(p);
    1984         pNext(p) = pOne();
    1985         pIter(p);
    1986         pSetComp(p,1+j+k);
    1987         pSetmComp(p);
    1988         p = pNeg(p);
    1989       }
    1990     }
    1991   }
    1992   ideal s_result = kNF(s_h3,currQuotient,s_temp,k);
    1993   s_result->rank = s_h3->rank;
    1994   ideal s_rest = idInit(IDELEMS(s_result),k);
    1995   idDelete(&s_h3);
    1996   idDelete(&s_temp);
    1997 
    1998   for (j=0;j<IDELEMS(s_result);j++)
    1999   {
    2000     if (s_result->m[j]!=NULL)
    2001     {
    2002       if (pGetComp(s_result->m[j])<=k)
    2003       {
    2004         if (!divide)
    2005         {
    2006           if (isSB)
    2007           {
    2008             WarnS("first module not a standardbasis\n"
    2009               "// ** or second not a proper submodule");
    2010           }
    2011           else
    2012             WerrorS("2nd module does not lie in the first");
    2013           idDelete(&s_result);
    2014           idDelete(&s_rest);
    2015           s_result=idInit(IDELEMS(submod),submod->rank);
    2016           break;
    2017         }
    2018         else
    2019         {
    2020           p = s_rest->m[j] = s_result->m[j];
    2021           while ((pNext(p)!=NULL) && (pGetComp(pNext(p))<=k)) pIter(p);
    2022           s_result->m[j] = pNext(p);
    2023           pNext(p) = NULL;
    2024         }
    2025       }
    2026       pShift(&(s_result->m[j]),-k);
    2027       pNeg(s_result->m[j]);
    2028     }
    2029   }
    2030   if ((lsmod==0) && (!idIs0(s_rest)))
    2031   {
    2032     for (j=IDELEMS(s_rest);j>0;j--)
    2033     {
    2034       if (s_rest->m[j-1]!=NULL)
    2035       {
    2036         pShift(&(s_rest->m[j-1]),-1);
    2037         s_rest->m[j-1] = s_rest->m[j-1];
    2038       }
    2039     }
    2040   }
    2041   if(syz_ring!=orig_ring)
    2042   {
    2043     idDelete(&s_mod);
    2044     rChangeCurrRing(orig_ring);
    2045     s_result = idrMoveR_NoSort(s_result, syz_ring);
    2046     s_rest = idrMoveR_NoSort(s_rest, syz_ring);
    2047     rKill(syz_ring);
    2048   }
    2049   if (rest!=NULL)
    2050     *rest = s_rest;
    2051   else
    2052     idDelete(&s_rest);
    2053 //idPrint(s_result);
    2054   if (unit!=NULL)
    2055   {
    2056     *unit=mpNew(comps_to_add,comps_to_add);
    2057     int i;
    2058     for(i=0;i<IDELEMS(s_result);i++)
    2059     {
    2060       poly p=s_result->m[i];
    2061       poly q=NULL;
    2062       while(p!=NULL)
    2063       {
    2064         if(pGetComp(p)<=comps_to_add)
    2065         {
    2066           pSetComp(p,0);
    2067           if (q!=NULL)
    2068           {
    2069             pNext(q)=pNext(p);
    2070           }
    2071           else
    2072           {
    2073             pIter(s_result->m[i]);
    2074           }
    2075           pNext(p)=NULL;
    2076           MATELEM(*unit,i+1,i+1)=pAdd(MATELEM(*unit,i+1,i+1),p);
    2077           if(q!=NULL)   p=pNext(q);
    2078           else          p=s_result->m[i];
    2079         }
    2080         else
    2081         {
    2082           q=p;
    2083           pIter(p);
    2084         }
    2085       }
    2086       pShift(&s_result->m[i],-comps_to_add);
    2087     }
    2088   }
    2089   return s_result;
    2090 }
    2091 
    2092 /*2
    2093 *computes division of P by Q with remainder up to (w-weighted) degree n
    2094 *P, Q, and w are not changed
    2095 */
    2096 void idLiftW(ideal P,ideal Q,int n,matrix &T, ideal &R,short *w)
    2097 {
    2098   long N=0;
    2099   int i;
    2100   for(i=IDELEMS(Q)-1;i>=0;i--)
    2101     if(w==NULL)
    2102       N=si_max(N,pDeg(Q->m[i]));
    2103     else
    2104       N=si_max(N,pDegW(Q->m[i],w));
    2105   N+=n;
    2106 
    2107   T=mpNew(IDELEMS(Q),IDELEMS(P));
    2108   R=idInit(IDELEMS(P),P->rank);
    2109 
    2110   for(i=IDELEMS(P)-1;i>=0;i--)
    2111   {
    2112     poly p;
    2113     if(w==NULL)
    2114       p=ppJet(P->m[i],N);
    2115     else
    2116       p=ppJetW(P->m[i],N,w);
    2117 
    2118     int j=IDELEMS(Q)-1;
    2119     while(p!=NULL)
    2120     {
    2121       if(pDivisibleBy(Q->m[j],p))
    2122       {
    2123         poly p0=pDivideM(pHead(p),pHead(Q->m[j]));
    2124         if(w==NULL)
    2125           p=pJet(pSub(p,ppMult_mm(Q->m[j],p0)),N);
    2126         else
    2127           p=pJetW(pSub(p,ppMult_mm(Q->m[j],p0)),N,w);
    2128         pNormalize(p);
    2129         if((w==NULL)&&(pDeg(p0)>n)||(w!=NULL)&&(pDegW(p0,w)>n))
    2130           pDelete(&p0);
    2131         else
    2132           MATELEM(T,j+1,i+1)=pAdd(MATELEM(T,j+1,i+1),p0);
    2133         j=IDELEMS(Q)-1;
    2134       }
    2135       else
    2136       {
    2137         if(j==0)
    2138         {
    2139           poly p0=p;
    2140           pIter(p);
    2141           pNext(p0)=NULL;
    2142           if(((w==NULL)&&(pDeg(p0)>n))
    2143           ||((w!=NULL)&&(pDegW(p0,w)>n)))
    2144             pDelete(&p0);
    2145           else
    2146             R->m[i]=pAdd(R->m[i],p0);
    2147           j=IDELEMS(Q)-1;
    2148         }
    2149         else
    2150           j--;
    2151       }
    2152     }
    2153   }
    2154 }
    2155 
    2156 /*2
    2157 *computes the quotient of h1,h2 : internal routine for idQuot
    2158 *BEWARE: the returned ideals may contain incorrectly ordered polys !
    2159 *
    2160 */
    2161 static ideal idInitializeQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb,
    2162                                BOOLEAN *addOnlyOne, int *kkmax)
    2163 {
    2164   ideal temph1;
    2165   poly     p,q = NULL;
    2166   int i,l,ll,k,kkk,kmax;
    2167   int j = 0;
    2168   int k1 = idRankFreeModule(h1);
    2169   int k2 = idRankFreeModule(h2);
    2170   tHomog   hom=isNotHomog;
    2171 
    2172   k=si_max(k1,k2);
    2173   if (k==0)
    2174     k = 1;
    2175   if ((k2==0) && (k>1)) *addOnlyOne = FALSE;
    2176 
    2177   intvec * weights;
    2178   hom = (tHomog)idHomModule(h1,currQuotient,&weights);
    2179   if (/**addOnlyOne &&*/ (!h1IsStb))
    2180     temph1 = kStd(h1,currQuotient,hom,&weights,NULL);
    2181   else
    2182     temph1 = idCopy(h1);
    2183   if (weights!=NULL) delete weights;
    2184   idTest(temph1);
    2185 /*--- making a single vector from h2 ---------------------*/
    2186   for (i=0; i<IDELEMS(h2); i++)
    2187   {
    2188     if (h2->m[i] != NULL)
    2189     {
    2190       p = pCopy(h2->m[i]);
    2191       if (k2 == 0)
    2192         pShift(&p,j*k+1);
    2193       else
    2194         pShift(&p,j*k);
    2195       q = pAdd(q,p);
    2196       j++;
    2197     }
    2198   }
    2199   *kkmax = kmax = j*k+1;
    2200 /*--- adding a monomial for the result (syzygy) ----------*/
    2201   p = q;
    2202   while (pNext(p)!=NULL) pIter(p);
    2203   pNext(p) = pOne();
    2204   pIter(p);
    2205   pSetComp(p,kmax);
    2206   pSetmComp(p);
    2207 /*--- constructing the big matrix ------------------------*/
    2208   ideal h4 = idInit(16,kmax+k-1);
    2209   h4->m[0] = q;
    2210   if (k2 == 0)
    2211   {
    2212     if (k > IDELEMS(h4))
    2213     {
    2214       pEnlargeSet(&(h4->m),IDELEMS(h4),k-IDELEMS(h4));
    2215       IDELEMS(h4) = k;
    2216     }
    2217     for (i=1; i<k; i++)
    2218     {
    2219       if (h4->m[i-1]!=NULL)
    2220       {
    2221         p = pCopy_noCheck(h4->m[i-1]);
    2222         pShift(&p,1);
    2223         h4->m[i] = p;
    2224       }
    2225     }
    2226   }
    2227   idSkipZeroes(h4);
    2228   kkk = IDELEMS(h4);
    2229   i = IDELEMS(temph1);
    2230   for (l=0; l<i; l++)
    2231   {
    2232     if(temph1->m[l]!=NULL)
    2233     {
    2234       for (ll=0; ll<j; ll++)
    2235       {
    2236         p = pCopy(temph1->m[l]);
    2237         if (k1 == 0)
    2238           pShift(&p,ll*k+1);
    2239         else
    2240           pShift(&p,ll*k);
    2241         if (kkk >= IDELEMS(h4))
    2242         {
    2243           pEnlargeSet(&(h4->m),IDELEMS(h4),16);
    2244           IDELEMS(h4) += 16;
    2245         }
    2246         h4->m[kkk] = p;
    2247         kkk++;
    2248       }
    2249     }
    2250   }
    2251 /*--- if h2 goes in as single vector - the h1-part is just SB ---*/
    2252   if (*addOnlyOne)
    2253   {
    2254     idSkipZeroes(h4);
    2255     p = h4->m[0];
    2256     for (i=0;i<IDELEMS(h4)-1;i++)
    2257     {
    2258       h4->m[i] = h4->m[i+1];
    2259     }
    2260     h4->m[IDELEMS(h4)-1] = p;
    2261     test |= Sy_bit(OPT_SB_1);
    2262   }
    2263   idDelete(&temph1);
    2264   return h4;
    2265 }
    2266 /*2
    2267 *computes the quotient of h1,h2
    2268 */
    2269 ideal idQuot (ideal  h1, ideal h2, BOOLEAN h1IsStb, BOOLEAN resultIsIdeal)
    2270 {
    2271   // first check for special case h1:(0)
    2272   if (idIs0(h2))
    2273   {
    2274     ideal res;
    2275     if (resultIsIdeal)
    2276     {
    2277       res = idInit(1,1);
    2278       res->m[0] = pOne();
    2279     }
    2280     else
    2281       res = idFreeModule(h1->rank);
    2282     return res;
    2283   }
    2284   BITSET old_test=test;
    2285   int i,l,ll,k,kkk,kmax;
    2286   BOOLEAN  addOnlyOne=TRUE;
    2287   tHomog   hom=isNotHomog;
    2288   intvec * weights1;
    2289 
    2290   ideal s_h4 = idInitializeQuot (h1,h2,h1IsStb,&addOnlyOne,&kmax);
    2291 
    2292   hom = (tHomog)idHomModule(s_h4,currQuotient,&weights1);
    2293 
    2294   ring orig_ring=currRing;
    2295   ring syz_ring=rCurrRingAssure_SyzComp();
    2296   rSetSyzComp(kmax-1);
    2297   if (orig_ring!=syz_ring)
    2298   //  s_h4 = idrMoveR_NoSort(s_h4,orig_ring);
    2299     s_h4 = idrMoveR(s_h4,orig_ring);
    2300   idTest(s_h4);
    2301   #if 0
    2302   void ipPrint_MA0(matrix m, const char *name);
    2303   matrix m=idModule2Matrix(idCopy(s_h4));
    2304   PrintS("start:\n");
    2305   ipPrint_MA0(m,"Q");
    2306   idDelete((ideal *)&m);
    2307   PrintS("last elem:");wrp(s_h4->m[IDELEMS(s_h4)-1]);PrintLn();
    2308   #endif
    2309   ideal s_h3;
    2310   if (addOnlyOne)
    2311   {
    2312     s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,0/*kmax-1*/,IDELEMS(s_h4)-1);
    2313   }
    2314   else
    2315   {
    2316     s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1);
    2317   }
    2318   test = old_test;
    2319   #if 0
    2320   // only together with the above debug stuff
    2321   idSkipZeroes(s_h3);
    2322   m=idModule2Matrix(idCopy(s_h3));
    2323   Print("result, kmax=%d:\n",kmax);
    2324   ipPrint_MA0(m,"S");
    2325   idDelete((ideal *)&m);
    2326   #endif
    2327   idTest(s_h3);
    2328   if (weights1!=NULL) delete weights1;
    2329   idDelete(&s_h4);
    2330 
    2331   for (i=0;i<IDELEMS(s_h3);i++)
    2332   {
    2333     if ((s_h3->m[i]!=NULL) && (pGetComp(s_h3->m[i])>=kmax))
    2334     {
    2335       if (resultIsIdeal)
    2336         pShift(&s_h3->m[i],-kmax);
    2337       else
    2338         pShift(&s_h3->m[i],-kmax+1);
    2339     }
    2340     else
    2341       pDelete(&s_h3->m[i]);
    2342   }
    2343   if (resultIsIdeal)
    2344     s_h3->rank = 1;
    2345   else
    2346     s_h3->rank = h1->rank;
    2347   if(syz_ring!=orig_ring)
    2348   {
    2349     rChangeCurrRing(orig_ring);
    2350     s_h3 = idrMoveR_NoSort(s_h3, syz_ring);
    2351     rKill(syz_ring);
    2352   }
    2353   idSkipZeroes(s_h3);
    2354   idTest(s_h3);
    2355   return s_h3;
    23561088}
    23571089
     
    25521284  idSkipZeroes(result);
    25531285  return result;
    2554 }
    2555 
    2556 /*2
    2557 * eliminate delVar (product of vars) in h1
    2558 */
    2559 ideal idElimination (ideal h1,poly delVar,intvec *hilb)
    2560 {
    2561   int    i,j=0,k,l;
    2562   ideal  h,hh, h3;
    2563   int    *ord,*block0,*block1;
    2564   int    ordersize=2;
    2565   int    **wv;
    2566   tHomog hom;
    2567   intvec * w;
    2568   ring tmpR;
    2569   ring origR = currRing;
    2570 
    2571   if (delVar==NULL)
    2572   {
    2573     return idCopy(h1);
    2574   }
    2575   if ((currQuotient!=NULL) && rIsPluralRing(origR))
    2576   {
    2577     WerrorS("cannot eliminate in a qring");
    2578     return NULL;
    2579   }
    2580   if (idIs0(h1)) return idInit(1,h1->rank);
    2581 #ifdef HAVE_PLURAL
    2582   if (rIsPluralRing(origR))
    2583     /* in the NC case, we have to check the admissibility of */
    2584     /* the subalgebra to be intersected with */
    2585   {
    2586     if ((ncRingType(origR) != nc_skew) && (ncRingType(origR) != nc_exterior)) /* in (quasi)-commutative algebras every subalgebra is admissible */
    2587     {
    2588       if (nc_CheckSubalgebra(delVar,origR))
    2589       {
    2590         WerrorS("no elimination is possible: subalgebra is not admissible");
    2591         return NULL;
    2592       }
    2593     }
    2594   }
    2595 #endif
    2596   hom=(tHomog)idHomModule(h1,NULL,&w); //sets w to weight vector or NULL
    2597   h3=idInit(16,h1->rank);
    2598   for (k=0;; k++)
    2599   {
    2600     if (origR->order[k]!=0) ordersize++;
    2601     else break;
    2602   }
    2603 #if 0
    2604   if (rIsPluralRing(origR)) // we have too keep the odering: it may be needed
    2605                             // for G-algebra
    2606   {
    2607     for (k=0;k<ordersize-1; k++)
    2608     {
    2609       block0[k+1] = origR->block0[k];
    2610       block1[k+1] = origR->block1[k];
    2611       ord[k+1] = origR->order[k];
    2612       if (origR->wvhdl[k]!=NULL) wv[k+1] = (int*) omMemDup(origR->wvhdl[k]);
    2613     }
    2614   }
    2615   else
    2616   {
    2617     block0[1] = 1;
    2618     block1[1] = pVariables;
    2619     if (origR->OrdSgn==1) ord[1] = ringorder_wp;
    2620     else                  ord[1] = ringorder_ws;
    2621     wv[1]=(int*)omAlloc0(pVariables*sizeof(int));
    2622     double wNsqr = (double)2.0 / (double)pVariables;
    2623     wFunctional = wFunctionalBuch;
    2624     int  *x= (int * )omAlloc(2 * (pVariables + 1) * sizeof(int));
    2625     int sl=IDELEMS(h1) - 1;
    2626     wCall(h1->m, sl, x, wNsqr);
    2627     for (sl = pVariables; sl!=0; sl--)
    2628       wv[1][sl-1] = x[sl + pVariables + 1];
    2629     omFreeSize((ADDRESS)x, 2 * (pVariables + 1) * sizeof(int));
    2630 
    2631     ord[2]=ringorder_C;
    2632     ord[3]=0;
    2633   }
    2634 #else
    2635 #endif
    2636   if ((hom==TRUE) && (origR->OrdSgn==1) && (!rIsPluralRing(origR)))
    2637   {
    2638     #if 1
    2639     // we change to an ordering:
    2640     // aa(1,1,1,...,0,0,0),wp(...),C
    2641     // this seems to be better than version 2 below,
    2642     // according to Tst/../elimiate_[3568].tat (- 17 %)
    2643     ord=(int*)omAlloc0(4*sizeof(int));
    2644     block0=(int*)omAlloc0(4*sizeof(int));
    2645     block1=(int*)omAlloc0(4*sizeof(int));
    2646     wv=(int**) omAlloc0(4*sizeof(int**));
    2647     block0[0] = block0[1] = 1;
    2648     block1[0] = block1[1] = rVar(origR);
    2649     wv[0]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
    2650     // use this special ordering: like ringorder_a, except that pFDeg, pWeights
    2651     // ignore it
    2652     ord[0] = ringorder_aa;
    2653     for (j=0;j<rVar(origR);j++)
    2654       if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
    2655     BOOLEAN wp=FALSE;
    2656     for (j=0;j<rVar(origR);j++)
    2657       if (pWeight(j+1,origR)!=1) { wp=TRUE;break; }
    2658     if (wp)
    2659     {
    2660       wv[1]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
    2661       for (j=0;j<rVar(origR);j++)
    2662         wv[1][j]=pWeight(j+1,origR);
    2663       ord[1] = ringorder_wp;
    2664     }
    2665     else
    2666       ord[1] = ringorder_dp;
    2667     #else
    2668     // we change to an ordering:
    2669     // a(w1,...wn),wp(1,...0.....),C
    2670     ord=(int*)omAlloc0(4*sizeof(int));
    2671     block0=(int*)omAlloc0(4*sizeof(int));
    2672     block1=(int*)omAlloc0(4*sizeof(int));
    2673     wv=(int**) omAlloc0(4*sizeof(int**));
    2674     block0[0] = block0[1] = 1;
    2675     block1[0] = block1[1] = rVar(origR);
    2676     wv[0]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
    2677     wv[1]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
    2678     ord[0] = ringorder_a;
    2679     for (j=0;j<rVar(origR);j++)
    2680       wv[0][j]=pWeight(j+1,origR);
    2681     ord[1] = ringorder_wp;
    2682     for (j=0;j<rVar(origR);j++)
    2683       if (pGetExp(delVar,j+1)!=0) wv[1][j]=1;
    2684     #endif
    2685     ord[2] = ringorder_C;
    2686     ord[3] = 0;
    2687   }
    2688   else
    2689   {
    2690     // we change to an ordering:
    2691     // aa(....),orig_ordering
    2692     ord=(int*)omAlloc0(ordersize*sizeof(int));
    2693     block0=(int*)omAlloc0(ordersize*sizeof(int));
    2694     block1=(int*)omAlloc0(ordersize*sizeof(int));
    2695     wv=(int**) omAlloc0(ordersize*sizeof(int**));
    2696     for (k=0;k<ordersize-1; k++)
    2697     {
    2698       block0[k+1] = origR->block0[k];
    2699       block1[k+1] = origR->block1[k];
    2700       ord[k+1] = origR->order[k];
    2701       if (origR->wvhdl[k]!=NULL) wv[k+1] = (int*) omMemDup(origR->wvhdl[k]);
    2702     }
    2703     block0[0] = 1;
    2704     block1[0] = rVar(origR);
    2705     wv[0]=(int*)omAlloc0((rVar(origR) + 1)*sizeof(int));
    2706     for (j=0;j<rVar(origR);j++)
    2707       if (pGetExp(delVar,j+1)!=0) wv[0][j]=1;
    2708     // use this special ordering: like ringorder_a, except that pFDeg, pWeights
    2709     // ignore it
    2710     ord[0] = ringorder_aa;
    2711   }
    2712   // fill in tmp ring to get back the data later on
    2713   tmpR  = rCopy0(origR,FALSE,FALSE); // qring==NULL
    2714   //rUnComplete(tmpR);
    2715   tmpR->p_Procs=NULL;
    2716   tmpR->order = ord;
    2717   tmpR->block0 = block0;
    2718   tmpR->block1 = block1;
    2719   tmpR->wvhdl = wv;
    2720   rComplete(tmpR, 1);
    2721 
    2722 #ifdef HAVE_PLURAL
    2723   /* update nc structure on tmpR */
    2724   if (rIsPluralRing(origR))
    2725   {
    2726     if ( nc_rComplete(origR, tmpR, false) ) // no quotient ideal!
    2727     {
    2728       Werror("no elimination is possible: ordering condition is violated");
    2729       // cleanup
    2730       rDelete(tmpR);
    2731       if (w!=NULL)
    2732         delete w;
    2733       return NULL;
    2734     }
    2735   }
    2736 #endif
    2737   // change into the new ring
    2738   //pChangeRing(pVariables,currRing->OrdSgn,ord,block0,block1,wv);
    2739   rChangeCurrRing(tmpR);
    2740 
    2741   //h = idInit(IDELEMS(h1),h1->rank);
    2742   // fetch data from the old ring
    2743   //for (k=0;k<IDELEMS(h1);k++) h->m[k] = prCopyR( h1->m[k], origR);
    2744   h=idrCopyR(h1,origR,currRing);
    2745   if (origR->qideal!=NULL)
    2746   {
    2747     WarnS("eliminate in q-ring: experimental");
    2748     ideal q=idrCopyR(origR->qideal,origR,currRing);
    2749     ideal s=idSimpleAdd(h,q);
    2750     idDelete(&h);
    2751     idDelete(&q);
    2752     h=s;
    2753   }
    2754   // compute kStd
    2755 #if 1
    2756   //rWrite(tmpR);PrintLn();
    2757   BITSET save=test;
    2758   //test |=1;
    2759   //Print("h: %d gen, rk=%d\n",IDELEMS(h),h->rank);
    2760   //extern char * showOption();
    2761   //Print("%s\n",showOption());
    2762   hh = kStd(h,NULL,hom,&w,hilb);
    2763   test=save;
    2764   idDelete(&h);
    2765 #else
    2766   extern ideal kGroebner(ideal F, ideal Q);
    2767   hh=kGroebner(h,NULL);
    2768 #endif
    2769   // go back to the original ring
    2770   rChangeCurrRing(origR);
    2771   i = IDELEMS(hh)-1;
    2772   while ((i >= 0) && (hh->m[i] == NULL)) i--;
    2773   j = -1;
    2774   // fetch data from temp ring
    2775   for (k=0; k<=i; k++)
    2776   {
    2777     l=pVariables;
    2778     while ((l>0) && (p_GetExp( hh->m[k],l,tmpR)*pGetExp(delVar,l)==0)) l--;
    2779     if (l==0)
    2780     {
    2781       j++;
    2782       if (j >= IDELEMS(h3))
    2783       {
    2784         pEnlargeSet(&(h3->m),IDELEMS(h3),16);
    2785         IDELEMS(h3) += 16;
    2786       }
    2787       h3->m[j] = prMoveR( hh->m[k], tmpR);
    2788       hh->m[k] = NULL;
    2789     }
    2790   }
    2791   id_Delete(&hh, tmpR);
    2792   idSkipZeroes(h3);
    2793   rDelete(tmpR);
    2794   if (w!=NULL)
    2795     delete w;
    2796   return h3;
    27971286}
    27981287
     
    35162005}
    35172006
    3518 /*3
    3519 *handles for some ideal operations the ring/syzcomp managment
    3520 *returns all syzygies (componentwise-)shifted by -syzcomp
    3521 *or -syzcomp-1 (in case of ideals as input)
    3522 static ideal idHandleIdealOp(ideal arg,int syzcomp,int isIdeal=FALSE)
    3523 {
    3524   ring orig_ring=currRing;
    3525   ring syz_ring=rCurrRingAssure_SyzComp();
    3526   rSetSyzComp(length);
    3527 
    3528   ideal s_temp;
    3529   if (orig_ring!=syz_ring)
    3530     s_temp=idrMoveR_NoSort(arg,orig_ring);
    3531   else
    3532     s_temp=arg;
    3533 
    3534   ideal s_temp1 = kStd(s_temp,currQuotient,testHomog,&w,NULL,length);
    3535   if (w!=NULL) delete w;
    3536 
    3537   if (syz_ring!=orig_ring)
    3538   {
    3539     idDelete(&s_temp);
    3540     rChangeCurrRing(orig_ring);
    3541   }
    3542 
    3543   idDelete(&temp);
    3544   ideal temp1=idRingCopy(s_temp1,syz_ring);
    3545 
    3546   if (syz_ring!=orig_ring)
    3547   {
    3548     rChangeCurrRing(syz_ring);
    3549     idDelete(&s_temp1);
    3550     rChangeCurrRing(orig_ring);
    3551     rKill(syz_ring);
    3552   }
    3553 
    3554   for (i=0;i<IDELEMS(temp1);i++)
    3555   {
    3556     if ((temp1->m[i]!=NULL)
    3557     && (pGetComp(temp1->m[i])<=length))
    3558     {
    3559       pDelete(&(temp1->m[i]));
    3560     }
    3561     else
    3562     {
    3563       pShift(&(temp1->m[i]),-length);
    3564     }
    3565   }
    3566   temp1->rank = rk;
    3567   idSkipZeroes(temp1);
    3568 
    3569   return temp1;
    3570 }
    3571 */
    3572 /*2
    3573 * represents (h1+h2)/h2=h1/(h1 intersect h2)
    3574 */
    3575 //ideal idModulo (ideal h2,ideal h1)
    3576 ideal idModulo (ideal h2,ideal h1, tHomog hom, intvec ** w)
    3577 {
    3578   intvec *wtmp=NULL;
    3579 
    3580   int i,j,k,rk,flength=0,slength,length;
    3581   poly p,q;
    3582 
    3583   if (idIs0(h2))
    3584     return idFreeModule(si_max(1,h2->ncols));
    3585   if (!idIs0(h1))
    3586     flength = idRankFreeModule(h1);
    3587   slength = idRankFreeModule(h2);
    3588   length  = si_max(flength,slength);
    3589   if (length==0)
    3590   {
    3591     length = 1;
    3592   }
    3593   ideal temp = idInit(IDELEMS(h2),length+IDELEMS(h2));
    3594   if ((w!=NULL)&&((*w)!=NULL))
    3595   {
    3596     //Print("input weights:");(*w)->show(1);PrintLn();
    3597     int d;
    3598     int k;
    3599     wtmp=new intvec(length+IDELEMS(h2));
    3600     for (i=0;i<length;i++)
    3601       ((*wtmp)[i])=(**w)[i];
    3602     for (i=0;i<IDELEMS(h2);i++)
    3603     {
    3604       poly p=h2->m[i];
    3605       if (p!=NULL)
    3606       {
    3607         d = pDeg(p);
    3608         k= pGetComp(p);
    3609         if (slength>0) k--;
    3610         d +=((**w)[k]);
    3611         ((*wtmp)[i+length]) = d;
    3612       }
    3613     }
    3614     //Print("weights:");wtmp->show(1);PrintLn();
    3615   }
    3616   for (i=0;i<IDELEMS(h2);i++)
    3617   {
    3618     temp->m[i] = pCopy(h2->m[i]);
    3619     q = pOne();
    3620     pSetComp(q,i+1+length);
    3621     pSetmComp(q);
    3622     if(temp->m[i]!=NULL)
    3623     {
    3624       if (slength==0) pShift(&(temp->m[i]),1);
    3625       p = temp->m[i];
    3626       while (pNext(p)!=NULL) pIter(p);
    3627       pNext(p) = q;
    3628     }
    3629     else
    3630       temp->m[i]=q;
    3631   }
    3632   rk = k = IDELEMS(h2);
    3633   if (!idIs0(h1))
    3634   {
    3635     pEnlargeSet(&(temp->m),IDELEMS(temp),IDELEMS(h1));
    3636     IDELEMS(temp) += IDELEMS(h1);
    3637     for (i=0;i<IDELEMS(h1);i++)
    3638     {
    3639       if (h1->m[i]!=NULL)
    3640       {
    3641         temp->m[k] = pCopy(h1->m[i]);
    3642         if (flength==0) pShift(&(temp->m[k]),1);
    3643         k++;
    3644       }
    3645     }
    3646   }
    3647 
    3648   ring orig_ring=currRing;
    3649   ring syz_ring=rCurrRingAssure_SyzComp();
    3650   rSetSyzComp(length);
    3651   ideal s_temp;
    3652 
    3653   if (syz_ring != orig_ring)
    3654   {
    3655     s_temp = idrMoveR_NoSort(temp, orig_ring);
    3656   }
    3657   else
    3658   {
    3659     s_temp = temp;
    3660   }
    3661 
    3662   idTest(s_temp);
    3663   ideal s_temp1 = kStd(s_temp,currQuotient,hom,&wtmp,NULL,length);
    3664 
    3665   //if (wtmp!=NULL)  Print("output weights:");wtmp->show(1);PrintLn();
    3666   if ((w!=NULL) && (*w !=NULL) && (wtmp!=NULL))
    3667   {
    3668     delete *w;
    3669     *w=new intvec(IDELEMS(h2));
    3670     for (i=0;i<IDELEMS(h2);i++)
    3671       ((**w)[i])=(*wtmp)[i+length];
    3672   }
    3673   if (wtmp!=NULL) delete wtmp;
    3674 
    3675   for (i=0;i<IDELEMS(s_temp1);i++)
    3676   {
    3677     if ((s_temp1->m[i]!=NULL)
    3678     && (pGetComp(s_temp1->m[i])<=length))
    3679     {
    3680       pDelete(&(s_temp1->m[i]));
    3681     }
    3682     else
    3683     {
    3684       pShift(&(s_temp1->m[i]),-length);
    3685     }
    3686   }
    3687   s_temp1->rank = rk;
    3688   idSkipZeroes(s_temp1);
    3689 
    3690   if (syz_ring!=orig_ring)
    3691   {
    3692     rChangeCurrRing(orig_ring);
    3693     s_temp1 = idrMoveR_NoSort(s_temp1, syz_ring);
    3694     rKill(syz_ring);
    3695     // Hmm ... here seems to be a memory leak
    3696     // However, simply deleting it causes memory trouble
    3697     // idDelete(&s_temp);
    3698   }
    3699   else
    3700   {
    3701     idDelete(&temp);
    3702   }
    3703   idTest(s_temp1);
    3704   return s_temp1;
    3705 }
    3706 
    37072007int idElem(const ideal F)
    37082008{
Note: See TracChangeset for help on using the changeset viewer.