source: git/Singular/LIB/JMBTest.lib

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