Changeset 67bd4c in git for Singular/LIB/primdec.lib


Ignore:
Timestamp:
Nov 5, 1997, 7:16:57 PM (26 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '2fa36c576e6a4ddbb1093b43c7f8e9835e17e52a')
Children:
8e8aca339d1600a8eb489af59da08ac141f0477d
Parents:
67a4495135101f4651e8ccc9b7618eb02713c97d
Message:
* hannes/pfister/greuel: update primdec.lib, add normal.lib


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/primdec.lib

    r67a449 r67bd4c  
    1 // $Id: primdec.lib,v 1.5 1997-09-18 09:58:26 Singular Exp $
     1// $Id: primdec.lib,v 1.6 1997-11-05 18:16:57 Singular Exp $
    22///////////////////////////////////////////////////////
    33// primdec.lib
    … …  
    77//////////////////////////////////////////////////////
    88
    9 LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (G/T/Z)
    10 
    11   minAssPrimes (ideal I, list choose)
    12   // minimal associated primes
    13   // The list choose must be either emty (minAssPrimes(I)) or 1
    14   // (minAssPrimes(I,1))
    15   // In the second case the factorizing Buchberger Algorithm is used
    16   // which in most cases may considerably speed up the algorithm
    17 
    18   primdec (ideal I)
    19 
    20   // Computes a complete primary decomposition via
     9LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I)
     10
     11  minAssPrimes (ideal I)
     12  //computes the minimal associated primes of the ideal I
     13
     14  decomp (ideal I)
     15  // Computes a complete primary decomposition via Gianni,Trager,Zacharias
    2116
    2217  radical(ideal I)
    2318  //computes the radical of the ideal I
    2419
     20  equiRadical(ideal I)
     21  //computes the radical of the equidimensional part of the ideal I
     22
     23  prepareAss(ideal I)
     24  //computes the radicals of the equidimensional parts of I
     25
    2526LIB "random.lib";
    26 LIB "poly.lib";
     27LIB "elim.lib";
    2728///////////////////////////////////////////////////////////////////////////////
    2829
    … …  
    399400   int i;
    400401   ideal j=pr[1];
    401    for (i=2;i<=size(pr) div 2;i++)
     402   for (i=2;i<=size(pr)/2;i++)
    402403   {
    403404       j=intersect(j,pr[2*i-1]);
    … …  
    438439   }
    439440   int k;
    440    for (k=1;k<=size(l) div 2;k=k+1)
     441   for (k=1;k<=size(l)/2;k=k+1)
    441442   {
    442443      "                                            ";
    … …  
    485486
    486487proc primaryTest (ideal i, poly p)
    487 USAGE:   primaryTest(i,p); i ideal p poly
    488 RETURN:  ideal = radical of i, if i is primary in general position,
    489          zerodimensional and the radical of i intersected with K[z]
    490          is (p), z the smallest variable with respect to the lexico-
    491          graphical ordering, and 0 else
    492 NOTE:    It is necessary that i is a standardbasis with respect to
    493          the lexicographical ordering and the first element in i is
    494          a power of p.
    495 EXAMPLE: example primaryTest; shows an example
    496488{
    497489   int m=1;
    … …  
    501493   ideal h;
    502494
    503    //the first generator of the prim ideal for the result
    504495   ideal prm=p;
    505496   attrib(prm,"isSB",1);
    … …  
    527518      //if not (0) is returned, else var(n)+h is added to prm
    528519
    529       e=deg(lead(i[m]));
    530       t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
    531       i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
    532 
    533       if (reduce(i[m]-t^e,prm) !=0)
    534       {
    535         return(ideal(0));
    536       }
    537       h=interred(t);
    538       t=h[1];
     520         e=deg(lead(i[m]));
     521        // t=hilfe1(i,prm,m,n);
     522         t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
     523
     524         i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
     525
     526         if (reduce(i[m]-t^e,prm) !=0)
     527         {
     528           return(ideal(0));
     529         }
     530         h=interred(t);
     531         t=h[1];
    539532
    540533      prm = prm,t;
    … …  
    543536   return(prm);
    544537}
    545 example
    546 { "EXAMPLE:"; echo=2;
    547    ring  r = (0,a,b),(x,y,z),lp;
    548    poly  p = z^2+1;
    549    ideal i = p^2,(a*y-z^3)^3,(b*x-yz+z4)^4;
    550    ideal pr= primaryTest(i,p);
    551    pr;
    552    i = p^2,(y-z3)^3,(x-yz+z4)^4+1;
    553    pr= primaryTest(i,p);
    554    pr;
    555    ring s=(0,e),(d,c,b,a,y,x,g,f),lp;
    556    ideal i=f,g,x4,y,a,b3,c,d;
    557    poly p=f;
    558    ideal pr= primaryTest(i,p);
    559    pr;
    560 }
    561 
     538
     539proc hilfe(ideal i,ideal prm,int m)
     540{
     541   poly t;
     542   int e;
     543
     544   if(size(i[m])==1)
     545   {
     546      t=var(n);
     547   }
     548   else
     549   {
     550      e=deg(lead(i[m]));
     551
     552      if(deg(poly(e))>=0)
     553      {
     554           t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
     555           i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
     556      {
     557      else
     558      {
     559         i[m]=i[m]/leadcoef(i[m]);
     560         t=reduce(coef(i[m],var(n))[2,e+1],prm);
     561         t=var(n)+factorize(t,1)[1];
     562      }
     563   }
     564   return(t);
     565}
     566proc hilfe1(ideal i,ideal prm,int m,int n)
     567{
     568   poly t;
     569   int e;
     570   if(size(i[m])==1)
     571   {
     572      t=var(n);
     573   }
     574   else
     575   {
     576      e=deg(lead(i[m]));
     577      i[m]=i[m]/leadcoef(i[m]);
     578      t=reduce(coeffs(i[m],var(n))[1,1],prm);
     579      if(size(t)==0){return(var(n));}
     580      t=var(n)+factorize(t,1)[1];
     581   }
     582   return(t);
     583}
    562584
    563585///////////////////////////////////////////////////////////////////////////////
    … …  
    569591   int sl=size(l);
    570592
    571    for(i=1;i<=sl div 2;i++)
     593   for(i=1;i<=sl/2;i++)
    572594   {
    573595      if(sact[2][i]>1)
    … …  
    579601         keepprime[i]=l[2*i-1];
    580602      }
    581    }
     603  }
    582604   i=0;
    583    while(i<size(l) div 2)
     605   while(i<size(l)/2)
    584606   {
    585607      i++;
    … …  
    600622         }
    601623         j=0;
    602          if(i<=sl div 2)
     624         if(i<=sl/2)
    603625         {
    604626            j=1;
    … …  
    628650               for(k=2;k<=r;k++)
    629651               {
    630                   keepprime[size(l) div 2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
     652                  keepprime[size(l)/2+k-1]=interred(keepprime[i]+ideal(act[1][k]));
    631653               }
    632654               keepprime[i]=interred(keepprime[i]+ideal(act[1][1]));
    … …  
    655677               {
    656678                  l[s+2*k-1]=keepresult[k];
    657                   keepprime[s div 2+k]=interred(keepresult[k]+ideal(act[1][k]));
     679                  keepprime[s/2+k]=interred(keepresult[k]+ideal(act[1][k]));
    658680                  if(vdim(keepresult[k])==deg(act[1][k]))
    659681                  {
    … …  
    664686                     l[s+2*k]=ideal(0);
    665687                  }
    666                   if((homog(keepresult[k])==1)||(homog(keepprime[s div 2+k])==1))
     688                  if((homog(keepresult[k])==1)||(homog(keepprime[s/2+k])==1))
    667689                  {
    668690                    l[s+2*k]=maxideal(1);
    … …  
    688710                     l[s+2]=ideal(0);
    689711                  }
    690                   keepprime[s div 2+1]=interred(keepprime[i]+ideal(@f));
    691                   if(homog(keepprime[s div 2+1])==1)
     712                  keepprime[s/2+1]=interred(keepprime[i]+ideal(@f));
     713                  if(homog(keepprime[s/2+1])==1)
    692714                  {
    693715                     l[s+2]=maxideal(1);
    … …  
    712734      return(l);
    713735   }
    714    for(i=1;i<=size(l) div 2;i++)
     736   for(i=1;i<=size(l)/2;i++)
    715737   {
    716738      if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1))
    … …  
    900922
    901923  @k=0;
    902   while(@k<(size(primary) div 2))
     924  int zz;
     925  while(@k<(size(primary)/2))
    903926  {
    904927    @k++;
    905928    if (size(primary[2*@k])==0)
    906929    {
     930       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
     931       {
     932          if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
     933          {
     934             primary[2*@k]=primary[2*@k-1];
     935          }
     936       }
     937    }
     938  }
     939  @k=0;
     940  while(@k<(size(primary)/2))
     941  {
     942    @k++;
     943    if (size(primary[2*@k])==0)
     944    {
     945//      "the univariat polynomials";
     946//       int qwe=timer;
     947//       system("finduni",primary[2*@k-1]);
     948
    907949       jmap=randomLast(100);
     950//       timer-qwe;
     951//       primary[2*@k-1];
     952//       pause;
    908953       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
    909954       {
    … …  
    10021047       else
    10031048       {
    1004           for(@n=1;@n<=size(@lr) div 2;@n++)
     1049          for(@n=1;@n<=size(@lr)/2;@n++)
    10051050          {
    10061051             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
    … …  
    10371082       else
    10381083       {
    1039           for(@n=1;@n<=size(@lr) div 2;@n++)
     1084          for(@n=1;@n<=size(@lr)/2;@n++)
    10401085          {
    10411086             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
    … …  
    10671112       primary[2*@k-1]=lres[1];
    10681113       primary[2*@k]=lres[2];
    1069        @s=size(primary) div 2;
    1070        for(@n=1;@n<=size(lres) div 2-1;@n++)
     1114       @s=size(primary)/2;
     1115       for(@n=1;@n<=size(lres)/2-1;@n++)
    10711116       {
    10721117         primary[2*@s+2*@n-1]=lres[2*@n+1];
    … …  
    11111156        return(1)
    11121157     }
    1113      p=gcd(p,i[k]);
     1158     p=GCD(p,i[k]);
    11141159     if(deg(p)==0)
    11151160     {
    … …  
    11551200   if(size(@pa)==0)
    11561201   {
    1157       return(lcm(h));
     1202      return(lcmP(h));
    11581203   }
    11591204   def bsr= basering;
    … …  
    11611206   execute "ring @r=0,("+@pa+","+varstr(bsr)+"),dp;";
    11621207   execute "ideal @i="+@id+";";
    1163    poly @p=lcm(@i);
     1208   poly @p=lcmP(@i);
    11641209   string @ps=string(@p);
    11651210   setring bsr;
    … …  
    11751220   ideal l=p,q;
    11761221   poly  pr= coeffLcm(l);
     1222   pr;
     1223}
     1224
     1225///////////////////////////////////////////////////////////////////////////////
     1226
     1227proc lcmP(ideal i)
     1228USAGE:   lcm(i); i list of polynomials
     1229RETURN:  poly = lcm(i[1],...,i[size(i)])
     1230NOTE:
     1231EXAMPLE: example lcm; shows an example
     1232{
     1233  int k,j;
     1234   poly p,q;
     1235  i=simplify(i,10);
     1236  for(j=1;j<=size(i);j++)
     1237  {
     1238    if(deg(i[j])>0)
     1239    {
     1240      p=i[j];
     1241      break;
     1242    }
     1243  }
     1244  if(deg(p)==-1)
     1245  {
     1246    return(1);
     1247  }
     1248  for (k=j+1;k<=size(i);k++)
     1249  {
     1250     if(deg(i[k])!=0)
     1251     {
     1252        q=GCD(p,i[k]);
     1253        if(deg(q)==0)
     1254        {
     1255           p=p*i[k];
     1256        }
     1257        else
     1258        {
     1259           p=p/q;
     1260           p=p*i[k];
     1261        }
     1262     }
     1263   }
     1264  return(p);
     1265}
     1266example
     1267{ "EXAMPLE:"; echo = 2;
     1268   ring  r = 0,(x,y,z),lp;
     1269   poly  p = (x+y)*(y+z);
     1270   poly  q = (z4+2)*(y+z);
     1271   ideal l=p,q;
     1272   poly  pr= lcmP(l);
     1273   pr;
     1274   l=1,-1,p,1,-1,q,1;
     1275   pr=lcmP(l);
    11771276   pr;
    11781277}
    … …  
    14341533}
    14351534
     1535///////////////////////////////////////////////////////////////////////
     1536
     1537proc projdim(list l)
     1538{
     1539   int i=size(l)+1;
     1540
     1541   while(i>2)
     1542   {
     1543      i--;
     1544      if((size(l[i])>0)&&(deg(l[i][1])>0))
     1545      {
     1546         return(i);
     1547      }
     1548   }
     1549}
     1550
    14361551///////////////////////////////////////////////////////////////////////////////
    14371552proc cleanPrimary(list l)
    … …  
    14391554   int i,j;
    14401555   list lh;
    1441    for(i=1;i<=size(l) div 2;i++)
     1556   for(i=1;i<=size(l)/2;i++)
    14421557   {
    14431558      if(deg(l[2*i-1][1])>0)
    … …  
    14601575EXAMPLE: example minAssPrimes; shows an example
    14611576{
     1577   #[1]=1;
    14621578   def @P=basering;
    1463    execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
    1464 //   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),dp;";
     1579   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
     1580             +ordstr(basering)+");";
     1581
     1582
    14651583   ideal i=fetch(@P,i);
    14661584   if(size(#)==0)
    … …  
    14781596   }
    14791597   list @res,empty;
     1598   ideal ser;
    14801599   option(redSB);
    14811600   list @pr=facstd(i);
     1601   if(size(@pr)==1)
     1602   {
     1603      attrib(@pr[1],"isSB",1);
     1604      if((dim(@pr[1])==0)&&(homog(@pr[1])==1))
     1605      {
     1606         setring @P;
     1607         list @res=maxideal(1);
     1608         return(@res);
     1609      }
     1610      if(dim(@pr[1])>1)
     1611      {
     1612         setring @P;
     1613         kill gnir;
     1614         execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
     1615         ideal i=fetch(@P,i);
     1616         list @pr=facstd(i);
     1617         ideal ser;
     1618      }
     1619   }
    14821620   option(noredSB);
    14831621   int j,k,odim,ndim,count;
    … …  
    15211659   else
    15221660   {
     1661     ser=ideal(1);
    15231662     for(j=1;j<=size(@pr);j++)
    15241663     {
    1525        @res[j]=decomp(@pr[j],2);
     1664         @res[j]=decomp(@pr[j],2);
     1665 //      @res[j]=decomp(@pr[j],2,@pr[j],ser);
     1666 //      for(k=1;k<=size(@res[j]);k++)
     1667 //      {
     1668 //         ser=intersect1(ser,@res[j][k]);
     1669 //      }
    15261670     }
    15271671   }
    … …  
    15461690///////////////////////////////////////////////////////////////////////////////
    15471691
    1548 proc union(list l)
     1692proc union(list li)
    15491693{
    15501694  int i,j,k;
     1695
     1696  def P=basering;
     1697
     1698  execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
     1699  list l=fetch(P,li);
    15511700  list @erg;
    1552   i=0;
    15531701
    15541702  for(k=1;k<=size(l);k++)
    15551703  {
    1556      for(j=1;j<=size(l[k]) div 2;j++)
     1704     for(j=1;j<=size(l[k])/2;j++)
    15571705     {
    15581706        if(deg(l[k][2*j][1])!=0)
    … …  
    16111759     @wos[size(@wos)+1]=@erg[size(@erg)];
    16121760  }
    1613   return(@wos);
     1761  setring P;
     1762  list @ser=fetch(ir,@wos);
     1763  return(@ser);
    16141764}
    16151765///////////////////////////////////////////////////////////////////////////////
    1616 proc radical(ideal i)
     1766proc radicalOld(ideal i)
    16171767{
    16181768   list pr=minAssPrimes(i,1);
    … …  
    16241774   }
    16251775   return(k);
     1776}
     1777///////////////////////////////////////////////////////////////////////////////
     1778proc equiRadical(ideal i)
     1779{
     1780   return(radical(i,1));
    16261781}
    16271782///////////////////////////////////////////////////////////////////////////////
    … …  
    16361791{
    16371792  def  @P = basering;
    1638   list primary,indep;
     1793  list primary,indep,ltras;
    16391794  intvec @vh,isat;
    16401795  int @wr,@k,@n,@m,@n1,@n2,@n3,homo;
    … …  
    16651820  if(homo==1)
    16661821  {
    1667      tras=std(i);
     1822     ltras=mstd(i);
     1823     attrib(ltras[1],"isSB",1);
     1824     tras=ltras[1];
    16681825     if(dim(tras)==0)
    16691826     {
    1670         primary[1]=tras;
     1827        primary[1]=ltras[2];
    16711828        primary[2]=maxideal(1);
    16721829        if(@wr>0)
    … …  
    16791836        return(primary);
    16801837     }
    1681      intvec @hilb=hilb(tras,1);
     1838      intvec @hilb=hilb(tras,1);
    16821839  }
    16831840
    … …  
    17371894  if(nvars(basering)==1)
    17381895  {
     1896
    17391897     list fac=factor(@j[1]);
    17401898     list gprimary;
    … …  
    17761934      return(primary);
    17771935   }
    1778 
    1779 
    1780   //------------------------------------------------------------------
    1781   //search for a maximal independent set indep,i.e.
    1782   //look for subring such that the intersection with the ideal is zero
    1783   //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
    1784   //indep[1] is the new varstring and indep[2] the string for the block-ordering
    1785   //------------------------------------------------------------------
    17861936
    17871937  poly @gs,@gh,@p;
    … …  
    17921942  list fett;
    17931943  int lauf,di;
     1944  //------------------------------------------------------------------
     1945  //search for a maximal independent set indep,i.e.
     1946  //look for subring such that the intersection with the ideal is zero
     1947  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
     1948  //indep[1] is the new varstring and indep[2] the string for the block-ordering
     1949  //------------------------------------------------------------------
    17941950
    17951951  if(@wr!=1)
    … …  
    19592115     @n3=@n2;
    19602116
    1961      for(@n1=1;@n1<=size(collectprimary) div 2;@n1++)
     2117     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
    19622118     {
    19632119        if(deg(collectprimary[2*@n1][1])>0)
    … …  
    19752131     //mentioned above is really computed
    19762132
    1977     for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++)
     2133    for(@n=@n3/2+1;@n<=@n2/2;@n++)
    19782134     {
    19792135        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
    … …  
    20422198     quprimary[2]=ideal(1);
    20432199  }
     2200
    20442201  //---------------------------------------------------------------
    20452202  //notice that j=sat(j,gh) intersected with (j,gh^n)
    … …  
    20552212           ideal htest=quprimary[1];
    20562213
    2057            for (@n1=2;@n1<=size(quprimary) div 2;@n1++)
     2214           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
    20582215           {
    20592216              htest=intersect(htest,quprimary[2*@n1-1]);
    … …  
    20642221           ideal htest=quprimary[2];
    20652222
    2066            for (@n1=2;@n1<=size(quprimary) div 2;@n1++)
     2223           for (@n1=2;@n1<=size(quprimary)/2;@n1++)
    20672224           {
    20682225              htest=intersect(htest,quprimary[2*@n1]);
    … …  
    21932350           @n3=@n2;
    21942351
    2195            for(@n1=1;@n1<=size(collectprimary) div 2;@n1++)
     2352           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
    21962353           {
    21972354              if(deg(collectprimary[2*@n1][1])>0)
    … …  
    22092366           //mentioned above is really computed
    22102367
    2211            for(@n=@n3 div 2+1;@n<=@n2 div 2;@n++)
     2368           for(@n=@n3/2+1;@n<=@n2/2;@n++)
    22122369           {
    22132370              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
    … …  
    22792436   testPrimary( pr, i);
    22802437}
     2438
     2439///////////////////////////////////////////////////////////////////////////////
     2440proc radical(ideal i,list #)
     2441{
     2442   def @P=basering;
     2443   int j,il;
     2444   if(size(#)>0)
     2445   {
     2446     il=#[1];
     2447   }
     2448   ideal re=1;
     2449   option(redSB);
     2450   list pr=facstd(i);
     2451
     2452   if(size(pr)==1)
     2453   {
     2454      attrib(pr[1],"isSB",1);
     2455      if((dim(pr[1])==0)&&(homog(pr[1])==1))
     2456      {
     2457         ideal @res=maxideal(1);
     2458         return(@res);
     2459      }
     2460      if(dim(pr[1])>1)
     2461      {
     2462         execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;";
     2463         ideal i=fetch(@P,i);
     2464         list @pr=facstd(i);
     2465         setring @P;
     2466         pr=fetch(gnir,@pr);
     2467      }
     2468   }
     2469   option(noredSB);
     2470   int s=size(pr);
     2471   if(s==1)
     2472   {
     2473      return(radicalEHV(i,ideal(1),il));
     2474   }
     2475   intvec pos;
     2476   pos[s]=0;
     2477
     2478   if(il==1)
     2479   {
     2480     int ndim,k;
     2481     attrib(pr[1],"isSB",1);
     2482     int odim=dim(pr[1]);
     2483     int count=1;
     2484
     2485     for(j=2;j<=s;j++)
     2486     {
     2487        attrib(pr[j],"isSB",1);
     2488        ndim=dim(pr[j]);
     2489        if(ndim>odim)
     2490        {
     2491           for(k=count;k<=j-1;k++)
     2492           {
     2493              pos[k]=1;
     2494           }
     2495           count=j;
     2496           odim=ndim;
     2497        }
     2498        if(ndim<odim)
     2499        {
     2500           pos[j]=1;
     2501        }
     2502     }
     2503   }
     2504
     2505   for(j=1;j<=s;j++)
     2506   {
     2507      if(pos[j]==0)
     2508      {
     2509         re=intersect1(re,radicalEHV(pr[s+1-j],re,il));
     2510      }
     2511   }
     2512   return(re);
     2513}
     2514
     2515proc intersect1(ideal i,ideal j)
     2516{
     2517   def R=basering;
     2518   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+",@t),dp;";
     2519   ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
     2520   ideal j=eliminate(i,var(nvars(basering)));
     2521   setring R;
     2522   map phi=gnir,maxideal(1);
     2523   return(phi(j));
     2524}
     2525
     2526
     2527///////////////////////////////////////////////////////////////////////////////
     2528proc radicalKL (list m,ideal ser,list #)
     2529USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
     2530         decomp(i,1);        (for the minimal associated primes) )
     2531RETURN:  list = list of primary ideals and their associated primes
     2532         (at even positions in the list)
     2533         (resp. a list of the minimal associated primes)
     2534NOTE:    Algorithm of Gianni, Traeger, Zacharias
     2535EXAMPLE: example decomp; shows an example
     2536{
     2537   ideal i=m[2];
     2538  //----------------------------------------------------------------
     2539  //i is the zero-ideal
     2540  //----------------------------------------------------------------
     2541
     2542   if(size(i)==0)
     2543   {
     2544      return(ideal(0));
     2545   }
     2546
     2547   def  @P = basering;
     2548   list indep,allindep,restindep,fett,@mu;
     2549   intvec @vh,isat;
     2550   int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di;
     2551   ideal @j,@j1,fac,@h,collectrad,htrad,lsau;
     2552   ideal rad=ideal(1);
     2553   ideal te=ser;
     2554   poly @p,@q;
     2555   string @va,quotring;
     2556   int  homo=homog(i);
     2557
     2558   if(size(#)>0)
     2559   {
     2560      @wr=#[1];
     2561   }
     2562   @j=m[1];
     2563   @j1=m[2];
     2564   int jdim=dim(@j);
     2565   if(size(reduce(ser,@j))==0)
     2566   {
     2567      return(ser);
     2568   }
     2569   if(homo==1)
     2570   {
     2571      if(jdim==0)
     2572      {
     2573         option(noredSB);
     2574         return(maxideal(1));
     2575      }
     2576      intvec @hilb=hilb(@j,1);
     2577   }
     2578
     2579
     2580  //----------------------------------------------------------------
     2581  //j is the ring
     2582  //----------------------------------------------------------------
     2583
     2584   if (jdim==-1)
     2585   {
     2586      option(noredSB);
     2587      return(ideal(0));
     2588   }
     2589
     2590  //----------------------------------------------------------------
     2591  //  the case of one variable
     2592  //----------------------------------------------------------------
     2593
     2594   if(nvars(basering)==1)
     2595   {
     2596      fac=factorize(@j[1],1);
     2597      poly @p=1;
     2598      for(@k=1;@k<=size(fac);@k++)
     2599      {
     2600         @p=@p*fac[@k];
     2601      }
     2602      option(noredSB);
     2603
     2604      return(ideal(@p));
     2605   }
     2606   //------------------------------------------------------------------
     2607   //the case of a complete intersection
     2608   //------------------------------------------------------------------
     2609    if(jdim+size(@j1)==nvars(basering))
     2610    {
     2611      // ideal jac=minor(jacob(@j1),size(@j1));
     2612      // return(quotient(@j1,jac));
     2613    }
     2614
     2615   //------------------------------------------------------------------
     2616   //the zero-dimensional case
     2617   //------------------------------------------------------------------
     2618
     2619   if (jdim==0)
     2620   {
     2621      @j1=system("finduni",@j);
     2622      for(@k=1;@k<=size(@j1);@k++)
     2623      {
     2624         fac=factorize(cleardenom(@j1[@k]),1);
     2625         @p=fac[1];
     2626         for(@n=2;@n<=size(fac);@n++)
     2627         {
     2628            @p=@p*fac[@n];
     2629         }
     2630         @j=@j,@p;
     2631      }
     2632      @j=std(@j);
     2633      option(noredSB);
     2634      return(@j);
     2635   }
     2636
     2637   //------------------------------------------------------------------
     2638   //search for a maximal independent set indep,i.e.
     2639   //look for subring such that the intersection with the ideal is zero
     2640   //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
     2641   //indep[1] is the new varstring and indep[2] the string for the block-ordering
     2642   //------------------------------------------------------------------
     2643
     2644   indep=maxIndependSet(@j);
     2645
     2646   di=dim(@j);
     2647
     2648   for(@m=1;@m<=size(indep);@m++)
     2649   {
     2650      if((indep[@m][1]==varstr(basering))&&(@m==1))
     2651      //this is the good case, nothing to do, just to have the same notations
     2652      //change the ring
     2653      {
     2654         execute "ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
     2655                              +ordstr(basering)+");";
     2656         ideal @j=fetch(@P,@j);
     2657         attrib(@j,"isSB",1);
     2658      }
     2659      else
     2660      {
     2661         @va=string(maxideal(1));
     2662         execute "ring gnir1 = ("+charstr(basering)+"),("+indep[@m][1]+"),("
     2663                              +indep[@m][2]+");";
     2664         execute "map phi=@P,"+@va+";";
     2665         if(homo==1)
     2666         {
     2667            ideal @j=std(phi(@j),@hilb);
     2668         }
     2669         else
     2670         {
     2671           ideal @j=std(phi(@j));
     2672         }
     2673      }
     2674      if((deg(@j[1])==0)||(dim(@j)<jdim))
     2675      {
     2676         setring @P;
     2677         break;
     2678      }
     2679      for (lauf=1;lauf<=size(@j);lauf++)
     2680      {
     2681         fett[lauf]=size(@j[lauf]);
     2682      }
     2683     //------------------------------------------------------------------------------------
     2684     //we have now the following situation:
     2685     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     2686     //to this quotientring, j is their still a standardbasis, the
     2687     //leading coefficients of the polynomials  there (polynomials in
     2688     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     2689     //we need their ggt, gh, because of the following:
     2690     //let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
     2691     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     2692     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
     2693
     2694     //------------------------------------------------------------------------------------
     2695
     2696     //the arrangement for the quotientring K(var(nnp+1),..,var(nva))[..the rest..]
     2697     //and the map phi:K[var(1),...,var(nva)] ----->K(var(nnpr+1),..,var(nva))[..the rest..]
     2698     //-------------------------------------------------------------------------------------
     2699
     2700      quotring=prepareQuotientring(nvars(basering)-indep[@m][3]);
     2701
     2702     //---------------------------------------------------------------------
     2703     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     2704     //---------------------------------------------------------------------
     2705
     2706      execute quotring;
     2707
     2708    // @j considered in the quotientring
     2709      ideal @j=imap(gnir1,@j);
     2710
     2711      kill gnir1;
     2712
     2713     //j is a standardbasis in the quotientring but usually not minimal
     2714     //here it becomes minimal
     2715
     2716      @j=clearSB(@j,fett);
     2717      attrib(@j,"isSB",1);
     2718
     2719     //we need later ggt(h[1],...)=gh for saturation
     2720      ideal @h;
     2721      if(deg(@j[1])>0)
     2722      {
     2723         for(@n=1;@n<=size(@j);@n++)
     2724         {
     2725            @h[@n]=leadcoef(@j[@n]);
     2726         }
     2727        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
     2728         option(redSB);
     2729         @j=interred(@j);
     2730         attrib(@j,"isSB",1);
     2731         list  @mo=@j,@j;
     2732         ideal zero_rad= radicalKL(@mo,ideal(1));
     2733      }
     2734      else
     2735      {
     2736         ideal zero_rad=ideal(1);
     2737      }
     2738
     2739     //we need the intersection of the ideals in the list quprimary with the
     2740     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     2741     //but fi polynomials, then the intersection of q with the polynomialring
     2742     //is the saturation of the ideal generated by f1,...,fr with respect to
     2743     //h which is the lcm of the leading coefficients of the fi considered in the
     2744     //quotientring: this is coded in saturn
     2745
     2746      ideal hpl;
     2747
     2748      for(@n=1;@n<=size(zero_rad);@n++)
     2749      {
     2750         hpl=hpl,leadcoef(zero_rad[@n]);
     2751      }
     2752
     2753     //--------------------------------------------------------------------
     2754     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     2755     //back to the polynomialring
     2756     //---------------------------------------------------------------------
     2757      setring @P;
     2758
     2759      collectrad=imap(quring,zero_rad);
     2760      lsau=simplify(imap(quring,hpl),2);
     2761      @h=imap(quring,@h);
     2762
     2763      kill quring;
     2764
     2765
     2766     //here the intersection with the polynomialring
     2767     //mentioned above is really computed
     2768
     2769      collectrad=sat2(collectrad,lsau)[1];
     2770
     2771      if(deg(@h[1])>=0)
     2772      {
     2773         fac=ideal(0);
     2774         for(lauf=1;lauf<=ncols(@h);lauf++)
     2775         {
     2776            if(deg(@h[lauf])>0)
     2777            {
     2778                fac=fac+factorize(@h[lauf],1);
     2779            }
     2780         }
     2781         fac=simplify(fac,4);
     2782         @q=1;
     2783         for(lauf=1;lauf<=size(fac);lauf++)
     2784         {
     2785            @q=@q*fac[lauf];
     2786         }
     2787
     2788
     2789         @mu=mstd(quotient(@j+ideal(@q),rad));
     2790         @j=@mu[1];
     2791         attrib(@j,"isSB",1);
     2792
     2793      }
     2794      if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
     2795      {
     2796int xyz=timer;
     2797"bei collecterad";
     2798         rad=intersect(rad,collectrad);
     2799timer-xyz;
     2800      }
     2801      else
     2802      {
     2803         if(deg(collectrad[1])>0)
     2804         {
     2805            rad=collectrad;
     2806         }
     2807      }
     2808
     2809      te=simplify(reduce(te*rad,@j),2);
     2810
     2811      if((dim(@j)<di)||(size(te)==0))
     2812      {
     2813         break;
     2814      }
     2815      if(homo==1)
     2816      {
     2817         @hilb=hilb(@j,1);
     2818      }
     2819   }
     2820
     2821   if(((@wr==1)&&(dim(@j)<di))||(deg(@j[1])==0)||(size(te)==0))
     2822   {
     2823      return(rad);
     2824   }
     2825   ideal tes=radicalKL(@mu,rad,@wr);
     2826int sml=timer;
     2827"bei rad";
     2828   rad=intersect(rad,tes);
     2829timer-sml;
     2830  // rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
     2831
     2832
     2833   option(noredSB);
     2834   return(rad);
     2835}
     2836
     2837
     2838example
     2839{ "EXAMPLE:"; echo = 2;
     2840}
     2841
     2842
     2843proc radicalEHV(ideal i,ideal re,list #)
     2844{
     2845   ideal J,I,I0,radI0,L,radI1,I2,radI2;
     2846   int l,il;
     2847   if(size(#)>0)
     2848   {
     2849      il=#[1];
     2850   }
     2851   option(redSB);
     2852   list m=mstd(i);
     2853        I=m[2];
     2854   option(noredSB);
     2855   if(size(reduce(re,m[1]))==0)
     2856   {
     2857      return(re);
     2858   }
     2859   int cod=nvars(basering)-dim(m[1]);
     2860   if(nvars(basering)<9)
     2861   {
     2862      if(cod==size(m[2]))
     2863      {
     2864         J=minor(jacob(I),cod);
     2865         return(quotient(I,J));
     2866      }
     2867
     2868      for(l=1;l<=cod;l++)
     2869      {
     2870         I0[l]=I[l];
     2871      }
     2872      if(dim(std(I0))+cod==nvars(basering))
     2873      {
     2874         J=minor(jacob(I0),cod);
     2875         radI0=quotient(I0,J);
     2876         L=quotient(radI0,I);
     2877         radI1=quotient(radI0,L);
     2878
     2879         if(size(reduce(radI1,m[1]))==0)
     2880         {
     2881            return(I);
     2882         }
     2883         if(il==1)
     2884         {
     2885            return(radI1);
     2886         }
     2887
     2888         I2=sat(I,radI1)[1];
     2889
     2890         if(deg(I2[1])<=0)
     2891         {
     2892            return(radI1);
     2893         }
     2894         return(intersect(radI1,radicalEHV(I2,re,il)));
     2895      }
     2896   }
     2897   return(radicalKL(m,re,il));
     2898}
     2899
     2900proc Ann(module M)
     2901{
     2902  M=prune(M);  //to obtain a small embedding
     2903  return(quotient(M,freemodule(nrows(M))));
     2904}
     2905
     2906//computes the equidimensional part of the ideal i of codimension e
     2907proc int_ass_primary_e(ideal i, int e)
     2908{
     2909  if(homog(i)!=1)
     2910  {
     2911     i=std(i);
     2912  }
     2913  list re=sres(i,0);                   //the resolution
     2914  re=minres(re);                       //minimized resolution
     2915  ideal ann=AnnExt_R(e,re);
     2916  if(nvars(basering)-dim(std(ann))!=e)
     2917  {
     2918    return(ideal(1));
     2919  }
     2920  return(ann);
     2921}
     2922
     2923//computes all equidimensional parts of the ideal i
     2924proc prepareAss(ideal i)
     2925{
     2926  ideal j=std(i);
     2927  int cod=nvars(basering)-dim(j);
     2928  int e;
     2929  list er;
     2930  ideal ann;
     2931  if(homog(i)==1)
     2932  {
     2933     list re=sres(i,0);                   //the resolution
     2934     re=minres(re);                       //minimized resolution
     2935  }
     2936  else
     2937  {
     2938     list re=mres(i,0);             //fehler in sres
     2939  }
     2940  for(e=cod;e<=nvars(basering);e++)
     2941  {
     2942     ann=AnnExt_R(e,re);
     2943     if(nvars(basering)-dim(std(ann))==e)
     2944     {
     2945        er[size(er)+1]=equiRadical(ann);
     2946     }
     2947  }
     2948  return(er);
     2949}
     2950
     2951//computes the annihilator of Ext^n(R/i,R) with given resolution re
     2952//n is not necessarily the number of variables
     2953proc AnnExt_R(int n,list re)
     2954{
     2955  if(n<nvars(basering))
     2956  {
     2957     matrix f=transpose(re[n+1]);      //Hom(_,R)
     2958     module k=res(f,2)[2];             //the kernel
     2959     matrix g=transpose(re[n]);        //the image of Hom(_,R)
     2960     ideal ann=quotient(g,k);           //the anihilator
     2961  }
     2962  else
     2963  {
     2964     ideal ann=Ann(transpose(re[n]));
     2965  }
     2966  return(ann);
     2967}
     2968
Note: See TracChangeset for help on using the changeset viewer.