source: git/Singular/LIB/swalk.lib @ 4083fa

spielwiese
Last change on this file since 4083fa was 4083fa, checked in by Stephan Oberfranz <oberfran@…>, 9 years ago
--- Merge branch 'spielwiese' of github.com:Singular/Sources into spielwiese Conflicts: Singular/LIB/grwalk.lib Singular/LIB/modwalk.lib Singular/walk.cc
  • Property mode set to 100755
File size: 20.4 KB
Line 
1////////////////////////////////////////////////////////
2version="version swalk.lib 4.0.0.0 Jun_2013 ";
3category="Commutative Algebra";
4info="
5LIBRARY: swalk.lib               Sagbi Walk Conversion Algorithm
6AUTHOR:  Junaid Alam Khan        junaidalamkhan@gmail.com
7
8OVERVIEW:
9 A library for computing the Sagbi basis of subalgebra through Sagbi walk
10 algorithm.
11
12THEORY: The concept of SAGBI ( Subalgebra Analog to Groebner Basis for Ideals)
13 is defined in [L. Robbiano, M. Sweedler: Subalgebra Bases, volume 42, volume
14 1430 of Lectures Note in Mathematics series, Springer-Verlag (1988),61-87].
15 The Sagbi Walk algorithm is the subalgebra analogue to the Groebner Walk
16 algorithm which has been proposed in [S. Collart, M. Kalkbrener and D.Mall:
17 Converting bases with the Grobner Walk. J. Symbolic Computation 24 (1997),
18 465-469].
19
20PROCEDURES:
21 swalk(ideal[,intvec]);   Sagbi basis of subalgebra via Sagbi walk algorithm
22 rswalk(ideal,int,int[,intvec]); Sagbi basis of subalgebra via Random Sagbi Walk Algorithm
23";
24
25LIB "sagbi.lib";
26//////////////////////////////////////////////////////////////////////////////
27proc swalk(ideal Gox, list #)
28"USAGE:  swalk(i[,v,w]); i ideal, v,w int vectors
29RETURN: The sagbi basis of the subalgebra defined by the generators of i,
30        calculated via the Sagbi walk algorithm from the ordering dp to lp
31        if v,w are not given (resp. from the ordering (a(v),lp) to the
32        ordering (a(w),lp) if v and w are given).
33EXAMPLE: example swalk; shows an example
34"
35{
36  /* we use ring with ordering (a(...),lp,C) */
37  list OSCTW    = OrderStringalp_NP("al", #);//"dp"
38  string ord_str =   OSCTW[2];
39  intvec icurr_weight   =   OSCTW[3]; /* original weight vector */
40  intvec itarget_weight =   OSCTW[4]; /* terget weight vector */
41  kill OSCTW;
42  option(redSB);
43  def xR = basering;
44  list rl=ringlist(xR);
45  rl[3][1][1]="dp";
46  def ostR=ring(rl);
47  setring ostR;
48  def new_ring = basering;
49  ideal Gnew = fetch(xR, Gox);
50  Gnew=sagbi(Gnew,1);
51  Gnew=interreduceSd(Gnew);
52  vector curr_weight=changeTypeInt(icurr_weight);
53  vector target_weight=changeTypeInt(itarget_weight);
54  ideal Gold;
55  list l;
56  intvec v;
57  int n=0;
58  while(n==0)
59    {
60       Gold=Gnew;
61       def old_ring=new_ring;
62       setring old_ring;
63       number ulast;
64       kill new_ring;
65       if(curr_weight==target_weight){n=1;}
66       else {
67              l=collectDiffExpo(Gold);
68              ulast=last(curr_weight, target_weight, l);
69              vector new_weight=(1-ulast)*curr_weight+ulast*target_weight;
70              vector w=cleardenom(new_weight);
71              v=changeType(w);
72              list p= ringlist(old_ring);
73              p[3][1][2]= v;
74              def new_ring=ring(p);
75              setring new_ring;
76              ideal Gold=fetch(old_ring,Gold);
77              vector curr_weight=fetch(old_ring,new_weight);
78              vector target_weight=fetch(old_ring,target_weight);
79              kill old_ring;
80              ideal Gnew=Convert(Gold);
81              Gnew=interreduceSd(Gnew);
82           }
83    }
84   setring xR;
85   ideal result = fetch(old_ring, Gnew);
86   attrib(result,"isSB",1);
87   return (result);
88}
89example
90{
91  "EXAMPLE:";echo = 2;
92  ring r = 0,(x,y), lp;
93  ideal I =x2,y2,xy+y,2xy2+y3;
94  swalk(I);
95}
96//////////////////////////////////////////////////////////////////////////////
97proc rswalk(ideal Gox, int weight_rad, int pdeg, list #)
98"USAGE:  rswalk(i,weight_rad,p_deg[,v,w]); i ideal, v,w int vectors
99RETURN: The sagbi basis of the subalgebra defined by the generators of i,
100        calculated via the Sagbi walk algorithm from the ordering dp to lp
101        if v,w are not given (resp. from the ordering (a(v),lp) to the
102        ordering (a(w),lp) if v and w are given).
103EXAMPLE: example swalk; shows an example
104"
105{
106  /* we use ring with ordering (a(...),lp,C) */
107  list OSCTW    = OrderStringalp_NP("al", #);//"dp"
108  string ord_str =   OSCTW[2];
109  intvec icurr_weight   =   OSCTW[3]; /* original weight vector */
110  intvec itarget_weight =   OSCTW[4]; /* terget weight vector */
111  kill OSCTW;
112  option(redSB);
113  def xR = basering;
114  list rl=ringlist(xR);
115  rl[3][1][1]="dp";
116  def ostR=ring(rl);
117  setring ostR;
118  def new_ring = basering;
119  ideal Gnew = fetch(xR, Gox);
120  Gnew=sagbi(Gnew,1);
121  Gnew=interreduceSd(Gnew);
122  vector curr_weight=changeTypeInt(icurr_weight);
123  vector target_weight=changeTypeInt(itarget_weight);
124  ideal Gold;
125  list l;
126  intvec v;
127  int n=0;
128  while(n==0)
129    {
130       Gold=Gnew;
131       def old_ring=new_ring;
132       setring old_ring;
133
134       kill new_ring;
135       if(curr_weight==target_weight){n=1;}
136       else {
137              l=collectDiffExpo(Gold);
138              vector new_weight=RandomNextWeight(Gold, l, curr_weight, target_weight, weight_rad, pdeg);
139              vector w=cleardenom(new_weight);
140              v=changeType(w);
141              list p= ringlist(old_ring);
142              p[3][1][2]= v;
143              def new_ring=ring(p);
144              setring new_ring;
145              ideal Gold=fetch(old_ring,Gold);
146              vector curr_weight=fetch(old_ring,new_weight);
147              vector target_weight=fetch(old_ring,target_weight);
148              kill old_ring;
149              ideal Gnew=Convert(Gold);
150              Gnew=interreduceSd(Gnew);
151           }
152    }
153   setring xR;
154   ideal result = fetch(old_ring, Gnew);
155   attrib(result,"isSB",1);
156   return (result);
157}
158example
159{
160  "EXAMPLE:";echo = 2;
161  ring r = 0,(x,y), lp;
162  ideal I =x2,y2,xy+y,2xy2+y3;
163  swalk(I,2,2);
164}
165//////////////////////////////////////////////////////////////////////////////
166static proc inprod(vector v,vector w)
167"USAGE:  inprod(v,w); v,w vectors
168RETURN:  inner product of vector v and w
169EXAMPLE: example inprod; shows an example
170"
171{
172  poly a;
173  int i;
174  for(i=1;i<=nvars(basering);i++)
175    {
176      a=a+v[i]*w[i] ;
177    }
178  return(a);
179}
180example
181{ "EXAMPLE:"; echo = 2;
182   ring r=0,(x,y,z),lp;
183   vector v =[1,-1,2];
184   vector w = [1,0,3];
185   inprod(v,w);
186}
187
188//////////////////////////////////////////////////////////////////////////////
189static proc diffExpo(poly f)
190"USAGE:  diffExpo(f); f polynomial
191RETURN:  a list of integers vectors which are the difference of exponent
192         vector of leading monomial of f with the exponent vector of of other
193         monmials in f.
194EXAMPLE: example diffExpo; shows an example
195"
196{
197  list l;
198  int i;
199  intvec v;
200  for(i=size(f);i>=2;i--)
201   {
202     v=leadexp(f[1])-leadexp(f[i]);
203     l[i-1]=v;
204   }
205 return(l);
206}
207example
208{ "EXAMPLE:"; echo = 2;
209   ring r=0,(x,y,z),lp;
210   poly f = xy+z2 ;
211   diffExpo(f);
212}
213
214//////////////////////////////////////////////////////////////////////////////
215static proc collectDiffExpo( ideal i)
216"USAGE:  collectDiffExpo(i); i ideal
217RETURN:  a list which contains diffExpo(f), for all generators f of ideal i
218EXAMPLE: example collectDiffExpo; shows an example
219"
220{
221 list l;
222 int j;
223 for(j=ncols(i); j>=1;j--)
224  {
225   l[j]=diffExpo(i[j]);
226  }
227  return(l);
228}
229example
230{ "EXAMPLE:"; echo = 2;
231   ring r=0,(x,y,z),lp;
232   ideal I = xy+z2,y3+x2y2;
233   collectDiffExpo(I);
234}
235
236//////////////////////////////////////////////////////////////////////////////
237static proc changeType(vector v)
238"USAGE:  changeType(v); v  vector
239RETURN:  change the type of  vector
240         v into integer vector.
241
242EXAMPLE: example changeType; shows an example
243"
244{
245  intvec w ;
246  int j ;
247  for(j=1;j<=nvars(basering);j++)
248   {
249     w[j]=int(leadcoef(v[j]));
250   }
251  return(w);
252}
253example
254{ "EXAMPLE:"; echo = 2;
255   ring r=0,(x,y,z),lp;
256   vector v = [2,1,3];
257   changeType(v);
258}
259//////////////////////////////////////////////////////////////////////////////
260static proc changeTypeInt( intvec v)
261"USAGE:  changeTypeInt(v); v integer vector
262RETURN:  change the type of integer vector v into vector.
263EXAMPLE: example changeTypeInt; shows an example
264"
265{
266   vector w;
267   int j ;
268   for(j=1;j<=size(v);j++)
269   {
270     w=w+v[j]*gen(j);
271   }
272   return(w);
273}
274example
275{ "EXAMPLE:"; echo = 2;
276   ring r=0,(x,y,z),lp;
277   intvec v = 4,2,3;
278   changeTypeInt(v);
279}
280
281//////////////////////////////////////////////////////////////////////////////
282static proc last( vector c, vector t,list l)
283"USAGE: last(c,t,l); c,t vectors, l list
284RETURN: a  parametric value which corresponds to vector lies along the path
285        between c and t using list l of integer vectors. This vector is the
286        last vector on old Sagbi cone
287EXAMPLE: example last; shows an example
288"
289{
290 number ul=1;
291 int i,j,k;
292 number u;
293 vector v;
294 for(i=1;i<=size(l);i++)
295 {
296    for(j=1;j<=size(l[i]);j++)
297    {
298        v=0;
299        for(k=1;k<=size(l[i][j]);k++)
300        {
301            v=v+l[i][j][k]*gen(k);
302        }
303        poly n= inprod(c,v);
304        poly q= inprod(t,v);
305        number a=leadcoef(n);
306        number b=leadcoef(q);
307        number z=a-b;
308        if(b<0)
309        {
310            u=a/z;
311            if(u<ul) {ul=u;}
312        }
313        kill a,b,z,n,q ;
314    }
315 }
316 return(ul);
317}
318example
319{ "EXAMPLE:"; echo = 2;
320   ring r=0,(x,y,z),(a(0,0,1),lp);
321   vector v= [0,0,1];
322   vector w=[1,0,0];
323   ideal i=z2+xy,x2y2+y3;
324    list l=collectDiffExpo(i);
325    last(v,w,l)
326}
327
328//////////////////////////////////////////////////////////////////////////////
329static proc initialForm(poly P)
330"USAGE:  initialForm(P); P polynomial
331RETURN:  sum of monomials of P with maximum w-degree
332         where w is first row of matrix of a given monomial ordering
333EXAMPLE: example initialForm; shows an example
334"
335{
336 poly q;
337 int i=1;
338 while(deg(P[i])==deg(P[1]))
339 {
340     q=q+P[i];
341     i++;
342     if(i>size(P)) {break;}
343 }
344 return(q);
345}
346example
347{ "EXAMPLE:"; echo = 2;
348   ring r=0,(x,y,z),dp;
349   poly f = x2+yz+z;
350   initialForm(f);
351}
352
353//////////////////////////////////////////////////////////////////////////////
354static proc Initial(ideal I)
355"USAGE:  Initial(I); I ideal
356RETURN:  an ideal which is generate by the InitialForm
357         of the generators of I.
358EXAMPLE: example Initial; shows an example
359"
360{
361 ideal J;
362 int i;
363 for(i=1;i<=ncols(I);i++)
364 {
365     J[i]=initialForm(I[i]);
366 }
367 return(J);
368}
369example
370{ "EXAMPLE:"; echo = 2;
371   ring r=0,(x,y,z),dp;
372   ideal I = x+1,x+y+1;
373   Initial(I);
374}
375
376//////////////////////////////////////////////////////////////////////////////
377static proc Lift(ideal In,ideal InG,ideal Gold)
378"USAGE:  Lift(In, InG, Gold); In, InG, Gold ideals;
379         Gold given by Sagbi basis {g_1,...,g_t},
380         In given by tne initial forms In(g_1),...,In(g_t),
381         InG = {h_1,...,h_s} a Sagbi basis of In
382RETURN:  P_j, a polynomial in K[y_1,..,y_t] such that h_j =
383         P_j(In(g_1),...,In_(g_t))
384         and return f_j = P_j(g_1,...,g_t)
385EXAMPLE: example Lift; shows an example
386"
387{
388  int i;
389  ideal J;
390  for(i=1;i<=ncols(InG);i++)
391  {
392    poly f=InG[i];
393    list l=algebra_containment(f,In,1);
394    def s=l[2];
395    map F=s,maxideal(1),Gold ;
396    poly g=F(check);
397    ideal k=g;
398    J=J+k;
399    kill g,l,s,F,f,k;
400  }
401  return(J);
402 }
403example
404{ "EXAMPLE:"; echo = 2;
405   ring r=0,(x,y,z),(a(2,0,3),lp);
406   ideal In = xy+z2,x2y2;
407   ideal InG=xy+z2,x2y2,xyz2+1/2z4;
408   ideal Gold=xy+z2,y3+x2y2;
409   Lift(In,InG,Gold);
410}
411
412//////////////////////////////////////////////////////////////////////////////
413static proc Convert( ideal Gold )
414"USAGE: Convert(I); I ideal
415RETURN: Convert old Sagbi basis into new Sagbi basis
416EXAMPLE: example Convert; shows an example
417"
418{
419 ideal In=Initial(Gold);
420 ideal InG=sagbi(In,1)+In;
421 ideal Gnew=Lift(In,InG,Gold);
422 return(Gnew);
423}
424example
425{ "EXAMPLE:"; echo = 2;
426   ring r=0,(x,y,z),lp;
427   ideal I=xy+z2, y3+x2y2;
428   Convert(I);
429}
430
431//////////////////////////////////////////////////////////////////////////////
432static proc interreduceSd(ideal I)
433"USAGE:  interreduceSd(I); I ideal
434RETURN:  interreduceSd the set of generators if I with
435         respect to a given term ordering
436EXAMPLE: example interreduceSd; shows an example
437"
438{
439  list l,M;
440  ideal J,B;
441  int i,j,k;
442  poly f;
443  for(k=1;k<=ncols(I);k++)
444  {l[k]=I[k];}
445  for(j=1;j<=size(l);j++)
446  {
447     f=l[j];
448     M=delete(l,j);
449     for(i=1;i<=size(M);i++)
450     { B[i]=M[i];}
451     f=sagbiNF(f,B,1);
452     J=J+f;
453  }
454  return(J);
455}
456example
457{ "EXAMPLE:"; echo = 2;
458   ring r=0,(x,y,z),lp;
459   ideal I = xy+z2,x2y2+y3;
460   interreduceSd(I);
461}
462
463//////////////////////////////////////////////////////////////////////////////
464static proc OrderStringalp(string Wpal,list #)
465{
466  int n= nvars(basering);
467  string order_str;
468  intvec curr_weight, target_weight;
469  curr_weight = system("Mivdp",n);
470  target_weight = system("Mivlp",n);
471
472   if(size(#) != 0)
473   {
474     if(size(#) == 1)
475     {
476       if(typeof(#[1]) == "intvec")
477       {
478         if(Wpal == "al"){
479           order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
480         }
481         else {
482           order_str = "(Wp("+string(#[1])+"),C)";
483         }
484         curr_weight = #[1];
485       }
486       else
487       {
488        if(typeof(#[1]) == "string")
489        {
490          if(#[1] == "Dp") {
491            order_str = "Dp";
492          }
493          else {
494            order_str = "dp";
495          }
496        }
497        else {
498          order_str = "dp";
499        }
500     }
501    }
502    else
503    {
504     if(size(#) == 2)
505     {
506       if(typeof(#[2]) == "intvec")
507       {
508         target_weight = #[2];
509       }
510       if(typeof(#[1]) == "intvec")
511       {
512         if(Wpal == "al"){
513           order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
514         }
515         else {
516           order_str = "(Wp("+string(#[1])+"),C)";
517         }
518         curr_weight = #[1];
519       }
520       else
521       {
522        if(typeof(#[1]) == "string")
523        {
524          if(#[1] == "Dp") {
525            order_str = "Dp";
526           }
527           else {
528              order_str = "dp";
529           }
530        }
531        else {
532           order_str = "dp";
533        }
534      }
535     }
536    }
537   }
538   else {
539     order_str = "dp";
540   }
541   list result;
542   result[1] = order_str;
543   result[2] = curr_weight;
544   result[3] = target_weight;
545   return(result);
546}
547
548//////////////////////////////////////////////////////////////////////////////
549static proc OrderStringalp_NP(string Wpal,list #)
550{
551  int n= nvars(basering);
552  string order_str = "dp";
553  int nP = 1;// call LatsGB to compute the wanted GB  by pwalk
554  intvec curr_weight = system("Mivdp",n); //define (1,1,...,1)
555  intvec target_weight = system("Mivlp",n); //define (1,0,...,0)
556  if(size(#) != 0)
557  {
558    if(size(#) == 1)
559    {
560      if(typeof(#[1]) == "intvec")
561      {
562        curr_weight = #[1];
563
564        if(Wpal == "al")
565        {
566          order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
567        }
568        else
569        {
570          order_str = "(Wp("+string(#[1])+"),C)";
571        }
572      }
573      else {
574        if(typeof(#[1]) == "int")
575        {
576          nP = #[1];
577        }
578        else
579        {
580          print("// ** the input must be \"(ideal, int)\" or ");
581          print("// **                   \"(ideal, intvec)\"");
582          print("// ** a lex. GB will be computed from \"dp\" to \"lp\"");
583        }
584      }
585    }
586    else
587    {
588     if(size(#) == 2)
589     {
590       if(typeof(#[1]) == "intvec" and typeof(#[2]) == "int")
591       {
592         curr_weight = #[1];
593
594         if(Wpal == "al")
595         {
596           order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
597         }
598         else
599         {
600           order_str = "(Wp("+string(#[1])+"),C)";
601         }
602       }
603       else
604       {
605         if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec")
606         {
607           curr_weight = #[1];
608           target_weight = #[2];
609
610           if(Wpal == "al")
611           {
612             order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
613           }
614           else
615           {
616             order_str = "(Wp("+string(#[1])+"),C)";
617           }
618         }
619         else
620         {
621           print("// ** the input  must be \"(ideal,intvec,int)\" or ");
622           print("// **                    \"(ideal,intvec,intvec)\"");
623           print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");
624         }
625       }
626     }
627     else {
628       if(size(#) == 3)
629       {
630         if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec" and
631            typeof(#[3]) == "int")
632         {
633           curr_weight = #[1];
634           target_weight = #[2];
635           nP = #[3];
636           if(Wpal == "al")
637           {
638             order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
639           }
640           else
641           {
642             order_str = "(Wp("+string(#[1])+"),C)";
643           }
644         }
645         else
646         {
647           print("// ** the input must be \"(ideal,intvec,intvec,int)\"");
648           print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");
649
650         }
651       }
652       else
653       {
654         print("// ** The given input is wrong");
655         print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");
656       }
657     }
658    }
659  }
660  list result;
661  result[1] = nP;
662  result[2] = order_str;
663  result[3] = curr_weight;
664  result[4] = target_weight;
665  return(result);
666}
667//////////////////////////////////////////////////////////////////////////////
668static proc test_in_cone(vector w, list l)
669{
670 int i,j,k;
671 vector v;
672 poly n;
673 number a;
674 for(i=1;i<=size(l);i++)
675 {
676    for(j=1;j<=size(l[i]);j++)
677    {
678        v=0;
679        for(k=1;k<=size(l[i][j]);k++)
680        {
681            v=v+l[i][j][k]*gen(k);
682        }
683        n = inprod(w,v);
684        a = leadcoef(n);
685        if(a<0)
686        {
687           return(0);
688        }
689    }
690 }
691 return(1);
692}
693//////////////////////////////////////////////////////////////////////////////
694static proc PertVectors(ideal Gold, vector target_weight, int pdeg)
695{
696int nV = nvars(basering);
697int nG = size(Gold);
698int i;
699number ntemp, maxAi, maxA;
700if(pdeg > nV || pdeg <= 0)
701  {
702    intvec v_null=0;
703    return v_null;
704  }
705if(pdeg == 1)
706  {
707    return target_weight;
708  }
709maxAi=0;
710for(i=1; i<=nV; i++)
711  {
712    ntemp = leadcoef(inprod(target_weight,gen(i)));
713    if(ntemp < 0)
714      {
715        ntemp = -ntemp;
716      }
717    if(maxAi < ntemp)
718      {
719        maxAi = ntemp;
720      }
721  }
722maxA = maxAi+pdeg-1;
723number epsilon = maxA*deg(Gold)+1;
724vector pert_weight = epsilon^(pdeg-1)*target_weight;
725for(i=2; i<=pdeg; i++)
726  {
727    pert_weight = pert_weight + epsilon^(pdeg-i)*gen(i);
728  }
729return(pert_weight);
730}
731
732
733//////////////////////////////////////////////////////////////////////////////
734static proc RandomNextWeight(ideal Gold, list L, vector curr_weight,
735                      vector target_weight,int weight_rad, int pdeg)
736"USAGE: RandomNextWeight(Gold, L, curr_weight, target_weight);
737RETURN: Intermediate next weight vector
738EXAMPLE: example RandomNextWeight; shows an example
739"
740{
741  int i,n1,n2,n3;
742  number norm, weight_norm;
743  def Rold = basering;
744  int nV = nvars(basering);
745  number ulast=last(curr_weight, target_weight, L);
746  vector new_weight=(1-ulast)*curr_weight+ulast*target_weight;
747  vector w1=cleardenom(new_weight);
748  intvec v1=changeType(w1);
749  list p= ringlist(Rold);
750  p[3][1][2]= v1;
751  def new_ring=ring(p);
752  setring new_ring;
753  ideal Gold = fetch(Rold, Gold);
754  n1=size(Initial(Gold));
755  setring Rold;
756  intvec next_weight;
757  kill new_ring;
758  while(1)
759  {
760    weight_norm = 0;
761    while(weight_norm == 0)
762    {
763      for(i=1; i<=nV; i++)
764      {
765        next_weight[i] = random(1,10000)-5000;
766        weight_norm = weight_norm + next_weight[i]^2;
767      }
768      norm = 0;
769      while(norm^2 < weight_norm)
770      {
771         norm=norm+1;
772      }
773      weight_norm = 1+norm;
774    }
775    new_weight = 0;
776    for(i=1; i<=nV;i++)
777    {
778      if(next_weight[i] < 0)
779      {
780        new_weight = new_weight + (1 + round(weight_rad*leadcoef(next_weight[i])/weight_norm))*gen(i);
781      }
782      else
783      {
784        new_weight = new_weight + ( round(weight_rad*leadcoef(next_weight[i])/weight_norm))*gen(i);
785      }
786    }
787    new_weight = new_weight + curr_weight;
788    if(test_in_cone(new_weight, L)==1)
789    {
790      break;
791    }
792  }
793  kill next_weight;
794  kill norm;
795  vector w2=cleardenom(new_weight);
796  intvec v2=changeType(w2);
797  p[3][1][2]= v2;
798  def new_ring=ring(p);
799  setring new_ring;
800  ideal Gold = fetch(Rold, Gold);
801  n2=size(Initial(Gold));
802  setring Rold;
803  kill new_ring;
804
805  vector w3=cleardenom(PertVectors(Gold,target_weight,pdeg));
806  intvec v3=changeType(w3);
807  p[3][1][2]= v1;
808  def new_ring=ring(p);
809  setring new_ring;
810  ideal Gold = fetch(Rold, Gold);
811  n3=size(Initial(Gold));
812  setring Rold;
813  kill new_ring;
814  kill p;
815
816  if(n2<n1)
817  {
818    if(n3<n2)
819    {
820      // n3<n2<n1
821      return(w3);
822    }
823    else
824    {
825      // n2<n1 und n2<=n3
826      return(w2);
827    }
828  }
829  else
830  {
831    if(n3<n1)
832    {
833      //n3<n1<=n2
834      return(w3);
835    }
836    else
837    {
838      // n1<=n3 und n1<=n2
839      return(w1);
840    }
841  }
842}
843
844///////////////////////////////////////////////////////////////////////////////*
845Further examples
846ring r=0,(x,y,z),lp;
847
848ideal I=x2y4, y4z2, xy4z+y2z, xy6z2+y10z5;
849
850ideal Q=x2y4, y4z2, xy4z+y2z, xy6z2+y14z7;
851
852ideal J=x2y4, y4z2, xy4z+y2z, xy6z2+y18z9;
853
854ideal K=x2,y2,xy+y,2xy2+y5,z3+x;
855
856ideal L=x2+y,y2+z,x+z2;
857*/
858
859
Note: See TracBrowser for help on using the repository browser.