Changeset 12603f in git for Singular/LIB/swalk.lib


Ignore:
Timestamp:
Mar 8, 2010, 10:11:46 AM (14 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
56b0c82650ffab78c72a5629796b5fe751482448
Parents:
8e865cd8f45ac33fd409c1996172450700565c7b
Message:
rename signature.lib to surfsig.lib

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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/swalk.lib

    r8e865c r12603f  
    1 ///////////////////////////////////////////////////////////////
    2 
    31version="$Id$";
    42category="Commutative Algebra";
    5 
    63info="
    7 LIBRARY: swalk.lib   Sagbi Walk Conversion Algorithm
    8 AUTHOR: Junaid Alam Khan
     4LIBRARY: swalk.lib               Sagbi Walk Conversion Algorithm
     5AUTHOR:  Junaid Alam Khan        junaidalamkhan@gmail.com
     6
     7OVERVIEW:
     8 A library for computing the Sagbi basis of subalgebra through Sagbi walk
     9 algorithm.
     10
     11THEORY: The concept of SAGBI ( Subalgebra Analog to Groebner Basis for Ideals)
     12 is defined in [L. Robbiano, M. Sweedler: Subalgebra Bases, volume 42, volume
     13 1430 of Lectures Note in Mathematics series, Springer-Verlag (1988),61-87].
     14 The Sagbi Walk algorithm is the subalgebra analogue to the Groebner Walk
     15 algorithm which has been proposed in [S. Collart, M. Kalkbrener and D.Mall:
     16 Converting bases with the Grobner Walk. J. Symbolic Computation 24 (1997),
     17 465-469].
    918
    1019PROCEDURE:
    11   swalk(ideal[,intvec]);   Sagbi basis of subalgebra via Sagbiwalk algorithm
     20 swalk(ideal[,intvec]);   Sagbi basis of subalgebra via Sagbi walk algorithm
    1221";
    1322
    14 LIB "algebra.lib";
    1523LIB "sagbi.lib";
    1624//////////////////////////////////////////////////////////////////////////////
    1725proc swalk(ideal Gox, list #)
    18 "USAGE:  swalk(i); ideal;
    19          swalk( i,v,w); i ideal, v,w integer vectors;
    20 RETURN:  If list #= v,w (resp empty) then it compute the sagbi basis of
    21          the subalgebra defined by the generators of ideal,
    22          calculated via the Sagbi walk algorithm from the ordering
    23          (a(v),lp) (resp dp) tothe ordering (a(w),lp) ( resp lp).
     26"USAGE:  swalk(i[,v,w]); i ideal, v,w int vectors
     27RETURN: The sagbi basis of the subalgebra defined by the generators of i,
     28        calculated via the Sagbi walk algorithm from the ordering dp to lp
     29        if v,w are not given (resp. from the ordering (a(v),lp) to the
     30        ordering (a(w),lp) if v and w are given).
    2431EXAMPLE: example swalk; shows an example
    2532"
     
    3340  option(redSB);
    3441  def xR = basering;
    35   execute("ring ostR = ("+charstr(xR)+"),("+varstr(xR)+"),"+ord_str+";");
     42  list rl=ringlist(xR);
     43  rl[3][1][1]="dp";
     44  def ostR=ring(rl);
     45  setring ostR;
    3646  def new_ring = basering;
    3747  ideal Gnew = fetch(xR, Gox);
    3848  Gnew=sagbi(Gnew,1);
    3949  Gnew=interreduceSd(Gnew);
    40   vector curr_weight=Ctv(icurr_weight);
    41   vector target_weight=Ctv(itarget_weight);
     50  vector curr_weight=changeTypeInt(icurr_weight);
     51  vector target_weight=changeTypeInt(itarget_weight);
    4252  ideal Gold;
    4353  list l;
     
    5363       if(curr_weight==target_weight){n=1;}
    5464       else {
    55               l=Df(Gold);
     65              l=collectDiffExpo(Gold);
    5666              ulast=last(curr_weight, target_weight, l);
    57               vector new_weight=newvec(curr_weight,target_weight,ulast);
     67              vector new_weight=(1-ulast)*curr_weight+ulast*target_weight;
    5868              vector w=cleardenom(new_weight);
    59               v=CT(w);
     69              v=changeType(w);
    6070              list p= ringlist(old_ring);
    6171              p[3][1][2]= v;
     
    7787example
    7888{
    79   "EXAMPLE:";
    80   //** compute a Sagbi basis of I w.r.t. lp.
     89  "EXAMPLE:";echo = 2;
    8190  ring r = 0,(x,y), lp;
    8291  ideal I =x2,y2,xy+y,2xy2+y3;
    83   intvec v=1,1;
    84   intvec w=1,0;
    85   swalk(I,v,w);
    86 }
    87 
    88 proc inprod(vector v,vector w)
     92  swalk(I);
     93}
     94
     95//////////////////////////////////////////////////////////////////////////////
     96static proc inprod(vector v,vector w)
    8997"USAGE:  inprod(v,w); v,w vectors
    9098RETURN:  inner product of vector v and w
     
    108116}
    109117
    110 
    111 proc df(poly f)
    112 "USAGE:  df(f); f polynomial
    113 RETURN:  a list of integers vectors which are the differnce of exponent vector
    114          of leading monomial of f with the exponent vector of
    115          of other monmials in f.
    116 EXAMPLE: example df; shows an example
     118//////////////////////////////////////////////////////////////////////////////
     119static proc diffExpo(poly f)
     120"USAGE:  diffExpo(f); f polynomial
     121RETURN:  a list of integers vectors which are the difference of exponent
     122         vector of leading monomial of f with the exponent vector of of other
     123         monmials in f.
     124EXAMPLE: example diffExpo; shows an example
    117125"
    118126{
     
    120128  int i;
    121129  intvec v;
    122   for(i=2;i<=size(f);i++)
     130  for(i=size(f);i>=2;i--)
    123131   {
    124132     v=leadexp(f[1])-leadexp(f[i]);
    125      l[size(l)+1]=v;
     133     l[i-1]=v;
    126134   }
    127135 return(l);
     
    131139   ring r=0,(x,y,z),lp;
    132140   poly f = xy+z2 ;
    133    df(f);
    134 }
    135 
    136 
    137 proc Df( ideal i)
    138 "USAGE:  Df(i); i ideal
    139 RETURN:  a list which contain df(f), for all generators f of ideal i
    140 EXAMPLE: example Df; shows an example
     141   diffExpo(f);
     142}
     143
     144//////////////////////////////////////////////////////////////////////////////
     145static proc collectDiffExpo( ideal i)
     146"USAGE:  collectDiffExpo(i); i ideal
     147RETURN:  a list which contains diffExpo(f), for all generators f of ideal i
     148EXAMPLE: example collectDiffExpo; shows an example
    141149"
    142150{
    143151 list l;
    144152 int j;
    145  for(j=1;j<=size(i);j++)
     153 for(j=ncols(i); j>=1;j--)
    146154  {
    147    l[size(l)+1]=df(i[j]);
     155   l[j]=diffExpo(i[j]);
    148156  }
    149 return(l);
    150   }
     157  return(l);
     158}
    151159example
    152160{ "EXAMPLE:"; echo = 2;
    153161   ring r=0,(x,y,z),lp;
    154162   ideal I = xy+z2,y3+x2y2;
    155    Df(I);
    156 }
    157 
    158 
    159 proc newvec(vector v, vector w, number u)
    160 "USAGE:  newvec(v,w,u); v,w vectors, u number
    161  RETURN: the vector (1-u)v+(u)w.
    162  EXAMPLE: example newvec; shows an example
    163 "
    164 {
    165   vector wnew ;
    166   wnew=(1-u)*v+(u)*w ;
    167   return(wnew);
    168 }
    169 example
    170 { "EXAMPLE:"; echo = 2;
    171    ring r=0,(x,y,z),lp;
    172    vector v=[0,0,1];
    173    vector w=[1,0,0];
    174    number u=2/3;
    175    newvec(v,w,u);
    176 }
    177 
    178 proc CT(vector v)
    179 "USAGE:  CT(v); v  vector
    180 RETURN:  change the type of  vector v into integer vector.
    181 EXAMPLE: example CT; shows an example
     163   collectDiffExpo(I);
     164}
     165
     166//////////////////////////////////////////////////////////////////////////////
     167static proc changeType(vector v)
     168"USAGE:  changeType(v); v  vector
     169RETURN:  change the type of  vector
     170         v into integer vector.
     171
     172EXAMPLE: example changeType; shows an example
    182173"
    183174{
     
    194185   ring r=0,(x,y,z),lp;
    195186   vector v = [2,1,3];
    196    CT(v);
    197 }
    198 
    199 
    200 proc Ctv( intvec v)
    201 "USAGE:  Ctv(v); v integer vector
     187   changeType(v);
     188}
     189//////////////////////////////////////////////////////////////////////////////
     190static proc changeTypeInt( intvec v)
     191"USAGE:  changeTypeInt(v); v integer vector
    202192RETURN:  change the type of integer vector v into vector.
    203 EXAMPLE: example Ctv; shows an example
    204 "
    205 {
    206 vector w;
    207 int j ;
    208 for(j=1;j<=size(v);j++)
    209  {
    210    w=w+v[j]*gen(j);
    211  }
    212 return(w);
     193EXAMPLE: example changeTypeInt; shows an example
     194"
     195{
     196   vector w;
     197   int j ;
     198   for(j=1;j<=size(v);j++)
     199   {
     200     w=w+v[j]*gen(j);
     201   }
     202   return(w);
    213203}
    214204example
     
    216206   ring r=0,(x,y,z),lp;
    217207   intvec v = 4,2,3;
    218    Ctv(v);
    219 }
    220 
    221 proc last( vector c, vector t,list l)
     208   changeTypeInt(v);
     209}
     210
     211//////////////////////////////////////////////////////////////////////////////
     212static proc last( vector c, vector t,list l)
    222213"USAGE: last(c,t,l); c,t vectors, l list
    223 RETURN: a  parametric value which corresponds to vector lies along
    224         the path between c and t using list l of integer
    225         vectors. This vector is the last vector on old Sagbi cone
     214RETURN: a  parametric value which corresponds to vector lies along the path
     215        between c and t using list l of integer vectors. This vector is the
     216        last vector on old Sagbi cone
    226217EXAMPLE: example last; shows an example
    227218"
     
    232223 vector v;
    233224 for(i=1;i<=size(l);i++)
    234   {
     225 {
    235226    for(j=1;j<=size(l[i]);j++)
    236       {
     227    {
    237228        v=0;
    238229        for(k=1;k<=size(l[i][j]);k++)
    239           {
     230        {
    240231            v=v+l[i][j][k]*gen(k);
    241           }
     232        }
    242233        poly n= inprod(c,v);
    243234        poly q= inprod(t,v);
     
    246237        number z=a-b;
    247238        if(b<0)
    248           {
     239        {
    249240            u=a/z;
    250241            if(u<ul) {ul=u;}
    251           }
     242        }
    252243        kill a,b,z,n,q ;
    253        }
    254    }
    255 return(ul);
     244    }
     245 }
     246 return(ul);
    256247}
    257248example
     
    261252   vector w=[1,0,0];
    262253   ideal i=z2+xy,x2y2+y3;
    263     list l=Df(i);
     254    list l=collectDiffExpo(i);
    264255    last(v,w,l)
    265256}
    266257
    267 proc initialForm(poly P)
     258//////////////////////////////////////////////////////////////////////////////
     259static proc initialForm(poly P)
    268260"USAGE:  initialForm(P); P polynomial
    269261RETURN:  sum of monomials of P with maximum w-degree
     
    275267 int i=1;
    276268 while(deg(P[i])==deg(P[1]))
    277    {
     269 {
    278270     q=q+P[i];
    279271     i++;
    280272     if(i>size(P)) {break;}
    281    }
    282   return(q);
     273 }
     274 return(q);
    283275}
    284276example
     
    289281}
    290282
    291 proc Initial(ideal I)
     283//////////////////////////////////////////////////////////////////////////////
     284static proc Initial(ideal I)
    292285"USAGE:  Initial(I); I ideal
    293 RETURN:  an ideal which is generate by the InitialForm of the generators of I.
     286RETURN:  an ideal which is generate by the InitialForm
     287         of the generators of I.
    294288EXAMPLE: example Initial; shows an example
    295289"
     
    297291 ideal J;
    298292 int i;
    299  for(i=1;i<=size(I);i++)
    300    {
     293 for(i=1;i<=ncols(I);i++)
     294 {
    301295     J[i]=initialForm(I[i]);
    302    }
     296 }
    303297 return(J);
    304298}
     
    310304}
    311305
    312 proc Lift(ideal In,ideal InG,ideal Gold)
    313 "USAGE:  Lift(In, InG, Gold); In,InG, Gold ideals
    314          Gold given by Sagbi basis {g_1,...,g_t}
    315          In given by tne initial forms In(g_1),...,In(g_t)
    316          InG={h_1,...,h_s} a Sagbi basis of In
    317 RETURN:  P_j, a polynomial in K[y_1,..,y_t] such that
    318          h_j=P_j(In(g_1),...,In_(g_t)) and return f_j=P_j(g_1,...,g_t)
     306//////////////////////////////////////////////////////////////////////////////
     307static proc Lift(ideal In,ideal InG,ideal Gold)
     308"USAGE:  Lift(In, InG, Gold); In, InG, Gold ideals;
     309         Gold given by Sagbi basis {g_1,...,g_t},
     310         In given by tne initial forms In(g_1),...,In(g_t),
     311         InG = {h_1,...,h_s} a Sagbi basis of In
     312RETURN:  P_j, a polynomial in K[y_1,..,y_t] such that h_j =
     313         P_j(In(g_1),...,In_(g_t))
     314         and return f_j = P_j(g_1,...,g_t)
    319315EXAMPLE: example Lift; shows an example
    320316"
     
    322318  int i;
    323319  ideal J;
    324   for(i=1;i<=size(InG);i++)
     320  for(i=1;i<=ncols(InG);i++)
    325321  {
    326322    poly f=InG[i];
     
    344340}
    345341
    346 
    347 proc Convert( ideal Gold )
     342//////////////////////////////////////////////////////////////////////////////
     343static proc Convert( ideal Gold )
    348344"USAGE: Convert(I); I ideal
    349345RETURN: Convert old Sagbi basis into new Sagbi basis
     
    363359}
    364360
    365 proc interreduceSd(ideal I)
     361//////////////////////////////////////////////////////////////////////////////
     362static proc interreduceSd(ideal I)
    366363"USAGE:  interreduceSd(I); I ideal
    367364RETURN:  interreduceSd the set of generators if I with
     
    374371  int i,j,k;
    375372  poly f;
    376   for(k=1;k<=size(I);k++)
    377     {l[k]=I[k];}
     373  for(k=1;k<=ncols(I);k++)
     374  {l[k]=I[k];}
    378375  for(j=1;j<=size(l);j++)
    379    {
     376  {
    380377     f=l[j];
    381378     M=delete(l,j);
    382379     for(i=1;i<=size(M);i++)
    383        { B[i]=M[i];}
     380     { B[i]=M[i];}
    384381     f=sagbiNF(f,B,1);
    385382     J=J+f;
    386    }
     383  }
    387384  return(J);
    388385}
     
    394391}
    395392
     393//////////////////////////////////////////////////////////////////////////////
    396394static proc OrderStringalp(string Wpal,list #)
    397395{
     
    478476}
    479477
     478//////////////////////////////////////////////////////////////////////////////
    480479static proc OrderStringalp_NP(string Wpal,list #)
    481480{
     
    489488    if(size(#) == 1)
    490489    {
    491       if(typeof(#[1]) == "intvec") {
     490      if(typeof(#[1]) == "intvec")
     491      {
    492492        curr_weight = #[1];
    493493
    494         if(Wpal == "al"){
     494        if(Wpal == "al")
     495        {
    495496          order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
    496497        }
    497         else {
     498        else
     499        {
    498500          order_str = "(Wp("+string(#[1])+"),C)";
    499501        }
    500502      }
    501503      else {
    502         if(typeof(#[1]) == "int"){
     504        if(typeof(#[1]) == "int")
     505        {
    503506          nP = #[1];
    504507        }
    505         else {
     508        else
     509        {
    506510          print("// ** the input must be \"(ideal, int)\" or ");
    507511          print("// **                   \"(ideal, intvec)\"");
     
    514518     if(size(#) == 2)
    515519     {
    516        if(typeof(#[1]) == "intvec" and typeof(#[2]) == "int"){
     520       if(typeof(#[1]) == "intvec" and typeof(#[2]) == "int")
     521       {
    517522         curr_weight = #[1];
    518523
    519          if(Wpal == "al"){
     524         if(Wpal == "al")
     525         {
    520526           order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
    521527         }
    522          else {
     528         else
     529         {
    523530           order_str = "(Wp("+string(#[1])+"),C)";
    524531         }
    525532       }
    526        else{
    527          if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec"){
     533       else
     534       {
     535         if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec")
     536         {
    528537           curr_weight = #[1];
    529538           target_weight = #[2];
    530539
    531            if(Wpal == "al"){
     540           if(Wpal == "al")
     541           {
    532542             order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
    533543           }
    534            else {
     544           else
     545           {
    535546             order_str = "(Wp("+string(#[1])+"),C)";
    536547           }
    537548         }
    538          else{
     549         else
     550         {
    539551           print("// ** the input  must be \"(ideal,intvec,int)\" or ");
    540552           print("// **                    \"(ideal,intvec,intvec)\"");
     
    544556     }
    545557     else {
    546        if(size(#) == 3) {
     558       if(size(#) == 3)
     559       {
    547560         if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec" and
    548561            typeof(#[3]) == "int")
     
    551564           target_weight = #[2];
    552565           nP = #[3];
    553            if(Wpal == "al"){
     566           if(Wpal == "al")
     567           {
    554568             order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
    555569           }
    556            else {
     570           else
     571           {
    557572             order_str = "(Wp("+string(#[1])+"),C)";
    558573           }
    559574         }
    560          else{
     575         else
     576         {
    561577           print("// ** the input must be \"(ideal,intvec,intvec,int)\"");
    562578           print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");
     
    564580         }
    565581       }
    566        else{
     582       else
     583       {
    567584         print("// ** The given input is wrong");
    568585         print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");
     
    579596}
    580597
     598///////////////////////////////////////////////////////////////////////////////*
     599Further examples
     600ring r=0,(x,y,z),lp;
     601
     602ideal I=x2y4, y4z2, xy4z+y2z, xy6z2+y10z5;
     603
     604ideal Q=x2y4, y4z2, xy4z+y2z, xy6z2+y14z7;
     605
     606ideal J=x2y4, y4z2, xy4z+y2z, xy6z2+y18z9;
     607
     608ideal K=x2,y2,xy+y,2xy2+y5,z3+x;
     609
     610ideal L=x2+y,y2+z,x+z2;
     611*/
     612
     613
Note: See TracChangeset for help on using the changeset viewer.