Changeset 80654d in git


Ignore:
Timestamp:
Jan 8, 2001, 2:27:30 PM (23 years ago)
Author:
Gerhard Pfister <pfister@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
03f11fad12098603bb486b2c5eba652d4d581332
Parents:
c612a88504b2a9145a21af5d28ec4c707f75213f
Message:
lgorithm of Kemper for radical implemented


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/primdec.lib

    rc612a8 r80654d  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: primdec.lib,v 1.82 2000-12-31 02:02:05 greuel Exp $";
     2version="$Id: primdec.lib,v 1.83 2001-01-08 13:27:30 pfister Exp $";
    33category="Commutative Algebra";
    44info="
     
    1414 Shimoyama/Yokoyama was written by Wolfram Decker and Hans Schoenemann.   
    1515 These procedures are implemented to be used in characteristic 0.
    16  @*They also work in positive characteristic >> 0.
    17  @*In small characteristic and for algebraic extensions, primdecGTZ,
    18  minAssGTZ, radical and equiRadical may not terminate and primdecSY and
     16 They also work in positive characteristic >> 0.
     17 In small characteristic and for algebraic extensions, primdecGTZ and
     18 minAssGTZ may not terminate and primdecSY and
    1919 minAssChar may not give a complete decomposition.
    2020
     
    27982798
    27992799///////////////////////////////////////////////////////////////////////////////
    2800 proc radicalKL (list m,ideal ser,list #)
    2801 {
    2802    ideal i=m[2];
     2800proc powerCoeffs(poly f,int e)
     2801//computes a polynomial with the same monomials as f but coefficients
     2802//the p^e th power of the coefficients of f
     2803{
     2804   int i;
     2805   poly g;
     2806   int ex=char(basering)^e;
     2807   for(i=1;i<=size(f);i++)
     2808   {
     2809      g=g+leadcoef(f[i])^ex*leadmonom(f[i]);
     2810   }
     2811   return(g);
     2812}
     2813///////////////////////////////////////////////////////////////////////////////
     2814
     2815proc sep(poly f,int i, list #)
     2816"USAGE:  input: a polynomial f depending on the i-th variable and optional
     2817         an integer k considering the polynomial f defined over Fp(t1,...,tm)
     2818         as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
     2819 RETURN: the separabel part of f as polynomial in Fp(t1,...,tm)
     2820        and an integer k to indicate that f should be considerd
     2821        as polynomial over Fp(t(1)^(p^-k),...,t(m)^(p^-k))
     2822 EXAMPLE: example sep; shows an example
     2823{
     2824   def R=basering;
     2825   int k;
     2826   if(size(#)>0){k=#[1];}
     2827
     2828   poly h=gcd(f,diff(f,var(i)));
     2829   poly g1=lift(h,f)[1][1];    //  f/h
     2830   poly h1;
     2831
     2832   while(h!=h1)
     2833   {
     2834      h1=h;
     2835      h=gcd(h,diff(h,var(i)));
     2836   }
     2837
     2838   if(deg(h1)==0){return(list(g1,k));} //in characteristic 0 we return here
     2839
     2840   k++;
     2841
     2842   ideal ma=maxideal(1);
     2843   ma[i]=var(i)^char(R);
     2844   map phi=R,ma;
     2845   ideal hh=h;    //this is technical because preimage works only for ideals
     2846
     2847   poly u=preimage(R,phi,hh)[1]; //h=u(x(i)^p)
     2848
     2849   list g2=sep(u,i,k);           //we consider u(t(1)^(p^-1),...,t(m)^(p^-1))
     2850   g1=powerCoeffs(g1,g2[2]-k+1); //to have g1 over the same field as g2[1]
     2851
     2852   list g3=sep(g1*g2[1],i,g2[2]);
     2853   return(g3);
     2854}
     2855example
     2856{ "EXAMPLE:"; echo = 2;
     2857   ring R=(5,t,s),(x,y,z),dp;
     2858   poly f=(x^25-t*x^5+t)*(x^3+s);
     2859   sep(f,1);
     2860}
     2861
     2862///////////////////////////////////////////////////////////////////////////////
     2863proc zeroRad(ideal I)
     2864"USAGE:  zeroRad(I) , I a zero-dimensional ideal
     2865 RETURN: the radical of I
     2866 NOTE:  Algorithm of Kemper
     2867 EXAMPLE: example zeroRad; shows an example
     2868{
     2869   
     2870   //I needs to be a reduced standard basis
     2871   def R=basering;
     2872   int m=npars(R);
     2873   int n=nvars(R);
     2874   int p=char(R); 
     2875   int i,k;
     2876   list l;
     2877
     2878   option(redSB); 
     2879   ideal F=finduni(I);//F[i] generates I intersected with K[var(i)]
     2880   option(noredSB);
     2881
     2882   for(i=1;i<=n;i++)
     2883   {
     2884      l[i]=sep(F[i],i);
     2885      F[i]=l[i][1];
     2886      if(l[i][2]>k){k=l[i][2];}  //computation of the maximal k
     2887   }
     2888
     2889   if(k==0){return(I+F);}        //the separable case
     2890
     2891   for(i=1;i<=n;i++)             //consider all polynomials over
     2892   {                             //Fp(t(1)^(p^-k),...,t(m)^(p^-k))
     2893      F[i]=powerCoeffs(F[i],k-l[i][2]);
     2894   }
     2895
     2896   string cR="ring @R="+string(p)+",("+parstr(R)+","+varstr(R)+"),dp;";
     2897   execute(cR);
     2898   ideal F=imap(R,F);
     2899   string nR="ring @S="+string(p)+",(y(1..m),"+varstr(R)+","+parstr(R)+"),dp;";
     2900   execute(nR);
     2901
     2902   ideal G=fetch(@R,F);    //G[i](t(1)^(p^-k),...,t(m)^(p^-k),x(i))=sep(F[i])
     2903   ideal I=imap(R,I);
     2904   ideal J=I+G;
     2905   k=p^k;
     2906   for(i=1;i<=m;i++)
     2907   {
     2908      J=J,var(i)^k-var(m+n+i);
     2909   }
     2910   J=eliminate(J,y(1..m));
     2911
     2912   setring R;
     2913   ideal J=imap(@S,J);
     2914   return(J);
     2915}
     2916example
     2917{ "EXAMPLE:"; echo = 2;
     2918   ring R=(5,t),(x,y),dp;
     2919   ideal I=x^5-t,y^5-t;
     2920   zeroradp(I);
     2921}
     2922
     2923///////////////////////////////////////////////////////////////////////////////
     2924proc radicalKL (ideal i,ideal ser,list #)
     2925{
     2926  // i needs to be a reduced standard basis
    28032927  //--------------------------------------------------------------------------
    28042928  //i is the zero-ideal
     
    28122936   def  @P = basering;
    28132937   list indep,allindep,restindep,fett,@mu;
    2814    intvec @vh,isat;
     2938   intvec @vh,isat,@w,@hilb;
    28152939   int @wr,@k,@n,@m,@n1,@n2,@n3,lauf,di;
    28162940   ideal @j,@j1,fac,@h,collectrad,htrad,lsau;
     
    28192943
    28202944   poly @p,@q;
    2821    string @va,quotring;
     2945   string @va,quotring,@ri;
    28222946   int  homo=homog(i);
    2823    if((find(ordstr(basering),"w")!=0)||(find(ordstr(basering),"W")!=0)||(find(ordstr(basering),"a")!=0))
    2824    {
    2825       homo=0;
    2826    }
     2947   attrib(i,"isSB",1);
    28272948   if(size(#)>0)
    28282949   {
    28292950      @wr=#[1];
    28302951   }
    2831    @j=m[1];
    2832    @j1=m[2];
    2833    int jdim=dim(@j);
    2834    if(size(reduce(ser,@j,1))==0)
    2835    {
    2836       return(ser);
    2837    }
     2952   int jdim=dim(i);
     2953
    28382954   if(homo==1)
    28392955   {
    28402956      if(jdim==0)
    28412957      {
    2842          option(noredSB);
    28432958         return(maxideal(1));
    28442959      }
    2845       intvec @hilb=hilb(@j,1);
    2846    }
    2847 
    2848 
     2960      for(@n=1;@n<=nvars(basering);@n++)
     2961      {
     2962         @w[@n]=ord(var(@n));
     2963      }
     2964      @hilb=hilb(i,1,@w);     
     2965   }
    28492966  //---------------------------------------------------------------------------
    28502967  //j is the ring
     
    28532970   if (jdim==-1)
    28542971   {
    2855       option(noredSB);
    2856       return(ideal(0));
    2857    }
    2858 
    2859   //---------------------------------------------------------------------------
    2860   //  the case of one variable
    2861   //---------------------------------------------------------------------------
    2862 
    2863    if(nvars(basering)==1)
    2864    {
    2865       fac=factorize(@j[1],1);
    2866       @p=1;
    2867       for(@k=1;@k<=size(fac);@k++)
    2868       {
    2869          @p=@p*fac[@k];
    2870       }
    2871       option(noredSB);
    2872 
    2873       return(ideal(@p));
    2874    }
    2875   //---------------------------------------------------------------------------
    2876   //the case of a complete intersection
    2877   //---------------------------------------------------------------------------
    2878     if(jdim+size(@j1)==nvars(basering))
    2879     {
    2880       // ideal jac=minor(jacob(@j1),size(@j1));
    2881       // return(quotient(@j1,jac));
    2882     }
     2972
     2973      return(ideal(1));
     2974   }
    28832975
    28842976  //---------------------------------------------------------------------------
     
    28882980   if (jdim==0)
    28892981   {
    2890       if((npars(basering)>0)&&(char(basering)>0)&&(char(basering)<=181))
    2891       {
    2892          "        !  Warning  !  ";
    2893          "    The field is not perfect.";
    2894          "    The result may be wrong.";
    2895       }
    2896       @j1=finduni(@j);
    2897       for(@k=1;@k<=size(@j1);@k++)
    2898       {
    2899          fac=factorize(cleardenom(@j1[@k]),1);
    2900          @p=fac[1];
    2901          for(@n=2;@n<=size(fac);@n++)
    2902          {
    2903             @p=@p*fac[@n];
    2904          }
    2905          @j=@j,@p;
    2906       }
    2907       @j=std(@j);
    2908       option(noredSB);
    2909       return(@j);
     2982      return(zeroRad(i));
    29102983   }
    29112984
     
    29172990   //-------------------------------------------------------------------------
    29182991
    2919    indep=maxIndependSet(@j);
    2920 
    2921    di=dim(@j);
     2992   indep=maxIndependSet(i);
    29222993
    29232994   for(@m=1;@m<=size(indep);@m++)
     
    29293000        execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
    29303001                              +ordstr(basering)+");");
    2931          ideal @j=fetch(@P,@j);
     3002         ideal @j=fetch(@P,i);
    29323003         attrib(@j,"isSB",1);
    29333004      }
     
    29403011         if(homo==1)
    29413012         {
    2942             ideal @j=std(phi(@j),@hilb);
     3013            ideal @j=std(phi(i),@hilb,@w);
    29433014         }
    29443015         else
    29453016         {
    2946            ideal @j=groebner(phi(@j));
     3017           ideal @j=groebner(phi(i));
    29473018         }
    29483019      }
     
    29903061
    29913062      @j=clearSB(@j,fett);
    2992       attrib(@j,"isSB",1);
    29933063
    29943064     //we need later ggt(h[1],...)=gh for saturation
     
    30003070            @h[@n]=leadcoef(@j[@n]);
    30013071         }
    3002         //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
     3072     
    30033073         option(redSB);
    3004 
    3005          @j=interred(@j);
    3006 
     3074         @j=interred(@j);  //to obtain a reduced standardbasis
    30073075         attrib(@j,"isSB",1);
    3008          list  @mo=@j,@j;
    3009          ideal zero_rad= radicalKL(@mo,ideal(1));
     3076         option(noredSB);
     3077         ideal zero_rad= zeroRad(@j);
    30103078      }
    30113079      else
     
    30213089     //the quotientring: this is coded in saturn
    30223090
     3091     
     3092       zero_rad=std(zero_rad);
     3093     
    30233094      ideal hpl;
    30243095
     
    30453116
    30463117      collectrad=sat2(collectrad,lsau)[1];
    3047 
    30483118      if(deg(@h[1])>=0)
    30493119      {
     
    30623132            @q=@q*fac[lauf];
    30633133         }
    3064 
    3065 
    3066          @mu=mstd(quotient(@j+ideal(@q),rad));
    3067          @j=@mu[1];
    3068          attrib(@j,"isSB",1);
     3134         option(returnSB);
     3135         option(redSB);
     3136         i=quotient(i+ideal(@q),rad);
     3137         attrib(i,"isSB",1);
     3138         option(noreturnSB);
     3139         option(noredSB);
    30693140
    30703141      }
     
    30803151         }
    30813152      }
    3082 
    3083       te=simplify(reduce(te*rad,@j),2);
    3084 
    3085       if((dim(@j)<di)||(size(te)==0))
     3153      te=simplify(reduce(intersect(te,rad),i,1),2);
     3154      if((dim(i)<jdim)||(size(te)==0))
    30863155      {
    30873156         break;
     
    30893158      if(homo==1)
    30903159      {
    3091          @hilb=hilb(@j,1);
    3092       }
    3093    }
    3094 
    3095    if(((@wr==1)&&(dim(@j)<di))||(deg(@j[1])==0)||(size(te)==0))
     3160         @hilb=hilb(i,1,@w);
     3161      }
     3162   }
     3163   if(((@wr==1)&&(dim(i)<jdim))||(deg(i[1])==0)||(size(te)==0))
    30963164   {
    30973165      return(rad);
    30983166   }
    3099   // rad=intersect(rad,radicalKL(@mu,rad,@wr));
    3100    rad=intersect(rad,radicalKL(@mu,ideal(1),@wr));
    3101 
    3102 
    3103    option(noredSB);
     3167   rad=intersect(rad,radicalKL(i,ideal(1),@wr));
    31043168   return(rad);
    31053169}
     
    31643228      }
    31653229   }
    3166    return(radicalKL(m,re,il));
     3230   return(radicalKL(m[1],re,il));
    31673231}
    31683232///////////////////////////////////////////////////////////////////////////////
     
    44654529   i=qr[1];
    44664530
     4531   i=std(i);
     4532   int di=dim(i);
     4533   if(di==0)
     4534   {
     4535     i=zeroRad(i);
     4536     return(phi(i));
     4537   }
    44674538   list pr=facstd(i);
    4468    if(size(pr)==1)
    4469    {
    4470       attrib(pr[1],"isSB",1);
    4471       if((dim(pr[1])==0)&&(homog(pr[1])==1))
    4472       {
    4473          ideal @res=maxideal(1);
    4474          return(phi(@res));
    4475       }
    4476       if(dim(pr[1])>1)
    4477       {
    4478          execute("ring gnir = ("+charstr(basering)+"),
    4479                               ("+varstr(basering)+"),(C,lp);");
    4480          ideal i=fetch(@P,i);
    4481          list @pr=facstd(i);
    4482          setring @P;
    4483          pr=fetch(gnir,@pr);
    4484       }
    4485    }
    44864539   option(noredSB);
    44874540   int s=size(pr);
    4488    if(s==1)
    4489    {
    4490      i=radicalEHV(pr[1],ideal(1),il);
    4491      return(phi(i));
    4492    }
    4493    intvec pos;
    4494    pos[s]=0;
    4495    if(il==1)
    4496    {
    4497      int ndim,k;
    4498      attrib(pr[1],"isSB",1);
    4499      int odim=dim(pr[1]);
    4500      int count=1;
    4501 
    4502      for(j=2;j<=s;j++)
    4503      {
    4504         attrib(pr[j],"isSB",1);
    4505         ndim=dim(pr[j]);
    4506         if(ndim>odim)
    4507         {
    4508            for(k=count;k<=j-1;k++)
    4509            {
    4510               pos[k]=1;
    4511            }
    4512            count=j;
    4513            odim=ndim;
    4514         }
    4515         if(ndim<odim)
    4516         {
    4517            pos[j]=1;
    4518         }
    4519      }
    4520    }
    45214541   for(j=1;j<=s;j++)
    45224542   {
    4523       if(pos[s+1-j]==0)
    4524       {
    4525          re=intersect(re,radicalEHV(pr[s+1-j],re,il));
     4543      attrib(pr[s+1-j],"isSB",1);
     4544      if((size(reduce(re,pr[s+1-j],1))!=0)&&((dim(pr[s+1-j])==di)||!il))
     4545      {
     4546         re=intersect(re,radicalKL(pr[s+1-j],re,il));
    45264547      }
    45274548   }
Note: See TracChangeset for help on using the changeset viewer.