Changeset f71e8c5 in git for libpolys/polys/simpleideals.cc
- Timestamp:
- Apr 12, 2011, 3:00:18 PM (13 years ago)
- 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
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
libpolys/polys/simpleideals.cc
r1377c9 rf71e8c5 11 11 #include "config.h" 12 12 #include <misc/auxiliary.h> 13 #include <misc/options.h> 14 #include <omalloc/omalloc.h> 15 #include <polys/monomials/p_polys.h> 13 16 14 17 #include "simpleideals.h" … … 44 47 { 45 48 Print("Module of rank %ld,real rank %ld and %d generators.\n", 46 id->rank,id RankFreeModule(id, lmRing, tailRing),IDELEMS(id));49 id->rank,id_RankFreeModule(id, lmRing, tailRing),IDELEMS(id)); 47 50 48 51 int j = (id->ncols*id->nrows) - 1; … … 58 61 /* index of generator with leading term in ground ring (if any); 59 62 otherwise -1 */ 60 int id PosConstant(ideal id)63 int id_PosConstant(ideal id, const ring r) 61 64 { 62 65 int k; 63 66 for (k = IDELEMS(id)-1; k>=0; k--) 64 67 { 65 if (p_LmIsConstantComp(id->m[k], currRing) == TRUE)68 if (p_LmIsConstantComp(id->m[k], r) == TRUE) 66 69 return k; 67 70 } … … 72 75 * initialise the maximal ideal (at 0) 73 76 */ 74 ideal id MaxIdeal (void)77 ideal id_MaxIdeal (const ring r) 75 78 { 76 79 int l; 77 80 ideal hh=NULL; 78 81 79 hh=idInit( pVariables,1);80 for (l=0; l< pVariables; l++)81 { 82 hh->m[l] = p One();83 p SetExp(hh->m[l],l+1,1);84 p Setm(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); 85 88 } 86 89 return hh; … … 174 177 * (Note that the copied polynomials may be zero.) 175 178 */ 176 ideal id CopyFirstK (const ideal ide, const int k)179 ideal id_CopyFirstK (const ideal ide, const int k,const ring r) 177 180 { 178 181 ideal newI = idInit(k, 0); 179 182 for (int i = 0; i < k; i++) 180 newI->m[i] = p Copy(ide->m[i]);183 newI->m[i] = p_Copy(ide->m[i],r); 181 184 return newI; 182 185 } … … 229 232 #else 230 233 if (pComparePolys(id->m[i], id->m[j])) pDelete(&id->m[j]); 231 #endif 234 #endif 232 235 } 233 236 } … … 331 334 #ifdef HAVE_RINGS 332 335 } 333 #endif 336 #endif 334 337 } 335 338 } … … 404 407 * for idSort: compare a and b revlex inclusive module comp. 405 408 */ 406 static int p Comp_RevLex(poly a, poly b,BOOLEAN nolex)409 static int p_Comp_RevLex(poly a, poly b,BOOLEAN nolex, const ring r) 407 410 { 408 411 if (b==NULL) return 1; 409 412 if (a==NULL) return -1; 410 413 411 if (nolex) 414 if (nolex) 412 415 { 413 416 int r=pLmCmp(a,b); … … 418 421 return r; 419 422 } 420 int l= pVariables;423 int l=rVar(r); 421 424 while ((l>0) && (pGetExp(a,l)==pGetExp(b,l))) l--; 422 425 if (l==0) … … 758 761 759 762 /*2 760 *returns a minimized set of generators of h1761 */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 else780 {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 /*2839 763 *the minimal index of used variables - 1 840 764 */ … … 1156 1080 1157 1081 ideal res=idElimination(h,t); 1158 // cleanup 1082 // cleanup 1159 1083 idDelete(&h); 1160 if (res!=NULL)res=idrMoveR(res,r,origRing);1084 res=idrMoveR(res,r,origRing); 1161 1085 rChangeCurrRing(origRing); 1162 1086 rKill(r); 1163 1087 return res; 1164 }1165 /*21166 * h3 := h1 intersect h21167 */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 else1185 {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 else1217 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 else1235 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 else1259 {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 else1285 {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 /*21304 * ideal/module intersection for a list of objects1305 * 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 else1326 {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 else1372 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 else1398 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 /*21418 *computes syzygies of h1,1419 *if quot != NULL it computes in the quotient ring modulo "quot"1420 *works always in a ring with ringorder_s1421 */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 MYTEST1455 #ifdef RDEBUG1456 Print("Prepare::h2: ");1457 idPrint(h2);1458 1459 for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);1460 1461 #endif1462 #endif1463 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 else1476 h2->m[j]=q;1477 }1478 1479 #ifdef PDEBUG1480 for(j=0;j<IDELEMS(h2);j++) pTest(h2->m[j]);1481 1482 #if MYTEST1483 #ifdef RDEBUG1484 Print("Prepare::Input: ");1485 idPrint(h2);1486 1487 Print("Prepare::currQuotient: ");1488 idPrint(currQuotient);1489 #endif1490 #endif1491 1492 #endif1493 1494 idTest(h2);1495 1496 h3 = kStd(h2,currQuotient,hom,w,NULL,syzcomp);1497 1498 #if MYTEST1499 #ifdef RDEBUG1500 Print("Prepare::Output: ");1501 idPrint(h3);1502 for(j=0;j<IDELEMS(h2);j++) pTest(h3->m[j]);1503 #endif1504 #endif1505 1506 1507 idDelete(&h2);1508 return h3;1509 }1510 1511 /*21512 * compute the syzygies of h1 in R/quot,1513 * weights of components are in w1514 * if setRegularity, return the regularity in deg1515 * do not change h1, w1516 */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 PDEBUG1530 for(ii=0;ii<idElemens_h1 /*IDELEMS(h1)*/;ii++) pTest(h1->m[ii]);1531 #endif1532 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 else1559 {1560 s_h1 = h1;1561 }1562 1563 idTest(s_h1);1564 1565 ideal s_h3=idPrepare(s_h1,h,k,w); // main (syz) GB computation1566 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 else1582 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_PLURAL1591 if (rIsPluralRing(currRing))1592 {1593 idDelMultiples(s_h3);1594 idSkipZeroes(s_h3);1595 }1596 #endif1597 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,®, *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 else1648 {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 /*21662 */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 else1680 {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 matrix1709 * in ma1710 */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 else1752 s_h1 = h1;1753 1754 ideal s_h3=idPrepare(s_h1,hi,k,&w); // main (syz) GB computation1755 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_h31764 // T: put in s_h21765 // syz: put in *syz1766 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 == currRing1772 {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 else1783 {1784 pIter(q);1785 }1786 }1787 if (!inputIsIdeal) pShift(&(s_h3->m[j]), -1);1788 }1789 else1790 {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 else1799 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 -11885 pSetComp(q,k+1+j);1886 pSetmComp(q);1887 while (pNext(p)) pIter(p);1888 pNext(p) = q;1889 }1890 }1891 }1892 1893 /*21894 *computes a representation of the generators of submod with respect to those1895 * of mod1896 */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 else1945 {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 else1956 {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 else2012 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 else2019 {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 else2052 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 else2072 {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 else2081 {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 /*22093 *computes division of P by Q with remainder up to (w-weighted) degree n2094 *P, Q, and w are not changed2095 */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 else2104 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 else2116 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 else2127 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 else2132 MATELEM(T,j+1,i+1)=pAdd(MATELEM(T,j+1,i+1),p0);2133 j=IDELEMS(Q)-1;2134 }2135 else2136 {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 else2146 R->m[i]=pAdd(R->m[i],p0);2147 j=IDELEMS(Q)-1;2148 }2149 else2150 j--;2151 }2152 }2153 }2154 }2155 2156 /*22157 *computes the quotient of h1,h2 : internal routine for idQuot2158 *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 else2182 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 else2194 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 else2240 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 /*22267 *computes the quotient of h1,h22268 */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 else2281 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 02302 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 #endif2309 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 else2315 {2316 s_h3 = kStd(s_h4,currQuotient,hom,&weights1,NULL,kmax-1);2317 }2318 test = old_test;2319 #if 02320 // only together with the above debug stuff2321 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 #endif2327 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 else2338 pShift(&s_h3->m[i],-kmax+1);2339 }2340 else2341 pDelete(&s_h3->m[i]);2342 }2343 if (resultIsIdeal)2344 s_h3->rank = 1;2345 else2346 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;2356 1088 } 2357 1089 … … 2552 1284 idSkipZeroes(result); 2553 1285 return result; 2554 }2555 2556 /*22557 * eliminate delVar (product of vars) in h12558 */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_PLURAL2582 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 #endif2596 hom=(tHomog)idHomModule(h1,NULL,&w); //sets w to weight vector or NULL2597 h3=idInit(16,h1->rank);2598 for (k=0;; k++)2599 {2600 if (origR->order[k]!=0) ordersize++;2601 else break;2602 }2603 #if 02604 if (rIsPluralRing(origR)) // we have too keep the odering: it may be needed2605 // for G-algebra2606 {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 else2616 {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 #else2635 #endif2636 if ((hom==TRUE) && (origR->OrdSgn==1) && (!rIsPluralRing(origR)))2637 {2638 #if 12639 // we change to an ordering:2640 // aa(1,1,1,...,0,0,0),wp(...),C2641 // 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, pWeights2651 // ignore it2652 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 else2666 ord[1] = ringorder_dp;2667 #else2668 // we change to an ordering:2669 // a(w1,...wn),wp(1,...0.....),C2670 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 #endif2685 ord[2] = ringorder_C;2686 ord[3] = 0;2687 }2688 else2689 {2690 // we change to an ordering:2691 // aa(....),orig_ordering2692 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, pWeights2709 // ignore it2710 ord[0] = ringorder_aa;2711 }2712 // fill in tmp ring to get back the data later on2713 tmpR = rCopy0(origR,FALSE,FALSE); // qring==NULL2714 //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_PLURAL2723 /* 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 // cleanup2730 rDelete(tmpR);2731 if (w!=NULL)2732 delete w;2733 return NULL;2734 }2735 }2736 #endif2737 // change into the new ring2738 //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 ring2743 //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 kStd2755 #if 12756 //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 #else2766 extern ideal kGroebner(ideal F, ideal Q);2767 hh=kGroebner(h,NULL);2768 #endif2769 // go back to the original ring2770 rChangeCurrRing(origR);2771 i = IDELEMS(hh)-1;2772 while ((i >= 0) && (hh->m[i] == NULL)) i--;2773 j = -1;2774 // fetch data from temp ring2775 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;2797 1286 } 2798 1287 … … 3516 2005 } 3517 2006 3518 /*33519 *handles for some ideal operations the ring/syzcomp managment3520 *returns all syzygies (componentwise-)shifted by -syzcomp3521 *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 else3532 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 else3562 {3563 pShift(&(temp1->m[i]),-length);3564 }3565 }3566 temp1->rank = rk;3567 idSkipZeroes(temp1);3568 3569 return temp1;3570 }3571 */3572 /*23573 * 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 else3630 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 else3658 {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 else3683 {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 leak3696 // However, simply deleting it causes memory trouble3697 // idDelete(&s_temp);3698 }3699 else3700 {3701 idDelete(&temp);3702 }3703 idTest(s_temp1);3704 return s_temp1;3705 }3706 3707 2007 int idElem(const ideal F) 3708 2008 {
Note: See TracChangeset
for help on using the changeset viewer.