source: git/Singular/LIB/JMBTest.lib @ fcc6e7

spielwiese
Last change on this file since fcc6e7 was fcc6e7, checked in by Hans Schoenemann <hannes@…>, 5 years ago
format
  • Property mode set to 100644
File size: 20.3 KB
Line 
1//////////////////////////////////////////////////////////////////////
2version="version JMBTest.lib 4.1.2.0 Feb_2019 "; // $Id$
3category="Algebraic Geometry";
4// summary description of the library
5info="
6LIBRARY:   JMBTest.lib  A library for Singular which performs JM basis test.
7AUTHOR:    Michela Ceria, email: michela.ceria@unito.it
8
9SEE ALSO:  JMSConst_lib
10KEYWORDS:  J-marked schemes
11
12OVERVIEW:
13The library performs the J-marked basis test, as described in [CR], [BCLR].
14Such a test is performed via the criterion explained in [BCLR],
15concerning Eliahou-Kervaire polynomials (EK from now on).
16We point out that all the polynomials are homogeneous
17and they must be arranged by degree.
18The fundamental steps are the following:@*
19-construct the Vm polynomials, via the algorithm VConstructor
20explained in [CR];@*
21-construct the Eliahou-Kervaire polynomials defined in [BCLR];@*
22-reduce the  Eliahou-Kervaire polynomials using the Vm's;@*
23-if it exist an Eliahou-Kervaire polynomial such that its reduction
24mod Vm is different from zero, the given one is not a J-Marked basis.
25
26The algorithm terminates only if the ordering is rp.
27Anyway, the number of reduction steps is bounded.
28
29REFERENCES:
30[CR] Francesca Cioffi, Margherita Roggero,Flat Families by Strongly
31Stable Ideals and a Generalization of Groebner Bases,
32J. Symbolic Comput. 46,  1070-1084, (2011).@*
33[BCLR] Cristina Bertone, Francesca Cioffi, Paolo Lella,
34Margherita Roggero, Upgraded methods for the effective
35computation of marked schemes on a strongly stable ideal,
36Journal of Symbolic Computation
37(2012), http://dx.doi.org/10.1016/j.jsc.2012.07.006 @*
38
39PROCEDURES:
40  Minimus(ideal)            minimal variable in an ideal
41  Maximus(ideal)               maximal variable in an ideal
42  StartOrderingV(list,list)     ordering of polynomials as in [BCLR]
43  TestJMark(list)           tests whether we have a J-marked basis
44";
45LIB "qhmoduli.lib";
46LIB "monomialideal.lib";
47LIB "ring.lib";
48////////////////////////////////////////////////////////////////////
49static proc mod_init()
50"USAGE:     mod_init();
51RETURN:   struct: jmp
52EXAMPLE:  example  mod_init; shows an example"
53{
54newstruct("jmp", "poly h, poly t");
55}
56example
57{ "EXAMPLE:"; echo = 2;
58  mod_init();
59}
60////////////////////////////////////////////////////////////////////
61proc Terns(list G, int c)
62"USAGE:    Terns(G,c); G list, c int
63RETURN:   list: T
64NOTE:     Input is a list of J-marked  polynomials
65(arranged by degree) and an integer.
66EXAMPLE:  example Terns; shows an example"
67{
68list T=list();
69int z;
70 for(int k=1; k<=size(G[c]);k=k+1)
71    {
72//Loop on G[c] making positions of polynomials in G[c]
73      z=size(T);
74      T=insert(T,list(1,c,k) ,size(T));
75    }
76return(T);
77}
78example
79{ "EXAMPLE:"; echo = 2;
80ring r=0, (x,y,z), rp;
81jmp r1;
82r1.h=z^3;
83r1.t=poly(0);
84jmp r2;
85r2.h=z^2*y;
86r2.t=poly(0);
87jmp r3;
88r3.h=z*y^2 ;
89r3.t=-x^2*y;
90jmp r4;
91r4.h=y^5;
92r4.t=poly(0);
93list G2F=list(list(r1,r2,r3),list(r4));
94Terns(G2F, 1);
95Terns(G2F, 2);
96}
97////////////////////////////////////////////////////////////////////
98proc VConst(list G, int c)
99"USAGE:    VConst(G, c); G list, c int
100RETURN:   list: V
101NOTES: this procedure computes the Vm polynomials following the
102algorithm in [CR],but it only keeps in memory the monomials by
103which the G's must be multplied and their positions.
104EXAMPLE:  example VConst; shows an example"
105{
106jmp f=G[1][1];
107int aJ=deg(f.h);
108// minimal degree of polynomials in G
109list V=list();
110V[1]=Terns(G,1);
111// V[1]=G[1] (keeping in memory only [head, position])
112//print(c-aJ+1);
113int i;
114int j;
115int m;
116list OO;
117jmp p;
118  for(m=2; m<=c-aJ+1; m=m+1)
119  {
120     //print("entro nel form");
121     if(m>size(G))
122        {V[m]=list();
123//If we have not G[m] we insert a list()
124        //print("vuota prima");
125        }
126     else
127        {V[m]=Terns(G,m);
128        //print("piena prima");
129        }
130     for(i=1; i<nvars(basering)+1; i=i+1)
131        {
132          //print("entrata fori");
133          //print(i);
134          for(j=1; j<=size(V[m-1]); j=j+1)
135             {
136                   p=G[V[m-1][j][2]][V[m-1][j][3]];
137                  //print(p.h);
138                  //print(p.t);
139                  //print(var(i));
140                  //print(Minimus(V[m-1][j][1]*p.h));
141               if(var(i)<=Minimus(variables(V[m-1][j][1]*p.h)))
142                  {
143//Can I multiply by the current variable?
144                    //print("minoremin");
145                    //print("fin qui ci sono");
146                     //print(V[m-1][j][1]);
147                     OO=list(var(i)*V[m-1][j][1],V[m-1][j][2],V[m-1][j][3]);
148                    V[m]=insert(V[m], OO ,size(V[m]));
149                  }
150             }
151        }
152  }
153return (V);}
154example
155{ "EXAMPLE:"; echo = 2;
156 ring r=0, (x,y,z), rp;
157jmp r1;
158r1.h=z^3;
159r1.t=poly(0);
160jmp r2;
161r2.h=z^2*y;
162r2.t=poly(0);
163jmp r3;
164r3.h=z*y^2 ;
165r3.t=-x^2*y;
166jmp r4;
167r4.h=y^5;
168r4.t=poly(0);
169list G2F=list(list(r1,r2,r3),list(r4));
170VConst(G2F,4,basering);}
171////////////////////////////////////////////////////////////////////
172proc Minimus(ideal L)
173"USAGE:    Minimus(L); G list, c int
174RETURN:   list: V
175NOTES: it returns the minimal variable generating the ideal L.@*
176The input must be an ideal generated by variables.
177EXAMPLE:  example Minimus; shows an example"
178{
179poly min=L[1];
180int i;
181for(i=2;i<=size(L); i++)
182  {
183    if(L[i]<min){min=L[i];}
184  }
185return(min);
186}
187example
188{ "EXAMPLE:"; echo = 2;
189 ring r=0, (x,y,z), rp;
190ideal I=y,x,z;
191Minimus(I);
192}
193////////////////////////////////////////////////////////////////////
194proc Maximus(ideal L)
195"USAGE:    Maximus(L); G list, c int
196RETURN:   list: V
197NOTES: it returns the maximal variable generating the ideal L.@*
198The input must be an ideal generated by variables.
199EXAMPLE:  example Maximus; shows an example"
200{
201poly max=L[1];
202int i;
203for(i=2;i<=size(L); i++)
204  {
205    if(L[i]>max){max=L[i];}
206  }
207return(max);
208}
209example
210{ "EXAMPLE:"; echo = 2;
211 ring r=0, (x,y,z), rp;
212ideal I=y,x,z;
213Maximus(I);
214}
215////////////////////////////////////////////////////////////////////
216proc GJmpMins(jmp P, jmp Q)
217"USAGE:    GJmpMins(P,Q); P jmp, Q jmp
218RETURN:   int: d
219EXAMPLE:  example GJmpMins; shows an example"
220{
221  int d=1;
222//-1=lower, 0=equal, 1=higher
223//At the beginning suppose Q is higher
224  if(deg(P.h)<deg(Q.h))
225     {
226//Compare degrees;
227      d=-1;
228     }
229  if(deg(P.h)==deg(Q.h))
230     {
231      if(P.h==Q.h)
232         {
233          if(P.t==Q.t)
234              {
235//head=tail
236               d=0;
237              }
238         }
239     else
240        {
241//print(Minimus(variables(P.h/gcdMon(P.h,Q.h))));
242//print(Minimus(variables(Q.h/gcdMon(P.h,Q.h))));
243          if(Minimus(variables(P.h/gcdMon(P.h,Q.h)))<Minimus(variables(Q.h/gcdMon(P.h,Q.h))))
244             {
245             d=-1;
246            //print("Per Indice");
247             }
248        }
249     }
250 return(d);
251}
252example
253{ "EXAMPLE:"; echo = 2;
254 ring r=0, (x,y,z), rp;
255jmp p1;
256p1.h=poly(1);
257p1.t=poly(1);
258jmp p2;
259p2.h=x^2;
260p2.t=poly(0);
261jmp p3;
262p3.h=x;
263p3.t=poly(0);
264 GJmpMins(p1, p2);
265 GJmpMins(p2, p3);
266GJmpMins(p1,p1);
267}
268////////////////////////////////////////////////////////////////////
269proc TernCompare(list A, list B, list G)
270"USAGE:    TernCompare(A,B,C); A list, B list, G list
271RETURN:   int: d
272NOTE:     A and B are terns, while G is the given list of
273J-marked polynomials.
274EXAMPLE:  example TernCompare; shows an example"
275{
276int d=-1;
277//Start: A<B
278if(A[1]==B[1])
279  {
280   if(A[2]==B[2]&& A[3]==B[3])
281     {
282      //print("Uguali");
283      d=0;
284     }
285   else
286     {
287jmp g1=G[A[2]][A[3]];
288jmp g2=G[B[2]][B[3]];
289       if(GJmpMins(g1, g2)==1)
290          {
291            //print("Maggiore per il G");
292            d=1;
293          }
294     }
295  }
296else
297  {
298    if(A[1]>B[1])
299      {
300//the ordering MUST be rp
301        //print("Maggiore per Lex");
302        d=1;
303      }
304  }
305return(d);
306}
307example
308{ "EXAMPLE:"; echo = 2;
309 ring r=0, (x,y,z), rp;
310jmp r1;
311r1.h=z^3;
312r1.t=poly(0);
313jmp r2;
314r2.h=z^2*y;
315r2.t=poly(0);
316jmp r3;
317r3.h=z*y^2 ;
318r3.t=-x^2*y;
319jmp r4;
320r4.h=y^5;
321r4.t=poly(0);
322list G2F=list(list(r1,r2,r3),list(r4));
323TernCompare([1,1,1],[x,1,1],G2F);
324}
325////////////////////////////////////////////////////////////////////
326proc MinOfV(list V, list G)
327"USAGE:    Minimal(V,G); V list, G list
328RETURN:   int: R
329NOTE:     Input=lista(terne), G.
330EXAMPLE:  example Minimal; shows an example"
331{
332//Minimal element for a given degree
333list R=list();
334list MIN=V[1];
335int h=1;
336int i;
337for(i=2; i<=size(V); i++)
338    {
339//I consider the first as minimum
340//If I find something smaller I change minimum
341     if(TernCompare(V[i],MIN,G)<=0)
342       {
343        MIN=V[i];
344        h=i;
345       }
346    }
347//Return: [minimum,position of the minimum]
348R=MIN,h;
349return(R);
350}
351example
352{ "EXAMPLE:"; echo = 2;
353 ring r=0, (x,y,z), rp;
354jmp r1;
355r1.h=z^3;
356r1.t=poly(0);
357jmp r2;
358r2.h=z^2*y;
359r2.t=poly(0);
360jmp r3;
361r3.h=z*y^2 ;
362r3.t=-x^2*y;
363jmp r4;
364r4.h=y^5;
365r4.t=poly(0);
366list G2F=list(list(r1,r2,r3),list(r4));
367MinOfV(VConst(G2F,4,basering)[1],G2F);
368}
369////////////////////////////////////////////////////////////////////
370proc OrderingV(list V,list G,list R)
371"USAGE:   OrderingV(V,G,R); V list, G list, R list
372RETURN:   list: R
373NOTE:     Input: Vm,G,emptylist
374EXAMPLE:  example OrderingV; shows an example"
375{
376//Order V[m]
377//R  will contain results but at the beginning it is empty
378list M=list();
379if(size(V)==1)
380  {
381    R=insert(R,V[1],size(R));
382  }
383else
384  {
385   M=MinOfV(V,G);
386   R=insert(R,M[1],size(R));
387   V=delete(V,M[2]);
388//recursive call
389   R=OrderingV(V,G,R);
390  }
391return(R);
392}
393example
394{ "EXAMPLE:"; echo = 2;
395 ring r=0, (x,y,z), rp;
396jmp r1;
397r1.h=z^3;
398r1.t=poly(0);
399jmp r2;
400r2.h=z^2*y;
401r2.t=poly(0);
402jmp r3;
403r3.h=z*y^2;
404r3.t=-x^2*y;
405jmp r4;
406r4.h=y^5;
407r4.t=poly(0);
408list G2F=list(list(r1,r2,r3),list(r4));
409OrderingV(VConst(G2F,4,basering)[1],G2F,list());
410}
411////////////////////////////////////////////////////////////////////
412proc StartOrderingV(list V,list G)
413"USAGE:   StartOrdina(V,G); V list, G list
414RETURN:   list: R
415NOTE:     Input Vm,G. This procedure uses OrderingV to get
416the ordered polynomials as in [BCLR].
417EXAMPLE:  example StartOrderingV; shows an example"
418{
419 return(OrderingV(V,G, list()));
420}
421example
422{ "EXAMPLE:"; echo = 2;
423 ring r=0, (x,y,z), rp;
424jmp r1;
425r1.h=z^3;
426r1.t=poly(0);
427jmp r2;
428r2.h=z^2*y;
429r2.t=poly(0);
430jmp r3;
431r3.h=z*y^2;
432r3.t=-x^2*y;
433jmp r4;
434r4.h=y^5;
435r4.t=poly(0);
436list G2F=list(list(r1,r2,r3),list(r4));
437StartOrderingV(VConst(G2F,4,basering)[1],G2F);
438}
439////////////////////////////////////////////////////////////////////
440proc Multiply(list L, list G)
441"USAGE:    moltiplica(L,G); L list, G list
442RETURN:   jmp: K
443NOTE:     Input: a 3-ple,G. It performs the product associated
444to the 3-uple.
445EXAMPLE:  example Multiply; shows an example"
446{
447jmp g=G[L[2]][L[3]];
448jmp K;
449K.h=L[1]*g.h;
450K.t=L[1]*g.t;
451 return(K);
452}
453example
454{ "EXAMPLE:"; echo = 2;
455 ring r=0, (x,y,z), rp;
456 list P=x^2,1,1;
457jmp r1;
458r1.h=z^3;
459r1.t=poly(0);
460jmp r2;
461r2.h=z^2*y;
462r2.t=poly(0);
463jmp r3;
464r3.h=z*y^2 ;
465r3.t=-x^2*y;
466jmp r4;
467r4.h=y^5;
468r4.t=poly(0);
469list G2F=list(list(r1,r2,r3),list(r4));
470Multiply(P,G2F);
471}
472////////////////////////////////////////////////////////////////////
473proc IdealOfV(list V)
474"USAGE:    IdealOfV(V); V list
475RETURN:   ideal: I
476NOTES: this procedure takes a list of Vm's of a certain degree
477and construct their ideal, multiplying the head by the weighted
478variable t.
479EXAMPLE:  example IdealOfV; shows an example"
480{
481ideal I=0;
482int i;
483if (size(V)!=0)
484  {
485   list M=list();
486jmp g;
487   for(i=1; i<= size(V); i++)
488       {
489         g=V[i];
490         g.h=t*g.h;
491         M[i]=g.h+g.t;
492       }
493   I=M[1..size(M)];
494  //I=std(I);
495  }
496return(I);
497}
498example
499{ "EXAMPLE:"; echo = 2;
500 ring r=0, (x,y,z,t), rp;
501jmp r1;
502r1.h=z^3;
503r1.t=poly(0);
504jmp r2;
505r2.h=z^2*y;
506r2.t=poly(0);
507jmp r3;
508r3.h=z*y^2 ;
509r3.t=-x^2*y;
510jmp r4;
511r4.h=y^5;
512r4.t=poly(0);
513list G2F=list(list(r1,r2,r3),list(r4));
514IdealOfV(G2F[1]);
515}
516////////////////////////////////////////////////////////////////////
517proc NewWeight(int n)
518"USAGE:    NewWeight(n); n int
519RETURN:   intvec: u
520EXAMPLE:  example NewWeight; shows an example"
521{
522intvec u=0;
523u[n]=1;
524return(u);
525}
526example
527{ "EXAMPLE:"; echo = 2;
528 NewWeight(3);
529}
530////////////////////////////////////////////////////////////////////
531proc FinalVm(list V1 , list G1 ,def r)
532"USAGE:     FinalVm(V1, G1, r);  V1 list,  G1 list , r
533RETURN:   intvec: u
534EXAMPLE:  example NewWeight; shows an example"
535{
536  //multiply and reduce, degree by degree
537  intvec u=NewWeight(nvars(r)+1);
538  list L=ringlist(r);
539  L[2]=insert(L[2],"t",size(L[2]));
540  list ordlist="a",u;
541  L[3]=insert(L[3],ordlist,0);
542  def H=ring(L);
543  list M=list();
544  jmp p;
545  list N;
546  poly q;
547  poly s;
548  int i;
549  int j;
550  for(i=1; i<=size(G1); i++)
551  {
552    N=list();
553    for(j=1; j<=size(G1[i]); j++)
554    {
555      p=G1[i][j];
556      q=p.h;
557      s=p.t;
558      N[j]=list(q,s);
559    }
560    M[i]=N;
561  }
562  p.h=poly(0);
563  p.t=poly(0);
564  setring H;
565  list R=list();
566  list S=list();
567  //print("anello definito");
568  def V=imap(r,V1);
569  //def G=imap(r,G1);
570  //print(V);
571  def MM=imap(r,M);
572  list G=list();
573  list N=list();
574  for(i=1; i<=size(MM); i++)
575  {
576    for(j=1; j<=size(MM[i]); j++)
577    {
578      p.h=MM[i][j][1];
579      p.t=MM[i][j][2];
580      N[j]=p;
581    }
582    G[i]=N;
583  }
584  ideal I=0;
585  jmp LL;
586  jmp UU;
587  for(i=1; i<=size(V);i++)
588  {
589    R[i]=list();
590    S[i]=list();
591    I=0;
592    for(j=1;j<=size(V[i]); j++)
593    {
594      LL=Multiply(V[i][j],G);
595      LL.t=reduce(t*LL.t,I);
596      //I only reduce the tail
597      LL.t=subst(LL.t,t,1);
598      S[i]=insert(S[i],LL,size(S[i]));
599      LL.h=t*LL.h;
600      R[i]=insert(R[i],LL,size(R[i]));
601      UU=R[i][j];
602      I=I+ideal(UU.h+UU.t);
603      attrib(I,"isSB",1);
604    }
605  }
606  list M=list();
607  poly q;
608  poly s;
609  for(i=1; i<=size(S); i++)
610  {
611    N=list();
612    for(j=1; j<=size(S[i]); j++)
613    {
614      p=S[i][j];
615      q=p.h;
616      s=p.t;
617      N[j]=list(q,s);
618    }
619    M[i]=N;
620  }
621  p.h=poly(0);
622  p.t=poly(0);
623  setring r;
624  def MM=imap(H,M);
625  list MMM=list();
626  for(i=1; i<=size(MM); i++)
627  {
628    N=list();
629    for(j=1; j<=size(MM[i]); j++)
630    {
631      p.h=MM[i][j][1];
632      p.t=MM[i][j][2];
633      N[j]=p;
634    }
635    MMM[i]=N;
636  }
637  return(MMM);
638}
639example
640{ "EXAMPLE:"; echo = 2;
641 ring r=0, (x,y,z), rp;
642 jmp r1;
643 r1.h=z^3;
644 r1.t=poly(0);
645 jmp r2;
646 r2.h=z^2*y;
647 r2.t=poly(0);
648 jmp r3;
649 r3.h=z*y^2 ;
650 r3.t=-x^2*y;
651 jmp r4;
652 r4.h=y^5;
653 r4.t=poly(0);
654 list G2F=list(list(r1,r2,r3),list(r4));
655 FinalVm(VConst(G2F,6,r) , G2F, r);
656}
657////////////////////////////////////////////////////////////////////
658proc ConstructorMain(list G, int c,def r)
659"USAGE:    Costruttore(G,c); G list, c int
660RETURN:   list: R
661NOTE:     At the end separated by degree.
662EXAMPLE:  example Costruttore; shows an example"
663{
664list V=list();
665V= VConst(G,c);
666//V non ordered
667list L=list();
668list R=list();
669int i;
670// head, position
671//order the different degrees
672for(i=1; i<=size(V); i++)
673   {
674    L[i]=StartOrderingV(V[i], G);
675   }
676//multiply and reduce
677//print("Ordinare");
678R=FinalVm(L, G, r);
679//print("FinalVm");
680return(R);
681}
682example
683{ "EXAMPLE:"; echo = 2;
684 ring r=0, (x,y,z), rp;
685jmp r1;
686r1.h=z^3;
687r1.t=poly(0);
688jmp r2;
689r2.h=z^2*y;
690r2.t=poly(0);
691jmp r3;
692r3.h=z*y^2 ;
693r3.t=-x^2*y;
694jmp r4;
695r4.h=y^5;
696r4.t=poly(0);
697list G2F=list(list(r1,r2,r3),list(r4));
698ConstructorMain(G2F,6,r);
699}
700////////////////////////////////////////////////////////////////////
701proc EKCouples(jmp A, jmp B)
702"USAGE:    CoppiaEK(A,B); A list, B list
703RETURN:   list: L
704NOTE:     At the end the monomials involved by EK.
705EXAMPLE:  example EKCouples; shows an example"
706{
707poly E;
708list L=0,0;
709string s=varstr(basering);
710list VVV=varstr(basering);
711//L will contain results
712poly h=Minimus(variables(A.h));
713int l=findvars(h,1)[2][1];
714if(l!=nvars(basering))
715 {
716  for(int j=l+1;j<=nvars(basering); j++)
717   {
718    E=var(j)*A.h/B.h;
719//Candidate for * product
720    if(E!=0)
721      {
722        //print("primo if passato");
723        if(Minimus(variables(B.h))>=Maximus(variables(E)))
724           {
725//Does it work with * ?
726            //print("secondo if passato");
727            L[1]=j;
728            L[2]=E;
729            break;
730           }
731      }
732   }
733  }
734return (L);
735}
736example
737{ "EXAMPLE:"; echo = 2;
738 ring r=0, (x,y,z), rp;
739jmp A;
740A.h=y*z^2;
741A.t=poly(0);
742jmp B;
743B.h=y^2*z;
744B.t=poly(0);
745EKCouples(A,B);
746EKCouples(B,A);
747}
748////////////////////////////////////////////////////////////////////
749proc EKPolys(list G)
750"USAGE:    PolysEK(G); G list
751RETURN:   list: EK, list: D
752NOTE:     At the end EK polynomials and their degrees
753
754EXAMPLE:  example PolysEK; shows an example"
755{
756list D=list();
757list C=list();
758list N=0,0;
759list EK=list();
760int i;
761int j;
762int k;
763int l;
764jmp p;
765for(i=1; i<=size(G); i++)
766   {
767     for(j=1; j<=size(G[i]); j++)
768        {
769         for(k=1; k<=size(G); k++)
770             {
771               for(l=1; l<=size(G[k]); l++)
772                  {
773                     if(i!=k||j!=l)
774                       {
775//Loop on polynomials
776                        C=EKCouples(G[i][j], G[k][l]);
777                        if(C[2]!=0)
778                          {
779                           C=insert(C,list(i,j,k,l),size(C));
780                           EK=insert(EK,C,size(EK));
781                           p=G[k][l];
782                           D=insert(D,deg(C[2]*p.h),size(D));
783                          }
784                       }
785                  }
786             }
787        }
788   }
789//Double Return
790return(EK, D);
791}
792example
793{ "EXAMPLE:"; echo = 2;
794  ring r=0, (x,y,z), rp;
795jmp r1;
796r1.h=z^3;
797r1.t=poly(0);
798jmp r2;
799r2.h=z^2*y;
800r2.t=poly(0);
801jmp r3;
802r3.h=z*y^2;
803r3.t=-x^2*y;
804jmp r4;
805r4.h=y^5;
806r4.t=poly(0);
807list G2F=list(list(r1,r2,r3),list(r4));
808EKPolys(G2F);
809}
810////////////////////////////////////////////////////////////////////
811proc EKPolynomials(list EK, list G)
812"USAGE:     EKPolynomials(EK,G); EK list, G list
813RETURN:   list: p
814NOTE:     At the end I obtain the EK polynomials and
815their degrees.
816EXAMPLE:  example SpolyEK; shows an example"
817{
818jmp u=G[EK[3][1]][EK[3][2]];
819jmp q=G[EK[3][3]][EK[3][4]];
820return(var(EK[1])*(u.h+u.t)-EK[2]*(q.h+q.t));
821}
822example
823{ "EXAMPLE:"; echo = 2;
824 ring r=0, (x,y,z), rp;
825jmp r1;
826r1.h=z^3;
827r1.t=poly(0);
828jmp r2;
829r2.h=z^2*y;
830r2.t=poly(0);
831jmp r3;
832r3.h=z*y^2;
833r3.t=-x^2*y;
834jmp r4;
835r4.h=y^5;
836r4.t=poly(0);
837list G2F=list(list(r1,r2,r3),list(r4));
838list EK,D=EKPolys(G2F);
839EKPolynomials(EK[1],G2F);
840}
841////////////////////////////////////////////////////////////////////
842proc TestJMark(list G1,def r)
843"USAGE:    TestJMark(G);  G list
844RETURN:    int: i
845NOTE:
846This procedure performs J-marked basis test.@*
847The input is a list of J-marked polynomials (jmp) arranged
848by degree, so G1 is a list of list.@*
849The output is a boolean evaluation:
850True=1/False=0
851EXAMPLE:  example TestJMark; shows an example"
852{
853  int flag=1;
854  if(size(G1)==1 && size(G1[1])==1)
855  {
856    //Hypersurface
857    print("Only One Polynomial");
858    flag=1;
859  }
860  else
861  {
862    int d=0;
863    list EK,D=EKPolys(G1);
864    //I found EK couples
865    int massimo=Max(D);
866    list V1=ConstructorMain(G1,massimo,r);
867    jmp mi=V1[1][1];
868    int minimo=Min(deg(mi.h));
869    intvec u=NewWeight(nvars(r)+1);
870    list L=ringlist(r);
871    L[2]=insert(L[2],"t",size(L[2]));
872    list ordlist="a",u;
873    L[3]=insert(L[3],ordlist,0);
874    def H=ring(L);
875    list JJ=list();
876    jmp pp;
877    jmp qq;
878    int i;
879    int j;
880    list NN;
881    for(i=size(V1);i>0;i--)
882    {
883      NN=list();
884      for(j=size(V1[i]);j>0;j--)
885      {
886        pp=V1[i][j];
887        NN[j]=list(pp.h,pp.t);
888      }
889      JJ[i]=NN;
890    }
891    list KK=list();
892    list UU=list();
893    //jmp qq;
894    for(i=size(G1);i>0;i--)
895    {
896      for(j=size(G1[i]);j>0;j--)
897      {
898        qq=G1[i][j];
899        UU[j]=list(qq.h,qq.t);
900      }
901      KK[i]=UU;
902    }
903    setring H;
904    //I defined the new ring with the weighted
905    //variable t
906    poly p;
907    //print("anello definito");
908    def JJJ=imap(r,JJ);
909    def EK=imap(r,EK);
910    list V=list();
911    jmp fp;
912    //int i;
913    //int j;
914    list N;
915    for(i=size(JJJ); i>0; i--)
916    {
917      N=list();
918      for(j=size(JJJ[i]); j>0; j--)
919      {
920        fp.h=JJJ[i][j][1];
921        fp.t=JJJ[i][j][2];
922        N[j]=fp;
923      }
924      V[i]=N;
925    }
926    def KKJ=imap(r,KK);
927    list G=list();
928    list U=list();
929    for(i=1; i<=size(KKJ); i++)
930    {
931      for(j=1; j<=size(KKJ[i]); j++)
932      {
933        fp.h=KKJ[i][j][1];
934        fp.t=KKJ[i][j][2];
935        U[j]=fp;
936      }
937      G[i]=U;
938    }
939    //I imported in H everithing I need
940    poly q;
941    ideal I;
942    for(j=1; j<=size(EK);j++)
943    {
944      d=D[j];
945      p=EKPolynomials(EK[j],G);
946      I=IdealOfV(V[d-minimo+1]);
947      attrib(I,"isSB",1);
948      q=reduce(t*p,I);
949      q=subst(q,t,1);
950      //I reduce all the EK polynomials
951      //       q=RiduzPoly(V[d-minimo+1], p);
952      if(q!=0)
953      {
954        //check whether reduction is 0
955        print("NOT A BASIS");
956        flag=0;
957        break;
958      }
959    }
960  }
961  setring r;
962  //typeof(flag);
963  return(flag);
964}
965example
966{ "EXAMPLE:"; echo = 2;
967 ring r=0, (x,y,z), rp;
968 jmp r1;
969 r1.h=z^3;
970 r1.t=poly(0);
971 jmp r2;
972 r2.h=z^2*y;
973 r2.t=poly(0);
974 jmp r3;
975 r3.h=z*y^2 ;
976 r3.t=-x^2*y;
977 jmp r4;
978 r4.h=y^5;
979 r4.t=poly(0);
980 list G2F=list(list(r1,r2,r3),list(r4));
981 TestJMark(G2F,r);
982}
Note: See TracBrowser for help on using the repository browser.