Changeset a36e78 in git


Ignore:
Timestamp:
Mar 21, 2007, 4:39:12 PM (17 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'a719bcf0b8dbc648b128303a49777a094b57592c')
Children:
3e5387dc706c674258eaa7f351c895212ad64c0e
Parents:
70ab7345177754da28df06fc80c572286ed3e7ee
Message:
*hannes: format, par2var


git-svn-id: file:///usr/local/Singular/svn/trunk@9954 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/primdec.lib

    r70ab73 ra36e78  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: primdec.lib,v 1.130 2007-03-20 17:31:37 Singular Exp $";
     2version="$Id: primdec.lib,v 1.131 2007-03-21 15:39:12 Singular Exp $";
    33category="Commutative Algebra";
    44info="
     
    226226      g=1;
    227227      verg=laedi;
    228 
    229228      for(j=1;j<=size(fac);j++)
    230229      {
     
    234233        }
    235234      }
    236 
    237235      verg=quotient(laedi,g);
    238236
     
    677675              l[2*i]=maxideal(1);
    678676            }
    679 
    680677            i--;
    681678            break;
     
    954951      if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
    955952      {
    956         primary[2*@k] = primary[2*@k-1];
     953        primary[2*@k]   = primary[2*@k-1];
    957954      }
    958955      else
     
    10761073            {
    10771074              jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
    1078                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
     1075                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
    10791076              jmap2[zz]=primary[2*@k-1][@n];
    10801077              @qht[@n]=var(zz);
     
    11641161      if(size(lres0)==2)
    11651162      {
    1166         lres0[2]=groebner(lres0[2],"par2var");
     1163        lres0[2]=groebner(lres0[2]);
    11671164      }
    11681165      else
     
    11721169          if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1)
    11731170          {
    1174             lres0[2*@n-1]=groebner(lres0[2*@n-1],"par2var");
     1171            lres0[2*@n-1]=groebner(lres0[2*@n-1]);
    11751172            lres0[2*@n]=lres0[2*@n-1];
    11761173            attrib(lres0[2*@n],"isSB",1);
     
    11781175          else
    11791176          {
    1180             lres0[2*@n-1]=groebner(lres0[2*@n-1],"par2var");
    1181             lres0[2*@n]=groebner(lres0[2*@n],"par2var");
     1177            lres0[2*@n-1]=groebner(lres0[2*@n-1]);
     1178            lres0[2*@n]=groebner(lres0[2*@n]);
     1179          }
     1180        }
     1181      }
     1182      primary=primary+lres0;
     1183
     1184//=============================================================
     1185//       if(npars(@P)>0)
     1186//       {
     1187//          @ri= "ring @Phelp ="
     1188//                  +string(char(@P))+",
     1189//                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
     1190//       }
     1191//       else
     1192//       {
     1193//          @ri= "ring @Phelp ="
     1194//                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
     1195//       }
     1196//       execute(@ri);
     1197//       list @lvec;
     1198//       list @lr=imap(@P,lres0);
     1199//       ideal @lr1;
     1200//
     1201//       if(size(@lr)==2)
     1202//       {
     1203//          @lr[2]=homog(@lr[2],@t);
     1204//          @lr1=std(@lr[2]);
     1205//          @lvec[2]=hilb(@lr1,1);
     1206//       }
     1207//       else
     1208//       {
     1209//          for(@n=1;@n<=size(@lr)/2;@n++)
     1210//          {
     1211//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
     1212//             {
     1213//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
     1214//                @lr1=std(@lr[2*@n-1]);
     1215//                @lvec[2*@n-1]=hilb(@lr1,1);
     1216//                @lvec[2*@n]=@lvec[2*@n-1];
     1217//             }
     1218//             else
     1219//             {
     1220//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
     1221//                @lr1=std(@lr[2*@n-1]);
     1222//                @lvec[2*@n-1]=hilb(@lr1,1);
     1223//                @lr[2*@n]=homog(@lr[2*@n],@t);
     1224//                @lr1=std(@lr[2*@n]);
     1225//                @lvec[2*@n]=hilb(@lr1,1);
     1226//
     1227//             }
     1228//         }
     1229//       }
     1230//       @ri= "ring @Phelp1 ="
     1231//                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
     1232//       execute(@ri);
     1233//       list @lr=imap(@Phelp,@lr);
     1234//
     1235//       kill @Phelp;
     1236//       if(size(@lr)==2)
     1237//      {
     1238//          @lr[2]=std(@lr[2],@lvec[2]);
     1239//          @lr[2]=subst(@lr[2],@t,1);
     1240//       }
     1241//       else
     1242//       {
     1243//          for(@n=1;@n<=size(@lr)/2;@n++)
     1244//          {
     1245//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
     1246//             {
     1247//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
     1248//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
     1249//                @lr[2*@n]=@lr[2*@n-1];
     1250//                attrib(@lr[2*@n],"isSB",1);
     1251//             }
     1252//             else
     1253//             {
     1254//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
     1255//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
     1256//                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
     1257//                @lr[2*@n]=subst(@lr[2*@n],@t,1);
     1258//             }
     1259//          }
     1260//       }
     1261//       kill @lvec;
     1262//       setring @P;
     1263//       lres0=imap(@Phelp1,@lr);
     1264//       kill @Phelp1;
     1265//       for(@n=1;@n<=size(lres0);@n++)
     1266//       {
     1267//          lres0[@n]=clearSB(lres0[@n]);
     1268//          attrib(lres0[@n],"isSB",1);
     1269//       }
     1270//
     1271//       primary[2*@k-1]=lres0[1];
     1272//       primary[2*@k]=lres0[2];
     1273//       @s=size(primary)/2;
     1274//       for(@n=1;@n<=size(lres0)/2-1;@n++)
     1275//       {
     1276//         primary[2*@s+2*@n-1]=lres0[2*@n+1];
     1277//         primary[2*@s+2*@n]=lres0[2*@n+2];
     1278//       }
     1279//       @k--;
     1280//=============================================================
     1281    }
     1282  }
     1283  return(primary);
     1284}
     1285example
     1286{ "EXAMPLE:"; echo = 2;
     1287   ring  r = 0,(x,y,z),lp;
     1288   poly  p = z2+1;
     1289   poly  q = z4+2;
     1290   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     1291   i=std(i);
     1292   list  pr= zero_decomp(i,ideal(0),0);
     1293   pr;
     1294}
     1295///////////////////////////////////////////////////////////////////////////////
     1296proc extF(list l,list #)
     1297{
     1298//zero_dimensional primary decomposition after finite field extension
     1299  def R=basering;
     1300  int p=char(R);
     1301
     1302  if((p==0)||(p>13)||(npars(R)>0)){return(l);}
     1303
     1304  int ex=3;
     1305  if(size(#)>0){ex=#[1];}
     1306
     1307  list peek,peek1;
     1308  while(size(l)>0)
     1309  {
     1310    if(size(l[2])==0)
     1311    {
     1312      peek[size(peek)+1]=l[1];
     1313    }
     1314    else
     1315    {
     1316      peek1[size(peek1)+1]=l[1];
     1317      peek1[size(peek1)+1]=l[2];
     1318    }
     1319    l=delete(l,1);
     1320    l=delete(l,1);
     1321  }
     1322  if(size(peek)==0){return(peek1);}
     1323
     1324  string gnir="ring RH=("+string(p)+"^"+string(ex)+",a),("+varstr(R)+"),lp;";
     1325  execute(gnir);
     1326  string mp="minpoly="+string(minpoly)+";";
     1327  gnir="ring RL=("+string(p)+",a),("+varstr(R)+"),lp;";
     1328  execute(gnir);
     1329  execute(mp);
     1330  list L=imap(R,peek);
     1331  list pr, keep;
     1332  int i;
     1333  for(i=1;i<=size(L);i++)
     1334  {
     1335    attrib(L[i],"isSB",1);
     1336    pr=zero_decomp(L[i],0,0);
     1337    keep=keep+pr;
     1338  }
     1339  for(i=1;i<=size(keep);i++)
     1340  {
     1341    keep[i]=simplify(keep[i],1);
     1342  }
     1343  mp="poly pp="+string(minpoly)+";";
     1344
     1345  string gnir1="ring RS="+string(p)+",("+varstr(R)+",a),lp;";
     1346  execute(gnir1);
     1347  execute(mp);
     1348  list L=imap(RL,keep);
     1349
     1350  for(i=1;i<=size(L);i++)
     1351  {
     1352    L[i]=eliminate(L[i]+ideal(pp),a);
     1353  }
     1354  i=0;
     1355  int j;
     1356  while(i<size(L)/2-1)
     1357  {
     1358    i++;
     1359    j=i;
     1360    while(j<size(L)/2)
     1361    {
     1362      j++;
     1363      if(idealsEqual(L[2*i-1],L[2*j-1]))
     1364      {
     1365        L=delete(L,2*j-1);
     1366        L=delete(L,2*j-1);
     1367        j--;
     1368      }
     1369    }
     1370  }
     1371  setring R;
     1372  list re=imap(RS,L);
     1373  re=re+peek1;
     1374
     1375  return(extF(re,ex+1));
     1376}
     1377
     1378///////////////////////////////////////////////////////////////////////////////
     1379proc zeroSp(ideal i)
     1380{
     1381//preparation for the separable closure
     1382//decomposition into ideals of special type
     1383//i.e. the minimal polynomials of every variable mod i are irreducible
     1384//returns a list of 2 lists: rr=pe,qe
     1385//the ideals in pe[l] are special, their special elements are in qe[l]
     1386//pe[l] is a dp-Groebnerbasis
     1387//the radical of the intersection of the pe[l] is equal to the radical of i
     1388
     1389  def R=basering;
     1390
     1391  //i has to be a reduced groebner basis
     1392  ideal F=finduni(i);
     1393
     1394  int j,k,l,ready;
     1395  list fa;
     1396  fa[1]=factorize(F[1],1);
     1397  poly te,ti;
     1398  ideal tj;
     1399  //avoid factorization of the same polynomial
     1400  for(j=2;j<=size(F);j++)
     1401  {
     1402    for(k=1;k<=j-1;k++)
     1403    {
     1404      ti=F[k];
     1405      te=subst(ti,var(k),var(j));
     1406      if(te==F[j])
     1407      {
     1408        tj=fa[k];
     1409        fa[j]=subst(tj,var(k),var(j));
     1410        ready=1;
     1411        break;
     1412      }
     1413    }
     1414    if(!ready)
     1415    {
     1416      fa[j]=factorize(F[j],1);
     1417    }
     1418    ready=0;
     1419  }
     1420  execute( "ring P=("+charstr(R)+"),("+varstr(R)+"),(C,dp);");
     1421  ideal i=imap(R,i);
     1422  if(npars(basering)==0)
     1423  {
     1424    ideal J=fglm(R,i);
     1425  }
     1426  else
     1427  {
     1428    ideal J=groebner(i);
     1429  }
     1430  list fa=imap(R,fa);
     1431  list qe=J;          //collects a dp-Groebnerbasis of the special ideals
     1432  list keep=ideal(0); //collects the special elements
     1433
     1434  list re,em,ke;
     1435  ideal K,L;
     1436
     1437  for(j=1;j<=nvars(basering);j++)
     1438  {
     1439    for(l=1;l<=size(qe);l++)
     1440    {
     1441      for(k=1;k<=size(fa[j]);k++)
     1442      {
     1443        L=std(qe[l],fa[j][k]);
     1444        K=keep[l],fa[j][k];
     1445        if(deg(L[1])>0)
     1446        {
     1447          re[size(re)+1]=L;
     1448          ke[size(ke)+1]=K;
     1449        }
     1450      }
     1451    }
     1452    qe=re;
     1453    re=em;
     1454    keep=ke;
     1455    ke=em;
     1456  }
     1457
     1458  setring R;
     1459  list qe=imap(P,keep);
     1460  list pe=imap(P,qe);
     1461  for(l=1;l<=size(qe);l++)
     1462  {
     1463    qe[l]=simplify(qe[l],2);
     1464  }
     1465  list rr=pe,qe;
     1466  return(rr);
     1467}
     1468///////////////////////////////////////////////////////////////////////////////
     1469
     1470proc zeroSepClos(ideal I,ideal F)
     1471{
     1472//computes the separable closure of the special ideal I
     1473//F is the set of special elements of I
     1474//returns the separable closure sc(I) of I and an intvec v
     1475//such that sc(I)=preimage(frobenius definde by v)
     1476//i.e. var(i)----->var(i)^(p^v[i])
     1477
     1478  if(homog(I)==1){return(maxideal(1));}
     1479
     1480  //assume F[i] irreducible in I and depending only on var(i)
     1481
     1482  def R=basering;
     1483  int n=nvars(R);
     1484  int p=char(R);
     1485  intvec v;
     1486  v[n]=0;
     1487  int i,k;
     1488  list l;
     1489
     1490  for(i=1;i<=n;i++)
     1491  {
     1492    l[i]=sep(F[i],i);
     1493    F[i]=l[i][1];
     1494    if(l[i][2]>k){k=l[i][2];}
     1495  }
     1496
     1497  if(k==0){return(list(I,v));}        //the separable case
     1498  ideal m;
     1499
     1500  for(i=1;i<=n;i++)
     1501  {
     1502    m[i]=var(i)^(p^l[i][2]);
     1503    v[i]=l[i][2];
     1504  }
     1505  map phi=R,m;
     1506  ideal J=preimage(R,phi,I);
     1507  return(list(J,v));
     1508}
     1509///////////////////////////////////////////////////////////////////////////////
     1510
     1511proc insepDecomp(ideal i)
     1512{
     1513//decomposes i into special ideals
     1514//computes the prime decomposition of the special ideals
     1515//and transforms it back to a decomposition of i
     1516
     1517  def R=basering;
     1518  list pr=zeroSp(i);
     1519  int l,k;
     1520  list re,wo,qr;
     1521  ideal m=maxideal(1);
     1522  ideal K;
     1523  map phi=R,m;
     1524  int p=char(R);
     1525  intvec op=option(get);
     1526
     1527  for(l=1;l<=size(pr[1]);l++)
     1528  {
     1529    wo=zeroSepClos(pr[1][l],pr[2][l]);
     1530    for(k=1;k<=nvars(basering);k++)
     1531    {
     1532      m[k]=var(k)^(p^wo[2][k]);
     1533    }
     1534    phi=R,m;
     1535    qr=decomp(wo[1],2);
     1536
     1537    option(redSB);
     1538    for(k=1;k<=size(qr)/2;k++)
     1539    {
     1540      K=qr[2*k];
     1541      K=phi(K);
     1542      K=groebner(K);
     1543      re[size(re)+1]=zeroRad(K);
     1544    }
     1545    option(noredSB);
     1546  }
     1547  option(set,op);
     1548  return(re);
     1549}
     1550
     1551
     1552///////////////////////////////////////////////////////////////////////////////
     1553
     1554static proc clearSB (ideal i,list #)
     1555"USAGE:   clearSB(i); i ideal which is SB ordered by monomial ordering
     1556RETURN:  ideal = minimal SB
     1557NOTE:
     1558EXAMPLE: example clearSB; shows an example
     1559"
     1560{
     1561  int k,j;
     1562  poly m;
     1563  int c=size(i);
     1564
     1565  if(size(#)==0)
     1566  {
     1567    for(j=1;j<c;j++)
     1568    {
     1569      if(deg(i[j])==0)
     1570      {
     1571        i=ideal(1);
     1572        return(i);
     1573      }
     1574      if(deg(i[j])>0)
     1575      {
     1576        m=lead(i[j]);
     1577        for(k=j+1;k<=c;k++)
     1578        {
     1579          if(size(lead(i[k])/m)>0)
     1580          {
     1581            i[k]=0;
     1582          }
     1583        }
     1584      }
     1585    }
     1586  }
     1587  else
     1588  {
     1589    j=0;
     1590    while(j<c-1)
     1591    {
     1592      j++;
     1593      if(deg(i[j])==0)
     1594      {
     1595        i=ideal(1);
     1596        return(i);
     1597      }
     1598      if(deg(i[j])>0)
     1599      {
     1600        m=lead(i[j]);
     1601        for(k=j+1;k<=c;k++)
     1602        {
     1603          if(size(lead(i[k])/m)>0)
     1604          {
     1605            if((leadexp(m)!=leadexp(i[k]))||(#[j]<=#[k]))
     1606            {
     1607              i[k]=0;
     1608            }
     1609            else
     1610            {
     1611              i[j]=0;
     1612              break;
     1613            }
     1614          }
     1615        }
     1616      }
     1617    }
     1618  }
     1619  return(simplify(i,2));
     1620}
     1621example
     1622{ "EXAMPLE:"; echo = 2;
     1623   ring  r = (0,a,b),(x,y,z),dp;
     1624   ideal i=ax2+y,a2x+y,bx;
     1625   list l=1,2,1;
     1626   ideal j=clearSB(i,l);
     1627   j;
     1628}
     1629
     1630///////////////////////////////////////////////////////////////////////////////
     1631static proc clearSBNeu (ideal i,list #)
     1632"USAGE:   clearSB(i); i ideal which is SB ordered by monomial ordering
     1633RETURN:  ideal = minimal SB
     1634NOTE:
     1635EXAMPLE: example clearSB; shows an example
     1636"
     1637{
     1638 int k,j;
     1639 intvec m,n,v,w;
     1640 int c=size(i);
     1641 w=leadexp(0);
     1642 v[size(i)]=0;
     1643
     1644 j=0;
     1645 while(j<c-1)
     1646 {
     1647   j++;
     1648   if(deg(i[j])>=0)
     1649   {
     1650      m=leadexp(i[j]);
     1651      for(k=j+1;k<=c;k++)
     1652      {
     1653        n=leadexp(i[k]);
     1654        if(n!=w)
     1655        {
     1656           if(((m==n)&&(#[j]>#[k]))||((teilt(n,m))&&(n!=m)))
     1657           {
     1658             i[j]=0;
     1659             v[j]=1;
     1660             break;
     1661           }
     1662           if(((m==n)&&(#[j]<=#[k]))||((teilt(m,n))&&(n!=m)))
     1663           {
     1664             i[k]=0;
     1665             v[k]=1;
     1666           }
     1667        }
     1668      }
     1669    }
     1670  }
     1671  return(v);
     1672}
     1673
     1674static proc teilt(intvec a, intvec b)
     1675{
     1676  int i;
     1677  for(i=1;i<=size(a);i++)
     1678  {
     1679    if(a[i]>b[i]){return(0);}
     1680  }
     1681  return(1);
     1682}
     1683///////////////////////////////////////////////////////////////////////////////
     1684
     1685static proc independSet (ideal j)
     1686"USAGE:   independentSet(i); i ideal
     1687RETURN:  list = new varstring with the independent set at the end,
     1688                ordstring with the corresponding block ordering,
     1689                the integer where the independent set starts in the varstring
     1690NOTE:
     1691EXAMPLE: example independentSet; shows an example
     1692"
     1693{
     1694  int n,k,di;
     1695  list resu,hilf;
     1696  string var1,var2;
     1697  list v=indepSet(j,1);
     1698
     1699  for(n=1;n<=size(v);n++)
     1700  {
     1701    di=0;
     1702    var1="";
     1703    var2="";
     1704    for(k=1;k<=size(v[n]);k++)
     1705    {
     1706      if(v[n][k]!=0)
     1707      {
     1708        di++;
     1709        var2=var2+"var("+string(k)+"),";
     1710      }
     1711      else
     1712      {
     1713        var1=var1+"var("+string(k)+"),";
     1714      }
     1715    }
     1716    if(di>0)
     1717    {
     1718      var1=var1+var2;
     1719      var1=var1[1..size(var1)-1];
     1720      hilf[1]=var1;
     1721      hilf[2]="lp";
     1722      //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
     1723      hilf[3]=di;
     1724      resu[n]=hilf;
     1725    }
     1726    else
     1727    {
     1728      resu[n]=varstr(basering),ordstr(basering),0;
     1729    }
     1730  }
     1731  return(resu);
     1732}
     1733example
     1734{ "EXAMPLE:"; echo = 2;
     1735   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
     1736   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
     1737   i=std(i);
     1738   list  l=independSet(i);
     1739   l;
     1740   i=i,g;
     1741   l=independSet(i);
     1742   l;
     1743
     1744   ring s=0,(x,y,z),lp;
     1745   ideal i=z,yx;
     1746   list l=independSet(i);
     1747   l;
     1748
     1749
     1750}
     1751///////////////////////////////////////////////////////////////////////////////
     1752
     1753static proc maxIndependSet (ideal j)
     1754"USAGE:   maxIndependentSet(i); i ideal
     1755RETURN:  list = new varstring with the maximal independent set at the end,
     1756                ordstring with the corresponding block ordering,
     1757                the integer where the independent set starts in the varstring
     1758NOTE:
     1759EXAMPLE: example maxIndependentSet; shows an example
     1760"
     1761{
     1762  int n,k,di;
     1763  list resu,hilf;
     1764  string var1,var2;
     1765  list v=indepSet(j,0);
     1766
     1767  for(n=1;n<=size(v);n++)
     1768  {
     1769    di=0;
     1770    var1="";
     1771    var2="";
     1772    for(k=1;k<=size(v[n]);k++)
     1773    {
     1774      if(v[n][k]!=0)
     1775      {
     1776        di++;
     1777        var2=var2+"var("+string(k)+"),";
     1778      }
     1779      else
     1780      {
     1781        var1=var1+"var("+string(k)+"),";
     1782      }
     1783    }
     1784    if(di>0)
     1785    {
     1786      var1=var1+var2;
     1787      var1=var1[1..size(var1)-1];
     1788      hilf[1]=var1;
     1789      hilf[2]="lp";
     1790      hilf[3]=di;
     1791      resu[n]=hilf;
     1792    }
     1793    else
     1794    {
     1795      resu[n]=varstr(basering),ordstr(basering),0;
     1796    }
     1797  }
     1798  return(resu);
     1799}
     1800example
     1801{ "EXAMPLE:"; echo = 2;
     1802   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
     1803   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
     1804   i=std(i);
     1805   list  l=maxIndependSet(i);
     1806   l;
     1807   i=i,g;
     1808   l=maxIndependSet(i);
     1809   l;
     1810
     1811   ring s=0,(x,y,z),lp;
     1812   ideal i=z,yx;
     1813   list l=maxIndependSet(i);
     1814   l;
     1815
     1816
     1817}
     1818
     1819///////////////////////////////////////////////////////////////////////////////
     1820
     1821static proc prepareQuotientring (int nnp)
     1822"USAGE:   prepareQuotientring(nnp); nnp int
     1823RETURN:  string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
     1824NOTE:
     1825EXAMPLE: example independentSet; shows an example
     1826"
     1827{
     1828  ideal @ih,@jh;
     1829  int npar=npars(basering);
     1830  int @n;
     1831
     1832  string quotring= "ring quring = ("+charstr(basering);
     1833  for(@n=nnp+1;@n<=nvars(basering);@n++)
     1834  {
     1835     quotring=quotring+",var("+string(@n)+")";
     1836     @ih=@ih+var(@n);
     1837  }
     1838
     1839  quotring=quotring+"),(var(1)";
     1840  @jh=@jh+var(1);
     1841  for(@n=2;@n<=nnp;@n++)
     1842  {
     1843    quotring=quotring+",var("+string(@n)+")";
     1844    @jh=@jh+var(@n);
     1845  }
     1846  quotring=quotring+"),(C,lp);";
     1847
     1848  return(quotring);
     1849
     1850}
     1851example
     1852{ "EXAMPLE:"; echo = 2;
     1853   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
     1854   def @Q=basering;
     1855   list l= prepareQuotientring(3);
     1856   l;
     1857   execute(l[1]);
     1858   execute(l[2]);
     1859   basering;
     1860   phi;
     1861   setring @Q;
     1862
     1863}
     1864
     1865///////////////////////////////////////////////////////////////////////////////
     1866static proc cleanPrimary(list l)
     1867{
     1868   int i,j;
     1869   list lh;
     1870   for(i=1;i<=size(l)/2;i++)
     1871   {
     1872      if(deg(l[2*i-1][1])>0)
     1873      {
     1874         j++;
     1875         lh[j]=l[2*i-1];
     1876         j++;
     1877         lh[j]=l[2*i];
     1878      }
     1879   }
     1880   return(lh);
     1881}
     1882///////////////////////////////////////////////////////////////////////////////
     1883
     1884
     1885proc minAssPrimesold(ideal i, list #)
     1886"USAGE:   minAssPrimes(i); i ideal
     1887         minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
     1888RETURN:  list = the minimal associated prime ideals of i
     1889EXAMPLE: example minAssPrimes; shows an example
     1890"
     1891{
     1892   def @P=basering;
     1893   if(size(i)==0){return(list(ideal(0)));}
     1894   list qr=simplifyIdeal(i);
     1895   map phi=@P,qr[2];
     1896   i=qr[1];
     1897
     1898   execute ("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
     1899             +ordstr(basering)+");");
     1900
     1901
     1902   ideal i=fetch(@P,i);
     1903   if(size(#)==0)
     1904   {
     1905      int @wr;
     1906      list tluser,@res;
     1907      list primary=decomp(i,2);
     1908
     1909      @res[1]=primary;
     1910
     1911      tluser=union(@res);
     1912      setring @P;
     1913      list @res=imap(gnir,tluser);
     1914      return(phi(@res));
     1915   }
     1916   list @res,empty;
     1917   ideal ser;
     1918   option(redSB);
     1919   list @pr=facstd(i);
     1920   //if(size(@pr)==1)
     1921//   {
     1922//      attrib(@pr[1],"isSB",1);
     1923//      if((dim(@pr[1])==0)&&(homog(@pr[1])==1))
     1924//      {
     1925//         setring @P;
     1926//         list @res=maxideal(1);
     1927//         return(phi(@res));
     1928//      }
     1929//      if(dim(@pr[1])>1)
     1930//      {
     1931//         setring @P;
     1932//        // kill gnir;
     1933//         execute ("ring gnir1 = ("+charstr(basering)+"),
     1934//                              ("+varstr(basering)+"),(C,lp);");
     1935//         ideal i=fetch(@P,i);
     1936//         list @pr=facstd(i);
     1937//        // ideal ser;
     1938//         setring gnir;
     1939//         @pr=fetch(gnir1,@pr);
     1940//         kill gnir1;
     1941//      }
     1942//   }
     1943    option(noredSB);
     1944   int j,k,odim,ndim,count;
     1945   attrib(@pr[1],"isSB",1);
     1946   if(#[1]==77)
     1947   {
     1948     odim=dim(@pr[1]);
     1949     count=1;
     1950     intvec pos;
     1951     pos[size(@pr)]=0;
     1952     for(j=2;j<=size(@pr);j++)
     1953     {
     1954        attrib(@pr[j],"isSB",1);
     1955        ndim=dim(@pr[j]);
     1956        if(ndim>odim)
     1957        {
     1958           for(k=count;k<=j-1;k++)
     1959           {
     1960              pos[k]=1;
     1961           }
     1962           count=j;
     1963           odim=ndim;
     1964        }
     1965        if(ndim<odim)
     1966        {
     1967           pos[j]=1;
     1968        }
     1969     }
     1970     for(j=1;j<=size(@pr);j++)
     1971     {
     1972        if(pos[j]!=1)
     1973        {
     1974            @res[j]=decomp(@pr[j],2);
     1975        }
     1976        else
     1977        {
     1978           @res[j]=empty;
     1979        }
     1980     }
     1981   }
     1982   else
     1983   {
     1984     ser=ideal(1);
     1985     for(j=1;j<=size(@pr);j++)
     1986     {
     1987//@pr[j];
     1988//pause();
     1989        @res[j]=decomp(@pr[j],2);
     1990//       @res[j]=decomp(@pr[j],2,@pr[j],ser);
     1991//       for(k=1;k<=size(@res[j]);k++)
     1992//       {
     1993//          ser=intersect(ser,@res[j][k]);
     1994//       }
     1995     }
     1996   }
     1997
     1998   @res=union(@res);
     1999   setring @P;
     2000   list @res=imap(gnir,@res);
     2001   return(phi(@res));
     2002}
     2003example
     2004{ "EXAMPLE:"; echo = 2;
     2005   ring  r = 32003,(x,y,z),lp;
     2006   poly  p = z2+1;
     2007   poly  q = z4+2;
     2008   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     2009   list pr= minAssPrimes(i);  pr;
     2010
     2011   minAssPrimes(i,1);
     2012}
     2013
     2014static proc primT(ideal i)
     2015{
     2016   //assumes that all generators of i are irreducible
     2017   //i is standard basis
     2018
     2019   attrib(i,"isSB",1);
     2020   int j=size(i);
     2021   int k;
     2022   while(j>0)
     2023   {
     2024     if(deg(i[j])>1){break;}
     2025     j--;
     2026   }
     2027   if(j==0){return(1);}
     2028   if(deg(i[j])==vdim(i)){return(1);}
     2029   return(0);
     2030}
     2031
     2032static proc minAssPrimes(ideal i, list #)
     2033"USAGE:   minAssPrimes(i); i ideal
     2034      Optional parameters in list #: (can be entered in any order)
     2035      0, "facstd"   ->   uses facstd to first decompose the ideal
     2036      1, "noFacstd" ->  does not use facstd (default)
     2037      "SL" ->     the new algorithm is used (default)
     2038      "GTZ" ->     the old algorithm is used
     2039RETURN:  list = the minimal associated prime ideals of i
     2040EXAMPLE: example minAssPrimes; shows an example
     2041"
     2042{
     2043  if(size(i) == 0){return(list(ideal(0)));}
     2044  string algorithm;    // Algorithm to be used
     2045  string facstdOption;    // To uses proc facstd
     2046  int j;          // Counter
     2047  def P0 = basering;
     2048  list Pl=ringlist(P0);
     2049  intvec dp_w;
     2050  for(j=nvars(P0);j>0;j--) {dp_w[j]=1;}
     2051  Pl[3]=list(list("dp",dp_w),list("C",0));
     2052  def P=ring(Pl);
     2053  setring P;
     2054  ideal i=imap(P0,i);
     2055
     2056  // Set input parameters
     2057  algorithm = "SL";         // Default: SL algorithm
     2058  facstdOption = "noFacstd";    // Default: facstd is not used
     2059  if(size(#) > 0)
     2060  {
     2061    int valid;
     2062    for(j = 1; j <= size(#); j++)
     2063    {
     2064      valid = 0;
     2065      if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
     2066      {
     2067        if (#[j] == 0) {facstdOption = "noFacstd"; valid = 1;}    // If #[j] == 0, facstd is not used.
     2068        if (#[j] == 1) {facstdOption = "facstd";   valid = 1;}    // If #[j] == 1, facstd is used.
     2069      }
     2070      if(typeof(#[j]) == "string")
     2071      {
     2072        if(#[j] == "GTZ" || #[j] == "SL")
     2073        {
     2074          algorithm = #[j];
     2075          valid = 1;
     2076        }
     2077        if(#[j] == "noFacstd" || #[j] == "facstd")
     2078        {
     2079          facstdOption = #[j];
     2080          valid = 1;
     2081        }
     2082      }
     2083      if(valid == 0)
     2084      {
     2085        dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
     2086      }
     2087    }
     2088  }
     2089
     2090  list q = simplifyIdeal(i);
     2091  list re = maxideal(1);
     2092  int a, k;
     2093  intvec op = option(get);
     2094  map phi = P,q[2];
     2095
     2096  list result;
     2097
     2098  if(npars(P) == 0){option(redSB);}
     2099
     2100  if(attrib(i,"isSB")!=1)
     2101  {
     2102    i=groebner(q[1]);
     2103  }
     2104  else
     2105  {
     2106    for(j=1;j<=nvars(basering);j++)
     2107    {
     2108      if(q[2][j]!=var(j)){k=1;break;}
     2109    }
     2110    if(k)
     2111    {
     2112      i=groebner(q[1]);
     2113    }
     2114  }
     2115
     2116  if(dim(i) == -1){setring P0;return(ideal(1));}
     2117  if((dim(i) == 0) && (npars(P) == 0))
     2118  {
     2119    int di = vdim(i);
     2120    execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
     2121    ideal J = std(imap(P,i));
     2122    attrib(J, "isSB", 1);
     2123    if(vdim(J) != di)
     2124    {
     2125      J = fglm(P, i);
     2126    }
     2127    list pr = triangMH(J, 2);
     2128    list qr, re;
     2129
     2130    for(k = 1; k <= size(pr); k++)
     2131    {
     2132      if(primT(pr[k]))
     2133      {
     2134        re[size(re) + 1] = pr[k];
     2135      }
     2136      else
     2137      {
     2138        attrib(pr[k], "isSB", 1);
     2139        // Lines changed
     2140        if (algorithm == "GTZ")
     2141        {
     2142          qr = decomp(pr[k], 2);
     2143        }
     2144        else
     2145        {
     2146          qr = minAssSL(pr[k]);
     2147        }
     2148        for(j = 1; j <= size(qr) / 2; j++)
     2149        {
     2150          re[size(re) + 1] = qr[2 * j];
     2151        }
     2152      }
     2153    }
     2154    setring P;
     2155    re = imap(gnir, re);
     2156    re=phi(re);
     2157    option(set, op);
     2158    setring(P0);
     2159    list re=imap(P,re);
     2160    return(re);
     2161  }
     2162
     2163  // Lines changed
     2164  if ((facstdOption == "noFacstd") || (dim(i) == 0))
     2165  {
     2166    if (algorithm == "GTZ")
     2167    {
     2168      re[1] = decomp(i, 2);
     2169    }
     2170    else
     2171    {
     2172      re[1] = minAssSL(i);
     2173    }
     2174    re = union(re);
     2175    option(set, op);
     2176    re=phi(re);
     2177    setring(P0);
     2178    list re=imap(P,re);
     2179    return(re);
     2180  }
     2181  q = facstd(i);
     2182
     2183/*
     2184  if((size(q) == 1) && (dim(i) > 1))
     2185  {
     2186    execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
     2187    list p = facstd(fetch(P, i));
     2188    if(size(p) > 1)
     2189    {
     2190      a = 1;
     2191      setring P;
     2192      q = fetch(gnir,p);
     2193    }
     2194    else
     2195    {
     2196      setring P;
     2197    }
     2198    kill gnir;
     2199  }
     2200*/
     2201  option(set,op);
     2202  // Debug
     2203  dbprint(printlevel - voice, "Components returned by facstd", size(q), q);
     2204  for(j = 1; j <= size(q); j++)
     2205  {
     2206    if(a == 0){attrib(q[j], "isSB", 1);}
     2207    // Debug
     2208    dbprint(printlevel - voice, "We compute the decomp of component", j);
     2209    // Lines changed
     2210    if (algorithm == "GTZ")
     2211    {
     2212      re[j] = decomp(q[j], 2);
     2213    }
     2214    else
     2215    {
     2216      re[j] = minAssSL(q[j]);
     2217    }
     2218    // Debug
     2219    dbprint(printlevel - voice, "Number of components obtained for this component:", size(re[j]) / 2);
     2220    dbprint(printlevel - voice, "re[j]:", re[j]);
     2221  }
     2222  re = union(re);
     2223  re=phi(re);
     2224  setring(P0);
     2225  list re=imap(P,re);
     2226  return(re);
     2227}
     2228example
     2229{ "EXAMPLE:"; echo = 2;
     2230   ring  r = 32003,(x,y,z),lp;
     2231   poly  p = z2+1;
     2232   poly  q = z4+2;
     2233   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     2234   list pr= minAssPrimes(i);  pr;
     2235
     2236   minAssPrimes(i,1);
     2237}
     2238
     2239static proc union(list li)
     2240{
     2241  int i,j,k;
     2242
     2243  def P=basering;
     2244
     2245  execute("ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);");
     2246  list l=fetch(P,li);
     2247  list @erg;
     2248
     2249  for(k=1;k<=size(l);k++)
     2250  {
     2251     for(j=1;j<=size(l[k])/2;j++)
     2252     {
     2253        if(deg(l[k][2*j][1])!=0)
     2254        {
     2255           i++;
     2256           @erg[i]=l[k][2*j];
     2257        }
     2258     }
     2259  }
     2260
     2261  list @wos;
     2262  i=0;
     2263  ideal i1,i2;
     2264  while(i<size(@erg)-1)
     2265  {
     2266     i++;
     2267     k=i+1;
     2268     i1=lead(@erg[i]);
     2269      attrib(i1,"isSB",1);
     2270      attrib(@erg[i],"isSB",1);
     2271
     2272     while(k<=size(@erg))
     2273     {
     2274        if(deg(@erg[i][1])==0)
     2275        {
     2276           break;
     2277        }
     2278        i2=lead(@erg[k]);
     2279        attrib(@erg[k],"isSB",1);
     2280        attrib(i2,"isSB",1);
     2281
     2282        if(size(reduce(i1,i2,1))==0)
     2283        {
     2284           if(size(reduce(@erg[i],@erg[k],1))==0)
     2285           {
     2286              @erg[k]=ideal(1);
     2287              i2=ideal(1);
     2288           }
     2289        }
     2290        if(size(reduce(i2,i1,1))==0)
     2291        {
     2292           if(size(reduce(@erg[k],@erg[i],1))==0)
     2293           {
     2294              break;
     2295           }
     2296        }
     2297        k++;
     2298        if(k>size(@erg))
     2299        {
     2300           @wos[size(@wos)+1]=@erg[i];
     2301        }
     2302     }
     2303  }
     2304  if(deg(@erg[size(@erg)][1])!=0)
     2305  {
     2306     @wos[size(@wos)+1]=@erg[size(@erg)];
     2307  }
     2308  setring P;
     2309  list @ser=fetch(ir,@wos);
     2310  return(@ser);
     2311}
     2312///////////////////////////////////////////////////////////////////////////////
     2313proc equidim(ideal i,list #)
     2314"USAGE:  equidim(i) or equidim(i,1) ; i ideal
     2315RETURN: list of equidimensional ideals a[1],...,a[s] with:
     2316        - a[s] the equidimensional locus of i, i.e. the intersection
     2317          of the primary ideals of dimension of i
     2318        - a[1],...,a[s-1] the lower dimensional equidimensional loci.
     2319NOTE:    An embedded component q (primary ideal) of i can be replaced in the
     2320         decomposition by a primary ideal q1 with the same radical as q. @*
     2321         @code{equidim(i,1)} uses the algorithm of Eisenbud/Huneke/Vasconcelos.
     2322
     2323EXAMPLE:example equidim; shows an example
     2324"
     2325{
     2326  if(attrib(basering,"global")!=1)
     2327  {
     2328      ERROR(
     2329      "// Not implemented for this ordering, please change to global ordering."
     2330      );
     2331  }
     2332  intvec op ;
     2333  def  P = basering;
     2334  list eq;
     2335  intvec w;
     2336  int n,m;
     2337  int g=size(i);
     2338  int a=attrib(i,"isSB");
     2339  int homo=homog(i);
     2340  if(size(#)!=0)
     2341  {
     2342     m=1;
     2343  }
     2344
     2345  if(((homo==1)||(a==1))&&(find(ordstr(basering),"l")==0)
     2346                                &&(find(ordstr(basering),"s")==0))
     2347  {
     2348     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
     2349                              +ordstr(basering)+");");
     2350     ideal i=imap(P,i);
     2351     ideal j=i;
     2352     if(a==1)
     2353     {
     2354       attrib(j,"isSB",1);
     2355     }
     2356     else
     2357     {
     2358       j=groebner(i);
     2359     }
     2360  }
     2361  else
     2362  {
     2363     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;");
     2364     ideal i=imap(P,i);
     2365     ideal j=groebner(i);
     2366  }
     2367  if(homo==1)
     2368  {
     2369     for(n=1;n<=nvars(basering);n++)
     2370     {
     2371        w[n]=ord(var(n));
     2372     }
     2373     intvec hil=hilb(j,1,w);
     2374  }
     2375
     2376  if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1)
     2377                  ||(dim(j)==0)||(dim(j)+g==nvars(basering)))
     2378  {
     2379    setring P;
     2380    eq[1]=i;
     2381    return(eq);
     2382  }
     2383
     2384  if(m==0)
     2385  {
     2386     ideal k=equidimMax(j);
     2387  }
     2388  else
     2389  {
     2390     ideal k=equidimMaxEHV(j);
     2391  }
     2392  if(size(reduce(k,j,1))==0)
     2393  {
     2394    setring P;
     2395    eq[1]=i;
     2396    kill gnir;
     2397    return(eq);
     2398  }
     2399  op=option(get);
     2400  option(returnSB);
     2401  j=quotient(j,k);
     2402  option(set,op);
     2403
     2404  list equi=equidim(j);
     2405  if(deg(equi[size(equi)][1])<=0)
     2406  {
     2407      equi[size(equi)]=k;
     2408  }
     2409  else
     2410  {
     2411    equi[size(equi)+1]=k;
     2412  }
     2413  setring P;
     2414  eq=imap(gnir,equi);
     2415  kill gnir;
     2416  return(eq);
     2417}
     2418example
     2419{ "EXAMPLE:"; echo = 2;
     2420   ring  r = 32003,(x,y,z),dp;
     2421   ideal i = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
     2422   equidim(i);
     2423}
     2424
     2425///////////////////////////////////////////////////////////////////////////////
     2426proc equidimMax(ideal i)
     2427"USAGE:  equidimMax(i); i ideal
     2428RETURN:  ideal of equidimensional locus (of maximal dimension) of i.
     2429EXAMPLE: example equidimMax; shows an example
     2430"
     2431{
     2432  if(attrib(basering,"global")!=1)
     2433  {
     2434      ERROR(
     2435      "// Not implemented for this ordering, please change to global ordering."
     2436      );
     2437  }
     2438  def  P = basering;
     2439  ideal eq;
     2440  intvec w;
     2441  int n;
     2442  int g=size(i);
     2443  int a=attrib(i,"isSB");
     2444  int homo=homog(i);
     2445
     2446  if(((homo==1)||(a==1))&&(find(ordstr(basering),"l")==0)
     2447                                &&(find(ordstr(basering),"s")==0))
     2448  {
     2449     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
     2450                              +ordstr(basering)+");");
     2451     ideal i=imap(P,i);
     2452     ideal j=i;
     2453     if(a==1)
     2454     {
     2455       attrib(j,"isSB",1);
     2456     }
     2457     else
     2458     {
     2459       j=groebner(i);
     2460     }
     2461  }
     2462  else
     2463  {
     2464     execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;");
     2465     ideal i=imap(P,i);
     2466     ideal j=groebner(i);
     2467  }
     2468  list indep;
     2469  ideal equ,equi;
     2470  if(homo==1)
     2471  {
     2472     for(n=1;n<=nvars(basering);n++)
     2473     {
     2474        w[n]=ord(var(n));
     2475     }
     2476     intvec hil=hilb(j,1,w);
     2477  }
     2478  if ((dim(j)==-1)||(size(j)==0)||(nvars(basering)==1)
     2479                  ||(dim(j)==0)||(dim(j)+g==nvars(basering)))
     2480  {
     2481    setring P;
     2482    return(i);
     2483  }
     2484
     2485  indep=maxIndependSet(j);
     2486
     2487  execute("ring gnir1 = ("+charstr(basering)+"),("+indep[1][1]+"),("
     2488                              +indep[1][2]+");");
     2489  if(homo==1)
     2490  {
     2491     ideal j=std(imap(gnir,j),hil,w);
     2492  }
     2493  else
     2494  {
     2495     ideal j=groebner(imap(gnir,j));
     2496  }
     2497  string quotring=prepareQuotientring(nvars(basering)-indep[1][3]);
     2498  execute(quotring);
     2499  ideal j=imap(gnir1,j);
     2500  kill gnir1;
     2501  j=clearSB(j);
     2502  ideal h;
     2503  for(n=1;n<=size(j);n++)
     2504  {
     2505     h[n]=leadcoef(j[n]);
     2506  }
     2507  setring gnir;
     2508  ideal h=imap(quring,h);
     2509  kill quring;
     2510
     2511  list l=minSat(j,h);
     2512
     2513  if(deg(l[2])>0)
     2514  {
     2515    equ=l[1];
     2516    attrib(equ,"isSB",1);
     2517    j=std(j,l[2]);
     2518
     2519    if(dim(equ)==dim(j))
     2520    {
     2521      equi=equidimMax(j);
     2522      equ=interred(intersect(equ,equi));
     2523    }
     2524  }
     2525  else
     2526  {
     2527    equ=i;
     2528  }
     2529
     2530  setring P;
     2531  eq=imap(gnir,equ);
     2532  kill gnir;
     2533  return(eq);
     2534}
     2535example
     2536{ "EXAMPLE:"; echo = 2;
     2537   ring  r = 32003,(x,y,z),dp;
     2538   ideal i = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
     2539   equidimMax(i);
     2540}
     2541///////////////////////////////////////////////////////////////////////////////
     2542static proc islp()
     2543{
     2544   string s=ordstr(basering);
     2545   int n=find(s,"lp");
     2546   if(!n){return(0);}
     2547   int k=find(s,",");
     2548   string t=s[k+1..size(s)];
     2549   int l=find(t,",");
     2550   t=s[1..k-1];
     2551   int m=find(t,",");
     2552   if(l+m){return(0);}
     2553   return(1);
     2554}
     2555///////////////////////////////////////////////////////////////////////////////
     2556
     2557proc algeDeco(ideal i, int w)
     2558{
     2559//reduces primery decomposition over algebraic extensions to
     2560//the other cases
     2561   def R=basering;
     2562   int n=nvars(R);
     2563
     2564//---Anfang Provisorium
     2565   if((size(i)==2) && (w==2))
     2566   {
     2567      option(redSB);
     2568      ideal J=std(i);
     2569      option(noredSB);
     2570      if((size(J)==2)&&(deg(J[1])==1))
     2571      {
     2572         ideal keep;
     2573         poly f;
     2574         int j;
     2575         for(j=1;j<=nvars(basering);j++)
     2576         {
     2577           f=J[2];
     2578           while((f/var(j))*var(j)-f==0)
     2579           {
     2580             f=f/var(j);
     2581             keep=keep,var(j);
     2582           }
     2583           J[2]=f;
     2584         }
     2585         ideal K=factorize(J[2],1);
     2586         if(deg(K[1])==0){K=0;}
     2587         K=K+std(keep);
     2588         ideal L;
     2589         list resu;
     2590         for(j=1;j<=size(K);j++)
     2591         {
     2592            L=J[1],K[j];
     2593            resu[j]=L;
     2594         }
     2595         return(resu);
     2596      }
     2597   }
     2598//---Ende Provisorium
     2599   string mp="poly p="+string(minpoly)+";";
     2600   string gnir="ring RH="+string(char(R))+",("+varstr(R)+","+string(par(1))
     2601                +"),dp;";
     2602   execute(gnir);
     2603   execute(mp);
     2604   ideal i=imap(R,i);
     2605   ideal I=subst(i,var(nvars(basering)),0);
     2606   int j;
     2607   for(j=1;j<=ncols(i);j++)
     2608   {
     2609     if(i[j]!=I[j]){break;}
     2610   }
     2611   if((j>ncols(i))&&(deg(p)==1))
     2612   {
     2613     setring R;
     2614     kill RH;
     2615     kill gnir;
     2616     string gnir="ring RH="+string(char(R))+",("+varstr(R)+"),dp;";
     2617     execute(gnir);
     2618     ideal i=imap(R,i);
     2619     ideal J;
     2620   }
     2621   else
     2622   {
     2623      i=i,p;
     2624   }
     2625   list pr;
     2626
     2627   if(w==0)
     2628   {
     2629      pr=decomp(i);
     2630   }
     2631   if(w==1)
     2632   {
     2633      pr=prim_dec(i,1);
     2634      pr=reconvList(pr);
     2635   }
     2636   if(w==2)
     2637   {
     2638      pr=minAssPrimes(i);
     2639   }
     2640   if(n<nvars(basering))
     2641   {
     2642      gnir="ring RS="+string(char(R))+",("+varstr(RH)
     2643                +"),(dp("+string(n)+"),lp);";
     2644      execute(gnir);
     2645      list pr=imap(RH,pr);
     2646      ideal K;
     2647      for(j=1;j<=size(pr);j++)
     2648      {
     2649         K=groebner(pr[j]);
     2650         K=K[2..size(K)];
     2651         pr[j]=K;
     2652      }
     2653      setring R;
     2654      list pr=imap(RS,pr);
     2655   }
     2656   else
     2657   {
     2658      setring R;
     2659      list pr=imap(RH,pr);
     2660   }
     2661   list re;
     2662   if(w==2)
     2663   {
     2664      re=pr;
     2665   }
     2666   else
     2667   {
     2668      re=convList(pr);
     2669   }
     2670   return(re);
     2671}
     2672
     2673///////////////////////////////////////////////////////////////////////////////
     2674static proc decomp(ideal i,list #)
     2675"USAGE:  decomp(i); i ideal  (for primary decomposition)   (resp.
     2676         decomp(i,1);        (for the associated primes of dimension of i) )
     2677         decomp(i,2);        (for the minimal associated primes) )
     2678         decomp(i,3);        (for the absolute primary decomposition) )
     2679RETURN:  list = list of primary ideals and their associated primes
     2680         (at even positions in the list)
     2681         (resp. a list of the minimal associated primes)
     2682NOTE:    Algorithm of Gianni/Trager/Zacharias
     2683EXAMPLE: example decomp; shows an example
     2684"
     2685{
     2686  intvec op,@vv;
     2687  def  @P = basering;
     2688  list primary,indep,ltras;
     2689  intvec @vh,isat,@w;
     2690  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn;
     2691  ideal peek=i;
     2692  ideal ser,tras;
     2693  int isS=(attrib(i,"isSB")==1);
     2694
     2695
     2696  if(size(#)>0)
     2697  {
     2698     if((#[1]==1)||(#[1]==2)||(#[1]==3))
     2699     {
     2700        @wr=#[1];
     2701        if(@wr==3){abspri=1;@wr=0;}
     2702        if(size(#)>1)
     2703        {
     2704          seri=1;
     2705          peek=#[2];
     2706          ser=#[3];
     2707        }
     2708      }
     2709      else
     2710      {
     2711        seri=1;
     2712        peek=#[1];
     2713        ser=#[2];
     2714      }
     2715  }
     2716  if(abspri)
     2717  {
     2718     list absprimary,abskeep,absprimarytmp,abskeeptmp;
     2719  }
     2720  homo=homog(i);
     2721  if(homo==1)
     2722  {
     2723    if(attrib(i,"isSB")!=1)
     2724    {
     2725      //ltras=mstd(i);
     2726      tras=groebner(i);
     2727      ltras=tras,tras;
     2728      attrib(ltras[1],"isSB",1);
     2729    }
     2730    else
     2731    {
     2732      ltras=i,i;
     2733      attrib(ltras[1],"isSB",1);
     2734    }
     2735    tras=ltras[1];
     2736    attrib(tras,"isSB",1);
     2737    if(dim(tras)==0)
     2738    {
     2739        primary[1]=ltras[2];
     2740        primary[2]=maxideal(1);
     2741        if(@wr>0)
     2742        {
     2743           list l;
     2744           l[1]=maxideal(1);
     2745           l[2]=maxideal(1);
     2746           return(l);
     2747        }
     2748        return(primary);
     2749     }
     2750     for(@n=1;@n<=nvars(basering);@n++)
     2751     {
     2752        @w[@n]=ord(var(@n));
     2753     }
     2754     intvec @hilb=hilb(tras,1,@w);
     2755     intvec keephilb=@hilb;
     2756  }
     2757
     2758  //----------------------------------------------------------------
     2759  //i is the zero-ideal
     2760  //----------------------------------------------------------------
     2761
     2762  if(size(i)==0)
     2763  {
     2764     primary=i,i;
     2765     return(primary);
     2766  }
     2767
     2768  //----------------------------------------------------------------
     2769  //pass to the lexicographical ordering and compute a standardbasis
     2770  //----------------------------------------------------------------
     2771
     2772  int lp=islp();
     2773
     2774  execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);");
     2775  op=option(get);
     2776  option(redSB);
     2777
     2778  ideal ser=fetch(@P,ser);
     2779
     2780  if(homo==1)
     2781  {
     2782     if(!lp)
     2783     {
     2784       ideal @j=std(fetch(@P,i),@hilb,@w);
     2785     }
     2786     else
     2787     {
     2788        ideal @j=fetch(@P,tras);
     2789        attrib(@j,"isSB",1);
     2790     }
     2791  }
     2792  else
     2793  {
     2794      if(lp&&isS)
     2795      {
     2796        ideal @j=fetch(@P,i);
     2797        attrib(@j,"isSB",1);
     2798      }
     2799      else
     2800      {
     2801        ideal @j=groebner(fetch(@P,i));
     2802      }
     2803  }
     2804  option(set,op);
     2805  if(seri==1)
     2806  {
     2807    ideal peek=fetch(@P,peek);
     2808    attrib(peek,"isSB",1);
     2809  }
     2810  else
     2811  {
     2812    ideal peek=@j;
     2813  }
     2814  if((size(ser)==0)&&(!abspri))
     2815  {
     2816    ideal fried;
     2817    @n=size(@j);
     2818    for(@k=1;@k<=@n;@k++)
     2819    {
     2820      if(deg(lead(@j[@k]))==1)
     2821      {
     2822        fried[size(fried)+1]=@j[@k];
     2823        @j[@k]=0;
     2824      }
     2825    }
     2826    if(size(fried)==nvars(basering))
     2827    {
     2828       setring @P;
     2829       primary[1]=i;
     2830       primary[2]=i;
     2831       return(primary);
     2832    }
     2833    if(size(fried)>0)
     2834    {
     2835       string newva;
     2836       string newma;
     2837       for(@k=1;@k<=nvars(basering);@k++)
     2838       {
     2839          @n1=0;
     2840          for(@n=1;@n<=size(fried);@n++)
     2841          {
     2842             if(leadmonom(fried[@n])==var(@k))
     2843             {
     2844                @n1=1;
     2845                break;
     2846             }
     2847          }
     2848          if(@n1==0)
     2849          {
     2850            newva=newva+string(var(@k))+",";
     2851            newma=newma+string(var(@k))+",";
     2852          }
     2853          else
     2854          {
     2855            newma=newma+string(0)+",";
     2856          }
     2857       }
     2858       newva[size(newva)]=")";
     2859       newma[size(newma)]=";";
     2860       execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;");
     2861       execute("map @kappa=gnir,"+newma);
     2862       ideal @j= @kappa(@j);
     2863       @j=simplify(@j,2);
     2864       attrib(@j,"isSB",1);
     2865       list pr=decomp(@j);
     2866       setring gnir;
     2867       list pr=imap(@deirf,pr);
     2868       for(@k=1;@k<=size(pr);@k++)
     2869       {
     2870         @j=pr[@k]+fried;
     2871         pr[@k]=@j;
     2872       }
     2873       setring @P;
     2874       return(imap(gnir,pr));
     2875    }
     2876  }
     2877  //----------------------------------------------------------------
     2878  //j is the ring
     2879  //----------------------------------------------------------------
     2880
     2881  if (dim(@j)==-1)
     2882  {
     2883    setring @P;
     2884    primary=ideal(1),ideal(1);
     2885    return(primary);
     2886  }
     2887
     2888  //----------------------------------------------------------------
     2889  //  the case of one variable
     2890  //----------------------------------------------------------------
     2891
     2892  if(nvars(basering)==1)
     2893  {
     2894
     2895     list fac=factor(@j[1]);
     2896     list gprimary;
     2897     for(@k=1;@k<=size(fac[1]);@k++)
     2898     {
     2899        if(@wr==0)
     2900        {
     2901           gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
     2902           gprimary[2*@k]=ideal(fac[1][@k]);
     2903        }
     2904        else
     2905        {
     2906           gprimary[2*@k-1]=ideal(fac[1][@k]);
     2907           gprimary[2*@k]=ideal(fac[1][@k]);
     2908        }
     2909     }
     2910     setring @P;
     2911     primary=fetch(gnir,gprimary);
     2912
     2913//HIER
     2914    if(abspri)
     2915    {
     2916       list resu,tempo;
     2917       string absotto;
     2918       for(ab=1;ab<=size(primary)/2;ab++)
     2919       {
     2920          absotto= absFactorize(primary[2*ab][1],77);
     2921          tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
     2922          resu[ab]=tempo;
     2923       }
     2924       primary=resu;
     2925     }
     2926     return(primary);
     2927  }
     2928
     2929 //------------------------------------------------------------------
     2930 //the zero-dimensional case
     2931 //------------------------------------------------------------------
     2932  if (dim(@j)==0)
     2933  {
     2934    op=option(get);
     2935    option(redSB);
     2936    list gprimary= zero_decomp(@j,ser,@wr);
     2937
     2938    setring @P;
     2939    primary=fetch(gnir,gprimary);
     2940
     2941    if(size(ser)>0)
     2942    {
     2943      primary=cleanPrimary(primary);
     2944    }
     2945//HIER
     2946    if(abspri)
     2947    {
     2948       list resu,tempo;
     2949       string absotto;
     2950       for(ab=1;ab<=size(primary)/2;ab++)
     2951       {
     2952          absotto= absFactorize(primary[2*ab][1],77);
     2953          tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
     2954          resu[ab]=tempo;
     2955       }
     2956       primary=resu;
     2957    }
     2958    return(primary);
     2959  }
     2960
     2961  poly @gs,@gh,@p;
     2962  string @va,quotring;
     2963  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
     2964  ideal @h;
     2965  int jdim=dim(@j);
     2966  list fett;
     2967  int lauf,di,newtest;
     2968  //------------------------------------------------------------------
     2969  //search for a maximal independent set indep,i.e.
     2970  //look for subring such that the intersection with the ideal is zero
     2971  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
     2972  //indep[1] is the new varstring and indep[2] the string for block-ordering
     2973  //------------------------------------------------------------------
     2974  if(@wr!=1)
     2975  {
     2976     allindep=independSet(@j);
     2977     for(@m=1;@m<=size(allindep);@m++)
     2978     {
     2979        if(allindep[@m][3]==jdim)
     2980        {
     2981           di++;
     2982           indep[di]=allindep[@m];
     2983        }
     2984        else
     2985        {
     2986           lauf++;
     2987           restindep[lauf]=allindep[@m];
     2988        }
     2989     }
     2990   }
     2991   else
     2992   {
     2993      indep=maxIndependSet(@j);
     2994   }
     2995
     2996  ideal jkeep=@j;
     2997  if(ordstr(@P)[1]=="w")
     2998  {
     2999     execute("ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");");
     3000  }
     3001  else
     3002  {
     3003     execute( "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);");
     3004  }
     3005
     3006  if(homo==1)
     3007  {
     3008    if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
     3009       ||(ordstr(@P)[3]=="w"))
     3010    {
     3011      ideal jwork=imap(@P,tras);
     3012      attrib(jwork,"isSB",1);
     3013    }
     3014    else
     3015    {
     3016      ideal jwork=std(imap(gnir,@j),@hilb,@w);
     3017    }
     3018
     3019  }
     3020  else
     3021  {
     3022    ideal jwork=groebner(imap(gnir,@j));
     3023  }
     3024  list hquprimary;
     3025  poly @p,@q;
     3026  ideal @h,fac,ser;
     3027  ideal @Ptest=1;
     3028  di=dim(jwork);
     3029  keepdi=di;
     3030
     3031  setring gnir;
     3032  for(@m=1;@m<=size(indep);@m++)
     3033  {
     3034     isat=0;
     3035     @n2=0;
     3036     if((indep[@m][1]==varstr(basering))&&(@m==1))
     3037     //this is the good case, nothing to do, just to have the same notations
     3038     //change the ring
     3039     {
     3040        execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
     3041                              +ordstr(basering)+");");
     3042        ideal @j=fetch(gnir,@j);
     3043        attrib(@j,"isSB",1);
     3044        ideal ser=fetch(gnir,ser);
     3045
     3046     }
     3047     else
     3048     {
     3049        @va=string(maxideal(1));
     3050        if(@m==1)
     3051        {
     3052           @j=fetch(@P,i);
     3053        }
     3054        execute("ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
     3055                              +indep[@m][2]+");");
     3056        execute("map phi=gnir,"+@va+";");
     3057        op=option(get);
     3058        option(redSB);
     3059        if(homo==1)
     3060        {
     3061           ideal @j=std(phi(@j),@hilb,@w);
     3062        }
     3063        else
     3064        {
     3065          ideal @j=groebner(phi(@j));
     3066        }
     3067        ideal ser=phi(ser);
     3068
     3069        option(set,op);
     3070     }
     3071     if((deg(@j[1])==0)||(dim(@j)<jdim))
     3072     {
     3073       setring gnir;
     3074       break;
     3075     }
     3076     for (lauf=1;lauf<=size(@j);lauf++)
     3077     {
     3078         fett[lauf]=size(@j[lauf]);
     3079     }
     3080     //------------------------------------------------------------------------
     3081     //we have now the following situation:
     3082     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     3083     //to this quotientring, j is their still a standardbasis, the
     3084     //leading coefficients of the polynomials  there (polynomials in
     3085     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     3086     //we need their ggt, gh, because of the following: let
     3087     //(j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
     3088     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     3089     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
     3090
     3091     //------------------------------------------------------------------------
     3092
     3093     //arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and
     3094     //map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..rest..]
     3095     //------------------------------------------------------------------------
     3096
     3097     quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
     3098
     3099     //---------------------------------------------------------------------
     3100     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     3101     //---------------------------------------------------------------------
     3102
     3103     ideal @jj=lead(@j);               //!! vorn vereinbaren
     3104     execute(quotring);
     3105
     3106     ideal @jj=imap(gnir1,@jj);
     3107     @vv=clearSBNeu(@jj,fett);  //!! vorn vereinbaren
     3108     setring gnir1;
     3109     @k=size(@j);
     3110     for (lauf=1;lauf<=@k;lauf++)
     3111     {
     3112         if(@vv[lauf]==1)
     3113         {
     3114            @j[lauf]=0;
     3115         }
     3116     }
     3117     @j=simplify(@j,2);
     3118     setring quring;
     3119      // @j considered in the quotientring
     3120     ideal @j=imap(gnir1,@j);
     3121
     3122     ideal ser=imap(gnir1,ser);
     3123
     3124     kill gnir1;
     3125
     3126     //j is a standardbasis in the quotientring but usually not minimal
     3127     //here it becomes minimal
     3128
     3129     attrib(@j,"isSB",1);
     3130
     3131     //we need later ggt(h[1],...)=gh for saturation
     3132     ideal @h;
     3133     if(deg(@j[1])>0)
     3134     {
     3135        for(@n=1;@n<=size(@j);@n++)
     3136        {
     3137           @h[@n]=leadcoef(@j[@n]);
     3138        }
     3139        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
     3140        op=option(get);
     3141        option(redSB);
     3142
     3143        list uprimary= zero_decomp(@j,ser,@wr);
     3144//HIER
     3145        if(abspri)
     3146        {
     3147           ideal II;
     3148           ideal jmap;
     3149           map sigma;
     3150           nn=nvars(basering);
     3151           map invsigma=basering,maxideal(1);
     3152           for(ab=1;ab<=size(uprimary)/2;ab++)
     3153           {
     3154              II=uprimary[2*ab];
     3155              attrib(II,"isSB",1);
     3156              if(deg(II[1])!=vdim(II))
     3157              {
     3158                 jmap=randomLast(50);
     3159                 sigma=basering,jmap;
     3160                 jmap[nn]=2*var(nn)-jmap[nn];
     3161                 invsigma=basering,jmap;
     3162                 II=groebner(sigma(II));
     3163               }
     3164               absprimarytmp[ab]= absFactorize(II[1],77);
     3165               II=var(nn);
     3166               abskeeptmp[ab]=string(invsigma(II));
     3167               invsigma=basering,maxideal(1);
     3168           }
     3169        }
     3170        option(set,op);
     3171     }
     3172     else
     3173     {
     3174       list uprimary;
     3175       uprimary[1]=ideal(1);
     3176       uprimary[2]=ideal(1);
     3177     }
     3178     //we need the intersection of the ideals in the list quprimary with the
     3179     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     3180     //but fi polynomials, then the intersection of q with the polynomialring
     3181     //is the saturation of the ideal generated by f1,...,fr with respect to
     3182     //h which is the lcm of the leading coefficients of the fi considered in
     3183     //in the quotientring: this is coded in saturn
     3184
     3185     list saturn;
     3186     ideal hpl;
     3187
     3188     for(@n=1;@n<=size(uprimary);@n++)
     3189     {
     3190        uprimary[@n]=interred(uprimary[@n]); // temporary fix
     3191        hpl=0;
     3192        for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
     3193        {
     3194           hpl=hpl,leadcoef(uprimary[@n][@n1]);
     3195        }
     3196        saturn[@n]=hpl;
     3197     }
     3198
     3199     //--------------------------------------------------------------------
     3200     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     3201     //back to the polynomialring
     3202     //---------------------------------------------------------------------
     3203     setring gnir;
     3204
     3205     collectprimary=imap(quring,uprimary);
     3206     lsau=imap(quring,saturn);
     3207     @h=imap(quring,@h);
     3208
     3209     kill quring;
     3210
     3211     @n2=size(quprimary);
     3212     @n3=@n2;
     3213
     3214     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
     3215     {
     3216        if(deg(collectprimary[2*@n1][1])>0)
     3217        {
     3218           @n2++;
     3219           quprimary[@n2]=collectprimary[2*@n1-1];
     3220           lnew[@n2]=lsau[2*@n1-1];
     3221           @n2++;
     3222           lnew[@n2]=lsau[2*@n1];
     3223           quprimary[@n2]=collectprimary[2*@n1];
     3224           if(abspri)
     3225           {
     3226              absprimary[@n2/2]=absprimarytmp[@n1];
     3227              abskeep[@n2/2]=abskeeptmp[@n1];
     3228           }
     3229        }
     3230     }
     3231     //here the intersection with the polynomialring
     3232     //mentioned above is really computed
     3233     for(@n=@n3/2+1;@n<=@n2/2;@n++)
     3234     {
     3235        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
     3236        {
     3237           quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     3238           quprimary[2*@n]=quprimary[2*@n-1];
     3239        }
     3240        else
     3241        {
     3242           if(@wr==0)
     3243           {
     3244              quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     3245           }
     3246           quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
     3247        }
     3248     }
     3249
     3250     if(size(@h)>0)
     3251     {
     3252           //---------------------------------------------------------------
     3253           //we change to @Phelp to have the ordering dp for saturation
     3254           //---------------------------------------------------------------
     3255        setring @Phelp;
     3256        @h=imap(gnir,@h);
     3257        if(@wr!=1)
     3258        {
     3259          if(defined(@LL)){kill @LL;}
     3260          list @LL=minSat(jwork,@h);
     3261          @Ptest=intersect(@Ptest,@LL[1]);
     3262          @q=@LL[2];
     3263        }
     3264        else
     3265        {
     3266            fac=ideal(0);
     3267            for(lauf=1;lauf<=ncols(@h);lauf++)
     3268           {
     3269              if(deg(@h[lauf])>0)
     3270              {
     3271                 fac=fac+factorize(@h[lauf],1);
     3272              }
     3273           }
     3274           fac=simplify(fac,6);
     3275           @q=1;
     3276           for(lauf=1;lauf<=size(fac);lauf++)
     3277           {
     3278             @q=@q*fac[lauf];
     3279           }
     3280        }
     3281        jwork=std(jwork,@q);
     3282        keepdi=dim(jwork);
     3283        if(keepdi<di)
     3284        {
     3285           setring gnir;
     3286           @j=imap(@Phelp,jwork);
     3287           break;
     3288        }
     3289        if(homo==1)
     3290        {
     3291              @hilb=hilb(jwork,1,@w);
     3292        }
     3293
     3294        setring gnir;
     3295        @j=imap(@Phelp,jwork);
     3296     }
     3297  }
     3298
     3299  if((size(quprimary)==0)&&(@wr==1))
     3300  {
     3301     @j=ideal(1);
     3302     quprimary[1]=ideal(1);
     3303     quprimary[2]=ideal(1);
     3304  }
     3305  if((size(quprimary)==0))
     3306  {
     3307    keepdi=di-1;
     3308    quprimary[1]=ideal(1);
     3309    quprimary[2]=ideal(1);
     3310  }
     3311  //---------------------------------------------------------------
     3312  //notice that j=sat(j,gh) intersected with (j,gh^n)
     3313  //we finished with sat(j,gh) and have to start with (j,gh^n)
     3314  //---------------------------------------------------------------
     3315  if((deg(@j[1])!=0)&&(@wr!=1))
     3316  {
     3317     if(size(quprimary)>0)
     3318     {
     3319        setring @Phelp;
     3320        ser=imap(gnir,ser);
     3321        hquprimary=imap(gnir,quprimary);
     3322        if(@wr==0)
     3323        {
     3324           //HIER STATT DURCHSCHNITT SATURIEREN!
     3325           ideal htest=@Ptest;
     3326        }
     3327        else
     3328        {
     3329           ideal htest=hquprimary[2];
     3330
     3331           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
     3332           {
     3333              htest=intersect(htest,hquprimary[2*@n1]);
     3334           }
     3335        }
     3336
     3337        if(size(ser)>0)
     3338        {
     3339           ser=intersect(htest,ser);
     3340        }
     3341        else
     3342        {
     3343          ser=htest;
     3344        }
     3345        setring gnir;
     3346        ser=imap(@Phelp,ser);
     3347     }
     3348     if(size(reduce(ser,peek,1))!=0)
     3349     {
     3350        for(@m=1;@m<=size(restindep);@m++)
     3351        {
     3352         // if(restindep[@m][3]>=keepdi)
     3353         // {
     3354           isat=0;
     3355           @n2=0;
     3356
     3357           if(restindep[@m][1]==varstr(basering))
     3358           //the good case, nothing to do, just to have the same notations
     3359           //change the ring
     3360           {
     3361              execute("ring gnir1 = ("+charstr(basering)+"),("+
     3362                varstr(basering)+"),("+ordstr(basering)+");");
     3363              ideal @j=fetch(gnir,jkeep);
     3364              attrib(@j,"isSB",1);
     3365           }
     3366           else
     3367           {
     3368              @va=string(maxideal(1));
     3369              execute("ring gnir1 = ("+charstr(basering)+"),("+
     3370                      restindep[@m][1]+"),(" +restindep[@m][2]+");");
     3371              execute("map phi=gnir,"+@va+";");
     3372              op=option(get);
     3373              option(redSB);
     3374              if(homo==1)
     3375              {
     3376                 ideal @j=std(phi(jkeep),keephilb,@w);
     3377              }
     3378              else
     3379              {
     3380                ideal @j=groebner(phi(jkeep));
     3381              }
     3382              ideal ser=phi(ser);
     3383              option(set,op);
     3384           }
     3385
     3386           for (lauf=1;lauf<=size(@j);lauf++)
     3387           {
     3388              fett[lauf]=size(@j[lauf]);
     3389           }
     3390           //------------------------------------------------------------------
     3391           //we have now the following situation:
     3392           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may
     3393           //pass to this quotientring, j is their still a standardbasis, the
     3394           //leading coefficients of the polynomials  there (polynomials in
     3395           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     3396           //we need their ggt, gh, because of the following:
     3397           //let (j:gh^n)=(j:gh^infinity) then
     3398           //j*K(var(nnp+1),..,var(nva))[..the rest..]
     3399           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     3400           //on the other hand j=(j,gh^n) intersected with (j:gh^n)
     3401
     3402           //------------------------------------------------------------------
     3403
     3404           //the arrangement for the quotientring
     3405           // K(var(nnp+1),..,var(nva))[..the rest..]
     3406           //and the map phi:K[var(1),...,var(nva)] ---->
     3407           //--->K(var(nnpr+1),..,var(nva))[..the rest..]
     3408           //------------------------------------------------------------------
     3409
     3410           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
     3411
     3412           //------------------------------------------------------------------
     3413           //we pass to the quotientring  K(var(nnp+1),..,var(nva))[..rest..]
     3414           //------------------------------------------------------------------
     3415
     3416           execute(quotring);
     3417
     3418           // @j considered in the quotientring
     3419           ideal @j=imap(gnir1,@j);
     3420           ideal ser=imap(gnir1,ser);
     3421
     3422           kill gnir1;
     3423
     3424           //j is a standardbasis in the quotientring but usually not minimal
     3425           //here it becomes minimal
     3426           @j=clearSB(@j,fett);
     3427           attrib(@j,"isSB",1);
     3428
     3429           //we need later ggt(h[1],...)=gh for saturation
     3430           ideal @h;
     3431
     3432           for(@n=1;@n<=size(@j);@n++)
     3433           {
     3434              @h[@n]=leadcoef(@j[@n]);
     3435           }
     3436           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
     3437
     3438           op=option(get);
     3439           option(redSB);
     3440           list uprimary= zero_decomp(@j,ser,@wr);
     3441//HIER
     3442           if(abspri)
     3443           {
     3444              ideal II;
     3445              ideal jmap;
     3446              map sigma;
     3447              nn=nvars(basering);
     3448              map invsigma=basering,maxideal(1);
     3449              for(ab=1;ab<=size(uprimary)/2;ab++)
     3450              {
     3451                 II=uprimary[2*ab];
     3452                 attrib(II,"isSB",1);
     3453                 if(deg(II[1])!=vdim(II))
     3454                 {
     3455                    jmap=randomLast(50);
     3456                    sigma=basering,jmap;
     3457                    jmap[nn]=2*var(nn)-jmap[nn];
     3458                    invsigma=basering,jmap;
     3459                    II=groebner(sigma(II));
     3460                  }
     3461                  absprimarytmp[ab]= absFactorize(II[1],77);
     3462                  II=var(nn);
     3463                  abskeeptmp[ab]=string(invsigma(II));
     3464                  invsigma=basering,maxideal(1);
     3465              }
     3466           }
     3467           option(set,op);
     3468
     3469           //we need the intersection of the ideals in the list quprimary with
     3470           //the polynomialring, i.e. let q=(f1,...,fr) in the quotientring
     3471           //such an ideal but fi polynomials, then the intersection of q with
     3472           //the polynomialring is the saturation of the ideal generated by
     3473           //f1,...,fr with respect toh which is the lcm of the leading
     3474           //coefficients of the fi considered in the quotientring:
     3475           //this is coded in saturn
     3476
     3477           list saturn;
     3478           ideal hpl;
     3479
     3480           for(@n=1;@n<=size(uprimary);@n++)
     3481           {
     3482              hpl=0;
     3483              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
     3484              {
     3485                 hpl=hpl,leadcoef(uprimary[@n][@n1]);
     3486              }
     3487               saturn[@n]=hpl;
     3488           }
     3489           //------------------------------------------------------------------
     3490           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..rest..]
     3491           //back to the polynomialring
     3492           //------------------------------------------------------------------
     3493           setring gnir;
     3494           collectprimary=imap(quring,uprimary);
     3495           lsau=imap(quring,saturn);
     3496           @h=imap(quring,@h);
     3497
     3498           kill quring;
     3499
     3500
     3501           @n2=size(quprimary);
     3502           @n3=@n2;
     3503
     3504           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
     3505           {
     3506              if(deg(collectprimary[2*@n1][1])>0)
     3507              {
     3508                 @n2++;
     3509                 quprimary[@n2]=collectprimary[2*@n1-1];
     3510                 lnew[@n2]=lsau[2*@n1-1];
     3511                 @n2++;
     3512                 lnew[@n2]=lsau[2*@n1];
     3513                 quprimary[@n2]=collectprimary[2*@n1];
     3514                 if(abspri)
     3515                 {
     3516                   absprimary[@n2/2]=absprimarytmp[@n1];
     3517                   abskeep[@n2/2]=abskeeptmp[@n1];
     3518                 }
     3519              }
     3520           }
     3521
     3522
     3523           //here the intersection with the polynomialring
     3524           //mentioned above is really computed
     3525
     3526           for(@n=@n3/2+1;@n<=@n2/2;@n++)
     3527           {
     3528              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
     3529              {
     3530                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     3531                 quprimary[2*@n]=quprimary[2*@n-1];
     3532              }
     3533              else
     3534              {
     3535                 if(@wr==0)
     3536                 {
     3537                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     3538                 }
     3539                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
     3540              }
     3541           }
     3542           if(@n2>=@n3+2)
     3543           {
     3544              setring @Phelp;
     3545              ser=imap(gnir,ser);
     3546              hquprimary=imap(gnir,quprimary);
     3547              for(@n=@n3/2+1;@n<=@n2/2;@n++)
     3548              {
     3549                if(@wr==0)
     3550                {
     3551                   ser=intersect(ser,hquprimary[2*@n-1]);
     3552                }
     3553                else
     3554                {
     3555                   ser=intersect(ser,hquprimary[2*@n]);
     3556                }
     3557              }
     3558              setring gnir;
     3559              ser=imap(@Phelp,ser);
     3560           }
     3561
     3562         // }
     3563        }
     3564//HIER
     3565        if(abspri)
     3566        {
     3567          list resu,tempo;
     3568          for(ab=1;ab<=size(quprimary)/2;ab++)
     3569          {
     3570             if (deg(quprimary[2*ab][1])!=0)
     3571             {
     3572               tempo=quprimary[2*ab-1],quprimary[2*ab],
     3573                         absprimary[ab],abskeep[ab];
     3574               resu[ab]=tempo;
     3575             }
     3576          }
     3577          quprimary=resu;
     3578          @wr=3;
     3579        }
     3580        if(size(reduce(ser,peek,1))!=0)
     3581        {
     3582           if(@wr>0)
     3583           {
     3584              htprimary=decomp(@j,@wr,peek,ser);
     3585           }
     3586           else
     3587           {
     3588              htprimary=decomp(@j,peek,ser);
     3589           }
     3590           // here we collect now both results primary(sat(j,gh))
     3591           // and primary(j,gh^n)
     3592           @n=size(quprimary);
     3593           for (@k=1;@k<=size(htprimary);@k++)
     3594           {
     3595              quprimary[@n+@k]=htprimary[@k];
     3596           }
     3597        }
     3598     }
     3599
     3600   }
     3601   else
     3602   {
     3603      if(abspri)
     3604      {
     3605        list resu,tempo;
     3606        for(ab=1;ab<=size(quprimary)/2;ab++)
     3607        {
     3608           tempo=quprimary[2*ab-1],quprimary[2*ab],
     3609                   absprimary[ab],abskeep[ab];
     3610           resu[ab]=tempo;
     3611        }
     3612        quprimary=resu;
     3613      }
     3614   }
     3615  //---------------------------------------------------------------------------
     3616  //back to the ring we started with
     3617  //the final result: primary
     3618  //---------------------------------------------------------------------------
     3619
     3620  setring @P;
     3621  primary=imap(gnir,quprimary);
     3622  if(!abspri)
     3623  {
     3624     primary=cleanPrimary(primary);
     3625  }
     3626  return(primary);
     3627}
     3628
     3629
     3630example
     3631{ "EXAMPLE:"; echo = 2;
     3632   ring  r = 32003,(x,y,z),lp;
     3633   poly  p = z2+1;
     3634   poly  q = z4+2;
     3635   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     3636   list pr= decomp(i);
     3637   pr;
     3638   testPrimary( pr, i);
     3639}
     3640
     3641///////////////////////////////////////////////////////////////////////////////
     3642static proc powerCoeffs(poly f,int e)
     3643//computes a polynomial with the same monomials as f but coefficients
     3644//the p^e th power of the coefficients of f
     3645{
     3646   int i;
     3647   poly g;
     3648   int ex=char(basering)^e;
     3649   for(i=1;i<=size(f);i++)
     3650   {
     3651      g=g+leadcoef(f[i])^ex*leadmonom(f[i]);
     3652   }
     3653   return(g);
     3654}
     3655///////////////////////////////////////////////////////////////////////////////
     3656
     3657proc sep(poly f,int i, list #)
     3658"USAGE:  input: a polynomial f depending on the i-th variable and optional
     3659         an integer k considering the polynomial f defined over Fp(t1,...,tm)
     3660         as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
     3661 RETURN: the separabel part of f as polynomial in Fp(t1,...,tm)
     3662        and an integer k to indicate that f should be considerd
     3663        as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
     3664 EXAMPLE: example sep; shows an example
     3665{
     3666   def R=basering;
     3667   int k;
     3668   if(size(#)>0){k=#[1];}
     3669
     3670
     3671   poly h=gcd(f,diff(f,var(i)));
     3672   if((reduce(f,std(h))!=0)||(reduce(diff(f,var(i)),std(h))!=0))
     3673   {
     3674      ERROR("FEHLER IN GCD");
     3675   }
     3676   poly g1=lift(h,f)[1][1];    //  f/h
     3677   poly h1;
     3678
     3679   while(h!=h1)
     3680   {
     3681      h1=h;
     3682      h=gcd(h,diff(h,var(i)));
     3683   }
     3684
     3685   if(deg(h1)==0){return(list(g1,k));} //in characteristic 0 we return here
     3686
     3687   k++;
     3688
     3689   ideal ma=maxideal(1);
     3690   ma[i]=var(i)^char(R);
     3691   map phi=R,ma;
     3692   ideal hh=h;    //this is technical because preimage works only for ideals
     3693
     3694   poly u=preimage(R,phi,hh)[1]; //h=u(x(i)^p)
     3695
     3696   list g2=sep(u,i,k);           //we consider u(t(1)^(p^-1),...,t(m)^(p^-1))
     3697   g1=powerCoeffs(g1,g2[2]-k+1); //to have g1 over the same field as g2[1]
     3698
     3699   list g3=sep(g1*g2[1],i,g2[2]);
     3700   return(g3);
     3701}
     3702example
     3703{ "EXAMPLE:"; echo = 2;
     3704   ring R=(5,t,s),(x,y,z),dp;
     3705   poly f=(x^25-t*x^5+t)*(x^3+s);
     3706   sep(f,1);
     3707}
     3708
     3709///////////////////////////////////////////////////////////////////////////////
     3710 proc zeroRad(ideal I,list #)
     3711"USAGE:  zeroRad(I) , I a zero-dimensional ideal
     3712 RETURN: the radical of I
     3713 NOTE:  Algorithm of Kemper
     3714 EXAMPLE: example zeroRad; shows an example
     3715{
     3716   if(homog(I)==1){return(maxideal(1));}
     3717   //I needs to be a reduced standard basis
     3718   def R=basering;
     3719   int m=npars(R);
     3720   int n=nvars(R);
     3721   int p=char(R);
     3722   int d=vdim(I);
     3723   int i,k;
     3724   list l;
     3725   if(((p==0)||(p>d))&&(d==deg(I[1])))
     3726   {
     3727     intvec e=leadexp(I[1]);
     3728     for(i=1;i<=nvars(basering);i++)
     3729     {
     3730       if(e[i]!=0) break;
     3731     }
     3732     I[1]=sep(I[1],i)[1];
     3733     return(interred(I));
     3734   }
     3735   intvec op=option(get);
     3736
     3737   option(redSB);
     3738   ideal F=finduni(I);//F[i] generates I intersected with K[var(i)]
     3739
     3740   option(set,op);
     3741   if(size(#)>0){I=#[1];}
     3742
     3743   for(i=1;i<=n;i++)
     3744   {
     3745      l[i]=sep(F[i],i);
     3746      F[i]=l[i][1];
     3747      if(l[i][2]>k){k=l[i][2];}  //computation of the maximal k
     3748   }
     3749
     3750   if((k==0)||(m==0)){return(interred(I+F));}        //the separable case
     3751
     3752   for(i=1;i<=n;i++)             //consider all polynomials over
     3753   {                             //Fp(t(1)^(p^-k),...,t(m)^(p^-k))
     3754      F[i]=powerCoeffs(F[i],k-l[i][2]);
     3755   }
     3756
     3757   string cR="ring @R="+string(p)+",("+parstr(R)+","+varstr(R)+"),dp;";
     3758   execute(cR);
     3759   ideal F=imap(R,F);
     3760
     3761   string nR="ring @S="+string(p)+",(y(1..m),"+varstr(R)+","+parstr(R)+"),dp;";
     3762   execute(nR);
     3763
     3764   ideal G=fetch(@R,F);    //G[i](t(1)^(p^-k),...,t(m)^(p^-k),x(i))=sep(F[i])
     3765
     3766   ideal I=imap(R,I);
     3767   ideal J=I+G;
     3768   poly el=1;
     3769   k=p^k;
     3770   for(i=1;i<=m;i++)
     3771   {
     3772      J=J,var(i)^k-var(m+n+i);
     3773      el=el*y(i);
     3774   }
     3775
     3776   J=eliminate(J,el);
     3777   setring R;
     3778   ideal J=imap(@S,J);
     3779   return(J);
     3780}
     3781example
     3782{ "EXAMPLE:"; echo = 2;
     3783   ring R=(5,t),(x,y),dp;
     3784   ideal I=x^5-t,y^5-t;
     3785   zeroRad(I);
     3786}
     3787
     3788///////////////////////////////////////////////////////////////////////////////
     3789
     3790proc radicalEHV(ideal i)
     3791"USAGE:   radicalEHV(i); i ideal.
     3792RETURN:  ideal, the radical of i.
     3793NOTE:    Uses the algorithm of Eisenbud/Huneke/Vasconcelos, which
     3794         reduces the computation to the complete intersection case,
     3795         by taking, in the general case, a generic linear combination
     3796         of the input.
     3797         Works only in characteristic 0 or p large.
     3798EXAMPLE: example radicalEHV; shows an example
     3799"
     3800{
     3801   if(attrib(basering,"global")!=1)
     3802   {
     3803      ERROR(
     3804      "// Not implemented for this ordering, please change to global ordering."
     3805      );
     3806   }
     3807   if((char(basering)<100)&&(char(basering)!=0))
     3808   {
     3809      "WARNING: The characteristic is too small, the result may be wrong";
     3810   }
     3811   ideal J,I,I0,radI0,L,radI1,I2,radI2;
     3812   int l,n;
     3813   intvec op=option(get);
     3814   matrix M;
     3815
     3816   option(redSB);
     3817   list m=mstd(i);
     3818        I=m[2];
     3819   option(set,op);
     3820
     3821   int cod=nvars(basering)-dim(m[1]);
     3822   //-------------------complete intersection case:----------------------
     3823   if(cod==size(m[2]))
     3824   {
     3825     J=minor(jacob(I),cod);
     3826     return(quotient(I,J));
     3827   }
     3828   //-----first codim elements of I are a complete intersection:---------
     3829   for(l=1;l<=cod;l++)
     3830   {
     3831      I0[l]=I[l];
     3832   }
     3833   n=dim(std(I0))+cod-nvars(basering);
     3834   //-----last codim elements of I are a complete intersection:----------
     3835   if(n!=0)
     3836   {
     3837      for(l=1;l<=cod;l++)
     3838      {
     3839         I0[l]=I[size(I)-l+1];
     3840      }
     3841      n=dim(std(I0))+cod-nvars(basering);
     3842   }
     3843   //-----taking a generic linear combination of the input:--------------
     3844   if(n!=0)
     3845   {
     3846      M=transpose(sparsetriag(size(m[2]),cod,95,1));
     3847      I0=ideal(M*transpose(I));
     3848      n=dim(std(I0))+cod-nvars(basering);
     3849   }
     3850   //-----taking a more generic linear combination of the input:---------
     3851   if(n!=0)
     3852   {
     3853      M=transpose(sparsetriag(size(m[2]),cod,0,100));
     3854      I0=ideal(M*transpose(I));
     3855      n=dim(std(I0))+cod-nvars(basering);
     3856   }
     3857   if(n==0)
     3858   {
     3859      J=minor(jacob(I0),cod);
     3860      radI0=quotient(I0,J);
     3861      L=quotient(radI0,I);
     3862      radI1=quotient(radI0,L);
     3863
     3864      if(size(reduce(radI1,m[1],1))==0)
     3865      {
     3866         return(I);
     3867      }
     3868
     3869      I2=sat(I,radI1)[1];
     3870
     3871      if(deg(I2[1])<=0)
     3872      {
     3873         return(radI1);
     3874      }
     3875      return(intersect(radI1,radicalEHV(I2)));
     3876   }
     3877   //---------------------general case-------------------------------------
     3878   return(radical(I));
     3879}
     3880example
     3881{ "EXAMPLE:";  echo = 2;
     3882   ring  r = 0,(x,y,z),dp;
     3883   poly  p = z2+1;
     3884   poly  q = z3+2;
     3885   ideal i = p*q^2,y-z2;
     3886   ideal pr= radicalEHV(i);
     3887   pr;
     3888}
     3889
     3890///////////////////////////////////////////////////////////////////////////////
     3891
     3892proc Ann(module M)
     3893"USAGE:   Ann(M);  M module
     3894RETURN:  ideal, the annihilator of coker(M)
     3895NOTE:    The output is the ideal of all elements a of the basering R such that
     3896         a * R^m is contained in M  (m=number of rows of M).
     3897EXAMPLE: example Ann; shows an example
     3898"
     3899{
     3900  M=prune(M);  //to obtain a small embedding
     3901  ideal ann=quotient1(M,freemodule(nrows(M)));
     3902  return(ann);
     3903}
     3904example
     3905{ "EXAMPLE:"; echo = 2;
     3906   ring  r = 0,(x,y,z),lp;
     3907   module M = x2-y2,z3;
     3908   Ann(M);
     3909   M = [1,x2],[y,x];
     3910   Ann(M);
     3911   qring Q=std(xy-1);
     3912   module M=imap(r,M);
     3913   Ann(M);
     3914}
     3915
     3916///////////////////////////////////////////////////////////////////////////////
     3917
     3918//computes the equidimensional part of the ideal i of codimension e
     3919static proc int_ass_primary_e(ideal i, int e)
     3920{
     3921  if(homog(i)!=1)
     3922  {
     3923     i=std(i);
     3924  }
     3925  list re=sres(i,0);                   //the resolution
     3926  re=minres(re);                       //minimized resolution
     3927  ideal ann=AnnExt_R(e,re);
     3928  if(nvars(basering)-dim(std(ann))!=e)
     3929  {
     3930    return(ideal(1));
     3931  }
     3932  return(ann);
     3933}
     3934
     3935///////////////////////////////////////////////////////////////////////////////
     3936
     3937//computes the annihilator of Ext^n(R/i,R) with given resolution re
     3938//n is not necessarily the number of variables
     3939static proc AnnExt_R(int n,list re)
     3940{
     3941  if(n<nvars(basering))
     3942  {
     3943     matrix f=transpose(re[n+1]);      //Hom(_,R)
     3944     module k=nres(f,2)[2];            //the kernel
     3945     matrix g=transpose(re[n]);        //the image of Hom(_,R)
     3946
     3947     ideal ann=quotient1(g,k);           //the anihilator
     3948  }
     3949  else
     3950  {
     3951     ideal ann=Ann(transpose(re[n]));
     3952  }
     3953  return(ann);
     3954}
     3955///////////////////////////////////////////////////////////////////////////////
     3956
     3957static proc analyze(list pr)
     3958{
     3959   int ii,jj;
     3960   for(ii=1;ii<=size(pr)/2;ii++)
     3961   {
     3962      dim(std(pr[2*ii]));
     3963      idealsEqual(pr[2*ii-1],pr[2*ii]);
     3964      "===========================";
     3965   }
     3966
     3967   for(ii=size(pr)/2;ii>1;ii--)
     3968   {
     3969      for(jj=1;jj<ii;jj++)
     3970      {
     3971         if(size(reduce(pr[2*jj],std(pr[2*ii],1)))==0)
     3972         {
     3973            "eingebette Komponente";
     3974            jj;
     3975            ii;
     3976         }
     3977      }
     3978   }
     3979}
     3980
     3981///////////////////////////////////////////////////////////////////////////////
     3982//
     3983//                  Shimoyama-Yokoyama
     3984//
     3985///////////////////////////////////////////////////////////////////////////////
     3986
     3987static proc simplifyIdeal(ideal i)
     3988{
     3989  def r=basering;
     3990
     3991  int j,k;
     3992  map phi;
     3993  poly p;
     3994
     3995  ideal iwork=i;
     3996  ideal imap1=maxideal(1);
     3997  ideal imap2=maxideal(1);
     3998
     3999
     4000  for(j=1;j<=nvars(basering);j++)
     4001  {
     4002    for(k=1;k<=size(i);k++)
     4003    {
     4004      if(deg(iwork[k]/var(j))==0)
     4005      {
     4006        p=-1/leadcoef(iwork[k]/var(j))*iwork[k];
     4007        imap1[j]=p+2*var(j);
     4008        phi=r,imap1;
     4009        iwork=phi(iwork);
     4010        iwork=subst(iwork,var(j),0);
     4011        iwork[k]=var(j);
     4012        imap1=maxideal(1);
     4013        imap2[j]=-p;
     4014        break;
     4015      }
     4016    }
     4017  }
     4018  return(iwork,imap2);
     4019}
     4020
     4021
     4022///////////////////////////////////////////////////////
     4023// ini_mod
     4024// input: a polynomial p
     4025// output: the initial term of p as needed
     4026// in the context of characteristic sets
     4027//////////////////////////////////////////////////////
     4028
     4029static proc ini_mod(poly p)
     4030{
     4031  if (p==0)
     4032  {
     4033    return(0);
     4034  }
     4035  int n; matrix m;
     4036  for( n=nvars(basering); n>0; n--)
     4037  {
     4038    m=coef(p,var(n));
     4039    if(m[1,1]!=1)
     4040    {
     4041      p=m[2,1];
     4042      break;
     4043    }
     4044  }
     4045  if(deg(p)==0)
     4046  {
     4047    p=0;
     4048  }
     4049  return(p);
     4050}
     4051///////////////////////////////////////////////////////
     4052// min_ass_prim_charsets
     4053// input: generators of an ideal PS and an integer cho
     4054// If cho=0, the given ordering of the variables is used.
     4055// Otherwise, the system tries to find an "optimal ordering",
     4056// which in some cases may considerably speed up the algorithm
     4057// output: the minimal associated primes of PS
     4058// algorithm: via characteriostic sets
     4059//////////////////////////////////////////////////////
     4060
     4061
     4062static proc min_ass_prim_charsets (ideal PS, int cho)
     4063{
     4064  if((cho<0) and (cho>1))
     4065  {
     4066    ERROR("<int> must be 0 or 1");
     4067  }
     4068  option(notWarnSB);
     4069  if(cho==0)
     4070  {
     4071    return(min_ass_prim_charsets0(PS));
     4072  }
     4073  else
     4074  {
     4075     return(min_ass_prim_charsets1(PS));
     4076  }
     4077}
     4078///////////////////////////////////////////////////////
     4079// min_ass_prim_charsets0
     4080// input: generators of an ideal PS
     4081// output: the minimal associated primes of PS
     4082// algorithm: via characteristic sets
     4083// the given ordering of the variables is used
     4084//////////////////////////////////////////////////////
     4085
     4086
     4087static proc min_ass_prim_charsets0 (ideal PS)
     4088{
     4089  intvec op;
     4090  matrix m=char_series(PS);  // We compute an irreducible
     4091                             // characteristic series
     4092  int i,j,k;
     4093  list PSI;
     4094  list PHI;  // the ideals given by the characteristic series
     4095  for(i=nrows(m);i>=1; i--)
     4096  {
     4097    PHI[i]=ideal(m[i,1..ncols(m)]);
     4098  }
     4099  // We compute the radical of each ideal in PHI
     4100  ideal I,JS,II;
     4101  int sizeJS, sizeII;
     4102  for(i=size(PHI);i>=1; i--)
     4103  {
     4104    I=0;
     4105    for(j=size(PHI[i]);j>0;j--)
     4106    {
     4107      I=I+ini_mod(PHI[i][j]);
     4108    }
     4109    JS=std(PHI[i]);
     4110    sizeJS=size(JS);
     4111    for(j=size(I);j>0;j--)
     4112    {
     4113      II=0;
     4114      sizeII=0;
     4115      k=0;
     4116      while(k<=sizeII)                  // successive saturation
     4117      {
     4118        op=option(get);
     4119        option(returnSB);
     4120        II=quotient(JS,I[j]);
     4121        option(set,op);
     4122        sizeII=size(II);
     4123        if(sizeII==sizeJS)
     4124        {
     4125          for(k=1;k<=sizeII;k++)
     4126          {
     4127            if(leadexp(II[k])!=leadexp(JS[k])) break;
     4128          }
     4129        }
     4130        JS=II;
     4131        sizeJS=sizeII;
     4132      }
     4133    }
     4134    PSI=insert(PSI,JS);
     4135  }
     4136  int sizePSI=size(PSI);
     4137  // We eliminate redundant ideals
     4138  for(i=1;i<sizePSI;i++)
     4139  {
     4140    for(j=i+1;j<=sizePSI;j++)
     4141    {
     4142      if(size(PSI[i])!=0)
     4143      {
     4144        if(size(PSI[j])!=0)
     4145        {
     4146          if(size(NF(PSI[i],PSI[j],1))==0)
     4147          {
     4148            PSI[j]=ideal(0);
     4149          }
     4150          else
     4151          {
     4152            if(size(NF(PSI[j],PSI[i],1))==0)
     4153            {
     4154              PSI[i]=ideal(0);
     4155            }
     4156          }
     4157        }
     4158      }
     4159    }
     4160  }
     4161  for(i=sizePSI;i>=1;i--)
     4162  {
     4163    if(size(PSI[i])==0)
     4164    {
     4165      PSI=delete(PSI,i);
     4166    }
     4167  }
     4168  return (PSI);
     4169}
     4170
     4171///////////////////////////////////////////////////////
     4172// min_ass_prim_charsets1
     4173// input: generators of an ideal PS
     4174// output: the minimal associated primes of PS
     4175// algorithm: via characteristic sets
     4176// input: generators of an ideal PS and an integer i
     4177// The system tries to find an "optimal ordering" of
     4178// the variables
     4179//////////////////////////////////////////////////////
     4180
     4181
     4182static proc min_ass_prim_charsets1 (ideal PS)
     4183{
     4184  intvec op;
     4185  def oldring=basering;
     4186  string n=system("neworder",PS);
     4187  execute("ring r=("+charstr(oldring)+"),("+n+"),dp;");
     4188  ideal PS=imap(oldring,PS);
     4189  matrix m=char_series(PS);  // We compute an irreducible
     4190                             // characteristic series
     4191  int i,j,k;
     4192  ideal I;
     4193  list PSI;
     4194  list PHI;    // the ideals given by the characteristic series
     4195  list ITPHI;  // their initial terms
     4196  for(i=nrows(m);i>=1; i--)
     4197  {
     4198    PHI[i]=ideal(m[i,1..ncols(m)]);
     4199    I=0;
     4200    for(j=size(PHI[i]);j>0;j=j-1)
     4201    {
     4202      I=I,ini_mod(PHI[i][j]);
     4203    }
     4204    I=I[2..ncols(I)];
     4205    ITPHI[i]=I;
     4206  }
     4207  setring oldring;
     4208  matrix m=imap(r,m);
     4209  list PHI=imap(r,PHI);
     4210  list ITPHI=imap(r,ITPHI);
     4211  // We compute the radical of each ideal in PHI
     4212  ideal I,JS,II;
     4213  int sizeJS, sizeII;
     4214  for(i=size(PHI);i>=1; i--)
     4215  {
     4216    I=0;
     4217    for(j=size(PHI[i]);j>0;j--)
     4218    {
     4219      I=I+ITPHI[i][j];
     4220    }
     4221    JS=std(PHI[i]);
     4222    sizeJS=size(JS);
     4223    for(j=size(I);j>0;j--)
     4224    {
     4225      II=0;
     4226      sizeII=0;
     4227      k=0;
     4228      while(k<=sizeII)                  // successive iteration
     4229      {
     4230        op=option(get);
     4231        option(returnSB);
     4232        II=quotient(JS,I[j]);
     4233        option(set,op);
     4234//std
     4235//         II=std(II);
     4236        sizeII=size(II);
     4237        if(sizeII==sizeJS)
     4238        {
     4239          for(k=1;k<=sizeII;k++)
     4240          {
     4241            if(leadexp(II[k])!=leadexp(JS[k])) break;
     4242          }
     4243        }
     4244        JS=II;
     4245        sizeJS=sizeII;
     4246      }
     4247    }
     4248    PSI=insert(PSI,JS);
     4249  }
     4250  int sizePSI=size(PSI);
     4251  // We eliminate redundant ideals
     4252  for(i=1;i<sizePSI;i++)
     4253  {
     4254    for(j=i+1;j<=sizePSI;j++)
     4255    {
     4256      if(size(PSI[i])!=0)
     4257      {
     4258        if(size(PSI[j])!=0)
     4259        {
     4260          if(size(NF(PSI[i],PSI[j],1))==0)
     4261          {
     4262            PSI[j]=ideal(0);
     4263          }
     4264          else
     4265          {
     4266            if(size(NF(PSI[j],PSI[i],1))==0)
     4267            {
     4268              PSI[i]=ideal(0);
     4269            }
     4270          }
     4271        }
     4272      }
     4273    }
     4274  }
     4275  for(i=sizePSI;i>=1;i--)
     4276  {
     4277    if(size(PSI[i])==0)
     4278    {
     4279      PSI=delete(PSI,i);
     4280    }
     4281  }
     4282  return (PSI);
     4283}
     4284
     4285
     4286/////////////////////////////////////////////////////
     4287// proc prim_dec
     4288// input:  generators of an ideal I and an integer choose
     4289// If choose=0, min_ass_prim_charsets with the given
     4290// ordering of the variables is used.
     4291// If choose=1, min_ass_prim_charsets with the "optimized"
     4292// ordering of the variables is used.
     4293// If choose=2, minAssPrimes from primdec.lib is used
     4294// If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
     4295// output: a primary decomposition of I, i.e., a list
     4296// of pairs consisting of a standard basis of a primary component
     4297// of I and a standard basis of the corresponding associated prime.
     4298// To compute the minimal associated primes of a given ideal
     4299// min_ass_prim_l is called, i.e., the minimal associated primes
     4300// are computed via characteristic sets.
     4301// In the homogeneous case, the performance of the procedure
     4302// will be improved if I is already given by a minimal set of
     4303// generators. Apply minbase if necessary.
     4304//////////////////////////////////////////////////////////
     4305
     4306
     4307static proc prim_dec(ideal I, int choose)
     4308{
     4309  if((choose<0) or (choose>3))
     4310  {
     4311    ERROR("ERROR: <int> must be 0 or 1 or 2 or 3");
     4312  }
     4313  option(notWarnSB);
     4314  ideal H=1; // The intersection of the primary components
     4315  list U;    // the leaves of the decomposition tree, i.e.,
     4316             // pairs consisting of a primary component of I
     4317             // and the corresponding associated prime
     4318  list W;    // the non-leaf vertices in the decomposition tree.
     4319             // every entry has 6 components:
     4320                // 1- the vertex itself , i.e., a standard bais of the
     4321                //    given ideal I (type 1), or a standard basis of a
     4322                //    pseudo-primary component arising from
     4323                //    pseudo-primary decomposition (type 2), or a
     4324                //    standard basis of a remaining component arising from
     4325                //    pseudo-primary decomposition or extraction (type 3)
     4326                // 2- the type of the vertex as indicated above
     4327                // 3- the weighted_tree_depth of the vertex
     4328                // 4- the tester of the vertex
     4329                // 5- a standard basis of the associated prime
     4330                //    of a vertex of type 2, or 0 otherwise
     4331                // 6- a list of pairs consisting of a standard
     4332                //    basis of a minimal associated prime ideal
     4333                //    of the father of the vertex and the
     4334                //    irreducible factors of the "minimal
     4335                //    divisor" of the seperator or extractor
     4336                //    corresponding to the prime ideal
     4337                //    as computed by the procedure minsat,
     4338                //    if the vertex is of type 3, or
     4339                //    the empty list otherwise
     4340  ideal SI=std(I);
     4341  if(SI[1]==1)  // primdecSY(ideal(1))
     4342  {
     4343    return(list());
     4344  }
     4345  int ncolsSI=ncols(SI);
     4346  int ncolsH=1;
     4347  W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree
     4348  int weighted_tree_depth;
     4349  int i,j;
     4350  int check;
     4351  list V;  // current vertex
     4352  list VV; // new vertex
     4353  list QQ;
     4354  list WI;
     4355  ideal Qi,SQ,SRest,fac;
     4356  poly tester;
     4357
     4358  while(1)
     4359  {
     4360    i=1;
     4361    while(1)
     4362    {
     4363      while(i<=size(W)) // find vertex V of smallest weighted tree-depth
     4364      {
     4365        if (W[i][3]<=weighted_tree_depth) break;
     4366        i++;
     4367      }
     4368      if (i<=size(W)) break;
     4369      i=1;
     4370      weighted_tree_depth++;
     4371    }
     4372    V=W[i];
     4373    W=delete(W,i); // delete V from W
     4374
     4375    // now proceed by type of vertex V
     4376
     4377    if (V[2]==2)  // extraction needed
     4378    {
     4379      SQ,SRest,fac=extraction(V[1],V[5]);
     4380                        // standard basis of primary component,
     4381                        // standard basis of remaining component,
     4382                        // irreducible factors of
     4383                        // the "minimal divisor" of the extractor
     4384                        // as computed by the procedure minsat,
     4385      check=0;
     4386      for(j=1;j<=ncolsH;j++)
     4387      {
     4388        if (NF(H[j],SQ,1)!=0) // Q is not redundant
     4389        {
     4390          check=1;
     4391          break;
     4392        }
     4393      }
     4394      if(check==1)             // Q is not redundant
     4395      {
     4396        QQ=list();
     4397        QQ[1]=list(SQ,V[5]);  // primary component, associated prime,
     4398                              // i.e., standard bases thereof
     4399        U=U+QQ;
     4400        H=intersect(H,SQ);
     4401        H=std(H);
     4402        ncolsH=ncols(H);
     4403        check=0;
     4404        if(ncolsH==ncolsSI)
     4405        {
     4406          for(j=1;j<=ncolsSI;j++)
     4407          {
     4408            if(leadexp(H[j])!=leadexp(SI[j]))
     4409            {
     4410              check=1;
     4411              break;
     4412            }
     4413          }
     4414        }
     4415        else
     4416        {
     4417          check=1;
     4418        }
     4419        if(check==0) // H==I => U is a primary decomposition
     4420        {
     4421          return(U);
     4422        }
     4423      }
     4424      if (SRest[1]!=1)        // the remaining component is not
     4425                              // the whole ring
     4426      {
     4427        if (rad_con(V[4],SRest)==0) // the new vertex is not the
     4428                                    // root of a redundant subtree
     4429        {
     4430          VV[1]=SRest;     // remaining component
     4431          VV[2]=3;         // pseudoprimdec_special
     4432          VV[3]=V[3]+1;    // weighted depth
     4433          VV[4]=V[4];      // the tester did not change
     4434          VV[5]=ideal(0);
     4435          VV[6]=list(list(V[5],fac));
     4436          W=insert(W,VV,size(W));
     4437        }
     4438      }
     4439    }
     4440    else
     4441    {
     4442      if (V[2]==3) // pseudo_prim_dec_special is needed
     4443      {
     4444        QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose);
     4445                         // QQ = quadruples:
     4446                         // standard basis of pseudo-primary component,
     4447                         // standard basis of corresponding prime,
     4448                         // seperator, irreducible factors of
     4449                         // the "minimal divisor" of the seperator
     4450                         // as computed by the procedure minsat,
     4451                         // SRest=standard basis of remaining component
     4452      }
     4453      else     // V is the root, pseudo_prim_dec is needed
     4454      {
     4455        QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose);
     4456                         // QQ = quadruples:
     4457                         // standard basis of pseudo-primary component,
     4458                         // standard basis of corresponding prime,
     4459                         // seperator, irreducible factors of
     4460                         // the "minimal divisor" of the seperator
     4461                         // as computed by the procedure minsat,
     4462                         // SRest=standard basis of remaining component
     4463
     4464      }
     4465      //check
     4466      for(i=size(QQ);i>=1;i--)
     4467      //for(i=1;i<=size(QQ);i++)
     4468      {
     4469        tester=QQ[i][3]*V[4];
     4470        Qi=QQ[i][2];
     4471        if(NF(tester,Qi,1)!=0)  // the new vertex is not the
     4472                                // root of a redundant subtree
     4473        {
     4474          VV[1]=QQ[i][1];
     4475          VV[2]=2;
     4476          VV[3]=V[3]+1;
     4477          VV[4]=tester;      // the new tester as computed above
     4478          VV[5]=Qi;          // QQ[i][2];
     4479          VV[6]=list();
     4480          W=insert(W,VV,size(W));
     4481        }
     4482      }
     4483      if (SRest[1]!=1)        // the remaining component is not
     4484                              // the whole ring
     4485      {
     4486        if (rad_con(V[4],SRest)==0) // the vertex is not the root
     4487                                    // of a redundant subtree
     4488        {
     4489          VV[1]=SRest;
     4490          VV[2]=3;
     4491          VV[3]=V[3]+2;
     4492          VV[4]=V[4];      // the tester did not change
     4493          VV[5]=ideal(0);
     4494          WI=list();
     4495          for(i=1;i<=size(QQ);i++)
     4496          {
     4497            WI=insert(WI,list(QQ[i][2],QQ[i][4]));
     4498          }
     4499          VV[6]=WI;
     4500          W=insert(W,VV,size(W));
     4501        }
     4502      }
     4503    }
     4504  }
     4505}
     4506
     4507//////////////////////////////////////////////////////////////////////////
     4508// proc pseudo_prim_dec_charsets
     4509// input: Generators of an arbitrary ideal I, a standard basis SI of I,
     4510// and an integer choo
     4511// If choo=0, min_ass_prim_charsets with the given
     4512// ordering of the variables is used.
     4513// If choo=1, min_ass_prim_charsets with the "optimized"
     4514// ordering of the variables is used.
     4515// If choo=2, minAssPrimes from primdec.lib is used
     4516// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
     4517// output: a pseudo primary decomposition of I, i.e., a list
     4518// of pseudo primary components together with a standard basis of the
     4519// remaining component. Each pseudo primary component is
     4520// represented by a quadrupel: A standard basis of the component,
     4521// a standard basis of the corresponding associated prime, the
     4522// seperator of the component, and the irreducible factors of the
     4523// "minimal divisor" of the seperator as computed by the procedure minsat,
     4524// calls  proc pseudo_prim_dec_i
     4525//////////////////////////////////////////////////////////////////////////
     4526
     4527
     4528static proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo)
     4529{
     4530  list L;          // The list of minimal associated primes,
     4531                   // each one given by a standard basis
     4532  if((choo==0) or (choo==1))
     4533  {
     4534    L=min_ass_prim_charsets(I,choo);
     4535  }
     4536  else
     4537  {
     4538    if(choo==2)
     4539    {
     4540      L=minAssPrimes(I);
     4541    }
     4542    else
     4543    {
     4544      L=minAssPrimes(I,1);
     4545    }
     4546    for(int i=size(L);i>=1;i--)
     4547    {
     4548      L[i]=std(L[i]);
     4549    }
     4550  }
     4551  return (pseudo_prim_dec_i(SI,L));
     4552}
     4553
     4554////////////////////////////////////////////////////////////////
     4555// proc pseudo_prim_dec_special_charsets
     4556// input: a standard basis of an ideal I whose radical is the
     4557// intersection of the radicals of ideals generated by one prime ideal
     4558// P_i together with one polynomial f_i, the list V6 must be the list of
     4559// pairs (standard basis of P_i, irreducible factors of f_i),
     4560// and an integer choo
     4561// If choo=0, min_ass_prim_charsets with the given
     4562// ordering of the variables is used.
     4563// If choo=1, min_ass_prim_charsets with the "optimized"
     4564// ordering of the variables is used.
     4565// If choo=2, minAssPrimes from primdec.lib is used
     4566// If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
     4567// output: a pseudo primary decomposition of I, i.e., a list
     4568// of pseudo primary components together with a standard basis of the
     4569// remaining component. Each pseudo primary component is
     4570// represented by a quadrupel: A standard basis of the component,
     4571// a standard basis of the corresponding associated prime, the
     4572// seperator of the component, and the irreducible factors of the
     4573// "minimal divisor" of the seperator as computed by the procedure minsat,
     4574// calls  proc pseudo_prim_dec_i
     4575////////////////////////////////////////////////////////////////
     4576
     4577
     4578static proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo)
     4579{
     4580  int i,j,l;
     4581  list m;
     4582  list L;
     4583  int sizeL;
     4584  ideal P,SP; ideal fac;
     4585  int dimSP;
     4586  for(l=size(V6);l>=1;l--)   // creates a list of associated primes
     4587                             // of I, possibly redundant
     4588  {
     4589    P=V6[l][1];
     4590    fac=V6[l][2];
     4591    for(i=ncols(fac);i>=1;i--)
     4592    {
     4593      SP=P+fac[i];
     4594      SP=std(SP);
     4595      if(SP[1]!=1)
     4596      {
     4597        if((choo==0) or (choo==1))
     4598        {
     4599          m=min_ass_prim_charsets(SP,choo);  // a list of SB
     4600        }
     4601        else
     4602        {
     4603          if(choo==2)
     4604          {
     4605            m=minAssPrimes(SP);
     4606          }
     4607          else
     4608          {
     4609            m=minAssPrimes(SP,1);
     4610          }
     4611          for(j=size(m);j>=1;j=j-1)
     4612            {
     4613              m[j]=std(m[j]);
     4614            }
     4615        }
     4616        dimSP=dim(SP);
     4617        for(j=size(m);j>=1; j--)
     4618        {
     4619          if(dim(m[j])==dimSP)
     4620          {
     4621            L=insert(L,m[j],size(L));
     4622          }
     4623        }
     4624      }
     4625    }
     4626  }
     4627  sizeL=size(L);
     4628  for(i=1;i<sizeL;i++)     // get rid of redundant primes
     4629  {
     4630    for(j=i+1;j<=sizeL;j++)
     4631    {
     4632      if(size(L[i])!=0)
     4633      {
     4634        if(size(L[j])!=0)
     4635        {
     4636          if(size(NF(L[i],L[j],1))==0)
     4637          {
     4638            L[j]=ideal(0);
     4639          }
     4640          else
     4641          {
     4642            if(size(NF(L[j],L[i],1))==0)
     4643            {
     4644              L[i]=ideal(0);
     4645            }
     4646          }
     4647        }
     4648      }
     4649    }
     4650  }
     4651  for(i=sizeL;i>=1;i--)
     4652  {
     4653    if(size(L[i])==0)
     4654    {
     4655      L=delete(L,i);
     4656    }
     4657  }
     4658  return (pseudo_prim_dec_i(SI,L));
     4659}
     4660
     4661
     4662////////////////////////////////////////////////////////////////
     4663// proc pseudo_prim_dec_i
     4664// input: A standard basis of an arbitrary ideal I, and standard bases
     4665// of the minimal associated primes of I
     4666// output: a pseudo primary decomposition of I, i.e., a list
     4667// of pseudo primary components together with a standard basis of the
     4668// remaining component. Each pseudo primary component is
     4669// represented by a quadrupel: A standard basis of the component Q_i,
     4670// a standard basis of the corresponding associated prime P_i, the
     4671// seperator of the component, and the irreducible factors of the
     4672// "minimal divisor" of the seperator as computed by the procedure minsat,
     4673////////////////////////////////////////////////////////////////
     4674
     4675
     4676static proc pseudo_prim_dec_i (ideal SI, list L)
     4677{
     4678  list Q;
     4679  if (size(L)==1)               // one minimal associated prime only
     4680                                // the ideal is already pseudo primary
     4681  {
     4682    Q=SI,L[1],1;
     4683    list QQ;
     4684    QQ[1]=Q;
     4685    return (QQ,ideal(1));
     4686  }
     4687
     4688  poly f0,f,g;
     4689  ideal fac;
     4690  int i,j,k,l;
     4691  ideal SQi;
     4692  ideal I'=SI;
     4693  list QP;
     4694  int sizeL=size(L);
     4695  for(i=1;i<=sizeL;i++)
     4696  {
     4697    fac=0;
     4698    for(j=1;j<=sizeL;j++)           // compute the seperator sep_i
     4699                                    // of the i-th component
     4700    {
     4701      if (i!=j)                       // search g not in L[i], but L[j]
     4702      {
     4703        for(k=1;k<=ncols(L[j]);k++)
     4704        {
     4705          if(NF(L[j][k],L[i],1)!=0)
     4706          {
     4707            break;
     4708          }
     4709        }
     4710        fac=fac+L[j][k];
     4711      }
     4712    }
     4713    // delete superfluous polynomials
     4714    fac=simplify(fac,8+2);
     4715    // saturation
     4716    SQi,f0,f,fac=minsat_ppd(SI,fac);
     4717    I'=I',f;
     4718    QP=SQi,L[i],f0,fac;
     4719             // the quadrupel:
     4720             // a standard basis of Q_i,
     4721             // a standard basis of P_i,
     4722             // sep_i,
     4723             // irreducible factors of
     4724             // the "minimal divisor" of the seperator
     4725             //  as computed by the procedure minsat,
     4726    Q[i]=QP;
     4727  }
     4728  I'=std(I');
     4729  return (Q, I');
     4730                   // I' = remaining component
     4731}
     4732
     4733
     4734////////////////////////////////////////////////////////////////
     4735// proc extraction
     4736// input: A standard basis of a pseudo primary ideal I, and a standard
     4737// basis of the unique minimal associated prime P of I
     4738// output: an extraction of I, i.e., a standard basis of the primary
     4739// component Q of I with associated prime P, a standard basis of the
     4740// remaining component, and the irreducible factors of the
     4741// "minimal divisor" of the extractor as computed by the procedure minsat
     4742////////////////////////////////////////////////////////////////
     4743
     4744
     4745static proc extraction (ideal SI, ideal SP)
     4746{
     4747  list indsets=indepSet(SP,0);
     4748  poly f;
     4749  if(size(indsets)!=0)      //check, whether dim P != 0
     4750  {
     4751    intvec v;               // a maximal independent set of variables
     4752                            // modulo P
     4753    string U;               // the independent variables
     4754    string A;               // the dependent variables
     4755    int j,k;
     4756    int a;                  //  the size of A
     4757    int degf;
     4758    ideal g;
     4759    list polys;
     4760    int sizepolys;
     4761    list newpoly;
     4762    def R=basering;
     4763    //intvec hv=hilb(SI,1);
     4764    for (k=1;k<=size(indsets);k++)
     4765    {
     4766      v=indsets[k];
     4767      for (j=1;j<=nvars(R);j++)
     4768      {
     4769        if (v[j]==1)
     4770        {
     4771          U=U+varstr(j)+",";
     4772        }
     4773        else
     4774        {
     4775          A=A+varstr(j)+",";
     4776          a++;
     4777        }
     4778      }
     4779
     4780      U[size(U)]=")";           // we compute the extractor of I (w.r.t. U)
     4781      execute("ring RAU=("+charstr(basering)+"),("+A+U+",(dp("+string(a)+"),dp);");
     4782      ideal I=imap(R,SI);
     4783      //I=std(I,hv);            // the standard basis in (R[U])[A]
     4784      I=std(I);            // the standard basis in (R[U])[A]
     4785      A[size(A)]=")";
     4786      execute("ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;");
     4787      ideal I=imap(RAU,I);
     4788      //"std in lokalisierung:"+newline,I;
     4789      ideal h;
     4790      for(j=ncols(I);j>=1;j--)
     4791      {
     4792        h[j]=leadcoef(I[j]);  // consider I in (R(U))[A]
     4793      }
     4794      setring R;
     4795      g=imap(Rloc,h);
     4796      kill RAU,Rloc;
     4797      U="";
     4798      A="";
     4799      a=0;
     4800      f=lcm(g);
     4801      newpoly[1]=f;
     4802      polys=polys+newpoly;
     4803      newpoly=list();
     4804    }
     4805    f=polys[1];
     4806    degf=deg(f);
     4807    sizepolys=size(polys);
     4808    for (k=2;k<=sizepolys;k++)
     4809    {
     4810      if (deg(polys[k])<degf)
     4811      {
     4812        f=polys[k];
     4813        degf=deg(f);
     4814      }
     4815    }
     4816  }
     4817  else
     4818  {
     4819    f=1;
     4820  }
     4821  poly f0,h0; ideal SQ; ideal fac;
     4822  if(f!=1)
     4823  {
     4824    SQ,f0,h0,fac=minsat(SI,f);
     4825    return(SQ,std(SI+h0),fac);
     4826             // the tripel
     4827             // a standard basis of Q,
     4828             // a standard basis of remaining component,
     4829             // irreducible factors of
     4830             // the "minimal divisor" of the extractor
     4831             // as computed by the procedure minsat
     4832  }
     4833  else
     4834  {
     4835    return(SI,ideal(1),ideal(1));
     4836  }
     4837}
     4838
     4839/////////////////////////////////////////////////////
     4840// proc minsat
     4841// input:  a standard basis of an ideal I and a polynomial p
     4842// output: a standard basis IS of the saturation of I w.r. to p,
     4843// the maximal squarefree factor f0 of p,
     4844// the "minimal divisor" f of f0 such that the saturation of
     4845// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
     4846// the irreducible factors of f
     4847//////////////////////////////////////////////////////////
     4848
     4849
     4850static proc minsat(ideal SI, poly p)
     4851{
     4852  ideal fac=factorize(p,1);       //the irreducible factors of p
     4853  fac=sort(fac)[1];
     4854  int i,k;
     4855  poly f0=1;
     4856  for(i=ncols(fac);i>=1;i--)
     4857  {
     4858    f0=f0*fac[i];
     4859  }
     4860  poly f=1;
     4861  ideal iold;
     4862  list quotM;
     4863  quotM[1]=SI;
     4864  quotM[2]=fac;
     4865  quotM[3]=f0;
     4866  // we deal seperately with the first quotient;
     4867  // factors, which do not contribute to this one,
     4868  // are omitted
     4869  iold=quotM[1];
     4870  quotM=minquot(quotM);
     4871  fac=quotM[2];
     4872  if(quotM[3]==1)
     4873    {
     4874      return(quotM[1],f0,f,fac);
     4875    }
     4876  while(special_ideals_equal(iold,quotM[1])==0)
     4877    {
     4878      f=f*quotM[3];
     4879      iold=quotM[1];
     4880      quotM=minquot(quotM);
     4881    }
     4882  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
     4883}
     4884
     4885/////////////////////////////////////////////////////
     4886// proc minsat_ppd
     4887// input:  a standard basis of an ideal I and a polynomial p
     4888// output: a standard basis IS of the saturation of I w.r. to p,
     4889// the maximal squarefree factor f0 of p,
     4890// the "minimal divisor" f of f0 such that the saturation of
     4891// I w.r. to f equals the saturation of I w.r. to f0 (which is IS),
     4892// the irreducible factors of f
     4893//////////////////////////////////////////////////////////
     4894
     4895
     4896static proc minsat_ppd(ideal SI, ideal fac)
     4897{
     4898  fac=sort(fac)[1];
     4899  int i,k;
     4900  poly f0=1;
     4901  for(i=ncols(fac);i>=1;i--)
     4902  {
     4903    f0=f0*fac[i];
     4904  }
     4905  poly f=1;
     4906  ideal iold;
     4907  list quotM;
     4908  quotM[1]=SI;
     4909  quotM[2]=fac;
     4910  quotM[3]=f0;
     4911  // we deal seperately with the first quotient;
     4912  // factors, which do not contribute to this one,
     4913  // are omitted
     4914  iold=quotM[1];
     4915  quotM=minquot(quotM);
     4916  fac=quotM[2];
     4917  if(quotM[3]==1)
     4918    {
     4919      return(quotM[1],f0,f,fac);
     4920    }
     4921  while(special_ideals_equal(iold,quotM[1])==0)
     4922  {
     4923    f=f*quotM[3];
     4924    iold=quotM[1];
     4925    quotM=minquot(quotM);
     4926    k++;
     4927  }
     4928  return(quotM[1],f0,f,fac);           // the quadrupel ((I:p),f0,f, irr. factors of f)
     4929}
     4930/////////////////////////////////////////////////////////////////
     4931// proc minquot
     4932// input: a list with 3 components: a standard basis
     4933// of an ideal I, a set of irreducible polynomials, and
     4934// there product f0
     4935// output: a standard basis of the ideal (I:f0), the irreducible
     4936// factors of the "minimal divisor" f of f0 with (I:f0) = (I:f),
     4937// the "minimal divisor" f
     4938/////////////////////////////////////////////////////////////////
     4939
     4940static proc minquot(list tsil)
     4941{
     4942   intvec op;
     4943   int i,j,k,action;
     4944   ideal verg;
     4945   list l;
     4946   poly g;
     4947   ideal laedi=tsil[1];
     4948   ideal fac=tsil[2];
     4949   poly f=tsil[3];
     4950
     4951//std
     4952//   ideal star=quotient(laedi,f);
     4953//   star=std(star);
     4954   op=option(get);
     4955   option(returnSB);
     4956   ideal star=quotient(laedi,f);
     4957   option(set,op);
     4958   if(special_ideals_equal(laedi,star)==1)
     4959     {
     4960       return(laedi,ideal(1),1);
     4961     }
     4962   action=1;
     4963   while(action==1)
     4964   {
     4965      if(size(fac)==1)
     4966      {
     4967         action=0;
     4968         break;
     4969      }
     4970      for(i=1;i<=size(fac);i++)
     4971      {
     4972        g=1;
     4973         for(j=1;j<=size(fac);j++)
     4974         {
     4975            if(i!=j)
     4976            {
     4977               g=g*fac[j];
     4978            }
     4979         }
     4980//std
     4981//         verg=quotient(laedi,g);
     4982//         verg=std(verg);
     4983         op=option(get);
     4984         option(returnSB);
     4985         verg=quotient(laedi,g);
     4986         option(set,op);
     4987         if(special_ideals_equal(verg,star)==1)
     4988         {
     4989            f=g;
     4990            fac[i]=0;
     4991            fac=simplify(fac,2);
     4992            break;
     4993         }
     4994         if(i==size(fac))
     4995         {
     4996            action=0;
     4997         }
     4998      }
     4999   }
     5000   l=star,fac,f;
     5001   return(l);
     5002}
     5003/////////////////////////////////////////////////
     5004// proc special_ideals_equal
     5005// input: standard bases of ideal k1 and k2 such that
     5006// k1 is contained in k2, or k2 is contained ink1
     5007// output: 1, if k1 equals k2, 0 otherwise
     5008//////////////////////////////////////////////////
     5009
     5010static proc special_ideals_equal( ideal k1, ideal k2)
     5011{
     5012   int j;
     5013   if(size(k1)==size(k2))
     5014   {
     5015      for(j=1;j<=size(k1);j++)
     5016      {
     5017         if(leadexp(k1[j])!=leadexp(k2[j]))
     5018         {
     5019            return(0);
     5020         }
     5021      }
     5022      return(1);
     5023   }
     5024   return(0);
     5025}
     5026
     5027
     5028///////////////////////////////////////////////////////////////////////////////
     5029
     5030static proc convList(list l)
     5031{
     5032   int i;
     5033   list re,he;
     5034   for(i=1;i<=size(l)/2;i++)
     5035   {
     5036      he=l[2*i-1],l[2*i];
     5037      re[i]=he;
     5038   }
     5039   return(re);
     5040}
     5041///////////////////////////////////////////////////////////////////////////////
     5042
     5043static proc reconvList(list l)
     5044{
     5045   int i;
     5046   list re;
     5047   for(i=1;i<=size(l);i++)
     5048   {
     5049      re[2*i-1]=l[i][1];
     5050      re[2*i]=l[i][2];
     5051   }
     5052   return(re);
     5053}
     5054
     5055///////////////////////////////////////////////////////////////////////////////
     5056//
     5057//     The main procedures
     5058//
     5059///////////////////////////////////////////////////////////////////////////////
     5060
     5061proc primdecGTZ(ideal i)
     5062"USAGE:   primdecGTZ(i); i ideal
     5063RETURN:  a list pr of primary ideals and their associated primes:
     5064@format
     5065   pr[i][1]   the i-th primary component,
     5066   pr[i][2]   the i-th prime component.
     5067@end format
     5068NOTE:    Algorithm of Gianni/Trager/Zacharias.
     5069         Designed for characteristic 0, works also in char k > 0, if it
     5070         terminates (may result in an infinite loop in small characteristic!)
     5071EXAMPLE: example primdecGTZ; shows an example
     5072"
     5073{
     5074   if(attrib(basering,"global")!=1)
     5075   {
     5076      ERROR(
     5077      "// Not implemented for this ordering, please change to global ordering."
     5078      );
     5079   }
     5080   if(minpoly!=0)
     5081   {
     5082      return(algeDeco(i,0));
     5083      ERROR(
     5084      "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
     5085      );
     5086   }
     5087  return(convList(decomp(i)));
     5088}
     5089example
     5090{ "EXAMPLE:";  echo = 2;
     5091   ring  r = 0,(x,y,z),lp;
     5092   poly  p = z2+1;
     5093   poly  q = z3+2;
     5094   ideal i = p*q^2,y-z2;
     5095   list pr = primdecGTZ(i);
     5096   pr;
     5097}
     5098///////////////////////////////////////////////////////////////////////////////
     5099
     5100proc absPrimdecGTZ(ideal I)
     5101"USAGE:   absPrimdecGTZ(I); I ideal
     5102ASSUME:  Ground field has characteristic 0.
     5103RETURN:  a ring containing two lists: @code{absolute_primes} (the absolute
     5104         prime components of I) and @code{primary_decomp} (the output of
     5105         @code{primdecGTZ(I)}).
     5106         The list absolute_primes has to be interpreted as follows:
     5107         each entry describes a class of conjugated absolute primes,
     5108@format
     5109   absolute_primes[i][1]   the absolute prime component,
     5110   absolute_primes[i][2]   the number of conjugates.
     5111@end format
     5112         The first entry of @code{absolute_primes[i][1]} is the minimal
     5113         polynomial of a minimal finite field extension over which the
     5114         absolute prime component is defined.
     5115NOTE:    Algorithm of Gianni/Trager/Zacharias combined with the
     5116         @code{absFactorize} command.
     5117SEE ALSO: primdecGTZ; absFactorize
     5118EXAMPLE: example absPrimdecGTZ; shows an example
     5119"
     5120{
     5121  if (char(basering) != 0)
     5122  {
     5123    ERROR("primdec.lib::absPrimdecGTZ is only implemented for "+
     5124           +"characteristic 0");
     5125  }
     5126
     5127  if(attrib(basering,"global")!=1)
     5128  {
     5129    ERROR(
     5130      "// Not implemented for this ordering, please change to global ordering."
     5131    );
     5132  }
     5133  if(minpoly!=0)
     5134  {
     5135    //return(algeDeco(i,0));
     5136    ERROR(
     5137      "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
     5138    );
     5139  }
     5140  def R=basering;
     5141  int n=nvars(R);
     5142  list L=decomp(I,3);
     5143  string newvar=L[1][3];
     5144  int k=find(newvar,",",find(newvar,",")+1);
     5145  newvar=newvar[k+1..size(newvar)];
     5146  list lR=ringlist(R);
     5147  int i,d;
     5148  intvec vv;
     5149  for(i=1;i<=n;i++){vv[i]=1;}
     5150
     5151  list orst;
     5152  orst[1]=list("dp",vv);
     5153  orst[2]=list("dp",intvec(1));
     5154  orst[3]=list("C",0);
     5155  lR[3]=orst;
     5156  lR[2][n+1] = newvar;
     5157  def Rz = ring(lR);
     5158  setring Rz;
     5159  list L=imap(R,L);
     5160  list absolute_primes,primary_decomp;
     5161  ideal I,M,N,K;
     5162  M=maxideal(1);
     5163  N=maxideal(1);
     5164  poly p,q,f,g;
     5165  map phi,psi;
     5166  for(i=1;i<=size(L);i++)
     5167  {
     5168    I=L[i][2];
     5169    execute("K="+L[i][3]+";");
     5170    p=K[1];
     5171    q=K[2];
     5172    execute("f="+L[i][4]+";");
     5173    g=2*var(n)-f;
     5174    M[n]=f;
     5175    N[n]=g;
     5176    d=deg(p);
     5177    phi=Rz,M;
     5178    psi=Rz,N;
     5179    I=phi(I),p,q;
     5180    I=std(I);
     5181    absolute_primes[i]=list(psi(I),d);
     5182    primary_decomp[i]=list(L[i][1],L[i][2]);
     5183  }
     5184  export(primary_decomp);
     5185  export(absolute_primes);
     5186  setring R;
     5187  dbprint( printlevel-voice+3,"
     5188// 'absPrimdecGTZ' created a ring, in which two lists absolute_primes (the
     5189// absolute prime components) and primary_decomp (the primary and prime
     5190// components over the current basering) are stored.
     5191// To access the list of absolute prime components, type (if the name S was
     5192// assigned to the return value):
     5193        setring S; absolute_primes; ");
     5194
     5195  return(Rz);
     5196}
     5197example
     5198{ "EXAMPLE:";  echo = 2;
     5199   ring  r = 0,(x,y,z),lp;
     5200   poly  p = z2+1;
     5201   poly  q = z3+2;
     5202   ideal i = p*q^2,y-z2;
     5203   def S = absPrimdecGTZ(i);
     5204   setring S;
     5205   absolute_primes;
     5206}
     5207///////////////////////////////////////////////////////////////////////////////
     5208
     5209proc primdecSY(ideal i, list #)
     5210"USAGE:   primdecSY(I, c); I ideal, c int (optional)
     5211RETURN:  a list pr of primary ideals and their associated primes:
     5212@format
     5213   pr[i][1]   the i-th primary component,
     5214   pr[i][2]   the i-th prime component.
     5215@end format
     5216NOTE:    Algorithm of Shimoyama/Yokoyama.
     5217@format
     5218   if c=0,  the given ordering of the variables is used,
     5219   if c=1,  minAssChar tries to use an optimal ordering (default),
     5220   if c=2,  minAssGTZ is used,
     5221   if c=3,  minAssGTZ and facstd are used.
     5222@end format
     5223EXAMPLE: example primdecSY; shows an example
     5224"
     5225{
     5226   if(attrib(basering,"global")!=1)
     5227   {
     5228      ERROR(
     5229      "// Not implemented for this ordering, please change to global ordering."
     5230      );
     5231   }
     5232   i=simplify(i,2);
     5233   if ((i[1]==0)||(i[1]==1))
     5234   {
     5235     list L=list(ideal(i[1]),ideal(i[1]));
     5236     return(list(L));
     5237   }
     5238   if(minpoly!=0)
     5239   {
     5240      return(algeDeco(i,1));
     5241   }
     5242   if (size(#)==1)
     5243   { return(prim_dec(i,#[1])); }
     5244   else
     5245   { return(prim_dec(i,1)); }
     5246}
     5247example
     5248{ "EXAMPLE:";  echo = 2;
     5249   ring  r = 0,(x,y,z),lp;
     5250   poly  p = z2+1;
     5251   poly  q = z3+2;
     5252   ideal i = p*q^2,y-z2;
     5253   list pr = primdecSY(i);
     5254   pr;
     5255}
     5256///////////////////////////////////////////////////////////////////////////////
     5257proc minAssGTZ(ideal i,list #)
     5258"USAGE:    minAssGTZ(I[, l]); I ideal, l list (optional)
     5259   @* Optional parameters in list l (can be entered in any order):
     5260   @* 0, \"facstd\" -> uses facstd to first decompose the ideal (default)
     5261   @* 1, \"noFacstd\" -> does not use facstd
     5262   @* \"GTZ\" -> the original algorithm by Gianni, Trager and Zacharias is used
     5263   @* \"SL\" -> GTZ algorithm with modificiations by Laplagne is used (default)
     5264
     5265RETURN:  a list, the minimal associated prime ideals of I.
     5266NOTE:    Designed for characteristic 0, works also in char k > 0 based
     5267         on an algorithm of Yokoyama
     5268EXAMPLE: example minAssGTZ; shows an example
     5269"
     5270{
     5271  int j;
     5272  string algorithm;
     5273  string facstdOption;
     5274  int useFac;
     5275
     5276  // Set input parameters
     5277  algorithm = "SL";         // Default: SL algorithm
     5278  facstdOption = "facstd";
     5279  if(size(#) > 0)
     5280  {
     5281    int valid;
     5282    for(j = 1; j <= size(#); j++)
     5283    {
     5284      valid = 0;
     5285      if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
     5286      {
     5287        if (#[j] == 1) {facstdOption = "noFacstd"; valid = 1;}    // If #[j] == 1, facstd is not used.
     5288        if (#[j] == 0) {facstdOption = "facstd";   valid = 1;}    // If #[j] == 0, facstd is used.
     5289      }
     5290      if(typeof(#[j]) == "string")
     5291      {
     5292        if((#[j] == "GTZ") || (#[j] == "SL"))
     5293        {
     5294          algorithm = #[j];
     5295          valid = 1;
     5296        }
     5297        if((#[j] == "noFacstd") || (#[j] == "facstd"))
     5298        {
     5299          facstdOption = #[j];
     5300          valid = 1;
     5301        }
     5302      }
     5303      if(valid == 0)
     5304      {
     5305        dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
     5306      }
     5307    }
     5308  }
     5309
     5310  if(attrib(basering,"global")!=1)
     5311  {
     5312    ERROR(
     5313      "// Not implemented for this ordering, please change to global ordering."
     5314    );
     5315  }
     5316  if(minpoly!=0)
     5317  {
     5318    return(algeDeco(i,2));
     5319  }
     5320
     5321  list result = minAssPrimes(i, facstdOption, algorithm);
     5322  return(result);
     5323}
     5324example
     5325{ "EXAMPLE:";  echo = 2;
     5326   ring  r = 0,(x,y,z),dp;
     5327   poly  p = z2+1;
     5328   poly  q = z3+2;
     5329   ideal i = p*q^2,y-z2;
     5330   list pr = minAssGTZ(i);
     5331   pr;
     5332}
     5333
     5334///////////////////////////////////////////////////////////////////////////////
     5335proc minAssChar(ideal i, list #)
     5336"USAGE:   minAssChar(I[,c]); i ideal, c int (optional).
     5337RETURN:  list, the minimal associated prime ideals of i.
     5338NOTE:    If c=0, the given ordering of the variables is used. @*
     5339         Otherwise, the system tries to find an optimal ordering,
     5340         which in some cases may considerably speed up the algorithm. @*
     5341         Due to a bug in the factorization, the result may be not completely
     5342         decomposed in small characteristic.
     5343EXAMPLE: example minAssChar; shows an example
     5344"
     5345{
     5346   if(attrib(basering,"global")!=1)
     5347   {
     5348      ERROR(
     5349      "// Not implemented for this ordering, please change to global ordering."
     5350      );
     5351   }
     5352   if (size(#)==1)
     5353   { return(min_ass_prim_charsets(i,#[1])); }
     5354   else
     5355   { return(min_ass_prim_charsets(i,1)); }
     5356}
     5357example
     5358{ "EXAMPLE:";  echo = 2;
     5359   ring  r = 0,(x,y,z),dp;
     5360   poly  p = z2+1;
     5361   poly  q = z3+2;
     5362   ideal i = p*q^2,y-z2;
     5363   list pr = minAssChar(i);
     5364   pr;
     5365}
     5366///////////////////////////////////////////////////////////////////////////////
     5367proc equiRadical(ideal i)
     5368"USAGE:   equiRadical(I); I ideal
     5369RETURN:  ideal, intersection of associated primes of I of maximal dimension.
     5370NOTE:    A combination of the algorithms of Krick/Logar (with modifications by Laplagne) and Kemper is used.
     5371         Works also in positive characteristic (Kempers algorithm).
     5372EXAMPLE: example equiRadical; shows an example
     5373"
     5374{
     5375  if(attrib(basering,"global")!=1)
     5376  {
     5377     ERROR(
     5378     "// Not implemented for this ordering, please change to global ordering."
     5379     );
     5380  }
     5381  return(radical(i, 1));
     5382}
     5383example
     5384{ "EXAMPLE:";  echo = 2;
     5385   ring  r = 0,(x,y,z),dp;
     5386   poly  p = z2+1;
     5387   poly  q = z3+2;
     5388   ideal i = p*q^2,y-z2;
     5389   ideal pr= equiRadical(i);
     5390   pr;
     5391}
     5392
     5393///////////////////////////////////////////////////////////////////////////////
     5394proc radical(ideal i, list #)
     5395"USAGE: radical(I[, l]); I ideal, l list (optional)
     5396 @*  Optional parameters in list l (can be entered in any order):
     5397 @*  0, \"fullRad\" -> full radical is computed (default)
     5398 @*  1, \"equiRad\" -> equiRadical is computed
     5399 @*  \"KL\" -> Krick/Logar algorithm is used
     5400 @*  \"SL\" -> modifications by Laplagne are used (default)
     5401 @*  \"facstd\" -> uses facstd to first decompose the ideal (default for non homogeneous ideals)
     5402 @*  \"noFacstd\" -> does not use facstd (default for homogeneous ideals)
     5403RETURN:  ideal, the radical of I (or the equiradical if required in the input parameters)
     5404NOTE:    A combination of the algorithms of Krick/Logar (with modifications by Laplagne) and Kemper is used.
     5405         Works also in positive characteristic (Kempers algorithm).
     5406EXAMPLE: example radical; shows an example
     5407"
     5408{
     5409  dbprint(printlevel - voice, "Radical, version 2006.05.08");
     5410  if(attrib(basering,"global")!=1)
     5411  {
     5412    ERROR(
     5413     "// Not implemented for this ordering, please change to global ordering."
     5414    );
     5415  }
     5416  if(size(i) == 0){return(ideal(0));}
     5417  int j;
     5418  def P0 = basering;
     5419  list Pl=ringlist(P0);
     5420  intvec dp_w;
     5421  for(j=nvars(P0);j>0;j--) {dp_w[j]=1;}
     5422  Pl[3]=list(list("dp",dp_w),list("C",0));
     5423  def @P=ring(Pl);
     5424  setring @P;
     5425  ideal i=imap(P0,i);
     5426
     5427  int il;
     5428  string algorithm;
     5429  int useFac;
     5430
     5431  // Set input parameters
     5432  algorithm = "SL";                                 // Default: SL algorithm
     5433  il = 0;                                           // Default: Full radical (not only equiRadical)
     5434  if (homog(i) == 1)
     5435  {   // Default: facStd is used, except if the ideal is homogeneous.
     5436    useFac = 0;
     5437  }
     5438  else
     5439  {
     5440    useFac = 1;
     5441  }
     5442  if(size(#) > 0)
     5443  {
     5444    int valid;
     5445    for(j = 1; j <= size(#); j++)
     5446    {
     5447      valid = 0;
     5448      if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
     5449      {
     5450        il = #[j];          // If il == 1, equiRadical is computed
     5451        valid = 1;
     5452      }
     5453      if(typeof(#[j]) == "string")
     5454      {
     5455        if(#[j] == "KL")
     5456        {
     5457          algorithm = "KL";
     5458          valid = 1;
     5459        }
     5460        if(#[j] == "SL")
     5461        {
     5462          algorithm = "SL";
     5463          valid = 1;
     5464        }
     5465        if(#[j] == "noFacstd")
     5466        {
     5467          useFac = 0;
     5468          valid = 1;
     5469        }
     5470        if(#[j] == "facstd")
     5471        {
     5472          useFac = 1;
     5473          valid = 1;
     5474        }
     5475        if(#[j] == "equiRad")
     5476        {
     5477          il = 1;
     5478          valid = 1;
     5479        }
     5480        if(#[j] == "fullRad")
     5481        {
     5482          il = 0;
     5483          valid = 1;
     5484        }
     5485      }
     5486      if(valid == 0)
     5487      {
     5488        dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
     5489      }
     5490    }
     5491  }
     5492
     5493  ideal rad = 1;
     5494  intvec op = option(get);
     5495  list qr = simplifyIdeal(i);
     5496  map phi = @P, qr[2];
     5497
     5498  option(redSB);
     5499  i = groebner(qr[1]);
     5500  option(set, op);
     5501  int di = dim(i);
     5502
     5503  if(di == 0)
     5504  {
     5505    i = zeroRad(i, qr[1]);
     5506    i=interred(phi(i));
     5507    setring(P0);
     5508    i=imap(@P,i);
     5509    return(i);
     5510  }
     5511
     5512  option(redSB);
     5513  list pr;
     5514  if(useFac == 1)
     5515  {
     5516    pr = facstd(i);
     5517  }
     5518  else
     5519  {
     5520    pr = i;
     5521  }
     5522  option(set, op);
     5523  int s = size(pr);
     5524  if(useFac == 1)
     5525  {
     5526    dbprint(printlevel - voice, "Number of components returned by facstd: ", s);
     5527  }
     5528  for(j = 1; j <= s; j++)
     5529  {
     5530    attrib(pr[s + 1 - j], "isSB", 1);
     5531    if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il))
     5532    {
     5533      // SL Debug messages
     5534      dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
     5535      dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));
     5536
     5537      if(algorithm == "KL")
     5538      {
     5539        rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
     5540      }
     5541      if(algorithm == "SL")
     5542      {
     5543        rad = intersect(rad, radicalSL(pr[s + 1 - j], il));
     5544      }
     5545    }
     5546    else
     5547    {
     5548      // SL Debug
     5549      dbprint(printlevel-voice, "The radical of this component is not needed.");
     5550      dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))",
     5551              size(reduce(rad, pr[s + 1 - j], 1)));
     5552      dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
     5553      dbprint(printlevel-voice, "il", il);
     5554    }
     5555  }
     5556  rad=interred(phi(rad));
     5557  setring(P0);
     5558  i=imap(@P,rad);
     5559  return(i);
     5560}
     5561example
     5562{ "EXAMPLE:";  echo = 2;
     5563   ring  r = 0,(x,y,z),dp;
     5564   poly  p = z2+1;
     5565   poly  q = z3+2;
     5566   ideal i = p*q^2,y-z2;
     5567   ideal pr = radical(i);
     5568   pr;
     5569}
     5570
     5571///////////////////////////////////////////////////////////////////////////////
     5572//
     5573// Computes the radical of I using KL algorithm.
     5574// The only difference with the the previous implementation of KL algorithm is
     5575// that now it uses block dp instead of lp ordering for the reduction to the
     5576// zerodimensional case.
     5577// The reduction step has been moved to the new routine radicalReduction, so that it can be
     5578// used also by radicalSL procedure.
     5579//
     5580static proc radicalKL(ideal I, ideal ser, list #)
     5581{
     5582// ideal I     The ideal for which the radical is computed
     5583// ideal ser   Used to reduce components already obtained
     5584// list #      If #[1] = 1, equiradical is computed.
     5585
     5586  // I needs to be a Groebner basis.
     5587  if (attrib(I, "isSB") != 1)
     5588  {
     5589    I = groebner(I);
     5590  }
     5591
     5592  ideal rad;                                // The radical
     5593  int allIndep = 1;                // All max independent sets are used
     5594
     5595  list result = radicalReduction(I, ser, allIndep, #);
     5596  int done = result[3];
     5597  rad = result[1];
     5598  if (done == 0)
     5599  {
     5600    rad = intersect(rad, radicalKL(result[2], ideal(1), #));
     5601  }
     5602  return(rad);
     5603}
     5604
     5605
     5606///////////////////////////////////////////////////////////////////////////////
     5607//
     5608// Computes the radical of I via Laplagne algorithm, using zerodimensional radical in
     5609// the zero dimensional case.
     5610// For the reduction to the zerodimensional case, it uses the procedure
     5611// radical, with some modifications to avoid the recursion.
     5612//
     5613static proc radicalSL(ideal I, list #)
     5614// Input = I, ideal
     5615//         #, list. If #[1] = 1, then computes only the equiradical.
     5616// Output = (P, primaryDec) where P = rad(I) and primaryDec is the list of the radicals
     5617// obtained in intermediate steps.
     5618{
     5619  ideal rad = 1;
     5620  ideal equiRad = 1;
     5621  list primes;
     5622  int k;                        // Counter
     5623  int il;                 // If il = 1, only the equiradical is required.
     5624  int iDim;                // The dimension of I
     5625  int stop = 0;   // Checks if the radical has been obtained
     5626
     5627  if (attrib(I, "isSB") != 1)
     5628  {
     5629    I = groebner(I);
     5630  }
     5631  iDim = dim(I);
     5632
     5633  // Checks if only equiradical is required
     5634  if (size(#) > 0)
     5635  {
     5636    il = #[1];
     5637  }
     5638
     5639  while(stop == 0)
     5640  {
     5641    dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals.");
     5642    primes = radicalSLIteration(I, rad);                         // A list of primes or intersections of primes, not included in P
     5643    dbprint (printlevel - voice, "// Output of Iteration Step:");
     5644    dbprint (printlevel - voice, primes);
     5645    if (size(primes) > 0)
     5646    {
     5647      dbprint (printlevel - voice, "// We intersect P with the ideal just obtained.");
     5648      for(k = 1; k <= size(primes); k++)
     5649      {
     5650        rad = intersect(rad, primes[k]);
     5651        if (il == 1)
     5652        {
     5653          if (attrib(primes[k], "isSB") != 1)
     5654          {
     5655            primes[k] = groebner(primes[k]);
     5656          }
     5657          if (iDim == dim(primes[k]))
     5658          {
     5659            equiRad = intersect(equiRad, primes[k]);
     5660          }
     5661        }
     5662      }
     5663    }
     5664    else
     5665    {
     5666      stop = 1;
     5667    }
     5668  }
     5669  if (il == 0)
     5670  {
     5671    return(rad);
     5672  }
     5673  else
     5674  {
     5675    return(equiRad);
     5676  }
     5677}
     5678
     5679//////////////////////////////////////////////////////////////////////////
     5680// Based on radicalKL.
     5681// It contains all of old version of proc radicalKL except the recursion call.
     5682//
     5683// Output:
     5684// #1 -> output ideal, the part of the radical that has been computed
     5685// #2 -> complementary ideal, the part of the ideal I whose radical remains to be computed
     5686//       = (I, h) in KL algorithm
     5687//       This is not used in the new algorithm. It is part of KL algorithm
     5688// #3 -> done, 1: output = radical, there is no need to continue
     5689//                   0: radical = output \cap \sqrt{complementary ideal}
     5690//       This is not used in the new algorithm. It is part of KL algorithm
     5691
     5692static proc radicalReduction(ideal I, ideal ser, int allIndep, list #)
     5693{
     5694// allMaximal                1 -> Indicates that the reduction to the zerodim case
     5695//                       must be done for all indep set of the leading terms ideal
     5696//                                        0 -> Otherwise
     5697// ideal ser                Only for radicalKL. (Same as in radicalKL)
     5698// list #                        Only for radicalKL (If #[1] = 1, only equiradical is required.
     5699//                  It is used to set the value of done.)
     5700
     5701  attrib(I, "isSB", 1);   // I needs to be a reduced standard basis
     5702  list indep, fett;
     5703  intvec @w, @hilb, op;
     5704  int @wr, @n, @m, lauf, di;
     5705  ideal fac, @h, collectrad, lsau;
     5706  poly @q;
     5707  string @va, quotring;
     5708
     5709  def @P = basering;
     5710  int jdim = dim(I);               // Computes the dimension of I
     5711  int  homo = homog(I);            // Finds out if I is homogeneous
     5712  ideal rad = ideal(1);            // The unit ideal
     5713  ideal te = ser;
     5714  if(size(#) > 0)
     5715  {
     5716    @wr = #[1];
     5717  }
     5718  if(homo == 1)
     5719  {
     5720    for(@n = 1; @n <= nvars(basering); @n++)
     5721    {
     5722      @w[@n] = ord(var(@n));
     5723    }
     5724    @hilb = hilb(I, 1, @w);
     5725  }
     5726
     5727  // SL 2006.04.11 1 Debug messages
     5728  dbprint(printlevel-voice, "//Computes the radical of the ideal:", I);
     5729  // SL 2006.04.11 2 Debug messages
     5730
     5731  //---------------------------------------------------------------------------
     5732  //j is the ring
     5733  //---------------------------------------------------------------------------
     5734
     5735  if (jdim==-1)
     5736  {
     5737    return(ideal(1), ideal(1), 1);
     5738  }
     5739
     5740  //---------------------------------------------------------------------------
     5741  //the zero-dimensional case
     5742  //---------------------------------------------------------------------------
     5743
     5744  if (jdim==0)
     5745  {
     5746    return(zeroRad(I), ideal(1), 1);
     5747  }
     5748
     5749  //-------------------------------------------------------------------------
     5750  //search for a maximal independent set indep,i.e.
     5751  //look for subring such that the intersection with the ideal is zero
     5752  //j intersected with K[var(indep[3]+1),...,var(nvar)] is zero,
     5753  //indep[1] is the new varstring, indep[2] the string for the block-ordering
     5754  //-------------------------------------------------------------------------
     5755
     5756  // SL 2006-04-24 1         If allIndep = 0, then it only computes one maximal
     5757  //                                        independent set.
     5758  //                                        This looks better for the new algorithm but not for KL
     5759  //                                        algorithm
     5760  list parameters = allIndep;
     5761  indep = newMaxIndependSetDp(I, parameters);
     5762  // SL 2006-04-24 2
     5763
     5764  for(@m = 1; @m <= size(indep); @m++)
     5765  {
     5766    if((indep[@m][1] == varstr(basering)) && (@m == 1))
     5767    //this is the good case, nothing to do, just to have the same notations
     5768    //change the ring
     5769    {
     5770      execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
     5771                              +ordstr(basering)+");");
     5772      ideal @j = fetch(@P, I);
     5773      attrib(@j, "isSB", 1);
     5774    }
     5775    else
     5776    {
     5777      @va = string(maxideal(1));
     5778
     5779      execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),("
     5780                              + indep[@m][2] + ");");
     5781      execute("map phi = @P," + @va + ";");
     5782      if(homo == 1)
     5783      {
     5784        ideal @j = std(phi(I), @hilb, @w);
     5785      }
     5786      else
     5787      {
     5788        ideal @j = groebner(phi(I));
     5789      }
     5790    }
     5791    if((deg(@j[1]) == 0) || (dim(@j) < jdim))
     5792    {
     5793      setring @P;
     5794      break;
     5795    }
     5796    for (lauf = 1; lauf <= size(@j); lauf++)
     5797    {
     5798      fett[lauf] = size(@j[lauf]);
     5799    }
     5800    //------------------------------------------------------------------------
     5801    // We have now the following situation:
     5802    // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     5803    // to this quotientring, j is there still a standardbasis, the
     5804    // leading coefficients of the polynomials there (polynomials in
     5805    // K[var(nnp+1),..,var(nva)]) are collected in the list h,
     5806    // we need their LCM, gh, because of the following:
     5807    // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
     5808    // intersected with K[var(1),...,var(nva)] is (j:gh^n)
     5809    // on the other hand j = ((j, gh^n) intersected with (j : gh^n))
     5810
     5811    //------------------------------------------------------------------------
     5812    // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
     5813    // and the map phi:K[var(1),...,var(nva)] ----->
     5814    // K(var(nnpr+1),..,var(nva))[..the rest..]
     5815    //------------------------------------------------------------------------
     5816    quotring = prepareQuotientRingDp(nvars(basering) - indep[@m][3]);
     5817    //------------------------------------------------------------------------
     5818    // We pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     5819    //------------------------------------------------------------------------
     5820
     5821    execute(quotring);
     5822
     5823    // @j considered in the quotientring
     5824    ideal @j = imap(gnir1, @j);
     5825
     5826    kill gnir1;
     5827
     5828    // j is a standardbasis in the quotientring but usually not minimal
     5829    // here it becomes minimal
     5830
     5831    @j = clearSB(@j, fett);
     5832
     5833    // We need later LCM(h[1],...) = gh for saturation
     5834    ideal @h;
     5835    if(deg(@j[1]) > 0)
     5836    {
     5837      for(@n = 1; @n <= size(@j); @n++)
     5838      {
     5839        @h[@n] = leadcoef(@j[@n]);
     5840      }
     5841      op = option(get);
     5842      option(redSB);
     5843      @j = interred(@j);  //to obtain a reduced standardbasis
     5844      attrib(@j, "isSB", 1);
     5845      option(set, op);
     5846
     5847      // SL 1 Debug messages
     5848      dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j)));
     5849      ideal zero_rad = zeroRad(@j);
     5850      dbprint(printlevel - voice, "zero_rad passed");
     5851      // SL 2
     5852    }
     5853    else
     5854    {
     5855      ideal zero_rad = ideal(1);
     5856    }
     5857
     5858    // We need the intersection of the ideals in the list quprimary with the
     5859    // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     5860    // but fi polynomials, then the intersection of q with the polynomialring
     5861    // is the saturation of the ideal generated by f1,...,fr with respect to
     5862    // h which is the lcm of the leading coefficients of the fi considered in
     5863    // the quotientring: this is coded in saturn
     5864
     5865    zero_rad = std(zero_rad);
     5866
     5867    ideal hpl;
     5868
     5869    for(@n = 1; @n <= size(zero_rad); @n++)
     5870    {
     5871      hpl = hpl, leadcoef(zero_rad[@n]);
     5872    }
     5873
     5874    //------------------------------------------------------------------------
     5875    // We leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     5876    // back to the polynomialring
     5877    //------------------------------------------------------------------------
     5878    setring @P;
     5879
     5880    collectrad = imap(quring, zero_rad);
     5881    lsau = simplify(imap(quring, hpl), 2);
     5882    @h = imap(quring, @h);
     5883
     5884    kill quring;
     5885
     5886    // Here the intersection with the polynomialring
     5887    // mentioned above is really computed
     5888
     5889    collectrad = sat2(collectrad, lsau)[1];
     5890    if(deg(@h[1])>=0)
     5891    {
     5892      fac = ideal(0);
     5893      for(lauf = 1; lauf <= ncols(@h); lauf++)
     5894      {
     5895        if(deg(@h[lauf]) > 0)
     5896        {
     5897          fac = fac + factorize(@h[lauf], 1);
     5898        }
     5899      }
     5900      fac = simplify(fac, 6);
     5901      @q = 1;
     5902      for(lauf = 1; lauf <= size(fac); lauf++)
     5903      {
     5904        @q = @q * fac[lauf];
     5905      }
     5906      op = option(get);
     5907      option(returnSB);
     5908      option(redSB);
     5909      I = quotient(I + ideal(@q), rad);
     5910      attrib(I, "isSB", 1);
     5911      option(set, op);
     5912    }
     5913    if((deg(rad[1]) > 0) && (deg(collectrad[1]) > 0))
     5914    {
     5915      rad = intersect(rad, collectrad);
     5916      te = intersect(te, collectrad);
     5917      te = simplify(reduce(te, I, 1), 2);
     5918    }
     5919    else
     5920    {
     5921      if(deg(collectrad[1]) > 0)
     5922      {
     5923        rad = collectrad;
     5924        te = intersect(te, collectrad);
     5925        te = simplify(reduce(te, I, 1), 2);
     5926      }
     5927    }
     5928
     5929    if((dim(I) < jdim)||(size(te) == 0))
     5930    {
     5931      break;
     5932    }
     5933    if(homo==1)
     5934    {
     5935      @hilb = hilb(I, 1, @w);
     5936    }
     5937  }
     5938
     5939  // SL 2006.04.11 1 Debug messages
     5940  dbprint (printlevel-voice, "// Part of the Radical already computed:", rad);
     5941  dbprint (printlevel-voice, "// Dimension:", dim(groebner(rad)));
     5942  // SL 2006.04.11 2 Debug messages
     5943
     5944  // SL 2006.04.21 1         New variable "done".
     5945  //                                        It tells if the radical is already computed or
     5946  //                                        if it still has to be computed the radical of the new ideal I
     5947  int done;
     5948  if(((@wr == 1) && (dim(I)<jdim)) || (deg(I[1])==0) || (size(te) == 0))
     5949  {
     5950    done = 1;
     5951  }
     5952  else
     5953  {
     5954    done = 0;
     5955  }
     5956  // SL 2006.04.21 2
     5957
     5958  // SL 2006.04.21 1         See details of the output at the beggining of this proc.
     5959  list result = rad, I, done;
     5960  return(result);
     5961  // SL 2006.04.21 2
     5962}
     5963
     5964
     5965
     5966///////////////////////////////////////////////////////////////////////////////
     5967// Given an ideal I and an ideal P (intersection of some minimal prime ideals
     5968// associated to I), it calculates the intersection of new minimal prime ideals
     5969// associated to I which where not used to calculate P.
     5970// This version uses ZD Radical in the zerodimensional case.
     5971static proc radicalSLIteration (ideal I, ideal P);
     5972// Input: I, ideal. The ideal from which new prime components will be obtained.
     5973//        P, ideal. Intersection of some prime ideals of I.
     5974// Output: ideal. Intersection of some primes of I different from the ones in P.
     5975{
     5976  int k = 1;                                // Counter
     5977  int good  = 0;                        // Checks if an element of P is in rad(I)
     5978
     5979  dbprint (printlevel-voice, "// We search for an element in P - sqrt(I).");
     5980  while ((k <= size(P)) and (good == 0))
     5981  {
     5982    dbprint (printlevel-voice, "// We try with:", P[k]);
     5983    good = 1 - rad_con(P[k], I);
     5984    k++;
     5985  }
     5986  k--;
     5987  if (good == 0)
     5988  {
     5989    dbprint (printlevel-voice, "// No element was found, P = sqrt(I).");
     5990    list emptyList = list();
     5991    return (emptyList);
     5992  }
     5993  dbprint(printlevel - voice, "// That one was good!");
     5994  dbprint(printlevel - voice, "// We saturate I with respect to this element.");
     5995  if (P[k] != 1)
     5996  {
     5997    ideal J = sat(I, P[k])[1];
     5998  }
     5999  else
     6000  {
     6001    dbprint(printlevel - voice, "// The polynomial is 1, the saturation in not actually computed.");
     6002    ideal J = I;
     6003  }
     6004
     6005  // We now call proc radicalNew;
     6006  dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via radical.");
     6007  dbprint(printlevel - voice, "// The ideal is ", J);
     6008  dbprint(printlevel - voice, "// The dimension is ", dim(groebner(J)));
     6009
     6010  int allMaximal = 0;                   // Compute the zerodim reduction for only one indep set.
     6011  ideal re = 1;                           // No reduction is need, there are not redundant components.
     6012  list emptyList = list();   // Look for primes of any dimension, not only of max dimension.
     6013  list result = radicalReduction(J, re, allMaximal, emptyList);
     6014
     6015  return(result[1]);
     6016}
     6017
     6018///////////////////////////////////////////////////////////////////////////////////
     6019// Based on maxIndependSet
     6020// Added list # as parameter
     6021// If the first element of # is 0, the output is only 1 max indep set.
     6022// If no list is specified or #[1] = 1, the output is all the max indep set of the
     6023// leading terms ideal. This is the original output of maxIndependSet
     6024
     6025// The ordering given in the output has been changed to block dp instead of lp.
     6026
     6027proc newMaxIndependSetDp(ideal j, list #)
     6028"USAGE:   newMaxIndependentSetDp(I); I ideal (returns all maximal independent sets of the corresponding leading terms ideal)
     6029          newMaxIndependentSetDp(I, 0); I ideal (returns only one maximal independent set)
     6030RETURN:  list = #1. new varstring with the maximal independent set at the end,
     6031                #2. ordstring with the corresponding dp block ordering,
     6032                #3. the number of independent variables
     6033NOTE:
     6034EXAMPLE: example newMaxIndependentSetDp; shows an example
     6035"
     6036{
     6037  int n, k, di;
     6038  list resu, hilf;
     6039  string var1, var2;
     6040  list v = indepSet(j, 0);
     6041
     6042  // SL 2006.04.21 1 Lines modified to use only one independent Set
     6043  int allMaximal;
     6044  if (size(#) > 0)
     6045  {
     6046    allMaximal = #[1];
     6047  }
     6048  else
     6049  {
     6050    allMaximal = 1;
     6051  }
     6052
     6053  int nMax;
     6054  if (allMaximal == 1)
     6055  {
     6056    nMax = size(v);
     6057  }
     6058  else
     6059  {
     6060    nMax = 1;
     6061  }
     6062
     6063  for(n = 1; n <= nMax; n++)
     6064  // SL 2006.04.21 2
     6065  {
     6066    di = 0;
     6067    var1 = "";
     6068    var2 = "";
     6069    for(k = 1; k <= size(v[n]); k++)
     6070    {
     6071     if(v[n][k] != 0)
     6072      {
     6073        di++;
     6074        var2 = var2 + "var(" + string(k) + "), ";
     6075      }
     6076      else
     6077      {
     6078        var1 = var1 + "var(" + string(k) + "), ";
     6079      }
     6080    }
     6081    if(di > 0)
     6082    {
     6083      var1 = var1 + var2;
     6084      var1 = var1[1..size(var1) - 2];                         // The "- 2" removes the trailer comma
     6085      hilf[1] = var1;
     6086      // SL 2006.21.04 1 The order is now block dp instead of lp
     6087      hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")";
     6088      // SL 2006.21.04 2
     6089      hilf[3] = di;
     6090      resu[n] = hilf;
     6091    }
     6092    else
     6093    {
     6094      resu[n] = varstr(basering), ordstr(basering), 0;
     6095    }
     6096  }
     6097  return(resu);
     6098}
     6099example
     6100{ "EXAMPLE:"; echo = 2;
     6101   ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp;
     6102   ideal i = ea - fbg, fa + be, ec - fdg, fc + de;
     6103   i = std(i);
     6104   list l = newMaxIndependSetDp(i);
     6105   l;
     6106   i = i, g;
     6107   l = newMaxIndependSetDp(i);
     6108   l;
     6109
     6110   ring s = 0, (x, y, z), lp;
     6111   ideal i = z, yx;
     6112   list l = newMaxIndependSetDp(i);
     6113   l;
     6114}
     6115
     6116
     6117///////////////////////////////////////////////////////////////////////////////
     6118// based on prepareQuotientring
     6119// The order returned is now (C, dp) instead of (C, lp)
     6120
     6121static proc prepareQuotientRingDp (int nnp)
     6122"USAGE:   prepareQuotientRingDp(nnp); nnp int
     6123RETURN:  string = to define Kvar(nnp+1),...,var(nvars)[..rest ]
     6124NOTE:
     6125EXAMPLE: example prepareQuotientRingDp; shows an example
     6126"
     6127{
     6128  ideal @ih,@jh;
     6129  int npar=npars(basering);
     6130  int @n;
     6131
     6132  string quotring= "ring quring = ("+charstr(basering);
     6133  for(@n = nnp + 1; @n <= nvars(basering); @n++)
     6134  {
     6135     quotring = quotring + ", var(" + string(@n) + ")";
     6136     @ih = @ih + var(@n);
     6137  }
     6138
     6139  quotring = quotring+"),(var(1)";
     6140  @jh = @jh + var(1);
     6141  for(@n = 2; @n <= nnp; @n++)
     6142  {
     6143    quotring = quotring + ", var(" + string(@n) + ")";
     6144    @jh = @jh + var(@n);
     6145  }
     6146  // SL 2006-04-21 1 The order returned is now (C, dp) instead of (C, lp)
     6147  quotring = quotring + "), (C, dp);";
     6148  // SL 2006-04-21 2
     6149
     6150  return(quotring);
     6151}
     6152example
     6153{ "EXAMPLE:"; echo = 2;
     6154   ring s1=(0,x),(a,b,c,d,e,f,g),lp;
     6155   def @Q=basering;
     6156   list l= prepareQuotientRingDp(3);
     6157   l;
     6158   execute(l[1]);
     6159   execute(l[2]);
     6160   basering;
     6161   phi;
     6162   setring @Q;
     6163
     6164}
     6165
     6166///////////////////////////////////////////////////////////////////////////////
     6167proc prepareAss(ideal i)
     6168"USAGE:   prepareAss(I); I ideal
     6169RETURN:  list, the radicals of the maximal dimensional components of I.
     6170NOTE:    Uses algorithm of Eisenbud/Huneke/Vasconcelos.
     6171EXAMPLE: example prepareAss; shows an example
     6172"
     6173{
     6174  if(attrib(basering,"global")!=1)
     6175  {
     6176      ERROR(
     6177      "// Not implemented for this ordering, please change to global ordering."
     6178      );
     6179  }
     6180  ideal j=std(i);
     6181  int cod=nvars(basering)-dim(j);
     6182  int e;
     6183  list er;
     6184  ideal ann;
     6185  if(homog(i)==1)
     6186  {
     6187     list re=sres(j,0);                   //the resolution
     6188     re=minres(re);                       //minimized resolution
     6189  }
     6190  else
     6191  {
     6192    list re=mres(i,0);
     6193  }
     6194  for(e=cod;e<=nvars(basering);e++)
     6195  {
     6196     ann=AnnExt_R(e,re);
     6197
     6198     if(nvars(basering)-dim(std(ann))==e)
     6199     {
     6200        er[size(er)+1]=equiRadical(ann);
     6201     }
     6202  }
     6203  return(er);
     6204}
     6205example
     6206{ "EXAMPLE:";  echo = 2;
     6207   ring  r = 0,(x,y,z),dp;
     6208   poly  p = z2+1;
     6209   poly  q = z3+2;
     6210   ideal i = p*q^2,y-z2;
     6211   list pr = prepareAss(i);
     6212   pr;
     6213}
     6214///////////////////////////////////////////////////////////////////////////////
     6215proc equidimMaxEHV(ideal i)
     6216"USAGE:  equidimMaxEHV(I); I ideal
     6217RETURN:  ideal, the equidimensional component (of maximal dimension) of I.
     6218NOTE:    Uses algorithm of Eisenbud, Huneke and Vasconcelos.
     6219EXAMPLE: example equidimMaxEHV; shows an example
     6220"
     6221{
     6222  if(attrib(basering,"global")!=1)
     6223  {
     6224      ERROR(
     6225      "// Not implemented for this ordering, please change to global ordering."
     6226      );
     6227  }
     6228  ideal j=groebner(i);
     6229  int cod=nvars(basering)-dim(j);
     6230  int e;
     6231  ideal ann;
     6232  if(homog(i)==1)
     6233  {
     6234     list re=sres(j,0);                   //the resolution
     6235     re=minres(re);                       //minimized resolution
     6236  }
     6237  else
     6238  {
     6239    list re=mres(i,0);
     6240  }
     6241  ann=AnnExt_R(cod,re);
     6242  return(ann);
     6243}
     6244example
     6245{ "EXAMPLE:";  echo = 2;
     6246   ring  r = 0,(x,y,z),dp;
     6247   ideal i=intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
     6248   equidimMaxEHV(i);
     6249}
     6250
     6251proc testPrimary(list pr, ideal k)
     6252"USAGE:   testPrimary(pr,k); pr a list, k an ideal.
     6253ASSUME:  pr is the result of primdecGTZ(k) or primdecSY(k).
     6254RETURN:  int, 1 if the intersection of the ideals in pr is k, 0 if not
     6255EXAMPLE: example testPrimary; shows an example
     6256"
     6257{
     6258   int i;
     6259   pr=reconvList(pr);
     6260   ideal j=pr[1];
     6261   for (i=2;i<=size(pr)/2;i++)
     6262   {
     6263       j=intersect(j,pr[2*i-1]);
     6264   }
     6265   return(idealsEqual(j,k));
     6266}
     6267example
     6268{ "EXAMPLE:";  echo = 2;
     6269   ring  r = 32003,(x,y,z),dp;
     6270   poly  p = z2+1;
     6271   poly  q = z4+2;
     6272   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     6273   list pr = primdecGTZ(i);
     6274   testPrimary(pr,i);
     6275}
     6276
     6277///////////////////////////////////////////////////////////////////////////////
     6278proc zerodec(ideal I)
     6279"USAGE:   zerodec(I); I ideal
     6280ASSUME:  I is zero-dimensional, the characteristic of the ground field is 0
     6281RETURN:  list of primary ideals, the zero-dimensional decomposition of I
     6282NOTE:    The algorithm (of Monico), works well only for a small total number
     6283         of solutions (@code{vdim(std(I))} should be < 100) and without
     6284         parameters. In practice, it works also in large characteristic p>0
     6285         but may fail for small p.
     6286@*       If printlevel > 0 (default = 0) additional information is displayed.
     6287EXAMPLE: example zerodec; shows an example
     6288"
     6289{
     6290  if(attrib(basering,"global")!=1)
     6291  {
     6292    ERROR(
     6293    "// Not implemented for this ordering, please change to global ordering."
     6294    );
     6295  }
     6296  def R=basering;
     6297  poly q;
     6298  int j,time;
     6299  matrix m;
     6300  list re;
     6301  poly va=var(1);
     6302  ideal J=groebner(I);
     6303  ideal ba=kbase(J);
     6304  int d=vdim(J);
     6305  dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d));
     6306//------ compute matrix of multiplication on R/I with generic element p -----
     6307  int e=nvars(basering);
     6308  poly p=randomLast(100)[e]+random(-50,50);     //the generic element
     6309  matrix n[d][d];
     6310  time = timer;
     6311  for(j=2;j<=e;j++)
     6312  {
     6313    va=va*var(j);
     6314  }
     6315  for(j=1;j<=d;j++)
     6316  {
     6317    q=reduce(p*ba[j],J);
     6318    m=coeffs(q,ba,va);
     6319    n[j,1..d]=m[1..d,1];
     6320  }
     6321  dbprint(printlevel-voice+2,
     6322    "// time for computing multiplication matrix (with generic element) : "+
     6323    string(timer-time));
     6324//---------------- compute characteristic polynomial of matrix --------------
     6325  execute("ring P1=("+charstr(R)+"),T,dp;");
     6326  matrix n=imap(R,n);
     6327  time = timer;
     6328  poly charpol=det(n-T*freemodule(d));
     6329  dbprint(printlevel-voice+2,"// time for computing char poly: "+
     6330         string(timer-time));
     6331//------------------- factorize characteristic polynomial -------------------
     6332//check first if constant term of charpoly is != 0 (which is true for
     6333//sufficiently generic element)
     6334  if(charpol[size(charpol)]!=0)
     6335  {
     6336    time = timer;
     6337    list fac=factor(charpol);
     6338    testFactor(fac,charpol);
     6339    dbprint(printlevel-voice+2,"// time for factorizing char poly: "+
     6340            string(timer-time));
     6341    int f=size(fac[1]);
     6342//--------------------------- the irreducible case --------------------------
     6343    if(f==1)
     6344    {
     6345      setring R;
     6346      re=I;
     6347      return(re);
     6348    }
     6349//---------------------------- the reducible case ---------------------------
     6350//if f_i are the irreducible factors of charpoly, mult=ri, then <I,g_i^ri>
     6351//are the primary components where g_i = f_i(p). However, substituting p in
     6352//f_i may result in a huge object although the final result may be small.
     6353//Hence it is better to simultaneously reduce with I. For this we need a new
     6354//ring.
     6355    execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);");
     6356    list rfac=imap(P1,fac);
     6357    intvec ov=option(get);;
     6358    option(redSB);
     6359    list re1;
     6360    ideal new = T-imap(R,p),imap(R,J);
     6361    attrib(new, "isSB",1);    //we know that new is a standard basis
     6362    for(j=1;j<=f;j++)
     6363    {
     6364      re1[j]=reduce(rfac[1][j]^rfac[2][j],new);
     6365    }
     6366    setring R;
     6367    re = imap(P,re1);
     6368    for(j=1;j<=f;j++)
     6369    {
     6370      J=I,re[j];
     6371      re[j]=interred(J);
     6372    }
     6373    option(set,ov);
     6374    return(re);
     6375  }
     6376  else
     6377//------------------- choice of generic element failed -------------------
     6378  {
     6379    dbprint(printlevel-voice+2,"// try new generic element!");
     6380    setring R;
     6381    return(zerodec(I));
     6382  }
     6383}
     6384example
     6385{ "EXAMPLE:";  echo = 2;
     6386   ring r  = 0,(x,y),dp;
     6387   ideal i = x2-2,y2-2;
     6388   list pr = zerodec(i);
     6389   pr;
     6390}
     6391///////////////////////////////////////////////////////////////////////////////
     6392static proc newDecompStep(ideal i, list #)
     6393"USAGE:  newDecompStep(i); i ideal  (for primary decomposition)
     6394         newDecompStep(i,1);        (for the associated primes of dimension of i)
     6395         newDecompStep(i,2);        (for the minimal associated primes)
     6396         newDecompStep(i,3);        (for the absolute primary decomposition (not tested!))
     6397         "oneIndep";        (for using only one max indep set)
     6398         "intersect";        (returns alse the intersection of the components founded)
     6399
     6400RETURN:  list = list of primary ideals and their associated primes
     6401         (at even positions in the list)
     6402         (resp. a list of the minimal associated primes)
     6403NOTE:    Algorithm of Gianni/Trager/Zacharias
     6404EXAMPLE: example newDecompStep; shows an example
     6405"
     6406{
     6407  intvec op,@vv;
     6408  def  @P = basering;
     6409  list primary,indep,ltras;
     6410  intvec @vh,isat,@w;
     6411  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn;
     6412  ideal peek=i;
     6413  ideal ser,tras;
     6414  list data;
     6415  list result;
     6416  intvec @hilb;
     6417  int isS=(attrib(i,"isSB")==1);
     6418
     6419  // Debug
     6420  dbprint(printlevel - voice, "newDecompStep, v2.0");
     6421
     6422  string indepOption = "allIndep";
     6423  string intersectOption = "noIntersect";
     6424
     6425  if(size(#)>0)
     6426  {
     6427    int count = 1;
     6428    if(typeof(#[count]) == "string")
     6429    {
     6430      if ((#[count] == "oneIndep") or (#[count] == "allIndep"))
     6431      {
     6432        indepOption = #[count];
     6433        count++;
     6434      }
     6435    }
     6436    if(typeof(#[count]) == "string")
     6437    {
     6438      if ((#[count] == "intersect") or (#[count] == "noIntersect"))
     6439      {
     6440        intersectOption = #[count];
     6441        count++;
     6442      }
     6443    }
     6444    if((typeof(#[count]) == "int") or (typeof(#[count]) == "number"))
     6445    {
     6446      if ((#[count]==1)||(#[count]==2)||(#[count]==3))
     6447      {
     6448        @wr=#[count];
     6449        if(@wr==3){abspri = 1; @wr = 0;}
     6450        count++;
     6451      }
     6452    }
     6453    if(size(#)>count)
     6454    {
     6455      seri=1;
     6456      peek=#[count + 1];
     6457      ser=#[count + 2];
     6458    }
     6459  }
     6460  if(abspri)
     6461  {
     6462    list absprimary,abskeep,absprimarytmp,abskeeptmp;
     6463  }
     6464  homo=homog(i);
     6465  if(homo==1)
     6466  {
     6467    if(attrib(i,"isSB")!=1)
     6468    {
     6469      //ltras=mstd(i);
     6470      tras=groebner(i);
     6471      ltras=tras,tras;
     6472      attrib(ltras[1],"isSB",1);
     6473    }
     6474    else
     6475    {
     6476      ltras=i,i;
     6477      attrib(ltras[1],"isSB",1);
     6478    }
     6479    tras = ltras[1];
     6480    attrib(tras,"isSB",1);
     6481    if(dim(tras)==0)
     6482    {
     6483      primary[1]=ltras[2];
     6484      primary[2]=maxideal(1);
     6485      if(@wr>0)
     6486      {
     6487        list l;
     6488        l[2]=maxideal(1);
     6489        l[1]=maxideal(1);
     6490        if (intersectOption == "intersect")
     6491        {
     6492          return(list(l, maxideal(1)));
     6493        }
     6494        else
     6495        {
     6496          return(l);
     6497        }
     6498      }
     6499      if (intersectOption == "intersect")
     6500      {
     6501        return(list(primary, primary[1]));
     6502      }
     6503      else
     6504      {
     6505        return(primary);
     6506      }
     6507    }
     6508    for(@n=1;@n<=nvars(basering);@n++)
     6509    {
     6510      @w[@n]=ord(var(@n));
     6511    }
     6512    @hilb=hilb(tras,1,@w);
     6513    intvec keephilb=@hilb;
     6514  }
     6515
     6516  //----------------------------------------------------------------
     6517  //i is the zero-ideal
     6518  //----------------------------------------------------------------
     6519
     6520  if(size(i)==0)
     6521  {
     6522    primary=i,i;
     6523    if (intersectOption == "intersect")
     6524    {
     6525      return(list(primary, i));
     6526    }
     6527    else
     6528    {
     6529      return(primary);
     6530    }
     6531  }
     6532
     6533  //----------------------------------------------------------------
     6534  //pass to the lexicographical ordering and compute a standardbasis
     6535  //----------------------------------------------------------------
     6536
     6537  int lp=islp();
     6538
     6539  execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);");
     6540  op=option(get);
     6541  option(redSB);
     6542
     6543  ideal ser=fetch(@P,ser);
     6544  if(homo==1)
     6545  {
     6546    if(!lp)
     6547    {
     6548      ideal @j=std(fetch(@P,i),@hilb,@w);
     6549    }
     6550    else
     6551    {
     6552      ideal @j=fetch(@P,tras);
     6553      attrib(@j,"isSB",1);
     6554    }
     6555  }
     6556  else
     6557  {
     6558    if(lp&&isS)
     6559    {
     6560      ideal @j=fetch(@P,i);
     6561      attrib(@j,"isSB",1);
     6562    }
     6563    else
     6564    {
     6565      ideal @j=groebner(fetch(@P,i));
     6566    }
     6567  }
     6568  option(set,op);
     6569  if(seri==1)
     6570  {
     6571    ideal peek=fetch(@P,peek);
     6572    attrib(peek,"isSB",1);
     6573  }
     6574  else
     6575  {
     6576    ideal peek=@j;
     6577  }
     6578  if((size(ser)==0)&&(!abspri))
     6579  {
     6580    ideal fried;
     6581    @n=size(@j);
     6582    for(@k=1;@k<=@n;@k++)
     6583    {
     6584      if(deg(lead(@j[@k]))==1)
     6585      {
     6586        fried[size(fried)+1]=@j[@k];
     6587        @j[@k]=0;
     6588      }
     6589    }
     6590    if(size(fried)==nvars(basering))
     6591    {
     6592      setring @P;
     6593      primary[1]=i;
     6594      primary[2]=i;
     6595      if (intersectOption == "intersect")
     6596      {
     6597        return(list(primary, i));
     6598      }
     6599      else
     6600      {
     6601        return(primary);
     6602      }
     6603    }
     6604    if(size(fried)>0)
     6605    {
     6606      string newva;
     6607      string newma;
     6608      for(@k=1;@k<=nvars(basering);@k++)
     6609      {
     6610        @n1=0;
     6611        for(@n=1;@n<=size(fried);@n++)
     6612        {
     6613          if(leadmonom(fried[@n])==var(@k))
     6614          {
     6615            @n1=1;
     6616            break;
     6617          }
     6618        }
     6619        if(@n1==0)
     6620        {
     6621          newva=newva+string(var(@k))+",";
     6622          newma=newma+string(var(@k))+",";
     6623        }
     6624        else
     6625        {
     6626          newma=newma+string(0)+",";
     6627        }
     6628      }
     6629      newva[size(newva)]=")";
     6630      newma[size(newma)]=";";
     6631      execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;");
     6632      execute("map @kappa=gnir,"+newma);
     6633      ideal @j= @kappa(@j);
     6634      @j=simplify(@j, 2);
     6635      attrib(@j,"isSB",1);
     6636      result = newDecompStep(@j, indepOption, intersectOption, @wr);
     6637      if (intersectOption == "intersect")
     6638      {
     6639       list pr = result[1];
     6640       ideal intersection = result[2];
     6641      }
     6642      else
     6643      {
     6644        list pr = result;
     6645      }
     6646
     6647      setring gnir;
     6648      list pr=imap(@deirf,pr);
     6649      for(@k=1;@k<=size(pr);@k++)
     6650      {
     6651        @j=pr[@k]+fried;
     6652        pr[@k]=@j;
     6653      }
     6654      if (intersectOption == "intersect")
     6655      {
     6656        ideal intersection = imap(@deirf, intersection);
     6657        @j = intersection + fried;
     6658        intersection = @j;
     6659      }
     6660      setring @P;
     6661      if (intersectOption == "intersect")
     6662      {
     6663        return(list(imap(gnir,pr), imap(gnir,intersection)));
     6664      }
     6665      else
     6666      {
     6667        return(imap(gnir,pr));
     6668      }
     6669    }
     6670  }
     6671  //----------------------------------------------------------------
     6672  //j is the ring
     6673  //----------------------------------------------------------------
     6674
     6675  if (dim(@j)==-1)
     6676  {
     6677    setring @P;
     6678    primary=ideal(1),ideal(1);
     6679    if (intersectOption == "intersect")
     6680    {
     6681      return(list(primary, ideal(1)));
     6682    }
     6683    else
     6684    {
     6685      return(primary);
     6686    }
     6687  }
     6688
     6689  //----------------------------------------------------------------
     6690  //  the case of one variable
     6691  //----------------------------------------------------------------
     6692
     6693  if(nvars(basering)==1)
     6694  {
     6695    list fac=factor(@j[1]);
     6696    list gprimary;
     6697    poly generator;
     6698    ideal gIntersection;
     6699    for(@k=1;@k<=size(fac[1]);@k++)
     6700    {
     6701      if(@wr==0)
     6702      {
     6703        gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
     6704        gprimary[2*@k]=ideal(fac[1][@k]);
     6705      }
     6706      else
     6707      {
     6708        gprimary[2*@k-1]=ideal(fac[1][@k]);
     6709        gprimary[2*@k]=ideal(fac[1][@k]);
     6710      }
     6711      if (intersectOption == "intersect")
     6712      {
     6713        generator = generator * fac[1][@k];
     6714      }
     6715    }
     6716    if (intersectOption == "intersect")
     6717    {
     6718      gIntersection = generator;
     6719    }
     6720    setring @P;
     6721    primary=fetch(gnir,gprimary);
     6722    if (intersectOption == "intersect")
     6723    {
     6724      ideal intersection = fetch(gnir,gIntersection);
     6725    }
     6726
     6727//HIER
     6728    if(abspri)
     6729    {
     6730      list resu,tempo;
     6731      string absotto;
     6732      for(ab=1;ab<=size(primary)/2;ab++)
     6733      {
     6734        absotto= absFactorize(primary[2*ab][1],77);
     6735        tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
     6736        resu[ab]=tempo;
     6737      }
     6738      primary=resu;
     6739      intersection = 1;
     6740      for(ab=1;ab<=size(primary);ab++)
     6741      {
     6742        intersection = intersect(intersection, primary[ab][2]);
     6743      }
     6744    }
     6745    if (intersectOption == "intersect")
     6746    {
     6747      return(list(primary, intersection));
     6748    }
     6749    else
     6750    {
     6751      return(primary);
     6752    }
     6753  }
     6754
     6755 //------------------------------------------------------------------
     6756 //the zero-dimensional case
     6757 //------------------------------------------------------------------
     6758  if (dim(@j)==0)
     6759  {
     6760    op=option(get);
     6761    option(redSB);
     6762    list gprimary= newZero_decomp(@j,ser,@wr);
     6763
     6764    setring @P;
     6765    primary=fetch(gnir,gprimary);
     6766
     6767    if(size(ser)>0)
     6768    {
     6769      primary=cleanPrimary(primary);
     6770    }
     6771//HIER
     6772    if(abspri)
     6773    {
     6774      list resu,tempo;
     6775      string absotto;
     6776      for(ab=1;ab<=size(primary)/2;ab++)
     6777      {
     6778        absotto= absFactorize(primary[2*ab][1],77);
     6779        tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
     6780        resu[ab]=tempo;
     6781      }
     6782      primary=resu;
     6783    }
     6784    if (intersectOption == "intersect")
     6785    {
     6786      return(list(primary, fetch(gnir,@j)));
     6787    }
     6788    else
     6789    {
     6790      return(primary);
     6791    }
     6792  }
     6793
     6794  poly @gs,@gh,@p;
     6795  string @va,quotring;
     6796  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
     6797  ideal @h;
     6798  int jdim=dim(@j);
     6799  list fett;
     6800  int lauf,di,newtest;
     6801  //------------------------------------------------------------------
     6802  //search for a maximal independent set indep,i.e.
     6803  //look for subring such that the intersection with the ideal is zero
     6804  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
     6805  //indep[1] is the new varstring and indep[2] the string for block-ordering
     6806  //------------------------------------------------------------------
     6807  if(@wr!=1)
     6808  {
     6809    allindep = newMaxIndependSetLp(@j, indepOption);
     6810    for(@m=1;@m<=size(allindep);@m++)
     6811    {
     6812      if(allindep[@m][3]==jdim)
     6813      {
     6814        di++;
     6815        indep[di]=allindep[@m];
     6816      }
     6817      else
     6818      {
     6819        lauf++;
     6820        restindep[lauf]=allindep[@m];
     6821      }
     6822    }
     6823  }
     6824  else
     6825  {
     6826    indep = newMaxIndependSetLp(@j, indepOption);
     6827  }
     6828
     6829  ideal jkeep=@j;
     6830  if(ordstr(@P)[1]=="w")
     6831  {
     6832    execute("ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");");
     6833  }
     6834  else
     6835  {
     6836    execute( "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);");
     6837  }
     6838
     6839  if(homo==1)
     6840  {
     6841    if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
     6842       ||(ordstr(@P)[3]=="w"))
     6843    {
     6844      ideal jwork=imap(@P,tras);
     6845      attrib(jwork,"isSB",1);
     6846    }
     6847    else
     6848    {
     6849      ideal jwork=std(imap(gnir,@j),@hilb,@w);
     6850    }
     6851  }
     6852  else
     6853  {
     6854    ideal jwork=groebner(imap(gnir,@j));
     6855  }
     6856  list hquprimary;
     6857  poly @p,@q;
     6858  ideal @h,fac,ser;
     6859//Aenderung================
     6860  ideal @Ptest=1;
     6861//=========================
     6862  di=dim(jwork);
     6863  keepdi=di;
     6864
     6865  ser = 1;
     6866
     6867  setring gnir;
     6868  for(@m=1; @m<=size(indep); @m++)
     6869  {
     6870    data[1] = indep[@m];
     6871    result = newReduction(@j, ser, @hilb, @w, jdim, abspri, @wr, data);
     6872    quprimary = quprimary + result[1];
     6873    if(abspri)
     6874    {
     6875      absprimary = absprimary + result[2];
     6876      abskeep = abskeep + result[3];
     6877    }
     6878    @h = result[5];
     6879    ser = result[4];
     6880    if(size(@h)>0)
     6881    {
     6882      //---------------------------------------------------------------
     6883      //we change to @Phelp to have the ordering dp for saturation
     6884      //---------------------------------------------------------------
     6885
     6886      setring @Phelp;
     6887      @h=imap(gnir,@h);
     6888//Aenderung==================================
     6889      if(defined(@LL)){kill @LL;}
     6890      list @LL=minSat(jwork,@h);
     6891      @Ptest=intersect(@Ptest,@LL[1]);
     6892      ser = intersect(ser, @LL[1]);
     6893//===========================================
     6894
     6895      if(@wr!=1)
     6896      {
     6897//Aenderung==================================
     6898        @q=@LL[2];
     6899//===========================================
     6900        //@q=minSat(jwork,@h)[2];
     6901      }
     6902      else
     6903      {
     6904        fac=ideal(0);
     6905        for(lauf=1;lauf<=ncols(@h);lauf++)
     6906        {
     6907          if(deg(@h[lauf])>0)
     6908          {
     6909            fac=fac+factorize(@h[lauf],1);
     6910          }
     6911        }
     6912        fac=simplify(fac,6);
     6913        @q=1;
     6914        for(lauf=1;lauf<=size(fac);lauf++)
     6915        {
     6916          @q=@q*fac[lauf];
     6917        }
     6918      }
     6919      jwork = std(jwork,@q);
     6920      keepdi = dim(jwork);
     6921      if(keepdi < di)
     6922      {
     6923        setring gnir;
     6924        @j = imap(@Phelp, jwork);
     6925        ser = imap(@Phelp, ser);
     6926        break;
     6927      }
     6928      if(homo == 1)
     6929      {
     6930        @hilb = hilb(jwork, 1, @w);
     6931      }
     6932
     6933      setring gnir;
     6934      ser = imap(@Phelp, ser);
     6935      @j = imap(@Phelp, jwork);
     6936    }
     6937  }
     6938
     6939  if((size(quprimary)==0)&&(@wr==1))
     6940  {
     6941     @j=ideal(1);
     6942     quprimary[1]=ideal(1);
     6943     quprimary[2]=ideal(1);
     6944  }
     6945  if((size(quprimary)==0))
     6946  {
     6947    keepdi = di - 1;
     6948    quprimary[1]=ideal(1);
     6949    quprimary[2]=ideal(1);
     6950  }
     6951  //---------------------------------------------------------------
     6952  //notice that j=sat(j,gh) intersected with (j,gh^n)
     6953  //we finished with sat(j,gh) and have to start with (j,gh^n)
     6954  //---------------------------------------------------------------
     6955  if((deg(@j[1])!=0)&&(@wr!=1))
     6956  {
     6957     if(size(quprimary)>0)
     6958     {
     6959        setring @Phelp;
     6960        ser=imap(gnir,ser);
     6961
     6962        hquprimary=imap(gnir,quprimary);
     6963        if(@wr==0)
     6964        {
     6965//Aenderung====================================================
     6966//HIER STATT DURCHSCHNITT SATURIEREN!
     6967           ideal htest=@Ptest;
     6968/*
     6969           ideal htest=hquprimary[1];
     6970           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
     6971           {
     6972              htest=intersect(htest,hquprimary[2*@n1-1]);
     6973           }
     6974*/
     6975//=============================================================
     6976        }
     6977        else
     6978        {
     6979           ideal htest=hquprimary[2];
     6980
     6981           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
     6982           {
     6983              htest=intersect(htest,hquprimary[2*@n1]);
     6984           }
     6985        }
     6986
     6987        if(size(ser)>0)
     6988        {
     6989           ser=intersect(htest,ser);
     6990        }
     6991        else
     6992        {
     6993          ser=htest;
     6994        }
     6995        setring gnir;
     6996        ser=imap(@Phelp,ser);
     6997     }
     6998     if(size(reduce(ser,peek,1))!=0)
     6999     {
     7000        for(@m=1;@m<=size(restindep);@m++)
     7001        {
     7002         // if(restindep[@m][3]>=keepdi)
     7003         // {
     7004           isat=0;
     7005           @n2=0;
     7006
     7007           if(restindep[@m][1]==varstr(basering))
     7008           //the good case, nothing to do, just to have the same notations
     7009           //change the ring
     7010           {
     7011              execute("ring gnir1 = ("+charstr(basering)+"),("+
     7012                varstr(basering)+"),("+ordstr(basering)+");");
     7013              ideal @j=fetch(gnir,jkeep);
     7014              attrib(@j,"isSB",1);
     7015           }
     7016           else
     7017           {
     7018              @va=string(maxideal(1));
     7019              execute("ring gnir1 = ("+charstr(basering)+"),("+
     7020                      restindep[@m][1]+"),(" +restindep[@m][2]+");");
     7021              execute("map phi=gnir,"+@va+";");
     7022              op=option(get);
     7023              option(redSB);
     7024              if(homo==1)
     7025              {
     7026                 ideal @j=std(phi(jkeep),keephilb,@w);
     7027              }
     7028              else
     7029              {
     7030                ideal @j=groebner(phi(jkeep));
     7031              }
     7032              ideal ser=phi(ser);
     7033              option(set,op);
     7034           }
     7035
     7036           for (lauf=1;lauf<=size(@j);lauf++)
     7037           {
     7038              fett[lauf]=size(@j[lauf]);
     7039           }
     7040           //------------------------------------------------------------------
     7041           //we have now the following situation:
     7042           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may
     7043           //pass to this quotientring, j is their still a standardbasis, the
     7044           //leading coefficients of the polynomials  there (polynomials in
     7045           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     7046           //we need their ggt, gh, because of the following:
     7047           //let (j:gh^n)=(j:gh^infinity) then
     7048           //j*K(var(nnp+1),..,var(nva))[..the rest..]
     7049           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     7050           //on the other hand j=(j,gh^n) intersected with (j:gh^n)
     7051
     7052           //------------------------------------------------------------------
     7053
     7054           //the arrangement for the quotientring
     7055           // K(var(nnp+1),..,var(nva))[..the rest..]
     7056           //and the map phi:K[var(1),...,var(nva)] ---->
     7057           //--->K(var(nnpr+1),..,var(nva))[..the rest..]
     7058           //------------------------------------------------------------------
     7059
     7060           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
     7061
     7062           //------------------------------------------------------------------
     7063           //we pass to the quotientring  K(var(nnp+1),..,var(nva))[..rest..]
     7064           //------------------------------------------------------------------
     7065
     7066           execute(quotring);
     7067
     7068           // @j considered in the quotientring
     7069           ideal @j=imap(gnir1,@j);
     7070           ideal ser=imap(gnir1,ser);
     7071
     7072           kill gnir1;
     7073
     7074           //j is a standardbasis in the quotientring but usually not minimal
     7075           //here it becomes minimal
     7076           @j=clearSB(@j,fett);
     7077           attrib(@j,"isSB",1);
     7078
     7079           //we need later ggt(h[1],...)=gh for saturation
     7080           ideal @h;
     7081
     7082           for(@n=1;@n<=size(@j);@n++)
     7083           {
     7084              @h[@n]=leadcoef(@j[@n]);
     7085           }
     7086           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
     7087
     7088           op=option(get);
     7089           option(redSB);
     7090           list uprimary= newZero_decomp(@j,ser,@wr);
     7091//HIER
     7092           if(abspri)
     7093           {
     7094              ideal II;
     7095              ideal jmap;
     7096              map sigma;
     7097              nn=nvars(basering);
     7098              map invsigma=basering,maxideal(1);
     7099              for(ab=1;ab<=size(uprimary)/2;ab++)
     7100              {
     7101                 II=uprimary[2*ab];
     7102                 attrib(II,"isSB",1);
     7103                 if(deg(II[1])!=vdim(II))
     7104                 {
     7105                    jmap=randomLast(50);
     7106                    sigma=basering,jmap;
     7107                    jmap[nn]=2*var(nn)-jmap[nn];
     7108                    invsigma=basering,jmap;
     7109                    II=groebner(sigma(II));
     7110                  }
     7111                  absprimarytmp[ab]= absFactorize(II[1],77);
     7112                  II=var(nn);
     7113                  abskeeptmp[ab]=string(invsigma(II));
     7114                  invsigma=basering,maxideal(1);
     7115              }
     7116           }
     7117           option(set,op);
     7118
     7119           //we need the intersection of the ideals in the list quprimary with
     7120           //the polynomialring, i.e. let q=(f1,...,fr) in the quotientring
     7121           //such an ideal but fi polynomials, then the intersection of q with
     7122           //the polynomialring is the saturation of the ideal generated by
     7123           //f1,...,fr with respect toh which is the lcm of the leading
     7124           //coefficients of the fi considered in the quotientring:
     7125           //this is coded in saturn
     7126
     7127           list saturn;
     7128           ideal hpl;
     7129
     7130           for(@n=1;@n<=size(uprimary);@n++)
     7131           {
     7132              hpl=0;
     7133              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
     7134              {
     7135                 hpl=hpl,leadcoef(uprimary[@n][@n1]);
     7136              }
     7137              saturn[@n]=hpl;
     7138           }
     7139           //------------------------------------------------------------------
     7140           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..rest..]
     7141           //back to the polynomialring
     7142           //------------------------------------------------------------------
     7143           setring gnir;
     7144           collectprimary=imap(quring,uprimary);
     7145           lsau=imap(quring,saturn);
     7146           @h=imap(quring,@h);
     7147
     7148           kill quring;
     7149
     7150
     7151           @n2=size(quprimary);
     7152//================NEU=========================================
     7153           if(deg(quprimary[1][1])<=0){ @n2=0; }
     7154//============================================================
     7155
     7156           @n3=@n2;
     7157
     7158           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
     7159           {
     7160              if(deg(collectprimary[2*@n1][1])>0)
     7161              {
     7162                 @n2++;
     7163                 quprimary[@n2]=collectprimary[2*@n1-1];
     7164                 lnew[@n2]=lsau[2*@n1-1];
     7165                 @n2++;
     7166                 lnew[@n2]=lsau[2*@n1];
     7167                 quprimary[@n2]=collectprimary[2*@n1];
     7168                 if(abspri)
     7169                 {
     7170                   absprimary[@n2/2]=absprimarytmp[@n1];
     7171                   abskeep[@n2/2]=abskeeptmp[@n1];
     7172                 }
     7173              }
     7174           }
     7175
     7176
     7177           //here the intersection with the polynomialring
     7178           //mentioned above is really computed
     7179
     7180           for(@n=@n3/2+1;@n<=@n2/2;@n++)
     7181           {
     7182              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
     7183              {
     7184                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7185                 quprimary[2*@n]=quprimary[2*@n-1];
     7186              }
     7187              else
     7188              {
     7189                 if(@wr==0)
     7190                 {
     7191                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7192                 }
     7193                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
     7194              }
     7195           }
     7196           if(@n2>=@n3+2)
     7197           {
     7198              setring @Phelp;
     7199              ser=imap(gnir,ser);
     7200              hquprimary=imap(gnir,quprimary);
     7201              for(@n=@n3/2+1;@n<=@n2/2;@n++)
     7202              {
     7203                if(@wr==0)
     7204                {
     7205                   ser=intersect(ser,hquprimary[2*@n-1]);
     7206                }
     7207                else
     7208                {
     7209                   ser=intersect(ser,hquprimary[2*@n]);
     7210                }
     7211              }
     7212              setring gnir;
     7213              ser=imap(@Phelp,ser);
     7214           }
     7215
     7216         // }
     7217        }
     7218//HIER
     7219        if(abspri)
     7220        {
     7221          list resu,tempo;
     7222          for(ab=1;ab<=size(quprimary)/2;ab++)
     7223          {
     7224             if (deg(quprimary[2*ab][1])!=0)
     7225             {
     7226               tempo=quprimary[2*ab-1],quprimary[2*ab],
     7227                         absprimary[ab],abskeep[ab];
     7228               resu[ab]=tempo;
     7229             }
     7230          }
     7231          quprimary=resu;
     7232          @wr=3;
     7233        }
     7234        if(size(reduce(ser,peek,1))!=0)
     7235        {
     7236           if(@wr>0)
     7237           {
     7238              // The following line was dropped to avoid the recursion step:
     7239              //htprimary=newDecompStep(@j,@wr,peek,ser);
     7240              htprimary = list();
     7241           }
     7242           else
     7243           {
     7244              // The following line was dropped to avoid the recursion step:
     7245              //htprimary=newDecompStep(@j,peek,ser);
     7246              htprimary = list();
     7247           }
     7248           // here we collect now both results primary(sat(j,gh))
     7249           // and primary(j,gh^n)
     7250           @n=size(quprimary);
     7251           if (deg(quprimary[1][1])<=0) { @n=0; }
     7252           for (@k=1;@k<=size(htprimary);@k++)
     7253           {
     7254              quprimary[@n+@k]=htprimary[@k];
     7255           }
     7256        }
     7257     }
     7258   }
     7259   else
     7260   {
     7261      if(abspri)
     7262      {
     7263        list resu,tempo;
     7264        for(ab=1;ab<=size(quprimary)/2;ab++)
     7265        {
     7266           tempo=quprimary[2*ab-1],quprimary[2*ab],
     7267                   absprimary[ab],abskeep[ab];
     7268           resu[ab]=tempo;
     7269        }
     7270        quprimary=resu;
     7271      }
     7272   }
     7273  //---------------------------------------------------------------------------
     7274  //back to the ring we started with
     7275  //the final result: primary
     7276  //---------------------------------------------------------------------------
     7277
     7278  setring @P;
     7279  primary=imap(gnir,quprimary);
     7280
     7281  if (intersectOption == "intersect")
     7282  {
     7283     return(list(primary, imap(gnir, ser)));
     7284  }
     7285  else
     7286  {
     7287    return(primary);
     7288  }
     7289}
     7290example
     7291{ "EXAMPLE:"; echo = 2;
     7292   ring  r = 32003,(x,y,z),lp;
     7293   poly  p = z2+1;
     7294   poly  q = z4+2;
     7295   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     7296   list pr= newDecompStep(i);
     7297   pr;
     7298   testPrimary( pr, i);
     7299}
     7300
     7301// This was part of proc decomp.
     7302// In proc newDecompStep, used for the computation of the minimal associated primes,
     7303// this part was separated as a soubrutine to make the code more clear.
     7304// Also, since the reduction is performed twice in proc newDecompStep, it should use both times this routine.
     7305// This is not yet implemented, since the reduction is not exactly the same and some changes should be made.
     7306static proc newReduction(ideal @j, ideal ser, intvec @hilb, intvec @w, int jdim, int abspri, int @wr, list data)
     7307{
     7308   string @va;
     7309   string quotring;
     7310   intvec op;
     7311   intvec @vv;
     7312   def gnir = basering;
     7313   ideal isat=0;
     7314   int @n;
     7315   int @n1 = 0;
     7316   int @n2 = 0;
     7317   int @n3 = 0;
     7318   int homo = homog(@j);
     7319   int lauf;
     7320   int @k;
     7321   list fett;
     7322   int keepdi;
     7323   list collectprimary;
     7324   list lsau;
     7325   list lnew;
     7326   ideal @h;
     7327
     7328   list indepInfo = data[1];
     7329   list quprimary = list();
     7330
     7331   //if(abspri)
     7332   //{
     7333     int ab;
     7334     list absprimarytmp,abskeeptmp;
     7335     list absprimary, abskeep;
     7336   //}
     7337   // Debug
     7338   dbprint(printlevel - voice, "newReduction, v2.0");
     7339
     7340   if((indepInfo[1]==varstr(basering)))  // &&(@m==1)
     7341   //this is the good case, nothing to do, just to have the same notations
     7342   //change the ring
     7343   {
     7344     execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
     7345                              +ordstr(basering)+");");
     7346     ideal @j = fetch(gnir, @j);
     7347     attrib(@j,"isSB",1);
     7348     ideal ser = fetch(gnir, ser);
     7349   }
     7350   else
     7351   {
     7352     @va=string(maxideal(1));
     7353//Aenderung==============
     7354     //if(@m==1)
     7355     //{
     7356     //  @j=fetch(@P,i);
     7357     //}
     7358//=======================
     7359     execute("ring gnir1 = ("+charstr(basering)+"),("+indepInfo[1]+"),("
     7360                              +indepInfo[2]+");");
     7361     execute("map phi=gnir,"+@va+";");
     7362     op=option(get);
     7363     option(redSB);
     7364     if(homo==1)
     7365     {
     7366       ideal @j=std(phi(@j),@hilb,@w);
     7367     }
     7368     else
     7369     {
     7370       ideal @j=groebner(phi(@j));
     7371     }
     7372     ideal ser=phi(ser);
     7373
     7374     option(set,op);
     7375   }
     7376   if((deg(@j[1])==0)||(dim(@j)<jdim))
     7377   {
     7378     setring gnir;
     7379     break;
     7380   }
     7381   for (lauf=1;lauf<=size(@j);lauf++)
     7382   {
     7383     fett[lauf]=size(@j[lauf]);
     7384   }
     7385   //------------------------------------------------------------------------
     7386   //we have now the following situation:
     7387   //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     7388   //to this quotientring, j is their still a standardbasis, the
     7389   //leading coefficients of the polynomials  there (polynomials in
     7390   //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     7391   //we need their ggt, gh, because of the following: let
     7392   //(j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
     7393   //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     7394   //on the other hand j=(j,gh^n) intersected with (j:gh^n)
     7395
     7396   //------------------------------------------------------------------------
     7397
     7398   //arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and
     7399   //map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..rest..]
     7400   //------------------------------------------------------------------------
     7401
     7402   quotring=prepareQuotientring(nvars(basering)-indepInfo[3]);
     7403
     7404   //---------------------------------------------------------------------
     7405   //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     7406   //---------------------------------------------------------------------
     7407
     7408   ideal @jj=lead(@j);               //!! vorn vereinbaren
     7409   execute(quotring);
     7410
     7411   ideal @jj=imap(gnir1,@jj);
     7412   @vv=clearSBNeu(@jj,fett);  //!! vorn vereinbaren
     7413   setring gnir1;
     7414   @k=size(@j);
     7415   for (lauf=1;lauf<=@k;lauf++)
     7416   {
     7417     if(@vv[lauf]==1)
     7418     {
     7419       @j[lauf]=0;
     7420     }
     7421   }
     7422   @j=simplify(@j,2);
     7423   setring quring;
     7424   // @j considered in the quotientring
     7425   ideal @j=imap(gnir1,@j);
     7426
     7427   ideal ser=imap(gnir1,ser);
     7428
     7429   kill gnir1;
     7430
     7431   //j is a standardbasis in the quotientring but usually not minimal
     7432   //here it becomes minimal
     7433
     7434   attrib(@j,"isSB",1);
     7435
     7436   //we need later ggt(h[1],...)=gh for saturation
     7437   ideal @h;
     7438   if(deg(@j[1])>0)
     7439   {
     7440     for(@n=1;@n<=size(@j);@n++)
     7441     {
     7442       @h[@n]=leadcoef(@j[@n]);
     7443     }
     7444     //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
     7445     op=option(get);
     7446     option(redSB);
     7447
     7448     int zeroMinAss = @wr;
     7449     if (@wr == 2) {zeroMinAss = 1;}
     7450     list uprimary= newZero_decomp(@j, ser, zeroMinAss);
     7451
     7452//HIER
     7453     if(abspri)
     7454     {
     7455       ideal II;
     7456       ideal jmap;
     7457       map sigma;
     7458       nn=nvars(basering);
     7459       map invsigma=basering,maxideal(1);
     7460       for(ab=1;ab<=size(uprimary)/2;ab++)
     7461       {
     7462         II=uprimary[2*ab];
     7463         attrib(II,"isSB",1);
     7464         if(deg(II[1])!=vdim(II))
     7465         {
     7466           jmap=randomLast(50);
     7467           sigma=basering,jmap;
     7468           jmap[nn]=2*var(nn)-jmap[nn];
     7469           invsigma=basering,jmap;
     7470           II=groebner(sigma(II));
     7471         }
     7472         absprimarytmp[ab]= absFactorize(II[1],77);
     7473         II=var(nn);
     7474         abskeeptmp[ab]=string(invsigma(II));
     7475         invsigma=basering,maxideal(1);
     7476       }
     7477     }
     7478     option(set,op);
     7479   }
     7480   else
     7481   {
     7482     list uprimary;
     7483     uprimary[1]=ideal(1);
     7484     uprimary[2]=ideal(1);
     7485   }
     7486   //we need the intersection of the ideals in the list quprimary with the
     7487   //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     7488   //but fi polynomials, then the intersection of q with the polynomialring
     7489   //is the saturation of the ideal generated by f1,...,fr with respect to
     7490   //h which is the lcm of the leading coefficients of the fi considered in
     7491   //in the quotientring: this is coded in saturn
     7492
     7493   list saturn;
     7494   ideal hpl;
     7495
     7496   for(@n=1;@n<=size(uprimary);@n++)
     7497   {
     7498     uprimary[@n]=interred(uprimary[@n]); // temporary fix
     7499     hpl=0;
     7500     for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
     7501     {
     7502       hpl=hpl,leadcoef(uprimary[@n][@n1]);
     7503     }
     7504     saturn[@n]=hpl;
     7505   }
     7506
     7507   //--------------------------------------------------------------------
     7508   //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     7509   //back to the polynomialring
     7510   //---------------------------------------------------------------------
     7511   setring gnir;
     7512
     7513   collectprimary=imap(quring,uprimary);
     7514   lsau=imap(quring,saturn);
     7515   @h=imap(quring,@h);
     7516
     7517   kill quring;
     7518
     7519   @n2=size(quprimary);
     7520   @n3=@n2;
     7521
     7522   for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
     7523   {
     7524     if(deg(collectprimary[2*@n1][1])>0)
     7525     {
     7526       @n2++;
     7527       quprimary[@n2]=collectprimary[2*@n1-1];
     7528       lnew[@n2]=lsau[2*@n1-1];
     7529       @n2++;
     7530       lnew[@n2]=lsau[2*@n1];
     7531       quprimary[@n2]=collectprimary[2*@n1];
     7532       if(abspri)
     7533       {
     7534         absprimary[@n2/2]=absprimarytmp[@n1];
     7535         abskeep[@n2/2]=abskeeptmp[@n1];
     7536       }
     7537     }
     7538   }
     7539
     7540   //here the intersection with the polynomialring
     7541   //mentioned above is really computed
     7542   for(@n=@n3/2+1;@n<=@n2/2;@n++)
     7543   {
     7544     if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
     7545     {
     7546       quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7547       quprimary[2*@n]=quprimary[2*@n-1];
     7548     }
     7549     else
     7550     {
     7551       if(@wr==0)
     7552       {
     7553         quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7554       }
     7555       quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
     7556     }
     7557   }
     7558
     7559   return(quprimary, absprimary, abskeep, ser, @h);
     7560}
     7561
     7562
     7563////////////////////////////////////////////////////////////////////////////
     7564
     7565
     7566
     7567
     7568///////////////////////////////////////////////////////////////////////////////
     7569// Based on minAssGTZ
     7570
     7571proc minAss(ideal i,list #)
     7572"USAGE:   minAss(I[, l]); i ideal, l list (optional) of parameters, same as minAssGTZ
     7573RETURN:  a list, the minimal associated prime ideals of I.
     7574NOTE:    Designed for characteristic 0, works also in char k > 0 based
     7575         on an algorithm of Yokoyama
     7576EXAMPLE: example minAss; shows an example
     7577"
     7578{
     7579  return(minAssGTZ(i,#));
     7580}
     7581example
     7582{ "EXAMPLE:";  echo = 2;
     7583   ring  r = 0, (x, y, z), dp;
     7584   poly  p = z2 + 1;
     7585   poly  q = z3 + 2;
     7586   ideal i = p * q^2, y - z2;
     7587   list pr = minAss(i);
     7588   pr;
     7589}
     7590
     7591
     7592///////////////////////////////////////////////////////////////////////////////
     7593//
     7594// Computes the minimal associated primes of I via Laplagne algorithm,
     7595// using primary decomposition in the zero dimensional case.
     7596// For reduction to the zerodimensional case, it uses the procedure
     7597// decomp, with some modifications to avoid the recursion.
     7598//
     7599
     7600static proc minAssSL(ideal I)
     7601// Input = I, ideal
     7602// Output = primaryDec where primaryDec is the list of the minimal
     7603// associated primes and the primary components corresponding to these primes.
     7604{
     7605  ideal P = 1;
     7606  list pd = list();
     7607  int k;
     7608  int stop = 0;
     7609  list primaryDec = list();
     7610
     7611  while (stop == 0)
     7612  {
     7613    // Debug
     7614    dbprint(printlevel - voice, "// We call minAssSLIteration to find new prime ideals!");
     7615    pd = minAssSLIteration(I, P);
     7616    // Debug
     7617    dbprint(printlevel - voice, "// Output of minAssSLIteration:");
     7618    dbprint(printlevel - voice, pd);
     7619    if (size(pd[1]) > 0)
     7620    {
     7621      primaryDec = primaryDec + pd[1];
     7622      // Debug
     7623      dbprint(printlevel - voice, "// We intersect the prime ideals obtained.");
     7624      P = intersect(P, pd[2]);
     7625      // Debug
     7626      dbprint(printlevel - voice, "// Intersection finished.");
     7627    }
     7628    else
     7629    {
     7630      stop = 1;
     7631    }
     7632  }
     7633
     7634  // Returns only the primary components, not the radical.
     7635  return(primaryDec);
     7636}
     7637
     7638///////////////////////////////////////////////////////////////////////////////
     7639// Given an ideal I and an ideal P (intersection of some minimal prime ideals
     7640// associated to I), it calculates new minimal prime ideals associated to I
     7641// which were not used to calculate P.
     7642// This version uses Primary Decomposition in the zerodimensional case.
     7643static proc minAssSLIteration(ideal I, ideal P);
     7644{
     7645  int k = 1;
     7646  int good  = 0;
     7647  list primaryDec = list();
     7648  // Debug
     7649  dbprint (printlevel-voice, "// We search for an element in P - sqrt(I).");
     7650  while ((k <= size(P)) and (good == 0))
     7651  {
     7652    good = 1 - rad_con(P[k], I);
     7653    k++;
     7654  }
     7655  k--;
     7656  if (good == 0)
     7657  {
     7658    // Debug
     7659    dbprint (printlevel - voice, "// No element was found, P = sqrt(I).");
     7660    return (list(primaryDec, ideal(0)));
     7661  }
     7662  // Debug
     7663  dbprint (printlevel - voice, "// We found h = ", P[k]);
     7664  dbprint (printlevel - voice, "// We calculate the saturation of I with respect to the element just founded.");
     7665  ideal J = sat(I, P[k])[1];
     7666
     7667  // Uses decomp from primdec, modified to avoid the recursion.
     7668  // Debug
     7669  dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via decomp.");
     7670
     7671  primaryDec = newDecompStep(J, "oneIndep", "intersect", 2);
     7672  // Debug
     7673  dbprint(printlevel - voice, "// Proc decomp has found", size(primaryDec) / 2, "new primary components.");
     7674
     7675  return(primaryDec);
     7676}
     7677
     7678
     7679
     7680///////////////////////////////////////////////////////////////////////////////////
     7681// Based on maxIndependSet
     7682// Added list # as parameter
     7683// If the first element of # is 0, the output is only 1 max indep set.
     7684// If no list is specified or #[1] = 1, the output is all the max indep set of the
     7685// leading terms ideal. This is the original output of maxIndependSet
     7686
     7687proc newMaxIndependSetLp(ideal j, list #)
     7688"USAGE:   newMaxIndependentSetLp(i); i ideal (returns all maximal independent sets of the corresponding leading terms ideal)
     7689          newMaxIndependentSetLp(i, 0); i ideal (returns only one maximal independent set)
     7690RETURN:  list = #1. new varstring with the maximal independent set at the end,
     7691                #2. ordstring with the lp ordering,
     7692                #3. the number of independent variables
     7693NOTE:
     7694EXAMPLE: example newMaxIndependentSetLp; shows an example
     7695"
     7696{
     7697  int n, k, di;
     7698  list resu, hilf;
     7699  string var1, var2;
     7700  list v = indepSet(j, 0);
     7701
     7702  // SL 2006.04.21 1 Lines modified to use only one independent Set
     7703  string indepOption;
     7704  if (size(#) > 0)
     7705  {
     7706    indepOption = #[1];
     7707  }
     7708  else
     7709  {
     7710    indepOption = "allIndep";
     7711  }
     7712
     7713  int nMax;
     7714  if (indepOption == "allIndep")
     7715  {
     7716    nMax = size(v);
     7717  }
     7718  else
     7719  {
     7720    nMax = 1;
     7721  }
     7722
     7723  for(n = 1; n <= nMax; n++)
     7724  // SL 2006.04.21 2
     7725  {
     7726    di = 0;
     7727    var1 = "";
     7728    var2 = "";
     7729    for(k = 1; k <= size(v[n]); k++)
     7730    {
     7731      if(v[n][k] != 0)
     7732      {
     7733        di++;
     7734        var2 = var2 + "var(" + string(k) + "), ";
     7735      }
     7736      else
     7737      {
     7738        var1 = var1 + "var(" + string(k) + "), ";
     7739      }
     7740    }
     7741    if(di > 0)
     7742    {
     7743      var1 = var1 + var2;
     7744      var1 = var1[1..size(var1) - 2];       // The "- 2" removes the trailer comma
     7745      hilf[1] = var1;
     7746      // SL 2006.21.04 1 The order is now block dp instead of lp
     7747      //hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")";
     7748      // SL 2006.21.04 2
     7749      // For decomp, lp ordering is needed. Nothing is changed.
     7750      hilf[2] = "lp";
     7751      hilf[3] = di;
     7752      resu[n] = hilf;
     7753    }
     7754    else
     7755    {
     7756      resu[n] = varstr(basering), ordstr(basering), 0;
     7757    }
     7758  }
     7759  return(resu);
     7760}
     7761example
     7762{ "EXAMPLE:"; echo = 2;
     7763   ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp;
     7764   ideal i = ea - fbg, fa + be, ec - fdg, fc + de;
     7765   i = std(i);
     7766   list l = newMaxIndependSetLp(i);
     7767   l;
     7768   i = i, g;
     7769   l = newMaxIndependSetLp(i);
     7770   l;
     7771
     7772   ring s = 0, (x, y, z), lp;
     7773   ideal i = z, yx;
     7774   list l = newMaxIndependSetLp(i);
     7775   l;
     7776}
     7777
     7778
     7779///////////////////////////////////////////////////////////////////////////////
     7780
     7781proc newZero_decomp (ideal j, ideal ser, int @wr, list #)
     7782"USAGE:   newZero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
     7783         (@wr=0 for primary decomposition, @wr=1 for computation of associated
     7784         primes)
     7785         if #[1] = "nest", then #[2] indicates the nest level (number of recursive calls)
     7786         When the nest level is high it indicates that the computation is difficult,
     7787         and different methods are applied.
     7788RETURN:  list = list of primary ideals and their radicals (at even positions
     7789         in the list) if the input is zero-dimensional and a standardbases
     7790         with respect to lex-ordering
     7791         If ser!=(0) and ser is contained in j or if j is not zero-dimen-
     7792         sional then ideal(1),ideal(1) is returned
     7793NOTE:    Algorithm of Gianni/Trager/Zacharias
     7794EXAMPLE: example newZero_decomp; shows an example
     7795"
     7796{
     7797  def   @P = basering;
     7798  int uytrewq;
     7799  int nva = nvars(basering);
     7800  int @k,@s,@n,@k1,zz;
     7801  list primary,lres0,lres1,act,@lh,@wh;
     7802  map phi,psi,phi1,psi1;
     7803  ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1;
     7804  intvec @vh,@hilb;
     7805  string @ri;
     7806  poly @f;
     7807
     7808  // Debug
     7809  dbprint(printlevel - voice, "proc newZero_decomp");
     7810
     7811  if (dim(j)>0)
     7812  {
     7813    primary[1]=ideal(1);
     7814    primary[2]=ideal(1);
     7815    return(primary);
     7816  }
     7817  j=interred(j);
     7818
     7819  attrib(j,"isSB",1);
     7820
     7821  int nestLevel = 0;
     7822  if (size(#) > 0)
     7823  {
     7824    if (typeof(#[1]) == "string")
     7825    {
     7826      if (#[1] == "nest")
     7827      {
     7828        nestLevel = #[2];
     7829      }
     7830      # = list();
     7831    }
     7832  }
     7833
     7834  if(vdim(j)==deg(j[1]))
     7835  {
     7836    act=factor(j[1]);
     7837    for(@k=1;@k<=size(act[1]);@k++)
     7838    {
     7839      @qh=j;
     7840      if(@wr==0)
     7841      {
     7842        @qh[1]=act[1][@k]^act[2][@k];
     7843      }
     7844      else
     7845      {
     7846        @qh[1]=act[1][@k];
     7847      }
     7848      primary[2*@k-1]=interred(@qh);
     7849      @qh=j;
     7850      @qh[1]=act[1][@k];
     7851      primary[2*@k]=interred(@qh);
     7852      attrib( primary[2*@k-1],"isSB",1);
     7853
     7854      if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
     7855      {
     7856        primary[2*@k-1]=ideal(1);
     7857        primary[2*@k]=ideal(1);
     7858      }
     7859    }
     7860    return(primary);
     7861  }
     7862
     7863  if(homog(j)==1)
     7864  {
     7865    primary[1]=j;
     7866    if((size(ser)>0)&&(size(reduce(ser,j,1))==0))
     7867    {
     7868      primary[1]=ideal(1);
     7869      primary[2]=ideal(1);
     7870      return(primary);
     7871    }
     7872    if(dim(j)==-1)
     7873    {
     7874      primary[1]=ideal(1);
     7875      primary[2]=ideal(1);
     7876    }
     7877    else
     7878    {
     7879      primary[2]=maxideal(1);
     7880    }
     7881    return(primary);
     7882  }
     7883
     7884//the first element in the standardbase is factorized
     7885  if(deg(j[1])>0)
     7886  {
     7887    act=factor(j[1]);
     7888    testFactor(act,j[1]);
     7889  }
     7890  else
     7891  {
     7892    primary[1]=ideal(1);
     7893    primary[2]=ideal(1);
     7894    return(primary);
     7895  }
     7896
     7897//with the factors new ideals (hopefully the primary decomposition)
     7898//are created
     7899  if(size(act[1])>1)
     7900  {
     7901    if(size(#)>1)
     7902    {
     7903      primary[1]=ideal(1);
     7904      primary[2]=ideal(1);
     7905      primary[3]=ideal(1);
     7906      primary[4]=ideal(1);
     7907      return(primary);
     7908    }
     7909    for(@k=1;@k<=size(act[1]);@k++)
     7910    {
     7911      if(@wr==0)
     7912      {
     7913        primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
     7914      }
     7915      else
     7916      {
     7917        primary[2*@k-1]=std(j,act[1][@k]);
     7918      }
     7919      if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
     7920      {
     7921        primary[2*@k]   = primary[2*@k-1];
     7922      }
     7923      else
     7924      {
     7925        primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
     7926      }
     7927    }
     7928  }
     7929  else
     7930  {
     7931    primary[1]=j;
     7932    if((size(#)>0)&&(act[2][1]>1))
     7933    {
     7934      act[2]=1;
     7935      primary[1]=std(primary[1],act[1][1]);
     7936    }
     7937    if(@wr!=0)
     7938    {
     7939      primary[1]=std(j,act[1][1]);
     7940    }
     7941    if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
     7942    {
     7943      primary[2]=primary[1];
     7944    }
     7945    else
     7946    {
     7947      primary[2]=primaryTest(primary[1],act[1][1]);
     7948    }
     7949  }
     7950
     7951  if(size(#)==0)
     7952  {
     7953    primary=splitPrimary(primary,ser,@wr,act);
     7954  }
     7955
     7956  if((voice>=6)&&(char(basering)<=181))
     7957  {
     7958    primary=splitCharp(primary);
     7959  }
     7960
     7961  if((@wr==2)&&(npars(basering)>0)&&(voice>=6)&&(char(basering)>0))
     7962  {
     7963  //the prime decomposition of Yokoyama in characteristic p
     7964&nbs