Changeset ebecf83 in git


Ignore:
Timestamp:
May 15, 1998, 6:25:48 PM (26 years ago)
Author:
Olaf Bachmann <obachman@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '38dfc5131670d387a89455159ed1e071997eec94')
Children:
b7792751599796c0e032684a87f5a80b650a5eb1
Parents:
ac6554110ea01f1f5b8bb9a4d561299931b8ba7d
Message:
* commited Gerhard's polished version


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/primdec.lib

    rac6554 rebecf83  
    1 // $Id: primdec.lib,v 1.18 1998-05-14 18:45:13 Singular Exp $
    2 ///////////////////////////////////////////////////////
    3 // primdec.lib
    4 // algorithms for primary decomposition based on
    5 // the ideas of Gianni,Trager,Zacharias
    6 // written by Gerhard Pfister
    7 //
    8 // algorithms for primary decomposition based on
    9 // the ideas of Shimoyama/Yokoyama
    10 // written by Wolfram Decker and Hans Schoenemann
    11 //////////////////////////////////////////////////////
    12 
    13 version="$Id: primdec.lib,v 1.18 1998-05-14 18:45:13 Singular Exp $";
     1// $Id: primdec.lib,v 1.19 1998-05-15 16:25:48 obachman Exp $
     2////////////////////////////////////////////////////////////////////////////////
     3// primdec.lib                                                                //
     4// algorithms for primary decomposition based on                              //
     5// the ideas of Gianni,Trager,Zacharias                                       //
     6// written by Gerhard Pfister                                                 //
     7//                                                                            //
     8// algorithms for primary decomposition based on                              //
     9// the ideas of Shimoyama/Yokoyama                                            //
     10// written by Wolfram Decker and Hans Schoenemann                             //
     11////////////////////////////////////////////////////////////////////////////////
     12
     13version="$Id: primdec.lib,v 1.19 1998-05-15 16:25:48 obachman Exp $";
    1414info="
    15 LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION (I)
    16 
    17   minAssPrimes (ideal I)
    18   //computes the minimal associated primes of the ideal I
    19 
    20   primdecGTZ (ideal I)
    21   // Computes a complete primary decomposition via Gianni,Trager,Zacharias
    22 
    23   radical(ideal I)
    24   //computes the radical of the ideal I
    25 
    26   equiRadical(ideal I)
    27   //computes the radical of the equidimensional part of the ideal I
    28 
    29   prepareAss(ideal I)
    30   //computes the radicals of the equidimensional parts of I
    31 
    32   min_ass_prim_charsets (ideal I, int choose)
    33   // minimal associated primes via the  characteristic set
    34   // package written by Michael Messollen.
    35   // The integer choose must be either 0 or 1.
    36   // If choose=0, the given ordering of the variables is used.
    37   // If choose=1, the system tries to find an \"optimal ordering\",
    38   // which in some cases may considerably speed up the algorithm
    39 
    40   // You may also may want to try one of the algorithms for
    41   // minimal associated primes in the the library
    42   // primdec.lib, written by Gerhard Pfister.
    43   // These algorithms are variants of the algorithm
    44   // of Gianni-Trager-Zacharias
    45 
    46   primdecSY (ideal I, int choose)
    47 
    48   // Computes a complete primary decomposition via
    49   // a variant of the pseudoprimary approach of
    50   // Shimoyama-Yokoyama.
    51   // The integer choose must be either 0, 1, 2 or 3.
    52   // If choose=0, min_ass_prim_charsets with the given
    53   // ordering of the variables is used.
    54   // If choose=1, min_ass_prim_charsets with the \"optimized\"
    55   // ordering of the variables is used.
    56   // If choose=2, minAssPrimes from primdec.lib is used
    57   // If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used
     15LIBRARY: primdec.lib: PROCEDURE FOR PRIMARY DECOMPOSITION
     16
     17  minAssGTZ(I);     computes the minimal associated primes via Gianni,Trager,Zacharias
     18
     19  minAssChar(I);    computes the minimal associated primes using characteristic sets
     20
     21  primdecGTZ(I);    computes a complete primary decomposition via Gianni,Trager,Zacharias
     22
     23  primdecSY(I);     computes a complete primary decomposition via Shimoyama-Yokoyama
     24
     25  testPrimary(L,k); tests whether the result of the primary decomposition is correct
     26
     27  radical(I);       computes the radical of the ideal I
     28
     29  equiRadical(I);   computes the radical of the equidimensional part of the ideal I
     30
     31  prepareAss(I);    computes the radicals of the equidimensional parts of I
     32
    5833";
    5934
    … …  
    6237LIB "poly.lib";
    6338LIB "random.lib";
     39///////////////////////////////////////////////////////////////////////////////
     40
     41///////////////////////////////////////////////////////////////////////////////
     42//
     43//   Gianni/Trager/Zacharias
     44//
     45//
    6446///////////////////////////////////////////////////////////////////////////////
    6547
    … …  
    213195   ideal fac=tsil[2];
    214196   poly f=tsil[3];
    215 
     197   
    216198   ideal star=quotient(laedi,f);
    217199   action=1;
    … …  
    228210        g=1;
    229211        verg=laedi;
    230 
     212       
    231213         for(j=1;j<=size(fac);j++)
    232214         {
    … …  
    237219         }
    238220         verg=quotient(laedi,g);
    239 
     221         
    240222         if(specialIdealsEqual(verg,star)==1)
    241223         {
    … …  
    251233      }
    252234   }
    253    l=star,fac,f;
     235   l=star,fac,f;   
    254236   return(l);
    255237}
    … …  
    259241{
    260242  poly keep=p;
    261 
     243 
    262244   int i;
    263245   poly q=act[1][1]^act[2][1];
    … …  
    272254      "ERROR IN FACTOR";
    273255      basering;
    274 
     256     
    275257      act;
    276258      keep;
    277259      pause;
    278 
     260     
    279261      p;
    280262      q;
    … …  
    438420
    439421
    440 ////////////////////////////////////////////////////////////////////////////////
    441 
    442 proc testPrimary(list pr, ideal k)
    443 "USAGE:   testPrimary(pr,k) pr list, k ideal;
    444 RETURN:  int = 1, if the intersection of the ideals in pr is k, 0 if not
    445 NOTE:
    446 EXAMPLE: example testPrimary ; shows an example
    447 "
    448 {
    449    int i;
    450    ideal j=pr[1];
    451    for (i=2;i<=size(pr)/2;i++)
    452    {
    453        j=intersect(j,pr[2*i-1]);
    454    }
    455    return(idealsEqual(j,k));
    456 }
    457 example
    458 { "EXAMPLE:"; echo = 2;
    459    ring  s = 0,(x,y,z),lp;
    460    ideal i=x3-x2-x+1,xy-y;
    461    ideal i1=x-1;
    462    ideal i2=x-1;
    463    ideal i3=y,x2-2x+1;
    464    ideal i4=y,x-1;
    465    ideal i5=y,x+1;
    466    ideal i6=y,x+1;
    467    list pr=i1,i2,i3,i4,i5,i6;
    468    testPrimary(pr,i);
    469    pr[5]=y+1,x+1;
    470    testPrimary(pr,i);
    471 }
    472 
    473 ////////////////////////////////////////////////////////////////////////////////
    474 proc printPrimary( list l, list #)
    475 "USAGE:   printPrimary(l) l list;
    476 RETURN:  nothing
    477 NOTE:
    478 EXAMPLE: example printPrimary; shows an example
    479 "
    480 {
    481    if(size(#)>0)
    482    {
    483       "                                            ";
    484       "  The primary decomposition of the ideal    ";
    485       #[1];
    486       "                                            ";
    487       "  is:                                       ";
    488       "                                            ";
    489    }
    490    int k;
    491    for (k=1;k<=size(l)/2;k=k+1)
    492    {
    493       "                                            ";
    494       "primary ideal:                              ";
    495       l[2*k-1];
    496       "                                            ";
    497       "associated prime ideal                      ";
    498       l[2*k];
    499       "                                            ";
    500    }
    501 }
    502 example
    503 { "EXAMPLE:"; echo = 2;
    504 }
    505 
    506 ////////////////////////////////////////////////////////////////////////////////
    507422
    508423
    … …  
    565480         m=m-1;
    566481      }
    567       //check whether i[m] =(c*var(n)+h)^e modulo prm for some
     482      //check whether i[m] =(c*var(n)+h)^e modulo prm for some 
    568483      //h in K[var(n+1),...,var(nvars(basering))], c in K
    569484      //if not (0) is returned, else var(n)+h is added to prm
    570485
    571486         e=deg(lead(i[m]));
    572         // t=hilfe1(i,prm,m,n);
    573487         t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
    574488
    … …  
    588502}
    589503
    590 proc hilfe(ideal i,ideal prm,int m)
    591 {
    592    poly t;
    593    int e;
    594 
    595    if(size(i[m])==1)
    596    {
    597       t=var(n);
    598    }
    599    else
    600    {
    601       e=deg(lead(i[m]));
    602 
    603       if(deg(poly(e))>=0)
    604       {
    605            t=leadcoef(i[m])*e*var(n)+(i[m]-lead(i[m]))/var(n)^(e-1);
    606            i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m];
    607       }
    608       else
    609       {
    610          i[m]=i[m]/leadcoef(i[m]);
    611          t=reduce(coef(i[m],var(n))[2,e+1],prm);
    612          t=var(n)+factorize(t,1)[1];
    613       }
    614    }
    615    return(t);
    616 }
    617 proc hilfe1(ideal i,ideal prm,int m,int n)
    618 {
    619    poly t;
    620    int e;
    621    if(size(i[m])==1)
    622    {
    623       t=var(n);
    624    }
    625    else
    626    {
    627       e=deg(lead(i[m]));
    628       i[m]=i[m]/leadcoef(i[m]);
    629       t=reduce(coeffs(i[m],var(n))[1,1],prm);
    630       if(size(t)==0){return(var(n));}
    631       t=var(n)+factorize(t,1)[1];
    632    }
    633    return(t);
    634 }
    635504
    636505///////////////////////////////////////////////////////////////////////////////
    … …  
    788657   {
    789658     attrib(l[2*i-1],"isSB",1);
    790 
     659     
    791660     if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0)&&(deg(l[2*i-1][1])>0))
    792661     {
    793662       "Achtung in split";
    794 
     663       
    795664         l[2*i-1]=ideal(1);
    796665         l[2*i]=ideal(1);
    … …  
    857726     return(primary);
    858727  }
    859 
    860   j=interred(j);
     728 
     729  j=interred(j); 
    861730  attrib(j,"isSB",1);
    862731  if(vdim(j)==deg(j[1]))
    863   {
     732  {   
    864733     act=factor(j[1]);
    865734     for(@k=1;@k<=size(act[1]);@k++)
    … …  
    879748       primary[2*@k]=interred(@qh);
    880749       attrib( primary[2*@k-1],"isSB",1);
    881 
     750       
    882751       if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
    883752       {
    884753          primary[2*@k-1]=ideal(1);
    885           primary[2*@k]=ideal(1);
     754          primary[2*@k]=ideal(1);         
    886755       }
    887756     }
    … …  
    957826  }
    958827  else
    959   {
     828  { 
    960829     primary[1]=j;
    961830     if((size(#)>0)&&(act[2][1]>1))
    … …  
    976845  if(size(#)==0)
    977846  {
    978 
     847   
    979848     primary=splitPrimary(primary,ser,@wr,act);
    980 
     849     
    981850  }
    982851
    … …  
    1001870    }
    1002871  }
    1003 
     872   
    1004873  @k=0;
    1005874  ideal keep;
    … …  
    1027896                    jmap2[zz]=primary[2*@k-1][@n];
    1028897                    @qht[@n]=var(zz);
    1029 
     898                     
    1030899                }
    1031900             }
    … …  
    1039908       phi1=@P,jmap1;
    1040909       phi=@P,jmap;
    1041 
     910 
    1042911       for(@n=1;@n<=nva;@n++)
    1043912       {
    … …  
    1048917
    1049918       @qh=phi(@qht);
    1050 
     919       
    1051920       if(npars(@P)>0)
    1052921       {
    … …  
    1076945       kill @Phelp1;
    1077946       @qh=clearSB(@qh);
    1078        attrib(@qh,"isSB",1);
     947       attrib(@qh,"isSB",1);       
    1079948       ser1=phi1(ser);
    1080949       @lh=zero_decomp (@qh,phi(ser1),@wr);
    1081950//       @lh=zero_decomp (@qh,psi(ser),@wr);
    1082 
    1083 
     951       
     952       
    1084953       kill lres0;
    1085954       list lres0;
    … …  
    16761545"
    16771546{
    1678    #[1]=1;
    16791547   def @P=basering;
    16801548   list qr=simplifyIdeal(i);
    16811549   map phi=@P,qr[2];
    16821550   i=qr[1];
    1683 
     1551   
    16841552   execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),("
    16851553             +ordstr(basering)+");";
    … …  
    17931661   poly  q = z4+2;
    17941662   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
    1795    LIB "primaryDecomposition.lib";
     1663   LIB "primdec.lib";
    17961664   list pr= minAssPrimes(i);
    17971665   pr;
    17981666   pr= minAssPrimes(i,1);
    17991667}
    1800 
    1801 ///////////////////////////////////////////////////////////////////////////////
    18021668
    18031669proc union(list li)
    … …  
    18851751   }
    18861752   return(k);
    1887 }
    1888 ///////////////////////////////////////////////////////////////////////////////
    1889 proc equiRadical(ideal i)
    1890 {
    1891    return(radical(i,1));
    18921753}
    18931754
    … …  
    19291790      }
    19301791  }
    1931 
     1792 
    19321793  homo=homog(i);
    1933 
     1794 
    19341795  if(homo==1)
    19351796  {
    … …  
    19781839  option(redSB);
    19791840
    1980   ideal ser=fetch(@P,ser);
     1841  ideal ser=fetch(@P,ser); 
    19811842
    19821843  if(homo==1)
    … …  
    20321893    }
    20331894  }
    2034 
     1895 
    20351896  //----------------------------------------------------------------
    20361897  //j is the ring
    … …  
    20701931     return(primary);
    20711932  }
    2072 
     1933 
    20731934 //------------------------------------------------------------------
    20741935 //the zero-dimensional case
    … …  
    21241985      indep=maxIndependSet(@j);
    21251986   }
    2126 
     1987 
    21271988  ideal jkeep=@j;
    21281989
    … …  
    21482009      ideal jwork=std(imap(gnir,@j),@hilb);
    21492010    }
    2150 
     2011   
    21512012  }
    21522013  else
    … …  
    21592020  di=dim(jwork);
    21602021  keepdi=di;
    2161 
     2022 
    21622023  setring gnir;
    21632024  for(@m=1;@m<=size(indep);@m++)
    … …  
    21752036        attrib(@j,"isSB",1);
    21762037        ideal ser=fetch(gnir,ser);
    2177 
     2038       
    21782039     }
    21792040     else
    … …  
    21922053        }
    21932054        ideal ser=phi(ser);
    2194 
    2195      }
    2196      option(noredSB);
     2055       
     2056     }
     2057     option(noredSB);     
    21972058     if((deg(@j[1])==0)||(dim(@j)<jdim))
    21982059     {
    … …  
    22502111        }
    22512112        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
    2252         option(redSB);
     2113        option(redSB);       
    22532114        list uprimary= zero_decomp(@j,ser,@wr);
    22542115        option(noredSB);
    … …  
    23122173     //mentioned above is really computed
    23132174     for(@n=@n3/2+1;@n<=@n2/2;@n++)
    2314      {
     2175     {       
    23152176        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
    23162177        {
    … …  
    23272188        }
    23282189     }
    2329 
     2190   
    23302191     if(size(@h)>0)
    23312192     {
    … …  
    23372198        if(@wr!=1)
    23382199        {
    2339           @q=minSat(jwork,@h)[2];
     2200          @q=minSat(jwork,@h)[2];   
    23402201        }
    23412202        else
    … …  
    23572218        }
    23582219        jwork=std(jwork,@q);
    2359         keepdi=dim(jwork);
     2220        keepdi=dim(jwork);       
    23602221        if(keepdi<di)
    23612222        {
    … …  
    23822243  {
    23832244    keepdi=di-1;
    2384   }
     2245  } 
    23852246  //---------------------------------------------------------------
    23862247  //notice that j=sat(j,gh) intersected with (j,gh^n)
    … …  
    24132274
    24142275        if(size(ser)>0)
    2415         {
    2416            ser=intersect(htest,ser);
     2276        { 
     2277           ser=intersect(htest,ser); 
    24172278        }
    24182279        else
    24192280        {
    24202281          ser=htest;
    2421         }
     2282        } 
    24222283        setring gnir;
    24232284        ser=imap(@Phelp,ser);
    24242285     }
    24252286     if(size(reduce(ser,peek,1))!=0)
    2426      {
     2287     {       
    24272288        for(@m=1;@m<=size(restindep);@m++)
    24282289        {
    24292290         // if(restindep[@m][3]>=keepdi)
    2430          // {
     2291         // { 
    24312292           isat=0;
    24322293           @n2=0;
    24332294           option(redSB);
    2434 
     2295           
    24352296           if(restindep[@m][1]==varstr(basering))
    24362297           //this is the good case, nothing to do, just to have the same notations
    … …  
    24592320           }
    24602321           option(noredSB);
    2461 
     2322           
    24622323           for (lauf=1;lauf<=size(@j);lauf++)
    24632324           {
    … …  
    25082369           }
    25092370           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
    2510 
     2371           
    25112372           option(redSB);
    25122373           list uprimary= zero_decomp(@j,ser,@wr);
    25132374           option(noredSB);
    2514 
    2515 
     2375           
     2376           
    25162377           //we need the intersection of the ideals in the list quprimary with the
    25172378           //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
    … …  
    25482409           @n2=size(quprimary);
    25492410           @n3=@n2;
    2550 
     2411           
    25512412           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
    25522413           {
    … …  
    25612422              }
    25622423           }
    2563 
    2564 
     2424           
     2425           
    25652426           //here the intersection with the polynomialring
    25662427           //mentioned above is really computed
    … …  
    26012462              ser=imap(@Phelp,ser);
    26022463           }
    2603 
     2464           
    26042465         // }
    2605         }
     2466        }             
    26062467        if(size(reduce(ser,peek,1))!=0)
    26072468        {
    … …  
    26232484        }
    26242485     }
    2625 
     2486     
    26262487   }
    26272488  //------------------------------------------------------------
    … …  
    26292490  //the final result: primary
    26302491  //------------------------------------------------------------
    2631 
     2492 
    26322493  setring @P;
    26332494  primary=imap(gnir,quprimary);
    … …  
    26492510
    26502511///////////////////////////////////////////////////////////////////////////////
    2651 proc primdecGTZ(ideal i)
    2652 {
    2653    return(decomp(i));
    2654 }
    2655 ///////////////////////////////////////////////////////////////////////////////
    2656 proc primdecSY(ideal i,int j)
    2657 {
    2658    return(prim_dec(i,j));
    2659 }
    2660 ///////////////////////////////////////////////////////////////////////////////
    2661 
    2662 proc radical(ideal i,list #)
    2663 {
    2664    def @P=basering;
    2665    int j,il;
    2666    if(size(#)>0)
    2667    {
    2668      il=#[1];
    2669    }
    2670    ideal re=1;
    2671    option(redSB);
    2672    list qr=simplifyIdeal(i);
    2673    map phi=@P,qr[2];
    2674    i=qr[1];
    2675 
    2676    list pr=facstd(i);
    2677 
    2678    if(size(pr)==1)
    2679    {
    2680       attrib(pr[1],"isSB",1);
    2681       if((dim(pr[1])==0)&&(homog(pr[1])==1))
    2682       {
    2683          ideal @res=maxideal(1);
    2684          return(phi(@res));
    2685       }
    2686       if(dim(pr[1])>1)
    2687       {
    2688          execute "ring gnir = ("+charstr(basering)+"),
    2689                               ("+varstr(basering)+"),(C,lp);";
    2690          ideal i=fetch(@P,i);
    2691          list @pr=facstd(i);
    2692          setring @P;
    2693          pr=fetch(gnir,@pr);
    2694       }
    2695    }
    2696    option(noredSB);
    2697    int s=size(pr);
    2698 
    2699    if(s==1)
    2700    {
    2701      i=radicalEHV(i,ideal(1),il);
    2702      return(phi(i));
    2703    }
    2704    intvec pos;
    2705    pos[s]=0;
    2706    if(il==1)
    2707    {
    2708      int ndim,k;
    2709      attrib(pr[1],"isSB",1);
    2710      int odim=dim(pr[1]);
    2711      int count=1;
    2712 
    2713      for(j=2;j<=s;j++)
    2714      {
    2715         attrib(pr[j],"isSB",1);
    2716         ndim=dim(pr[j]);
    2717         if(ndim>odim)
    2718         {
    2719            for(k=count;k<=j-1;k++)
    2720            {
    2721               pos[k]=1;
    2722            }
    2723            count=j;
    2724            odim=ndim;
    2725         }
    2726         if(ndim<odim)
    2727         {
    2728            pos[j]=1;
    2729         }
    2730      }
    2731    }
    2732    for(j=1;j<=s;j++)
    2733    {
    2734       if(pos[j]==0)
    2735       {
    2736          re=intersect(re,radicalEHV(pr[s+1-j],re,il));
    2737       }
    2738    }
    2739    return(phi(re));
    2740 }
    27412512
    27422513proc intersect1(ideal i,ideal j)
    … …  
    27512522   return(phi(j));
    27522523}
    2753 
     2524 
    27542525
    27552526///////////////////////////////////////////////////////////////////////////////
    27562527proc radicalKL (list m,ideal ser,list #)
    2757 "USAGE:   decomp(i); i ideal  (for primary decomposition)   (resp.
    2758          decomp(i,1);        (for the minimal associated primes) )
    2759 RETURN:  list = list of primary ideals and their associated primes
    2760          (at even positions in the list)
    2761          (resp. a list of the minimal associated primes)
    2762 NOTE:    Algorithm of Gianni, Traeger, Zacharias
    2763 EXAMPLE: example decomp; shows an example
    2764 "
    27652528{
    27662529   ideal i=m[2];
    … …  
    27732536      return(ideal(0));
    27742537   }
    2775 
     2538   
    27762539   def  @P = basering;
    27772540   list indep,allindep,restindep,fett,@mu;
    … …  
    28332596
    28342597      return(ideal(@p));
    2835    }
     2598   }   
    28362599   //------------------------------------------------------------------
    28372600   //the case of a complete intersection
    … …  
    28492612   if (jdim==0)
    28502613   {
    2851       @j1=finduni(@j);
     2614      @j1=finduni(@j); 
    28522615      for(@k=1;@k<=size(@j1);@k++)
    28532616      {
    … …  
    28642627      return(@j);
    28652628   }
    2866 
     2629 
    28672630   //------------------------------------------------------------------
    28682631   //search for a maximal independent set indep,i.e.
    … …  
    28732636
    28742637   indep=maxIndependSet(@j);
    2875 
     2638 
    28762639   di=dim(@j);
    28772640
    … …  
    30022765      {
    30032766         fac=ideal(0);
    3004          for(lauf=1;lauf<=ncols(@h);lauf++)
     2767         for(lauf=1;lauf<=ncols(@h);lauf++)   
    30052768         {
    30062769            if(deg(@h[lauf])>0)
    … …  
    30162779         }
    30172780
    3018 
     2781       
    30192782         @mu=mstd(quotient(@j+ideal(@q),rad));
    30202783         @j=@mu[1];
    30212784         attrib(@j,"isSB",1);
    3022 
     2785       
    30232786      }
    30242787      if((deg(rad[1])>0)&&(deg(collectrad[1])>0))
    … …  
    30352798
    30362799      te=simplify(reduce(te*rad,@j),2);
    3037 
     2800 
    30382801      if((dim(@j)<di)||(size(te)==0))
    30392802      {
    … …  
    30492812   {
    30502813      return(rad);
    3051    }
     2814   }   
    30522815  // rad=intersect(rad,radicalKL(@mu,rad,@wr));
    30532816   rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
    … …  
    30582821}
    30592822
    3060 
    3061 example
    3062 { "EXAMPLE:"; echo = 2;
    3063 }
    3064 
     2823///////////////////////////////////////////////////////////////////////////////
    30652824
    30662825proc radicalEHV(ideal i,ideal re,list #)
    … …  
    30722831      il=#[1];
    30732832   }
    3074 
     2833   
    30752834   option(redSB);
    30762835   list m=mstd(i);
    … …  
    30892848        return(quotient(I,J));
    30902849      }
    3091 
     2850   
    30922851      for(l=1;l<=cod;l++)
    30932852      {
    … …  
    31072866         if(il==1)
    31082867         {
    3109 
     2868     
    31102869            return(radI1);
    31112870         }
    … …  
    31172876            return(radI1);
    31182877         }
    3119          return(intersect(radI1,radicalEHV(I2,re,il)));
     2878         return(intersect(radI1,radicalEHV(I2,re,il)));       
    31202879      }
    31212880   }
    31222881   return(radicalKL(m,re,il));
    31232882}
     2883///////////////////////////////////////////////////////////////////////////////
    31242884
    31252885proc Ann(module M)
    … …  
    31292889  return(ann);
    31302890}
     2891///////////////////////////////////////////////////////////////////////////////
    31312892
    31322893//computes the equidimensional part of the ideal i of codimension e
    … …  
    31452906  }
    31462907  return(ann);
    3147 }
    3148 
    3149 //computes all equidimensional parts of the ideal i
    3150 proc prepareAss(ideal i)
    3151 {
    3152   ideal j=std(i);
    3153   int cod=nvars(basering)-dim(j);
    3154   int e;
    3155   list er;
    3156   ideal ann;
    3157   if(homog(i)==1)
    3158   {
    3159      list re=sres(i,0);                   //the resolution
    3160      re=minres(re);                       //minimized resolution
    3161   }
    3162   else
    3163   {
    3164      list re=mres(i,0);             //fehler in sres
    3165   }
    3166   for(e=cod;e<=nvars(basering);e++)
    3167   {
    3168      ann=AnnExt_R(e,re);
    3169 
    3170      if(nvars(basering)-dim(std(ann))==e)
    3171      {
    3172         er[size(er)+1]=equiRadical(ann);
    3173      }
    3174   }
    3175   return(er);
    3176 }
     2908}
     2909 
     2910///////////////////////////////////////////////////////////////////////////////
    31772911
    31782912//computes the annihilator of Ext^n(R/i,R) with given resolution re
    … …  
    31922926     ideal ann=Ann(transpose(re[n]));
    31932927  }
    3194   return(ann);
    3195 }
     2928  return(ann);           
     2929}
     2930///////////////////////////////////////////////////////////////////////////////
    31962931
    31972932proc quotient1(module a,module b)
    … …  
    32042939      re=intersect1(re,quotient(a,module(b[i])));
    32052940   }
    3206    return(re);
    3207 }
    3208 
    3209 
     2941   return(re);   
     2942}
     2943
     2944
     2945///////////////////////////////////////////////////////////////////////////////
    32102946
    32112947proc analyze(list pr)
    3212 {
     2948{ 
    32132949   int ii,jj;
    32142950   for(ii=1;ii<=size(pr)/2;ii++)
    … …  
    32302966         }
    32312967      }
    3232    }
    3233 }
    3234 
     2968   }   
     2969}
     2970
     2971///////////////////////////////////////////////////////////////////////////////
     2972//
     2973//                  Shimoyama-Yokoyama
     2974//
     2975///////////////////////////////////////////////////////////////////////////////
    32352976
    32362977proc simplifyIdeal(ideal i)
    32372978{
    32382979  def r=basering;
    3239 
     2980 
    32402981  int j,k;
    32412982  map phi;
    32422983  poly p;
    3243 
     2984 
    32442985  ideal iwork=i;
    32452986  ideal imap1=maxideal(1);
    32462987  ideal imap2=maxideal(1);
    3247 
     2988 
    32482989
    32492990  for(j=1;j<=nvars(basering);j++)
    … …  
    32603001        iwork[k]=var(j);
    32613002        imap1=maxideal(1);
    3262         imap2[j]=-p;
     3003        imap2[j]=-p;       
    32633004        break;
    32643005      }
    … …  
    32683009}
    32693010
    3270 
     3011       
    32713012///////////////////////////////////////////////////////
    32723013// ini_mod
    … …  
    38663607            }
    38673608        }
    3868         dimSP=dim(SP);
     3609        dimSP=dim(SP);
    38693610        for(j=size(m);j>=1; j--)
    38703611        {
    … …  
    40633804      {
    40643805        f=polys[k];
    4065         degf=deg(f);
     3806        degf=deg(f);
    40663807      }
    40673808    }
    … …  
    42734014   return(0);
    42744015}
    4275 
    4276 
    4277  proc quotient2(module a,module b)
     4016       
     4017       
     4018proc quotient2(module a,module b)
    42784019{
    42794020  string s="ring @newr=("+charstr(basering)+"),(@t,"+varstr(basering)+"),dp;";
    … …  
    42914032  setring @newP;
    42924033 ideal re=imap(@newr,re);
    4293  return(re);
    4294 }
    4295 //Im homogenen Fall system("LaScala",i) verwenden
     4034 return(re);   
     4035}
     4036///////////////////////////////////////////////////////////////////////////////
     4037
     4038proc convList(list l)
     4039{
     4040   int i;
     4041   list re,he;
     4042   for(i=1;i<=size(l)/2;i++)
     4043   {
     4044      he=l[2*i-1],l[2*i];
     4045      re[i]=he;
     4046   }
     4047   return(re);   
     4048}
     4049///////////////////////////////////////////////////////////////////////////////
     4050
     4051proc reconvList(list l)
     4052{
     4053   int i;
     4054   list re;
     4055   for(i=1;i<=size(l);i++)
     4056   {
     4057      re[2*i-1]=l[i][1];
     4058      re[2*i]=l[i][2];
     4059   }
     4060   return(re);   
     4061}
     4062
     4063///////////////////////////////////////////////////////////////////////////////
     4064//
     4065//     The main procedures
     4066//
     4067///////////////////////////////////////////////////////////////////////////////
     4068
     4069proc primdecGTZ(ideal i)
     4070"USAGE:  primdecGTZ(i); i ideal
     4071RETURN:  list = list of primary ideals
     4072                and their associated primes
     4073NOTE:    Algorithm of Gianni, Traeger, Zacharias
     4074EXAMPLE: example primdecGTZ; shows an example
     4075"
     4076{
     4077   return(convList(decomp(i)));
     4078}
     4079example
     4080{ "EXAMPLE:";  echo = 2;
     4081   ring  r = 32003,(x,y,z),lp;
     4082   poly  p = z2+1;
     4083   poly  q = z4+2;
     4084   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4085   list pr= primdecGTZ(i);
     4086   pr;
     4087}
     4088///////////////////////////////////////////////////////////////////////////////
     4089
     4090proc primdecSY(ideal i)
     4091"USAGE:  primdecSY(i); i ideal
     4092RETURN:  list = list of primary ideals
     4093                and their associated primes
     4094NOTE:    Algorithm of Shimoyama-Yokoyama
     4095EXAMPLE: example primdecSY; shows an example
     4096"
     4097{
     4098   return(prim_dec(i,1));
     4099}
     4100example
     4101{ "EXAMPLE:";  echo = 2;
     4102   ring  r = 32003,(x,y,z),lp;
     4103   poly  p = z2+1;
     4104   poly  q = z4+2;
     4105   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4106   list pr= primdecSY(i);
     4107   pr;
     4108}
     4109///////////////////////////////////////////////////////////////////////////////
     4110proc minAssGTZ(ideal i)
     4111"USAGE:  minAssGTZ(i); i ideal
     4112RETURN:  list = the minimal associated prime ideals of i
     4113NOTE:
     4114EXAMPLE: example minAssGTZ; shows an example
     4115"
     4116{
     4117   return(minAssPrimes(i,1));
     4118}
     4119example
     4120{ "EXAMPLE:";  echo = 2;
     4121   ring  r = 32003,(x,y,z),dp;
     4122   poly  p = z2+1;
     4123   poly  q = z4+2;
     4124   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4125   list pr= minAssGTZ(i);
     4126   pr;
     4127}
     4128
     4129///////////////////////////////////////////////////////////////////////////////
     4130proc minAssChar(ideal i)
     4131"USAGE:  minAssChar(i); i ideal
     4132RETURN:  list = the minimal associated prime ideals of i
     4133NOTE:
     4134EXAMPLE: example minAssChar; shows an example
     4135"
     4136
     4137{
     4138   return(min_ass_prim_charsets(i,1));
     4139}
     4140example
     4141{ "EXAMPLE:";  echo = 2;
     4142   ring  r = 32003,(x,y,z),dp;
     4143   poly  p = z2+1;
     4144   poly  q = z4+2;
     4145   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4146   list pr= minAssChar(i);
     4147   pr;
     4148}
     4149
     4150///////////////////////////////////////////////////////////////////////////////
     4151proc equiRadical(ideal i)
     4152"USAGE:  equiRadical(i); i ideal
     4153RETURN:  ideal = the intersection of associated primes 
     4154         of i of dimension of i
     4155NOTE:
     4156EXAMPLE: example equiRadical; shows an example
     4157"
     4158{
     4159   return(radical(i,1));
     4160}
     4161example
     4162{ "EXAMPLE:";  echo = 2;
     4163   ring  r = 32003,(x,y,z),dp;
     4164   poly  p = z2+1;
     4165   poly  q = z4+2;
     4166   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4167   ideal pr= equiRadical(i);
     4168   pr;
     4169}
     4170///////////////////////////////////////////////////////////////////////////////
     4171proc radical(ideal i,list #)
     4172"USAGE:  radical(i); i ideal
     4173RETURN:  ideal = the radical of i
     4174NOTE:    a combination of the algorithms of Krick/Logar
     4175         and Eisenbud/Huneke/Vasconcelos
     4176EXAMPLE: example radical; shows an example
     4177"
     4178{
     4179   def @P=basering;
     4180   int j,il;
     4181   if(size(#)>0)
     4182   {
     4183     il=#[1];
     4184   }
     4185   ideal re=1;
     4186   option(redSB);
     4187   list qr=simplifyIdeal(i);
     4188   map phi=@P,qr[2];
     4189   i=qr[1];
     4190   
     4191   list pr=facstd(i);
     4192
     4193   if(size(pr)==1)
     4194   {
     4195      attrib(pr[1],"isSB",1);
     4196      if((dim(pr[1])==0)&&(homog(pr[1])==1))
     4197      {
     4198         ideal @res=maxideal(1);
     4199         return(phi(@res));
     4200      }
     4201      if(dim(pr[1])>1)
     4202      {
     4203         execute "ring gnir = ("+charstr(basering)+"),
     4204                              ("+varstr(basering)+"),(C,lp);";
     4205         ideal i=fetch(@P,i);
     4206         list @pr=facstd(i);
     4207         setring @P;
     4208         pr=fetch(gnir,@pr);
     4209      }
     4210   }
     4211   option(noredSB);
     4212   int s=size(pr);
     4213
     4214   if(s==1)
     4215   {
     4216     i=radicalEHV(i,ideal(1),il);
     4217     return(phi(i));
     4218   }
     4219   intvec pos;
     4220   pos[s]=0;
     4221   if(il==1)
     4222   {
     4223     int ndim,k;
     4224     attrib(pr[1],"isSB",1);
     4225     int odim=dim(pr[1]);
     4226     int count=1;
     4227   
     4228     for(j=2;j<=s;j++)
     4229     {
     4230        attrib(pr[j],"isSB",1);
     4231        ndim=dim(pr[j]);
     4232        if(ndim>odim)
     4233        {
     4234           for(k=count;k<=j-1;k++)
     4235           {
     4236              pos[k]=1;
     4237           }
     4238           count=j;
     4239           odim=ndim;
     4240        }
     4241        if(ndim<odim)
     4242        {
     4243           pos[j]=1;
     4244        }
     4245     }
     4246   }
     4247   for(j=1;j<=s;j++)
     4248   {
     4249      if(pos[j]==0)
     4250      {
     4251         re=intersect(re,radicalEHV(pr[s+1-j],re,il));
     4252      }
     4253   }
     4254   return(phi(re));
     4255}
     4256example
     4257{ "EXAMPLE:";  echo = 2;
     4258   ring  r = 32003,(x,y,z),dp;
     4259   poly  p = z2+1;
     4260   poly  q = z4+2;
     4261   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4262   ideal pr= radical(i);
     4263   pr;
     4264}
     4265///////////////////////////////////////////////////////////////////////////////
     4266proc prepareAss(ideal i)
     4267"USAGE:  prepareAss(i); i ideal
     4268RETURN:  list = the radicals of the equidimensional parts of i
     4269NOTE:    algorithm of Eisenbud,Huneke and Vasconcelos
     4270EXAMPLE: example prepareAss; shows an example
     4271"
     4272{
     4273  ideal j=std(i);
     4274  int cod=nvars(basering)-dim(j);
     4275  int e;
     4276  list er;
     4277  ideal ann;
     4278  if(homog(i)==1)
     4279  {
     4280     list re=sres(i,0);                   //the resolution
     4281     re=minres(re);                       //minimized resolution
     4282  }
     4283  else
     4284  {
     4285     list re=mres(i,0);             
     4286  } 
     4287  for(e=cod;e<=nvars(basering);e++)
     4288  {
     4289     ann=AnnExt_R(e,re);
     4290
     4291     if(nvars(basering)-dim(std(ann))==e)
     4292     {
     4293        er[size(er)+1]=equiRadical(ann);
     4294     }
     4295  }
     4296  return(er);
     4297} 
     4298example
     4299{ "EXAMPLE:";  echo = 2;
     4300   ring  r = 32003,(x,y,z),dp;
     4301   poly  p = z2+1;
     4302   poly  q = z4+2;
     4303   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4304   list pr= prepareAss(i);
     4305   pr;
     4306}
     4307
     4308proc testPrimary(list pr, ideal k)
     4309"USAGE:   testPrimary(pr,k) pr=primdecGTZ(k) or primdecSY(k)
     4310RETURN:   int = 1, if the intersection of the primary ideals
     4311                in pr is k, 0 if not
     4312NOTE:
     4313EXAMPLE: example testPrimary ; shows an example
     4314"
     4315{
     4316   int i;
     4317   pr=reconvList(pr);
     4318   ideal j=pr[1];
     4319   for (i=2;i<=size(pr)/2;i++)
     4320   {
     4321       j=intersect(j,pr[2*i-1]);
     4322   }
     4323   return(idealsEqual(j,k));
     4324}
     4325example
     4326{ "EXAMPLE:";  echo = 2;
     4327   ring  r = 32003,(x,y,z),dp;
     4328   poly  p = z2+1;
     4329   poly  q = z4+2;
     4330   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     4331   list pr= primdecGTZ(i);
     4332   testPrimary(pr,i);
     4333}
     4334
Note: See TracChangeset for help on using the changeset viewer.