Changeset 7a68965 in git for Singular/LIB/sagbi.lib


Ignore:
Timestamp:
Jan 4, 2007, 1:49:27 PM (17 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '2a584933abf2a2d3082034c7586d38bb6de1a30a')
Children:
50d47e7bb261355ebac0cfa45c79d74f074d4a45
Parents:
dfd4fc6ea94830e53cf1a0e58c158d96a83f0b47
Message:
*hannes: format


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/sagbi.lib

    rdfd4fc r7a68965  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="$Id: sagbi.lib,v 1.5 2006-11-22 21:40:29 levandov Exp $";
     2version="$Id: sagbi.lib,v 1.6 2007-01-04 12:49:27 Singular Exp $";
    33category="Commutative Algebra";
    44info="
    … …  
    1515";
    1616
    17  LIB "algebra.lib";
    18  LIB "elim.lib";
     17LIB "algebra.lib";
     18LIB "elim.lib";
    1919
    2020///////////////////////////////////////////////////////////////////////////////
    … …  
    3232{
    3333  if(size(#)==0)
    34     {
    35       #[1]=0;
    36     }
     34  {
     35    #[1]=0;
     36  }
    3737  degBound=0;
    3838  def bsr=basering;
    … …  
    4949
    5050  if(b!=0)
    51     {
    52       id =reduce(id,groebner(0));
    53     }
     51  {
     52    id =reduce(id,groebner(0));
     53  }
    5454  int n,m=nvars(bsr),ncols(id);
    5555  int z;
    … …  
    5959
    6060  if(id==0)
    61     {
    62       if(#[1]==0)
    63         {
    64           return(P);
    65         }
    66       else
    67         {
    68           return(L);
    69         }
    70     }
    71 
     61  {
     62    if(#[1]==0)
     63    {
     64      return(P);
     65    }
     66    else
     67    {
     68      return(L);
     69    }
     70  }
    7271  else
    73 
    74     {
     72  {
    7573    //==================create anew ring with extra variables================
    7674
    77       execute("ring R1="+charstr(bsr)+",("+varstr(bsr)+",@y(1..m)),(dp(n),dp(m));");
    78       execute("minpoly=number("+mp+");");
    79       ideal id=imap(bsr,id);
    80       ideal A;
    81 
    82       for(z=1;z<=m;z++)
    83         {
    84           A[z]=lead(id[z])-@y(z);
    85         }
    86 
    87       A=groebner(A);
    88       ideal kern=nselect(A,1,n);// "kern" is the kernel of the ring map:
    89                            // R1----->bsr ;y(z)----> lead(id[z]).
    90                           //"kern" is the ideal of algebraic relations between
    91                          // lead(id[z]).
    92 
    93       export kern,A;// we export:* the ideal A  to avoid useless computations
    94                //                  betwee 2 steps in sagbi procedure.
    95               //                 *the ideal kern : some times we can get intresting
    96              //                    informations from this ideal.
    97 
    98       setring bsr;
    99       map phi=R1,vars,id;
    100 
    101       // the sagbiSPolynomials are the image by phi of the generators of kern
    102 
    103       P=simplify(phi(kern),1);
    104       if(#[1]==0)
    105         {
    106           return(P);
    107         }
    108       else
    109         {
    110           L=P,R1;
    111           kill phi,vars;
    112 
    113                  dbprint(printlevel-voice+3,"
     75    execute("ring R1="+charstr(bsr)+",("+varstr(bsr)+",@y(1..m)),(dp(n),dp(m));");
     76    execute("minpoly=number("+mp+");");
     77    ideal id=imap(bsr,id);
     78    ideal A;
     79
     80    for(z=1;z<=m;z++)
     81    {
     82      A[z]=lead(id[z])-@y(z);
     83    }
     84
     85    A=groebner(A);
     86    ideal kern=nselect(A,1,n);// "kern" is the kernel of the ring map:
     87                        // R1----->bsr ;y(z)----> lead(id[z]).
     88                        //"kern" is the ideal of algebraic relations between
     89                        // lead(id[z]).
     90
     91    export kern,A;// we export:
     92                  // * the ideal A  to avoid useless computations
     93                  //   between 2 steps in sagbi procedure.
     94                  // * the ideal kern : some times we can get intresting
     95                  //   informations from this ideal.
     96
     97    setring bsr;
     98    map phi=R1,vars,id;
     99
     100    // the sagbiSPolynomials are the image by phi of the generators of kern
     101
     102    P=simplify(phi(kern),1);
     103    if(#[1]==0)
     104    {
     105      return(P);
     106    }
     107    else
     108    {
     109      L=P,R1;
     110      kill phi,vars;
     111
     112      dbprint(printlevel-voice+3,"
    114113// 'sagbiSPoly' created a ring as 2nd element of the list.
    115114// The ring contains the ideal 'kern'  of algebraic relations between the
    … …  
    118117// e.g.:
    119118               def S = L[2]; setring S; kern;
    120         ");
    121                  return(L);
    122         }
    123 
    124     }
     119      ");
     120      return(L);
     121    }
     122  }
    125123}
    126124example
    … …  
    147145
    148146
    149   if(size(#)>0)
    150     {tt=#[1];}
    151 
    152   if(size(I)==0)
    153     {@result=groebner(J);}
     147  if(size(#)>0) {tt=#[1];}
     148
     149  if(size(I)==0) {@result=groebner(J);}
    154150
    155151  if((size(I)!=0) && (size(J)-size(I)>=1))
    156     {
    157 
    158       qring Q=I;
    159       ideal J=fetch(br,J);
    160       J=groebner(J);
    161       setring br;
    162       Res=fetch(Q,J);// Res contains the generators that we add to I
    163                       //to get the generators of std(J);
    164       @result=Res+I;
    165 
    166     }
    167 
    168   if(tt==0)
    169     {return(@result);}
    170   else
    171     {return(Res);}
     152  {
     153    qring Q=I;
     154    ideal J=fetch(br,J);
     155    J=groebner(J);
     156    setring br;
     157    Res=fetch(Q,J);// Res contains the generators that we add to I
     158                   // to get the generators of std(J);
     159    @result=Res+I;
     160  }
     161
     162  if(tt==0) { return(@result);}
     163  else      { return(Res);}
    172164}
    173165
    … …  
    187179
    188180  if(b!=0)
    189     {
    190       I=reduce(I,groebner(0));
    191       J=reduce(J,groebner(0));
    192     }
    193 
    194 
     181  {
     182    I=reduce(I,groebner(0));
     183    J=reduce(J,groebner(0));
     184  }
    195185  int n,ii,jj=nvars(br),ncols(I),ncols(J);
    196186  int z;
    … …  
    199189
    200190  if(size(J)==0)
    201     {
    202       @L =sagbiSPoly(I,1);
    203     }
     191  {
     192    @L =sagbiSPoly(I,1);
     193  }
    204194  else
    205     {
    206       ideal @sum=I+J;
    207       ideal P1;
    208       ideal P=l[1];//P is the ideal of spolynomials of I;
    209       def R=l[2];setring R;int kk=nvars(R);
    210       ideal J=fetch(br,J);
    211 
    212 
    213 
    214       //================create a new ring with extra variables==============
    215 
    216       execute("ring R1="+charstr(R)+",("+varstr(R)+",@y((ii+1)..(ii+jj))),(dp(n),dp(kk+jj-n));");
    217       ideal kern1;
    218       ideal A=fetch(R,A);
    219       attrib(A,"isSB",1);
    220       ideal J=fetch(R,J);
    221       ideal kern=fetch(R,kern);
    222       ideal A1;
    223       for(z=1;z<=jj;z++)
    224         {
    225           A1[z]=lead(J[z])-var(z+kk);
    226         }
    227       A1=A+A1;
    228       ideal @Res=std1(A1,A,1);//the generators of @Res are whose we have to add
    229                              //to A to get std(A1).
    230       A=A+@Res;
    231       kern1=nselect(@Res,1,n);
    232       kern=kern+kern1;
    233       export kern,kern1,A;
    234       setring br;
    235       map phi=R1,vars,@sum;
    236       P1=simplify(phi(kern1),1);//P1 is th ideal we add to P to get the ideal
    237                                 //of Spolynomials of @sum.
    238       P=P+P1;
    239 
    240       if (a==1)
    241         {
    242           @L=P,R1;
    243           kill phi,vars;
    244           dbprint(printlevel-voice+3,"
     195  {
     196    ideal @sum=I+J;
     197    ideal P1;
     198    ideal P=l[1];//P is the ideal of spolynomials of I;
     199    def R=l[2];setring R;int kk=nvars(R);
     200    ideal J=fetch(br,J);
     201
     202    //================create a new ring with extra variables==============
     203
     204    execute("ring R1="+charstr(R)+",("+varstr(R)+",@y((ii+1)..(ii+jj))),(dp(n),dp(kk+jj-n));");
     205    ideal kern1;
     206    ideal A=fetch(R,A);
     207    attrib(A,"isSB",1);
     208    ideal J=fetch(R,J);
     209    ideal kern=fetch(R,kern);
     210    ideal A1;
     211    for(z=1;z<=jj;z++)
     212    {
     213      A1[z]=lead(J[z])-var(z+kk);
     214    }
     215    A1=A+A1;
     216    ideal @Res=std1(A1,A,1);// the generators of @Res are whose we have to add
     217                            // to A to get std(A1).
     218    A=A+@Res;
     219    kern1=nselect(@Res,1,n);
     220    kern=kern+kern1;
     221    export kern,kern1,A;
     222    setring br;
     223    map phi=R1,vars,@sum;
     224    P1=simplify(phi(kern1),1);//P1 is th ideal we add to P to get the ideal
     225                              //of Spolynomials of @sum.
     226    P=P+P1;
     227
     228    if (a==1)
     229    {
     230      @L=P,R1;
     231      kill phi,vars;
     232      dbprint(printlevel-voice+3,"
    245233// 'Spoly1' created a ring as 2nd element of the list.
    246234// The ring contains the ideal 'kern'  of algebraic relations between the
    … …  
    249237// e.g.:
    250238               def @ring = L[2]; setring @ring ; kern;
    251         ");
    252         }
    253       if(a==2)
    254         {
    255           @L=P1,R1;
    256           kill phi,vars;
    257         }
    258     }
    259 
     239      ");
     240    }
     241    if(a==2)
     242    {
     243      @L=P1,R1;
     244      kill phi,vars;
     245    }
     246  }
    260247  return(@L);
    261 
    262 
    263 
    264248}
    265249///////////////////////////////////////////////////////////////////////////////
    … …  
    287271
    288272  if(b !=0) //means that the basering is a quotient ring
    289     {
    290       p=reduce(p,groebner(0));
    291       dom=reduce(groebner,std(0));
    292     }
     273  {
     274    p=reduce(p,groebner(0));
     275    dom=reduce(groebner,std(0));
     276  }
    293277
    294278  int i,choose;
    … …  
    296280
    297281  if((size(#)>0) && (typeof(#[1])=="int"))
    298     {
    299       choose = #[1];
    300     }
     282  {
     283    choose = #[1];
     284  }
    301285  if (size(#)>1)
    302     {
    303       choose =#[2];
    304     }
     286  {
     287    choose =#[2];
     288  }
    305289
    306290  //=======================first algorithm(default)=========================
    307    if ( choose == 0 )
    308   {
    309      list L = algebra_containment(lead(p),lead(dom),1);
    310      if( L[1]==1 )
    311      {
    312   // the ring L[2] = char(bsr),(x(1..nvars(bsr)),y(1..z)),(dp(n),dp(m)),
    313   // contains poly check s.t. LT(p) is of the form check(LT(f1),...,LT(fr))
    314        def s1 = L[2];
    315        map psi = s1,maxideal(1),dom;
    316        poly re = p - psi(check);
    317        // divide by the maximal power of #[1]
    318           if ( (size(#)>0) && (typeof(#[1])=="poly") )
    319 
    320             {  while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
    321              {
    322                re=re/#[1];
    323              }
    324             }
    325           return(re);
    326      }
    327   return(p);
    328   }
    329    //======================2end variant of algorithm=========================
    330    //It uses two different commands for elimaination.
    331    //if(choose==1):"elimainate"command.
    332    //if (choose==2):"nselect" command.
    333    else
     291  if ( choose == 0 )
     292  {
     293    list L = algebra_containment(lead(p),lead(dom),1);
     294    if( L[1]==1 )
     295    {
     296      // the ring L[2] = char(bsr),(x(1..nvars(bsr)),y(1..z)),(dp(n),dp(m)),
     297      // contains poly check s.t. LT(p) is of the form check(LT(f1),...,LT(fr))
     298      def s1 = L[2];
     299      map psi = s1,maxideal(1),dom;
     300      poly re = p - psi(check);
     301      // divide by the maximal power of #[1]
     302      if ( (size(#)>0) && (typeof(#[1])=="poly") )
    334303      {
    335         poly v=product(maxideal(1));
    336 
    337   //------------- change the basering bsr to bsr[@(0),...,@(z)] ----------
    338   execute("ring s="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;");
    339 // Ev hier die Reihenfolge der Vars aendern. Dazu muss unten aber entsprechend
    340 // geaendert werden:
    341 //  execute("ring s="+charstr(basering)+",(@(0..z),"+varstr(basering)+"),dp;");
    342 
    343   //constructs the leading ideal of dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z))
    344      ideal dom=imap(bsr,dom);
    345      for (i=1;i<=z;i++)
    346      {
    347         dom[i]=lead(dom[i])-var(nvars(bsr)+i+1);
    348      }
    349      dom=lead(imap(bsr,p))-@(0),dom;
    350 
    351       //---------- eliminate the variables of the basering bsr --------------
    352   //i.e. computes dom intersected with K[@(0),...,@(z)].
    353 
    354      if(choose==1)
    355        {
    356          ideal kern=eliminate(dom,imap(bsr,v));//eliminate does not need a
    357                                               //standard basis as input.
    358        }
    359      if(choose==2)
    360        {
    361          ideal kern= nselect(groebner(dom),1,n);//"nselect" is combinatorial command
     304        while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
     305        {
     306          re=re/#[1];
     307        }
     308      }
     309      return(re);
     310    }
     311    return(p);
     312  }
     313  //======================2end variant of algorithm=========================
     314  //It uses two different commands for elimaination.
     315  //if(choose==1):"elimainate"command.
     316  //if (choose==2):"nselect" command.
     317  else
     318  {
     319    poly v=product(maxideal(1));
     320
     321    //------------- change the basering bsr to bsr[@(0),...,@(z)] ----------
     322    execute("ring s="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;");
     323    // Ev hier die Reihenfolge der Vars aendern. Dazu muss unten aber entsprechend
     324    // geaendert werden:
     325    //  execute("ring s="+charstr(basering)+",(@(0..z),"+varstr(basering)+"),dp;");
     326
     327    //constructs the leading ideal of dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z))
     328    ideal dom=imap(bsr,dom);
     329    for (i=1;i<=z;i++)
     330    {
     331      dom[i]=lead(dom[i])-var(nvars(bsr)+i+1);
     332    }
     333    dom=lead(imap(bsr,p))-@(0),dom;
     334
     335    //---------- eliminate the variables of the basering bsr --------------
     336    //i.e. computes dom intersected with K[@(0),...,@(z)].
     337
     338    if(choose==1)
     339    {
     340      ideal kern=eliminate(dom,imap(bsr,v));//eliminate does not need a
     341                                            //standard basis as input.
     342    }
     343    if(choose==2)
     344    {
     345      ideal kern= nselect(groebner(dom),1,n);//"nselect" is combinatorial command
    362346                                         //which uses the internal command
    363                                         // "simplify"
    364        }
    365 
    366   //---------  test wether @(0)-h(@(1),...,@(z)) is in ker ---------------
    367   // for some poly h and divide by maximal power of q=#[1]
    368      poly h;
    369      z=size(kern);
    370      for (i=1;i<=z;i++)
    371      {
    372         h=kern[i]/@(0);
    373         if (deg(h)==0)
    374         {  h=(1/h)*kern[i];
     347                                         // "simplify"
     348    }
     349
     350    //---------  test wether @(0)-h(@(1),...,@(z)) is in ker ---------------
     351    // for some poly h and divide by maximal power of q=#[1]
     352    poly h;
     353    z=size(kern);
     354    for (i=1;i<=z;i++)
     355    {
     356      h=kern[i]/@(0);
     357      if (deg(h)==0)
     358      {
     359        h=(1/h)*kern[i];
    375360        // define the map psi : s ---> bsr defined by @(i) ---> p,dom[i]
    376            setring bsr;
    377            map psi=s,maxideal(1),p,dom;
    378            poly re=psi(h);
    379            // divide by the maximal power of #[1]
    380            if ((size(#)>0) && (typeof(#[1])== "poly") )
    381            {  while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
    382               {  re=re/#[1];
    383               }
    384            }
    385            return(re);
     361        setring bsr;
     362        map psi=s,maxideal(1),p,dom;
     363        poly re=psi(h);
     364        // divide by the maximal power of #[1]
     365        if ((size(#)>0) && (typeof(#[1])== "poly") )
     366        {
     367          while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
     368          {
     369            re=re/#[1];
     370          }
    386371        }
    387      }
    388      setring bsr;
    389      return(p);
    390   }
    391 
    392 
     372        return(re);
     373      }
     374    }
     375    setring bsr;
     376    return(p);
     377  }
    393378}
    394379example
    … …  
    421406  p1=p;
    422407  while(p1!=0)
    423     {
    424       p2=reduction(p1,dom,#);
    425       if(p2!=p1)
    426         {
    427           p1=p2;
    428         }
    429       else
    430         {
    431           re=re+lead(p2);
    432           p1=p2-lead(p2);
    433         }
    434     }
     408  {
     409    p2=reduction(p1,dom,#);
     410    if(p2!=p1)
     411    {
     412      p1=p2;
     413    }
     414    else
     415    {
     416      re=re+lead(p2);
     417      p1=p2-lead(p2);
     418    }
     419  }
    435420  return(re);
    436421}
    … …  
    459444  poly re;
    460445  if(typeof(id)=="ideal")
    461     {
    462       int i=ncols(id);
    463       for(z=1;z<=i;z++)
     446  {
     447    int i=ncols(id);
     448    for(z=1;z<=i;z++)
    464449    {
    465450      if(k==0)
    466         {
    467 
    468           id[z]=completeReduction(id[z],dom,#);
    469         }
     451      {
     452        id[z]=completeReduction(id[z],dom,#);
     453      }
    470454      else
    471         {
    472           id[z]=completeReduction1(id[z],dom,#);//tail reduction.
    473         }
    474     }
    475   Red=simplify(id,7);
    476   return(Red);
    477     }
     455      {
     456        id[z]=completeReduction1(id[z],dom,#);//tail reduction.
     457      }
     458    }
     459    Red=simplify(id,7);
     460    return(Red);
     461  }
    478462  if(typeof(id)=="poly")
    479     {
    480       if(k==0)
    481         {
    482           re=completeReduction(id,dom,#);
    483         }
    484       else
    485         {
    486           re=completeReduction1(id,dom,#);
    487         }
    488       return(re);
    489     }
     463  {
     464    if(k==0)
     465    {
     466      re=completeReduction(id,dom,#);
     467    }
     468    else
     469    {
     470      re=completeReduction1(id,dom,#);
     471    }
     472    return(re);
     473  }
    490474}
    491475example
    … …  
    509493  z=ncols(id);
    510494  for(i=1;i<=z;i++)
    511     {
    512       Rest=id;
    513       Rest[i]=0;
    514       Rest=simplify(Rest,2);
    515       if(k==0)
    516         {
    517           intRed[i]=completeReduction(id[i],Rest,#);
    518         }
    519       else
    520         {
    521           intRed[i]=completeReduction1(id[i],Rest,#);
    522         }
    523     }
     495  {
     496    Rest=id;
     497    Rest[i]=0;
     498    Rest=simplify(Rest,2);
     499    if(k==0)
     500    {
     501      intRed[i]=completeReduction(id[i],Rest,#);
     502    }
     503    else
     504    {
     505      intRed[i]=completeReduction1(id[i],Rest,#);
     506    }
     507  }
    524508  intRed=simplify(intRed,7);//1+2+4 in simplify command
    525509  return(intRed);
    526 
    527510}
    528511//////////////////////////////////////////////////////////////////////////////
    … …  
    549532  S=intRed(id,k,#);
    550533  while(size(S)!=size(oldS))
    551     {
    552       L=Spoly1(L,S,Red,2);
    553       Red=L[1];
    554       Red=sagbiNF(Red,S,k,#);
    555       oldS=S;
    556       S=S+Red;
    557     }
     534  {
     535    L=Spoly1(L,S,Red,2);
     536    Red=L[1];
     537    Red=sagbiNF(Red,S,k,#);
     538    oldS=S;
     539    S=S+Red;
     540  }
    558541  return(S);
    559542}
    … …  
    591574  S=intRed(id,k,#);
    592575  while((size(S)!=size(oldS))&&(counter<=c))
    593     {
    594       L=Spoly1(L,S,Red,2);
    595       Red=L[1];
    596       Red=sagbiNF(Red,S,k,#);
    597       oldS=S;
    598       S=S+Red;
    599       counter=counter+1;
    600     }
     576  {
     577    L=Spoly1(L,S,Red,2);
     578    Red=L[1];
     579    Red=sagbiNF(Red,S,k,#);
     580    oldS=S;
     581    S=S+Red;
     582    counter=counter+1;
     583  }
    601584  return(S);
    602585}
Note: See TracChangeset for help on using the changeset viewer.