source: git/Singular/LIB/classify_aeq.lib @ f9b7f98

spielwiese
Last change on this file since f9b7f98 was f9b7f98, checked in by Hans Schoenemann <hannes@…>, 8 years ago
ifix by Gerhard (classify_aeq.lib)
  • Property mode set to 100644
File size: 58.5 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version classify_aeq.lib 4.0.0.0 Jun_2013 "; // $Id$
3category="Singularities";
4info="
5LIBRARY: classifyAeq.lib         Simple Space Curve singularities in characteristic 0
6
7AUTHORS: Faira Kanwal Janjua     fairakanwaljanjua@gmail.com
8         Gerhard Pfister         pfister@mathematik.uni-kl.de
9         Khawar Mehmood          khawar1073@gmail.com                         NEU
10
11OVERVIEW: A library for classifying the simple singularities
12          with respect to A equivalence in characteristic 0.
13  Simple Surface singularities in characteristic O have been classified by Bruce and Gaffney [4] resp.
14  Gibson and Hobbs [1] with respect to A equivalence. If the input is one of the simple singularities in
15  [1] it returns a normal form otherwise a zero ideal(i.e not simple).
16
17REFERENCES:
18  [1] Gibson,C.G; Hobbs,C.A.:Simple SIngularities of Space Curves.
19  Math.Proc. Comb.Phil.Soc.(1993),113,297.
20  [2] Hefez,A;Hernandes,M.E.:Standard bases for local rings of branches and their modules of differentials.
21  Journal of Symbolic Computation 42(2007) 178-191.
22  [3] Hefez,A;Hernandes,M.E.:The Analytic Classification Of Plane Branches. Bull.Lond Math Soc.43.(2011) 2,289-298.
23  [4] Bruce, J.W.,Gaffney, T.J.: Simple singularities of mappings (C, 0) ->(C2,0).
24  J. London Math. Soc. (2) 26 (1982), 465-474.
25  [5] Ishikawa,G; Janeczko,S.: The Complex Symplectic Moduli Spaces of Unimodal Parametric Plane Curve  NEU Singularities. Insitute of Mathematics of the Polish Academy of Sciences,Preprint 664(2006)
26PROCEDURES:
27          sagbiAlg(G);    Compute the Sagbi-basis of the Algebra.
28          sagbiMod(I,A);  Compute the Sagbi- basis of the Module.
29          semiGroup(G);   Compute the Semi-Group of the Algebra provided the input is Sagbi Bases of the Algebra.
30          semiMod(I,A);   Compute the Semi-Module provided that the input are the Sagbi Bases of the Algebra resp.Module.
31          planeCur(I);    Compute the type of the Simple Plane Curve singularity.
32          spaceCur(I);    Compute the type of the simple Space Curve singularity.
33";
34LIB "algebra.lib";
35LIB "curvepar.lib";
36///////////////////////////////////////////////////////////////////////////////
37proc planeCur(ideal I)
38"USAGE":  planeCur(I);  I ideal
39RETURN: An ideal.Ideal is one of the singularity in the list of Bruce and Gaffney [4]
40EXAMPLE: example planeCur;  shows an example
41"
42{
43   def R=basering;
44   I=sortMOD(I);
45   list M;
46   list K;
47   if(I==0)
48   {return(I);}
49   ideal G=sagbiAlg(I);
50   list L=semiGroup(G);
51   ideal J=diff(G,var(1));
52   J=sagbiMod(J,G);
53   M=semiMod(J,G);
54   int C=L[2];
55   ideal Z=0,0;
56   if(L[1][1]>4)
57   {
58      return(Z);
59   }
60   if(L[1][1]==1)
61   {
62        ideal A=var(1);
63        K=Guess(A);
64        if(CompareList(M,K,6)!=0)
65        {
66              return(A);
67        }
68        else
69        {
70            return(Z);
71        }
72    }
73    if(L[1][1]==2)
74    {
75         ideal A=var(1)^2,var(1)^(L[2]+1);
76         K=Guess(A);
77         if(CompareList(M,K,6)!=0)
78         {
79            return(A);
80         }
81         else
82         {
83            return(Z);
84         }
85    }
86    if(L[1][1]==4)
87    {
88          if(L[1][2]==5)
89          {
90                intvec q=4,5;
91                if((L[1]==q)&&(L[2]==12)&&(size(L[3])==7))
92                {
93                   intvec q1=3,4; intvec q2=3,4,10;
94                   if((M[4]==q1)&&(M[5]==11)&&(size(M[6])==6))
95                   {
96                       ideal A=var(1)^4,var(1)^5;
97                        K=Guess(A);
98                       if(CompareList(M,K,6)!=0)
99                       {
100                           return(A);
101                       }
102                   }
103                   if((M[4]==q2)&&(M[5]==7)&&(size(M[6])==3))
104                   {
105                       ideal A=var(1)^4,var(1)^5+var(1)^7;
106                       K=Guess(A);
107                       if(CompareList(M,K,6)!=0)
108                       {
109                         return(A);
110                       }
111                   }
112                   else
113                   {
114                      return(Z);
115                   }
116                }
117                else
118                {
119                   return(Z);
120                }
121          }
122          if(L[1][2]==6)
123          {
124                ideal A=var(1)^4,var(1)^6+var(1)^(L[1][3]-6);
125                K=Guess(A);
126                if(L[1][3] mod 2 !=0)
127                {
128                   ideal S=t4,t6+t^(M[2]-9);
129                   if(CompareList(M,K,6)!=0)
130                   {
131                      return(S);
132                   }
133                   if(CompareList(M,K,6)==0)
134                   {
135                      int m=size(K[4])+1;
136                      if(size(M[4])==m)
137                      {
138                        return(S);
139                      }
140                      else{return(Z);}
141                   }
142                }
143                else
144                {
145                   return(Z);
146                }
147          }
148          if(L[1][2]==7)
149          {
150                intvec q=4,7;list K;
151                ideal A=var(1)^4,var(1)^7;
152                ideal B=var(1)^4,var(1)^7+var(1)^9;
153                ideal T=var(1)^4,var(1)^7+var(1)^10;
154                list Q=A,B,T;
155                for(int i=1;i<=3;i++)
156                {    K=Guess(Q[i]);
157                     if(CompareList(M,K,6)!=0)
158                     {
159                        if(i==1)
160                        {
161                            return(A);
162                            break;
163                        }
164                        if(i==2)
165                        {
166                          return(B);
167                          break;
168                        }
169                        if(i==3)
170                        {
171                           return(T);
172                           break;
173                        }
174                     }
175                }
176                else
177                {
178                     return(Z);
179                }
180        }
181        else
182        {
183             return(Z);
184        }
185   }
186   if(L[1][1]==3)
187   {
188          int k=L[1][2]-1;
189          int p=L[1][2]-2;
190          if(k mod 3 ==0)
191          {
192              if(size(M[4])==2)
193              {
194                  ideal A=var(1)^3,var(1)^L[1][2];
195                  ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5];
196                  list Q=A,B;
197                  for(int i=1;i<=2;i++)
198                  {  K=Guess(Q[i]);
199                     if(CompareList(M,K,6)!=0)
200                     {
201                        return(Q[i]);
202                     }
203                  }
204              }
205              if(size(M[4])==3)
206              {
207                  ideal A=var(1)^3,var(1)^L[1][2];
208                  ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5];
209                  list Q=A,B;
210                  for(int i=1;i<=2;i++)
211                  {  K=Guess(Q[i]);
212                     if(CompareList(M,K,6)!=0)
213                     {
214                          return(Q[i]);
215                     }
216                  }
217              }
218              else
219              {
220                  return(Z);
221              }
222          }
223          if(p mod 3 ==0)
224          {
225               if(size(M[4])==2)
226               {
227                    ideal A=var(1)^3,var(1)^L[1][2];
228                    ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5];
229                    list Q=A,B;
230                    for(int i=1;i<=2;i++)
231                    {    K=Guess(Q[i]);
232                         if(CompareList(M,K,6)!=0)
233                         {
234                            return(Q[i]);
235                         }
236                     }
237                }
238                if(size(M[4])==3)
239                {
240                     ideal A=var(1)^3,var(1)^L[1][2];
241                     ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5];
242                     list Q=A,B;
243                     for(int i=1;i<=2;i++)
244                     {    K=Guess(Q[i]);
245                          if(CompareList(M,K,6)!=0)
246                          {
247                              return(Q[i]);
248                          }
249                     }
250                }
251                else
252                {
253                      return(Z);
254                }
255           }
256           else
257           {
258                 return(Z)
259           }
260   }
261}
262example
263{
264"EXAMPLE:"; echo=2;
265  ring R=0,t,Ds;
266  ideal I=t4+4t5+6t6+8t7+13t8+12t9+10t10+12t11+6t12+4t13+4t14+t16,t7+7t8+22t9+51t10+113t11+219t12+366t13+589t14+876t15+1170t16+1514t17
267+1828t18+2011t19+2165t20+2163t21+1982t22+1806t23+1491t24+1141t25+889t26+588t27+379t28+252t29+120t30+72t31+36t32+9t33+9t34+t36;
268  planeCur(I);
269}
270
271////////////////////////////////////////////////////////////////////////////////
272proc spaceCur(ideal I)
273"USAGE":  spaceCur(I);  I ideal
274RETURN: an ideal. Ideal is one of the singularity  in the list of C.G.Gibson and C.A.Hobbs.
275EXAMPLE: example spaceCur;  shows an example
276"
277{
278   def R=basering;
279   I=sortMOD(I);
280   list M;
281   list K;
282   if(I==0)
283   {return(I);}
284   ideal G=sagbiAlg(I);
285   if(size(G)<=2){return(planeCur(G));}
286   list L=semiGroup(G);
287   ideal J=diff(G,var(1));
288   J=sagbiMod(J,G);
289   M=semiMod(J,G);
290   int C=L[2];
291   ideal Z=0,0,0;
292   if(L[1][1]>5)
293   {
294      return(Z);
295   }
296   if(L[1][1]==3)
297   {
298          int k=L[1][2]-1;
299          int p=L[1][2]-2;
300          if(k mod 3 ==0)
301          {
302               poly q=var(1)*(J[2])-G[2];
303               if(leadexp(q)!=leadexp(J[3]))
304               {
305                  if(size(M[4])!=3)
306                  {
307                     ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
308                     return(B);
309                  }
310                  if(size(M[4])==3)
311                  {
312                     ideal I1=G[1],G[2];
313                     I1=sortMOD(I1);
314                     ideal T=sagbiAlg(I1);
315                     ideal J1=diff(T,var(1));
316                     J1=sagbiMod(J1,T);
317                     K=semiMod(J1,T);
318                     if(size(K[4])!=2)
319                     {
320                           ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
321                          return(B);
322                     }
323                     if(size(K[4])==2)
324                     {
325                         ideal A=var(1)^3,var(1)^L[1][2],var(1)^L[1][3];
326                          return(A);
327                     }
328                  }
329               }
330               if(leadexp(q)==leadexp(J[3]))
331               {
332                  if(size(M[4])!=3)
333                  {
334                      ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
335                      return(B);
336                  }
337                  if(size(M[4])==3)
338                  {
339                      ideal I1=G[1],G[2];
340                      I1=sortMOD(I1);
341                      ideal T=sagbiAlg(I1);
342                      ideal J1=diff(T,var(1));
343                      J1=sagbiMod(J1,T);
344                      K=semiMod(J1,T);
345                      if(size(K[4])!=2)
346                      {
347                            ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
348                            return(B);
349                      }
350                      if(size(K[4])==2)
351                      {
352                           ideal A=var(1)^3,var(1)^L[1][2],var(1)^L[1][3];
353                           return(A);
354                      }
355                  }
356               }
357          }
358          if(p mod 3 ==0)
359          {
360                poly q=var(1)^3*(J[2])-var(1)^2*(G[2]);
361               if(leadexp(q)!=leadexp(J[3]))
362               {
363                  if(size(M[4])!=3)
364                  {
365                     ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
366                     return(B);
367                  }
368                  if(size(M[4])==3)
369                  {
370                       ideal I1=G[1],G[2];
371                       I1=sortMOD(I1);
372                       ideal T=sagbiAlg(I1);
373                       ideal J1=diff(T,var(1));
374                       J1=sagbiMod(J1,T);
375                       K=semiMod(J1,T);
376                       if(size(K[4])!=2)
377                       {
378                           ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
379                           return(B);
380                       }
381                       if(size(K[4])==2)
382                       {
383                           ideal A=var(1)^3,var(1)^L[1][2],var(1)^L[1][3];
384                           return(A);
385                       }
386                 }
387              }
388              if(leadexp(q)==leadexp(J[3]))
389              {
390                  if(size(M[4])!=3)
391                  {
392                        ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
393                        return(B);
394                  }
395                  if(size(M[4])==3)
396                  {
397                       ideal I1=G[1],G[2];
398                       ideal T=sagbiAlg(I1);
399                       ideal J1=diff(T,var(1));
400                       J1=sagbiMod(J1,T);
401                       K=semiMod(J1,T);
402                       if(size(K[4])!=2)
403                       {
404                           ideal B=var(1)^3,var(1)^L[1][2]+var(1)^M[5],var(1)^L[1][3];
405                           return(B);
406                       }
407                       if(size(K[4])==2)
408                       {
409                           ideal A=var(1)^3,var(1)^L[1][2],var(1)^L[1][3];
410                           return(A);
411                       }
412                  }
413             }
414           }
415           else
416           {
417                 return(Z);
418           }
419   }
420   if(L[1][1]==4)
421   {
422      if(L[1][2]==5)
423      {
424         if(L[1][3]==11)
425         {
426             ideal A=var(1)^4,var(1)^5,var(1)^11;
427             ideal B=var(1)^4,var(1)^5+var(1)^7,var(1)^11;
428             list Q=A,B;
429             ideal Ij=jet(I,10);
430             Ij=simplify(Ij,2);
431             ideal Gj=sagbiAlg(Ij);
432             list Lj=semiGroup(Gj);
433             ideal Jj=diff(Gj,var(1));
434             Jj=sagbiMod(Jj,Gj);
435             list Mj=semiMod(Jj,Gj);
436             if(size(Mj[4])==2)
437             {
438                 K=Guess(Q[1]);
439                 if(CompareList(M,K,6)!=0)
440                 {
441                    return(Q[1]);
442                 }
443             }
444             if(size(Mj[4])==3)
445             {
446                 K=Guess(Q[2]);
447                 if(CompareList(M,K,6)!=0)
448                 {
449                    return(Q[2]);
450                 }
451             }
452         }
453         if(L[1][3]!=11)
454         {
455              ideal A=var(1)^4,var(1)^5,var(1)^6;
456              ideal B=var(1)^4,var(1)^5,var(1)^7;
457              list Q=A,B;
458              for(int i=1;i<=2;i++)
459              {
460                 K=Guess(Q[i]);
461                 if(CompareList(M,K,6)!=0)
462                 {
463                    return(Q[i]);
464                    break;
465                 }
466              }
467         }
468         else
469         {return(Z);
470         }
471      }
472      if(L[1][2]==6)
473      {
474          if(size(L[1])==3)
475          {
476               if(size(M[4])==3)
477               {
478                   ideal A=var(1)^4,var(1)^6,var(1)^L[1][3];
479                   K=Guess(A);
480                   if(CompareList(M,K,6)!=0)
481                   {
482                       return(A);
483                   }
484                   else
485                   {
486                       return(Z);
487                   }
488               }
489               if(size(M[4])==4)
490               {
491                   ideal A=var(1)^4,var(1)^6+var(1)^(L[1][3]-2),var(1)^L[1][3];
492                   K=Guess(A);
493                   if(CompareList(M,K,6)!=0)
494                   {
495                          return(A);
496                   }
497                   else
498                   {
499                       return(Z);
500                   }
501               }
502          }
503          if(size(L[1])==4)
504          {
505               if(size(M[4])==4)
506               {
507                   ideal A=var(1)^4,var(1)^6+var(1)^(L[1][3]-4),var(1)^L[1][3];
508                   K=Guess(A);
509                   if(CompareList(M,K,6)!=0)
510                   {
511                       return(A);
512                   }
513                   else
514                   {
515                       return(Z);
516                   }
517               }
518               if(size(M[4])==5)
519               {
520                   ideal A=var(1)^4,var(1)^6+var(1)^(L[1][4]-8),var(1)^L[1][4];
521                   K=Guess(A);
522                   if(CompareList(M,K,6)!=0)
523                   {
524                          return(A);
525                   }
526                   else
527                   {
528                       return(Z);
529                   }
530               }
531          }
532          else
533          {
534              return(Z);
535          }
536      }
537      if(L[1][2]==7)
538      {
539          if(L[1][3]==9)
540          {
541               ideal A=var(1)^4,var(1)^7,var(1)^9+var(1)^10;
542               ideal B=var(1)^4,var(1)^7,var(1)^9;
543               list Q=A,B;
544               for(int i=1;i<=2;i++)
545               {
546                  K=Guess(Q[i]);
547                  if(CompareList(M,K,6)!=0)
548                  {
549                     return(Q[i]);
550                     break;
551                  }
552               }
553          }
554          if(L[1][3]==10)
555          {
556               ideal A=var(1)^4,var(1)^7,var(1)^10;
557               ideal B=var(1)^4,var(1)^7+var(1)^9,var(1)^10;
558               list Q=A,B;
559               for(int i=1;i<=2;i++)
560               {
561                    K=Guess(Q[i]);
562                   if(CompareList(M,K,6)!=0)
563                   {
564                      return(Q[i]);
565                      break;
566                   }
567               }
568          }
569          if(L[1][3]==13)
570          {
571               ideal A=var(1)^4,var(1)^7,var(1)^13;
572               ideal B=var(1)^4,var(1)^7+var(1)^9,var(1)^13;
573               list Q=A,B;
574               ideal Ij=jet(I,12);
575               Ij=simplify(Ij,2);
576               ideal Gj=sagbiAlg(Ij);
577               list Lj=semiGroup(Gj);
578               ideal Jj=diff(Gj,var(1));
579               Jj=sagbiMod(Jj,Gj);
580               Jj=jet(Jj,12);
581               Jj=simplify(Jj,2);
582               list Mj=semiMod(Jj,Gj);
583               if(size(Jj)==2)
584               {
585                    K=Guess(Q[1]);
586                    if(CompareList(M,K,6)!=0)
587                    {
588                       return(A);
589                       break;
590                    }
591               }
592               if(size(Jj)==3)
593               {
594                    K=Guess(Q[2]);
595                    if(CompareList(M,K,6)!=0)
596                    {
597                       return(B);
598                       break;
599                    }
600                }
601          }
602          if(L[1][3]==17)
603          {
604                ideal A=var(1)^4,var(1)^7,var(1)^17;
605                ideal B=var(1)^4,var(1)^7+var(1)^9,var(1)^17;
606                ideal T=var(1)^4,var(1)^7+var(1)^10,var(1)^17;
607                list Q=A,B,T;
608                for(int i=1;i<=3;i++)
609                {
610                    K=Guess(Q[i]);
611                    if(CompareList(M,K,6)!=0)
612                    {
613                        if(i==2)
614                        {
615                           return(Q[i]);
616                           break;
617                        }
618                        else
619                        {
620                             ideal Ij=jet(I,16);
621                             Ij=simplify(Ij,2);
622                             ideal Gj=sagbiAlg(Ij);
623                             list Lj=semiGroup(Gj);
624                             ideal Jj=diff(Gj,var(1));
625                             Jj=sagbiMod(Jj,Gj);
626                             Jj=jet(Jj,16);
627                             Jj=simplify(Jj,2);
628                             list Mj=semiMod(Jj,Gj);
629                             if(size(Jj)==2)
630                             {
631                                 if(CompareList(M,K,6)!=0)
632                                 {
633                                    return(A);
634                                    break;
635                                 }
636                             }
637                             if(size(Jj)==3)
638                             {
639                                 if(CompareList(M,K,6)!=0)
640                                 {
641                                    return(T);
642                                    break;
643                                 }
644                             }
645                        }
646                    }
647                }
648          }
649          else
650          {
651              return(Z);
652          }
653      }
654   }
655}
656example
657{
658  "EXAMPLE:"; echo=2;
659   ring R=0,t,Ds;
660ideal I=t3+3t4+3t5+t6,t13+14t14+92t15+377t16+1079t17+2288t18+3718t19+4719t20+4719t21+3718t22+2288t23+1079t24+377t25+92t26+14t27+t28,t17+17t18+136t19+680t20+2380t21+6188t22+12376t23+19448t24+24310t25+24310t26+19448t27+12376t28+6188t29+2380t30+680t31+136t32+17t33+t34;
661  spaceCur(I);
662}
663
664////////////////////////////////////////////////////////////////////////////////
665proc sagbiAlg(ideal I,list #)
666"USAGE":  sagbiAlg(I);  I ideal
667RETURN: An ideal.The sagbi bases of I.
668EXAMPLE: example sagbiAlg;  shows an example
669{
670  def R=basering;
671  def O=changeord(list(list("Ds",nvars(R))));
672  setring O;
673  ideal I=imap(R,I);
674  ideal L,S;
675  poly h;
676  int z,n;
677
678  if(size(I)==0){return(I);}
679  if(size(#)==0)
680  {
681    int b=ConductorBound(I);
682    S=interReduceSagbi(I,b) ;
683    b=ConductorBound(S);
684    // int b=200;
685    // b=correctBound(I,b);
686  }
687  else
688  {
689    int b=#[1];
690  }
691  S=interReduceSagbi(I,b) ;
692  // b=correctBound(S,b);
693  while(size(S)!=n)
694  {
695    n=size(S);
696    L=sagbiSP(S);
697    for (z=1;z<=size(L);z++)
698    {
699      h=sagbiNF(L[z],S,b);
700      if(h!=0)
701      {
702        S=insertOne(h,S,b);
703      }
704    }
705  }
706  setring R;
707  ideal S=imap(O,S);
708  return(S);
709}
710example
711{
712  "EXAMPLE:"; echo=2;
713  ring R=0,t,ds;
714  ideal I=t8,t10+t13,t12+t15;
715  sagbiAlg(I);
716  I=t8,t10+t13,t12+2t15;
717  sagbiAlg(I);
718}
719
720//////////////////////////////////////////////////////////////////////////////// NEU
721proc reducedSagbiAlg(ideal I,list #)
722{
723
724   I=sagbiAlg(I,#);
725   intvec L=semiGroup(I)[3];
726   if(size(#)==0)
727   {
728      int b=findConductor(L);
729   }
730   else
731   {
732      int b=#[1];
733   }
734   int i;
735   poly q;
736   for(i=1;i<=size(I);i++)
737   {
738      q=I[i]-lead(I[i]);
739      q=sagbiNF(q,I,b);
740      I[i]=lead(I[i])+q;
741   }
742   return(I);
743}
744example
745{
746  "EXAMPLE:"; echo=2;
747ring R=0,t,ds;
748ideal I=t4+2t9,t9+t10+19/18t11-3t12+t13-t14;
749reducedSagbiAlg(I);
750}
751////////////////////////////////////////////////////////////////////////////////   NEU
752proc classifyAEQunimodal(ideal I)
753"USAGE":  classifyAEQunimodal(I);  I ideal generated by 2 polynomials
754RETURN: An ideal.Ideal is one of the singularity in the list of Ishikawa and Jenczko [5]
755EXAMPLE: example classifyAEQunimodal;  shows an example
756"
757{
758    def R=basering;
759    ring @S=0,t,ds;
760    ideal I=fetch(R,I);
761    ideal J;
762    poly g;
763    if(size(I)>=3){ERROR("not a plane curve");}
764    I=simplify(I,2);   //deletes zero’s i I
765    I=simplify(I,1);   //creates monic generators in I
766    if(ord(I[1])>ord(I[2])){poly q=I[2];I[2]=I[1];I[1]=q;}
767    if((ord(I[1])>=6)||(ord(I[1])<=3)){return("not in the unimodal list");}
768    //compute estimate of the term with the modulus
769    int c=ord(sagbiNF(I[2],ideal(I[1]),10));
770    if(c==10)
771    {
772       if(ord(I[1])!=4){return("not in the unimodal list");}
773       c=ord(I[2][2])+2;
774    }
775    else
776    {
777       c=0;
778       intvec v=ord(I[1]),ord(I[2]);
779       if(v==intvec(5,6)){c=14;}
780       if(v==intvec(5,7)){c=18;}
781       if(v==intvec(5,8)){c=22;}
782       if(v==intvec(4,9)){c=19;}
783       if(v==intvec(4,11)){c=25;}
784       if(c==0){return("not in the unimodal list");}
785    }
786    while(size(I[1])>1)
787    {
788      I=jet(subst(I,t,t-number(1)/number(ord(I[1]))*leadcoef(I[1][2])*t^(ord(I[1][2])-ord(I[1])+1)),c);
789    }
790    ideal G=I;
791    G[2]=sagbiNF(G[2],ideal(G[1]),c);
792    ideal M=sagbiMod(diff(G,t),G);
793    list K=semiMod(M,G);
794
795    if(K[1]==intvec(4,9))
796    {
797       if(K[4]==intvec(3,8)){J=t4,t9;}
798       if(K[4]==intvec(3,8,22)){J=t4,t9+t19;}
799       if(K[4]==intvec(3,8,18)){J=t4,t9+t15;}
800       if(K[4]==intvec(3,8,14)){J=t4,t9+t11;}
801       if(K[4]==intvec(3,8,13))
802       {
803          G=reducedSagbiAlg(G,15);
804          if(ord(G[2][4])==14)
805          {
806             //kill the term t14 by some transformation
807             G=subst(G,t,t-leadcoef(G[2][4])/9*t^6);
808             G=jet(G,15);
809             G[1]=sagbiNF(G[1],ideal(G[2]),15);
810             //arrange the first element to be t4
811             while(size(G[1])>1)
812             {
813                G=subst(G,t,t-1/4*(G[1]-lead(G[1]))/t^3);
814                G=jet(G,15);
815             }
816          }
817          G[2]=sagbiNF(G[2],ideal(G[1]),15);
818          //arrange the coefficient of t10 to become 1
819          number m=leadcoef(G[2]-lead(G[2]));
820          G[2]=m^9*subst(G[2],t,1/m*t);
821          J=G;
822       }
823       if(K[4]==intvec(3,8,13,18))
824       {
825          G=reducedSagbiAlg(G,11);
826          number m=leadcoef(G[2]-lead(G[2]));
827          G[2]=m^9*subst(G[2],t,1/m*t);
828          J=G;
829       }
830    }
831    if(K[1]==intvec(4,11))
832    {
833       if(K[4]==intvec(3,10)){J=t4,t11;}
834       if(K[4]==intvec(3,10,28)){J=t4,t11+t25;}
835       if(K[4]==intvec(3,10,24)){J=t4,t11+t21;}
836       if(K[4]==intvec(3,10,20)){J=t4,t11+t17;}
837       if(K[4]==intvec(3,10,16))
838       {
839          G=reducedSagbiAlg(G,14);
840          number m=leadcoef(G[2]-lead(G[2]));
841          number l=leadcoef(G[2][3]);
842          //lambda^2=l^2/m^3
843          J=G;
844       }
845       if(K[4]==intvec(3,10,17))
846       {
847          G=reducedSagbiAlg(G,21);
848          if(ord(G[2][4])==18)
849          {
850             //kill the term t18 by some transformation
851             G=subst(G,t,t-leadcoef(G[2][4])/11*t^8);
852             G=jet(G,21);
853             G[1]=sagbiNF(G[1],ideal(G[2]),21);
854             //arrange the first element to be t4
855             while(size(G[1])>1)
856             {
857                G=subst(G,t,t-1/4*(G[1]-lead(G[1]))/t^3);
858                G=jet(G,21);
859             }
860          }
861          G[2]=sagbiNF(G[2],ideal(G[1]),21);
862          //arrange the coefficient of t14 to become 1
863          number m=leadcoef(G[2]-lead(G[2]));
864          number l=leadcoef(G[2][4]);
865          //lambda^2=l^3/m^10
866          J=G;
867
868
869       }
870       if(K[4]==intvec(3,10,17,24))
871       {
872          G=reducedSagbiAlg(G,18);
873          //arrange the coefficient of t14 to become 1
874          number m=leadcoef(G[2]-lead(G[2]));
875          G[2]=t11+t14+leadcoef(G[2][3])/m^2*t17;
876          J=G;
877       }
878    }
879    if((size(K[1])==3)&&(K[1][1]==4)&&(K[1][2]==10))
880    {
881       int l=(K[1][3]-19) div 2;
882       G=reducedSagbiAlg(G,2*l+12);
883       number m=leadcoef(G[2]-lead(G[2]));
884       number s=leadcoef(G[2][3]);
885       //lambda^(2l-1)=s^(2l-1)/m^(2l+1)
886       J=G;
887    }
888    if(K[1]==intvec(5,6))
889    {
890       if(K[4]==intvec(4,5)){J=t5,t6;}
891       if(K[4]==intvec(4,5,18)){J=t5,t6+t14;}
892       if(K[4]==intvec(4,5,13)){J=t5,t6+t9;}
893       if(K[4]==intvec(4,5,12))
894       {
895          G=reducedSagbiAlg(G,9);
896          if(ord(G[2][2])==7)
897          {
898             //kill the term t7 by some transformation
899             G=subst(G,t,t-leadcoef(G[2][2])/6*t^2);
900             G=jet(G,10);
901             G[1]=sagbiNF(G[1],ideal(G[2]),9);
902             //arrange the first element to be t4
903             while(size(G[1])>1)
904             {
905                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
906                G=jet(G,9);
907             }
908          }
909          G[2]=sagbiNF(G[2],ideal(G[1]),9);
910          //arrange the coefficient of t8 to become 1
911          number m=leadcoef(G[2]-lead(G[2]));
912          number l=leadcoef(G[2][3]);
913          //lambda^2=l^2/m^3
914          J=G;
915
916       }
917    }
918    if(K[1]==intvec(5,7))
919    {
920       if(K[4]==intvec(4,6)){J=t5,t7;}
921       if(K[4]==intvec(4,6,22)){J=t5,t7+t18;}
922       if(K[4]==intvec(4,6,17)){J=t5,t7+t13;}
923       if(K[4]==intvec(4,6,12))
924       {
925          G=reducedSagbiAlg(G,11);
926          if(ord(G[2][3])==9)
927          {
928             //kill the term t9 by some transformation
929             G=subst(G,t,t-leadcoef(G[2][3])/7*t^3);
930             G=jet(G,11);
931             G[1]=sagbiNF(G[1],ideal(G[2]),11);
932             //arrange the first element to be t4
933             while(size(G[1])>1)
934             {
935                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
936                G=jet(G,11);
937             }
938          }
939          G[2]=sagbiNF(G[2],ideal(G[1]),11);
940          //arrange the coefficient of t8 to become 1
941          number m=leadcoef(G[2]-lead(G[2]));
942          G[2]=m^7*subst(G[2],t,1/m*t);
943          J=G;
944       }
945       if(K[4]==intvec(4,6,15))
946       {
947          G=reducedSagbiAlg(G,14);
948          if(ord(G[2][2])==9)
949          {
950             //kill the term t9 by some transformation
951             G=subst(G,t,t-leadcoef(G[2][2])/7*t^3);
952             G=jet(G,14);
953             G[1]=sagbiNF(G[1],ideal(G[2]),14);
954             //arrange the first element to be t4
955             while(size(G[1])>1)
956             {
957                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
958                G=jet(G,14);
959             }
960          }
961          G[2]=sagbiNF(G[2],ideal(G[1]),14);
962          //arrange the coefficient of t11 to become 1
963          number m=leadcoef(G[2]-lead(G[2]));
964          number l=leadcoef(G[2][3]);
965          //lambda^2=l^2/m^3
966          J=G;
967       }
968
969    }
970    if(K[1]==intvec(5,8))
971    {
972       if(K[4]==intvec(4,7)){J=t5,t8;}
973       if(K[4]==intvec(4,7,26)){J=t5,t8+t22;}
974       if(K[4]==intvec(4,7,21)){J=t5,t8+t17;}
975       if(K[4]==intvec(4,7,13))
976       {
977          G=reducedSagbiAlg(G,12);
978          if(ord(G[2][3])==11)
979          {
980             //kill the term t11 by some transformation
981             G=subst(G,t,t-leadcoef(G[2][3])/8*t^4);
982             G=jet(G,12);
983             G[1]=sagbiNF(G[1],ideal(G[2]),12);
984             //arrange the first element to be t4
985             while(size(G[1])>1)
986             {
987                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
988                G=jet(G,12);
989             }
990          }
991          G[2]=sagbiNF(G[2],ideal(G[1]),12);
992          //arrange the coefficient of t9 to become 1
993          number m=leadcoef(G[2]-lead(G[2]));
994          G[2]=m^8*subst(G[2],t,1/m*t);
995          J=G;
996       }
997
998       if(K[4]==intvec(4,7,16))
999       {
1000          G=reducedSagbiAlg(G,14);
1001          if(ord(G[2][2])==11)
1002          {
1003             //kill the term t11 by some transformation
1004             G=subst(G,t,t-leadcoef(G[2][2])/8*t^4);
1005             G=jet(G,14);
1006             G[1]=sagbiNF(G[1],ideal(G[2]),14);
1007             //arrange the first element to be t4
1008             while(size(G[1])>1)
1009             {
1010                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
1011                G=jet(G,14);
1012             }
1013          }
1014          G[2]=sagbiNF(G[2],ideal(G[1]),14);
1015          //arrange the coefficient of t12 to become 1
1016          number m=leadcoef(G[2]-lead(G[2]));
1017          number l=leadcoef(G[2][3]);
1018          //lambda^2=l^2/m^3
1019          J=G;
1020
1021       }
1022       if(K[4]==intvec(4,7,18))
1023       {
1024          G=reducedSagbiAlg(G,17);
1025          if(ord(G[2][2])==11)
1026          {
1027             //kill the term t11 by some transformation
1028             G=subst(G,t,t-leadcoef(G[2][2])/8*t^4);
1029             G=jet(G,17);
1030             G[1]=sagbiNF(G[1],ideal(G[2]),17);
1031             //arrange the first element to be t4
1032             while(size(G[1])>1)
1033             {
1034                G=subst(G,t,t-1/5*(G[1]-lead(G[1]))/t^4);
1035                G=jet(G,17);
1036             }
1037          }
1038          G[2]=sagbiNF(G[2],ideal(G[1]),17);
1039          //arrange the coefficient of t12 to become 1
1040          number m=leadcoef(G[2]-lead(G[2]));
1041          number l=leadcoef(G[2][3]);
1042          //lambda^2=l^2/m^3
1043          J=G;
1044       }
1045    }
1046    setring R;
1047    ideal J=fetch(@S,J);
1048    if(size(J)==0)
1049    {
1050       return("not in the unimodal list");
1051    }
1052    return(J);
1053}
1054example
1055{
1056  "EXAMPLE:"; echo=2;
1057ring R=0,t,ds;
1058ideal I=t4,9t9+18t10+38t11-216t12+144t13-288t14;
1059classifyAEQunimodal(I);
1060I=t4,9t9+18t10+40t11-216t12+144t13-288t14;
1061classifyAEQunimodal(I);
1062I=t4,t11+t12+3t14+2t15+7t16+7t17;
1063classifyAEQunimodal(I);
1064I=t4,t11+t14+25/22t17+3t18+4t21;
1065classifyAEQunimodal(I);
1066I=t5,t6+2t7+t8+3t9;
1067classifyAEQunimodal(I);
1068I=t5,t7+3t8+3t9+5t10;
1069classifyAEQunimodal(I);
1070I=t5,t7+3t11+3t12+5t13;
1071classifyAEQunimodal(I);
1072I=t5,t8+3t9+5t10+2t11+3t12+5t13;
1073classifyAEQunimodal(I);
1074I=t5,t8+5t11+3t12+7t13+5t14;
1075classifyAEQunimodal(I);
1076I=t5,t8+5t11+7t13+5t14+7t15+2t16+8t17;
1077classifyAEQunimodal(I);
1078I=subst(I,t,t+t2);
1079classifyAEQunimodal(I);
1080I=t4+2t5+3t6+5t7+t8,9t9+18t10+40t11-216t12+144t13-288t14;
1081classifyAEQunimodal(I);
1082}
1083////////////////////////////////////////////////////////////////////////////////   NEU
1084proc computeModulus(poly p)
1085"USAGE":  computeModulus(p);  p monic poly with 3 or 4 monomials
1086RETURN: A polynomial with first and second coefficient 1
1087EXAMPLE: computeModulus;  shows an example
1088ASSUME: the basering has one vearable and one parameter
1089"
1090{
1091   def R=basering;
1092   int a1=ord(p);
1093   int a2=ord(p-lead(p));
1094   number m=leadcoef(p-lead(p));
1095
1096   poly q=par(1)^(a2-a1)-1/m;
1097   ring S=(0,a),t,ds;
1098   number m=fetch(R,m);
1099   minpoly=par(1)^(a2-a1)-1/m;
1100   poly p=fetch(R,p);
1101   p=1/par(1)^a1*subst(p,var(1),par(1)*var(1));
1102   setring R;
1103   p=imap(S,p);
1104   return(list(p,q));
1105}
1106example
1107{
1108  "EXAMPLE:"; echo=2;
1109ring R=(0,a),t,ds;
1110poly p=t8-395/16t14+4931/32t17;
1111computeModulus(p);
1112p=t8+3t12-395/16t14;
1113computeModulus(p);
1114p=t8-395/16t14+4931/32t17;
1115computeModulus(p);
1116
1117}
1118
1119////////////////////////////////////////////////////////////////////////////////   NEU
1120static proc n_thRoot(poly p,int n, int b)
1121{
1122   //computes the n-th root of 1+p up to order b
1123   //assumes that p(0)=0
1124   poly s=1;
1125   poly q=jet(p,b);
1126   if(q==0){return(s);}
1127   int i;
1128   for(i=1;i<=b;i++)
1129   {
1130      s=s+bino(n,i)*q;
1131      q=jet(q*p,b);
1132   }
1133   return(jet(s,b));
1134}
1135////////////////////////////////////////////////////////////////////////////////   NEU
1136static proc bino(number n, int i)
1137{
1138//computes the i-th binomial coefficient of 1/n
1139   if(i==0){return(1);}
1140   return(bino(n,i-1)*(1/n-i+1)/i);
1141}
1142////////////////////////////////////////////////////////////////////////////////
1143proc sagbiMod(ideal I,ideal G)
1144"USAGE":  sagbiMod(I,G);  I an ideal module and ideal G being the sagbi bases of the Algebra
1145RETURN: An ideal. the sagbi bases for the differential module.
1146EXAMPLE: example sagbiMod;  shows an example
1147{
1148  def R=basering;//up till now the ordering of the base ring is ds
1149  def O=changeord(list(list("Ds",nvars(R))));
1150  setring O;
1151  ideal I=imap(R,I);
1152  ideal G=imap(R,G);
1153  int n=ncols(G);poly h;
1154  if(I==0)
1155  { return(I);}
1156  ideal S,J,M;
1157  I=sortMOD(I);
1158  if(deg(lead(I[1]))<=1)
1159  { setring R;
1160    return(imap(O,I));}
1161  int b=ConductorBound(lead(G))+deg(lead(I[1]));
1162  list P;int i;
1163  P=createP(I);
1164  while(size(P)!=0)
1165  {
1166      J=P[1][1],P[1][2];
1167      P=delete(P,1);
1168      S=SpolyMOD(J,G);
1169      for(i=1;i<=size(S);i++)
1170      {
1171         h=sagbiNFMOD(S[i],G,I,b);
1172         if(h!=0)
1173         {
1174             h=simplify(h,1);
1175             P=enlargeP(h,P,I);
1176             I[size(I)+1]=h;
1177         }
1178      }
1179  }
1180  I=sortMOD(I);
1181  setring R;
1182  ideal K=imap(O,I);
1183  return(K);
1184}
1185example
1186{
1187  "EXAMPLE:"; echo=2;
1188  ring r=0,t,Ds;
1189  ideal G=t8,t10+t13,t12+t15,t23-t29,t27;
1190  ideal I=diff(G,t);
1191  sagbiMod(I,G);
1192}
1193
1194////////////////////////////////////////////////////////////////////////////////
1195proc semiGroup(ideal I)
1196"USAGE": semiGroup(I);  I ideal the sagbi bases of Algebra.
1197RETURN: list L; list with three entries associated to the algebra generated by
1198   the sagbi basis:
1199   generators of the semigroup
1200   the conductor
1201    the semigroup
1202EXAMPLE: example planeCur;  shows an example
1203{
1204   list M;
1205   if(deg(I[1])<=1)
1206   {
1207      M[1]=intvec(1);
1208      M[2]=1;
1209      M[3]=intvec(0,1);
1210   }
1211   else
1212   {
1213      ideal J=lead(I);
1214      int b=ConductorBound(J);
1215      int i;
1216      list L=J[1];
1217      for(i=2;i<=size(J);i++)
1218      {
1219         L[i]=J[i];
1220      }
1221      M=WSemigroup(L,b);
1222      intvec v=0,M[3];
1223      M[3]=cutAfterConductor(v);
1224      M[2]=findConductor(M[3]);
1225   }
1226   return(M);
1227}
1228
1229example
1230{
1231 "EXAMPLE:"; echo=2;
1232ring R=0,t,ds;
1233ideal I=t8,t10+t13,t12+t15,t23-t29,t27;
1234semiGroup(I);
1235I=t8,t10+t13,t12+2t15,t27-3t33,t29;
1236semiGroup(I);
1237}
1238
1239////////////////////////////////////////////////////////////////////////////////
1240proc semiMod(ideal I,ideal G)
1241"USAGE":  semiMod(I,G);  I ideal,G ideal;I and G are the sagbi bases of the differential module resp.Algebra.
1242RETURN: list K;
1243      K[1]min generators of the semialgebra.
1244      K[2]conductor of the algebra.
1245      K[3]genrators for the semialgebra.
1246      K[4]min generators of the module.
1247      K[5]conductor of the module.
1248      K[6]semigroup of the module.
1249EXAMPLE: example semiMod;  shows an example
1250{
1251     list L=semiGroup(G);
1252     intvec M;
1253     list C;intvec S;
1254     int j; int k; int b;
1255     for(int i=1;i<=size(I);i++)
1256     {
1257         M[size(M)+1]=ord(I[i]);
1258     }
1259     M=M[2..size(M)];
1260     for(i=1;i<=size(M);i++)
1261     {
1262         C[size(C)+1]=M[i]+L[3];
1263     }
1264     int a=M[1]+L[2];
1265     for(j=1;j<=size(M);j++)
1266     {
1267        for(i=0;i<=a;i++)
1268        {
1269           for(k=1;k<=size(L[3]);k++)
1270           {
1271               if(i==C[j][k])
1272               {
1273                  S[size(S)+1]=i;
1274               }
1275           }
1276        }
1277      }
1278      S=S[2..size(S)];
1279      list K;
1280      K[1]=L[1];//generators of the semialgebra.
1281      K[2]=L[2];//conductor of the algebra.
1282      K[3]=L[3];//semi group of the algebra.
1283      K[4]=M;// generators of the semimodule.
1284      K[5]=findConductor(sortIntvec(S)); //conductor of the module.
1285      K[6]=cutAfterConductor(sortIntvec(S));//semigroup of the module.
1286      return(K);
1287}
1288example
1289{
1290  "EXAMPLE:"; echo=2;
1291  ring r=0,t,Ds;
1292  ideal G=t4,t7+t10;
1293  ideal I=diff(G,t);
1294  ideal k=sagbiMod(I,G);
1295  semiMod(k,G);
1296}
1297////////////////////////////////////////////////////////////////////////////////
1298static proc sagbiNF(poly f,ideal I,int b)
1299{
1300//computes the Sagbi normal form
1301   list L=1;
1302   map psi;
1303   f=jet(f,b);
1304   if(f==0){return(f);}
1305   while((f!=0) && (L[1]!=0))
1306   {
1307     L= algebra_containment(lead(f),lead(I),1);
1308     if (L[1]==1)
1309     {
1310        def S= L[2];
1311        psi= S,maxideal(1),I;
1312        f=jet(f-psi(check),b);
1313        kill S;
1314     }
1315   }
1316   return (lead(f)+sagbiNF(f-lead(f),I,b));
1317}
1318
1319/*
1320ring R=0,t,ds;
1321
1322ideal I=t5+t7,t4;
1323
1324sagbiNF(t7+2t9+3t11+t14+t13+6t15+t17,I,20);
1325
1326*/
1327////////////////////////////////////////////////////////////////////////////////
1328static proc sagbiSP(ideal I)
1329{
1330//computes the set of Sagbi-s-polys
1331   if(I==0){ return(I); }
1332   list L=algDependent(lead(I));
1333
1334   def S= L[2];
1335   map phi= S,maxideal(1),I;
1336   return(simplify(phi(ker),2));
1337}
1338
1339/*
1340
1341ring R=0,t,ds;
1342
1343ideal I=t4+t5,t7+t11,t9+t20;
1344
1345sagbiSP(I);
1346
1347*/
1348
1349////////////////////////////////////////////////////////////////////////////////
1350static proc sortSagbi(ideal I)
1351{
1352    //sorts, makes input monic and removes zeros
1353    I=simplify(I,2+1);
1354    int i;
1355    int n=1;
1356    poly p;
1357    while(n)
1358    {
1359       n=0;
1360       for(i=1;i<size(I);i++)
1361       {
1362           if(deg(lead(I[i]))>deg(lead(I[i+1])))
1363           {
1364                n=1;
1365                p=I[i];
1366                I[i]=I[i+1];
1367               I[i+1]=p;
1368               break;
1369            }
1370       }
1371    }
1372   return(I);
1373}
1374
1375/*
1376
1377ring R=0,t,ds;
1378
1379ideal I=3t5,7t2+t7,6t3+t8,3t+t7;
1380
1381sortSagbi(I);
1382
1383*/
1384
1385////////////////////////////////////////////////////////////////////////////////
1386static proc insertOne(poly p, ideal I, int b)
1387{
1388   //assume I is sorted, inserts p at the correct place
1389     int i,j;
1390     poly q;
1391     for(i=1;i<=size(I);i++)
1392     {
1393              if(deg(lead(p))<deg(lead(I[i])))
1394              {
1395                  break;
1396              }
1397     }
1398     if(i==size(I)+1)
1399     {
1400        I=I,simplify(p,1);
1401     }
1402     else
1403     {
1404         for(j=size(I)+1;j>i;j--)
1405         {
1406              I[j]=I[j-1];
1407         }
1408         I[i]=simplify(p,1);
1409     }
1410     if(i<size(I))
1411     {
1412           I=interReduceSagbi(I,b);
1413     }
1414     return(I);
1415}
1416
1417/*
1418
1419ring R=0,t,ds;
1420
1421ideal I=t8,t10+t13,t12+t15;
1422
1423insertOne(t17,I,20);
1424
1425I=t8,t10+t13,t12+t15,t23-t29;
1426
1427insertOne(-2t27,I,40);
1428
1429*/
1430
1431////////////////////////////////////////////////////////////////////////////////
1432static proc interReduceSagbi(ideal I, int b)
1433{
1434// reduces the elements of the dial against each other
1435     I=sortSagbi(I);
1436     ideal J;
1437     int n=1;
1438   int i;
1439     poly h;
1440     while(n)
1441     {
1442          n=0;
1443          i=1;
1444          while(i<size(I))
1445          {
1446              i++;
1447              J=I[1..i-1];
1448              h=sagbiNF(I[i],J,b);
1449              h=simplify(h,1);
1450              if(h!=I[i])
1451              {
1452                   n=1;
1453                   I[i]=h;
1454                   I=sortSagbi(I);
1455                   break;
1456              }
1457           }
1458      }
1459      return(I);
1460}
1461
1462/*
1463
1464    ring R=0,t,ds;
1465
1466    ideal I=t8,t8+t10+t13,t8+t12+t15;
1467
1468    interReduceSagbi(I,20);
1469
1470*/
1471////////////////////////////////////////////////////////////////////////////////
1472
1473static proc correctBound(ideal I, int b)
1474{
1475  //computes the conductor c of the semigroup associated to K[I]
1476  //if b>=c
1477   list L;
1478   int i;
1479   for(i=1;i<=size(I);i++)
1480   {
1481       L[i]=I[i];
1482   }
1483   list M=WSemigroup(L,b);
1484   if(b>M[2])
1485   {b=M[2]+1;}
1486   return(b);
1487}
1488
1489/*
1490
1491ring R=0,t,ds;
1492
1493ideal I=t8,t10+t13,t12+t15;
1494
1495correctBound(I,40);
1496
1497I=t8,t10+t13,t12+2t15;
1498
1499correctBound(I,40);
1500
1501*/
1502////////////////////////////////////////////////////////////////////////////////
1503static proc sortMinord(ideal I)
1504{
1505    //input an ideal
1506    //output a list L[1]=minimal order,
1507    //              L[2]=poly having the minimal order,
1508    //              L[3]=the k suchthat I[k] has the minimal order,
1509    //              L[4]=ideal I sorted in a way that minimal degree polynomial
1510
1511    //appears as the last polynomial of the ideal.ie I[size(I)]=I[k].
1512    int i;
1513    int n=1;
1514    list L;
1515    poly p;
1516    while(n)
1517    {
1518       n=0;
1519       for(i=1;i<size(I);i++)
1520       {
1521           if(ord(I[i])<ord(I[i+1]))
1522           {
1523                n=1;
1524                p=I[i];
1525                I[i]=I[i+1];
1526                I[i+1]=p;
1527                break;
1528            }
1529       }
1530    }
1531    L[1]=ord(I[size(I)]);
1532    L[2]=I[size(I)];
1533    L[3]=size(I);
1534    L[4]=I;
1535    return(L);
1536}
1537/*
1538ring r=0,t,Ds;
1539ideal I=t3,t6,t8,t4,t5,t9,t11,t3;
1540sortMinord(I);
1541*/
1542
1543////////////////////////////////////////////////////////////////////////////////
1544static proc inversP(poly p,int b)
1545{
1546  //computes the inverse of p upto the bound b
1547  if(size(p)==1)
1548  {
1549   return(p);
1550  }
1551  number c=leadcoef(p);
1552  p=p/c;
1553  poly q=1;
1554  poly s=1;
1555  while(deg(lead(q))<b)
1556  {
1557      q=q*(1-p);
1558      s=s+q;
1559  }
1560  s=1/c*jet(s,b);
1561  return(s);
1562}
1563
1564////////////////////////////////////////////////////////////////////////////////
1565static proc ConductorBound(ideal I)
1566{
1567    //input an ideal
1568    // output an integer which gives the bound of the semigroup conductor
1569    list M,L;
1570    int c,i,b;
1571    ideal J;
1572    poly p;
1573    if(size(I)<=1)
1574    {return(2);}
1575    while(1)
1576    {
1577        b=b+5;
1578        J=I;
1579        L=sortMinord(J);
1580        M[size(M)+1]=L[1];
1581        while((M[size(M)]!=1)&&(size(L[4])>1))
1582        {
1583             p=L[2]/var(1)^L[1];
1584             J=L[4];
1585             for(i=1;i<=L[3]-1;i++)
1586             {
1587               J[i]=J[i]/var(1)^L[1]*inversP(p,b);
1588               if(deg(lead(J[i]))==0){J[i]=J[i]-lead(J[i]);}
1589             }
1590             J=simplify(J,2);
1591             L=sortMinord(J);
1592             M[size(M)+1]=L[1];
1593        }
1594        if(M[size(M)]==1){break;}
1595    }
1596    for(i=1;i<=size(M)-1;i++)
1597    {
1598       c=c+M[i]*(M[i]-1);
1599    }
1600    return(c+1);
1601}
1602/*
1603ring r=0,t,Ds;
1604ideal I=t3+3t7,t8+5t9;
1605ConductorBound(I);
1606*/
1607////////////////////////////////////////////////////////////////////////////////
1608static proc sortMOD(ideal I)
1609{
1610    //sorts, makes input monic and removes zeros
1611    I=simplify(I,2);
1612    I=simplify(I,1);
1613    int i;
1614    int n=1;
1615    poly p;
1616    while(n)
1617    {
1618       n=0;
1619       for(i=1;i<size(I);i++)
1620       {
1621           if(deg(lead(I[i]))>deg(lead(I[i+1])))
1622           {
1623                n=1;
1624                p=I[i];
1625                I[i]=I[i+1];
1626                I[i+1]=p;
1627                break;
1628            }
1629       }
1630    }
1631   return(I);
1632}
1633////////////////////////////////////////////////////////////////////////////////
1634static proc SpolyMOD(ideal S,ideal P)
1635{
1636     //Assume that the basering is a ring in one variable.
1637    //input two ideals ideal S=<s_1,s_2> generators of the module and ideal P=<p_1,p_2,..,p_n> the sagbi basis of the algebra
1638    //output is an ideal generated by Q[p_1,p_2,...p_n]s_1-R[p_1,p_2,...p_n]s_2 for generators of
1639    //Q[lead(p_1),lead(p_2),.,lead(p_n)]lead(s_1)-R[lead(p_1),lead(p_2),.,lead(p_n)]lead(s_2)=0 .
1640   def br=basering;
1641   int n=ncols(P);
1642   ideal P1=lead(P);
1643   ideal S1=lead(S);
1644   execute
1645    ("ring T=("+charstr(br)+",x(1),z(1..n)),(y(1..2)),dp;");
1646    poly q;
1647    execute
1648   ("ring R=("+charstr(br)+"),(x(1),y(1..2),z(1..n)),(lp(3),dp(n));");
1649   map phi=br,x(1);
1650   ideal G=phi(P1);
1651   ideal I=phi(S1);
1652   ideal K,J;
1653   int d,o,s,j;
1654   poly q=I[1];
1655   if(deg(I[1])>deg(I[2]))
1656   {
1657       o=1;
1658       q=I[2];
1659   }
1660   I=I/q;
1661   for(int i=1;i<=2;i++)
1662   {
1663     K[i]=I[i]-y(i);
1664   }
1665   for(i=1;i<=n;i++)
1666   {
1667     K[2+i]=G[i]-z(i);
1668   }
1669   option(redSB);
1670   K=std(K);
1671   for(i=1;i<=size(K);i++)
1672   {
1673         if((K[i]/x(1)==0)&&((diff(K[i],y(1))!=0)||(diff(K[i],y(2))!=0)))
1674         {
1675                q=K[i];
1676                for(j=1;j<=2;j++)
1677                {
1678                    q=subst(q,y(j),0);
1679                }
1680                K[i]=K[i]-q+q*y(o+1);
1681                q=K[i];
1682                setring T;
1683                q=imap(R,q);
1684                s=deg(q);
1685                setring R;
1686                if(s==1){J[size(J)+1]=simplify(q,1);}
1687         }
1688   }
1689   setring br;
1690   map phi=R,maxideal(1),S,P;
1691   return(phi(J));
1692}
1693/*
1694ring r=0,t,dp;
1695ideal I=4t3,7t6+10t9;
1696ideal J=t4,t7+t10;
1697sortSagbi(SpolyMOD(I,J));
1698*/
1699////////////////////////////////////////////////////////////////////////////////
1700static proc sagbiNFMODO(poly p, ideal G, ideal I,int b)
1701{
1702    //input a poly ideal G ideal I int b is a bound
1703    //output an ideal K such that in each K[i] generators of I appear in linear.
1704    def br=basering;
1705    p=jet(p,b);
1706    if(p==0){return(p);}
1707    int n=ncols(G);
1708    int m=ncols(I);
1709    ideal G1=lead(G);
1710    ideal I1=lead(I);
1711     poly p1=lead(p);
1712    //create new ring with extra variables -
1713    execute
1714    ("ring T=("+charstr(br)+",x(1),z(1..n)),(x(2),y(1..m)),dp;");
1715     execute
1716    ("ring R=("+charstr(br)+"),(x(1..2),y(1..m),z(1..n)),(lp(2),dp(m),dp(n));");
1717
1718     map phi = br,x(1);
1719     ideal P = phi(G1);
1720     ideal S = phi(I1);
1721     poly check = phi(p1);
1722     poly keep=S[1];
1723     S=S/keep;
1724
1725     check=check/keep;
1726     ideal M;
1727     poly q;
1728     for (int i=1;i<=m;i=i+1)
1729     {
1730         M[i]=S[i]-y(i);
1731     }
1732     for (i=1;i<=n;i=i+1)
1733     {
1734        M[m+i]=P[i]-z(i);
1735     }
1736     M[size(M)+1]=check-x(2);
1737     check=check*keep;
1738
1739     option(redSB);
1740     M=std(M);
1741     int j,s;
1742
1743     for(i=1;i<=size(M);i++)
1744     {
1745         if((deg(M[i]/x(2))==0)&&(M[i]/x(1)==0))
1746         {
1747               q=subst(M[i],x(2),0);
1748               for(j=1;j<=m;j++)
1749               {
1750                  q=subst(q,y(j),0);
1751               }
1752               M[i]=M[i]-q+q*y(1);
1753             q=M[i];
1754               setring T;
1755               poly q=imap(R,q);
1756             s=deg(q);
1757               setring R;
1758               if(s==1){check=simplify(q,1);break;}
1759         }
1760     }
1761     setring br;
1762     map psi=R,maxideal(1),p,I,G;
1763     return(psi(check));
1764}
1765////////////////////////////////////////////////////////////////////////////////
1766
1767static proc sagbiNFMOD(poly p, ideal G, ideal I, int b)
1768{
1769    poly f=jet(p,b);
1770    if(f==0){return(f);}
1771    poly h;
1772    while(f!=h)
1773    {
1774         h=f;
1775         f=sagbiNFMODO(f,G,I,b);
1776    }
1777    return(lead(f)+sagbiNFMOD(p-lead(p),G,I,b));
1778}
1779////////////////////////////////////////////////////////////////////////////////
1780static proc createP(ideal I)
1781{
1782   list P;
1783   int i=1;
1784   int j;
1785   while(i<=size(I)-1)
1786   {
1787      j=i+1;
1788      while(j<=size(I))
1789      {
1790          P[size(P)+1]=list(I[i],I[j]);
1791          j++;
1792      }
1793      i++;
1794   }
1795   return(P);
1796}
1797////////////////////////////////////////////////////////////////////////////////
1798static proc enlargeP(poly h,list P,ideal I)
1799{
1800    int i;
1801    for(i=1;i<=size(I);i++)
1802    {
1803         P[size(P)+1]=list(I[i],h);
1804    }
1805    return(P);
1806}
1807/*
1808ring r=0,t,Ds;
1809ideal I=4t3,7t6+10t9;
1810ideal G=t4,t7+t10;
1811sagbiMod(I,G,18);
1812*/
1813
1814////////////////////////////////////////////////////////////////////////////////
1815static proc sortIntvec(intvec L)
1816{
1817    //input: intvec L.
1818    //output: L sorted, multiple elements canceled.
1819    int i;
1820    int j;
1821    int n=1;
1822    intvec M;
1823    while(n)
1824    {
1825        for(i=1;i<=size(L);i++)
1826        {
1827            for(j=i+1;j<=size(L);j++)
1828            {
1829                 if(L[i]==L[j])
1830                 {
1831                    L[j]=0;
1832                 }
1833             }
1834         }
1835         n=0;
1836     }
1837     for(i=1;i<=size(L);i++)
1838     {
1839         if((L[i]!=0)||(i==1))
1840         {
1841            M[size(M)+1]=L[i];
1842         }
1843     }
1844     int m=1;int p;
1845     while(m)
1846     {
1847        m=0;
1848        for(i=1;i<size(M);i++)
1849        {
1850           if(M[i]>M[i+1])
1851           {
1852                m=1;
1853                p=M[i];
1854                M[i]=M[i+1];
1855               M[i+1]=p;
1856               break;
1857            }
1858         }
1859      }
1860     M=M[2..size(M)];
1861     return(M);
1862}
1863////////////////////////////////////////////////////////////////////////////////
1864static proc findConductor(intvec L)
1865{
1866     //input a intvec L
1867     //output is an integer which came before the gap from right to left.
1868     int i;int j; list K;
1869     int c;
1870     for(i=size(L);i>=2;i--)
1871     {
1872        if(L[i]!=L[i-1]+1)
1873        {
1874            c=L[i];
1875            break;
1876        }
1877     }
1878     if(c==0){c=1;}
1879     return(c);
1880}
1881////////////////////////////////////////////////////////////////////////////////
1882static proc cutAfterConductor(intvec L)
1883{
1884     //input an integer vector
1885     //output cut all the integers in the intvec which came after the conductor
1886     int i;int j; intvec K;
1887     int c=findConductor(L);
1888     for(i=1;i<=size(L);i++)
1889     {
1890         if(L[i]==c)
1891         {
1892            K[1..i]=L[1..i];
1893         }
1894     }
1895     return(K);
1896}
1897////////////////////////////////////////////////////////////////////////////////
1898static proc CompareList(list L,list M,int n)
1899{
1900      //input two list L,M with the same size n
1901      //out put 0 if not equal 1 if equal.
1902      for(int i=1;i<=n;i++)
1903      {
1904           if(L[i]!=M[i])
1905           {
1906               i=0;
1907               break;
1908           }
1909      }
1910      return(i);
1911}
1912////////////////////////////////////////////////////////////////////////////////
1913static proc Guess(ideal I)
1914{
1915  // comput the sagbi basis of the module
1916  //which we guess .
1917   I=sagbiAlg(I);
1918   ideal H=diff(I,var(1));
1919   H=sagbiMod(H,I);
1920   list K=semiMod(H,I);
1921   return(K);
1922}
1923////////////////////////////////////////////////////////////////////////////////
1924/*
1925===============================   Examples==========================================
1926ideal I=t4+4t5+6t6+8t7+13t8+12t9+10t10+12t11+6t12+4t13+4t14+t16,t7+7t8+22t9+51t10+113t11+219t12+366t13+589t14+876t15+1170t16+1514t17+1828t1
19278+2011t19+2165t20+2163t21+1982t22+1806t23+1491t24+1141t25+889t26+588t27+379t28+2
192852t29+120t30+72t31+36t32+9t33+9t34+t36;
1929planeCur(I);
1930//=============================
1931ideal I=t4+4t5+6t6+8t7+13t8+12t9+10t10+12t11+6t12+4t13+4t14+t16,t7+7t8+21t9+42t10+77t11+126t12+168t13+211t14+252t15+252t16+245t17+231t18+17
19325t19+140t20+105t21+56t22+42t23+21t24+7t25+7t26+t28
1933planeCur(I);
1934//===============================
1935ideal I=t4+4t5+6t6+8t7+13t8+12t9+10t10+12t11+6t12+4t13+4t14+t16,t5+5t6+11t7+22t8+46t9+73t10+107t11+161t12+198t13+231t14+272t15+262t16+250t1
19367+236t18+175t19+141t20+105t21+56t22+42t23+21t24+7t25+7t26+t28
1937planeCur(I);
1938//===============================
1939ideal I=t4+4t5+6t6+8t7+13t8+12t9+10t10+12t11+6t12+4t13+4t14+t16,t6+7t7+22t8+47t9+87t10+143t11+202t12+258t13+307t14+332t15+327t16+305t17+266
1940t18+205t19+155t20+111t21+62t22+42t23+22t24+7t25+7t26+t28
1941planeCur(I);
1942//===============================
1943ideal I=t2+2t3+t4+2t5+2t6+t8,t+t2+t4;
1944planeCur(I);
1945//===============================
1946ideal I=t2+2t3+t4+2t5+2t6+t8,t3+3t4+3t5+4t6+6t7+3t8+3t9+3t10+t12;
1947planeCur(I);
1948//===============================
1949ideal I=t2+2t3+t4+2t5+2t6+t8,t5+5t6+10t7+15t8+25t9+31t10+30t11+35t12+30t13+20t14+20t15+10t16+5t17+5t18+t
195020;
1951planeCur(I);
1952//================================================================
1953ideal I=t2+2t3+t4+2t5+2t6+t8,t11+11t12+55t13+176t14+440t15+957t16+1837t17+3135t18+4917t19+7150t20+9581t2
19541+12046t22+14300t23+15851t24+16665t25+16687t26+15642t27+14025t28+12012t29+9570t3
19550+7392t31+5412t32+3630t33+2442t34+1485t35+825t36+495t37+220t38+110t39+55t40+11t4
19561+11t42+t44
1957planeCur(I);
1958//===============================
1959ideal I=t2+2t3+t4+2t5+2t6+t8,t45+45t46+990t47+14235t48+150975t49+1264329t50+8742030t51+51530985t52+26531
19607525t53+1216052255t54+5037384726t55+19091253735t56+66860434260t57+218159032410t5
19618+667743178590t59+1928258130018t60+5278946615910t61+13758022145340t62+3425642198
19621760t63+81743054778990t64+187438301870193t65+413998043743845t66+882643063827960t
196367+1819834573178925t68+3634672399863945t69+7042671464388093t70+13256726980146210
1964t71+24271349963247255t72+43270648586469315t73+75192560924341905t74+1274795590273
196539134t75+211037186585880765t76+341404127193205395t77+540109313678250885t78+83615
19662328502076770t79+1267494306126371433t80+1882391473790147350t81+27403488768330021
196760t82+3912426884928977910t83+5480608823069934180t84+7535946071701345419t85+10175
1968247273088233765t86+13496177050168252770t87+17590776929351920305t88+2253760903474
19699950330t89+28392934993342165732t90+35181553858703840610t91+42888103580926417860t
197092+51449748796644626670t93+60751205041524651720t94+70622965899108523296t95+80843
1971398349265488310t96+91145062374529367655t97+101225220090613564275t98+110760068529
1972877638960t99+119421810187582522995t100+126897320456330125725t101+132906930278955
1973392505t102+137221752614812709130t103+139678059865381605315t104+14018746206071963
19745683t105+138742016728357115865t106+135413875517988518550t107+1303495836626693311
197525t108+123759636437037165840t109+115904304930914703126t110+107077029168089360280
1976t111+97586814544772570280t112+87741050370279892245t113+77830012377996062865t114+
197768114044171037561004t115+58814074232856531765t116+50105762317964865600t117+42117
1978223130580686220t118+34929979773602146200t119+28582581501297657240t120+2307618932
19799698326690t121+18381388272325750530t122+14445518786710710480t123+111999120315284
198053530t124+8566543884036576384t125+6463772035817658320t126+4810966835075093880t12
19817+3531977599087147320t128+2557482632962404180t129+1826346112628778972t130+128615
19821054039308160t131+893096793855988260t132+611445912380539110t133+4126879484894709
198390t134+274559737461674588t135+180030436220988810t136+116328756134241090t137+7406
19841684381355110t138+46450833440621940t139+28695217633493598t140+17456561066064945t
1985141+10455665532950385t142+6164429567615550t143+3576677924170795t144+204174682346
19868917t145+1146414046643415t146+632953124099190t147+343522434444255t148+1832093883
198747205t149+95981896978935t150+49375510221510t151+24930700142535t152+1234956944936
19880t153+5998779092790t154+2855797655022t155+1331635383390t156+607860009900t157+271
1989401068250t158+118455934740t159+50498441136t160+20999419155t161+8518084355t162+33
199061582620t163+1290701115t164+481780299t165+173664315t166+61087950t167+20511645t16
19918+6704775t169+2115729t170+610170t171+191565t172+42570t173+15180t174+1980t175+990
1992t176+45t177+45t178+t180
1993planeCur(I);
1994//===============================
1995ideal I=t3+3t4+3t5+4t6+6t7+3t8+3t9+3t10+t12,t2+2t3+t4+2t5+2t6+t8
1996planeCur(I);
1997//===============================
1998 ideal I=t3+3t4+3t5+4t6+6t7+3t8+3t9+3t10+t12,t5+5t6+10t7+15t8+25t9+31t10+30t11+35t12+30t13+20t14+20t15+10t16+5t17+5t18+t20
1999 planeCur(I);
2000//===============================
2001ideal I=t3+3t4+3t5+4t6+6t7+3t8+3t9+3t10+t12,t4+4t5+6t6+8t7+13t8+12t9+10t10+12t11+6t12+4t13+4t14+t16
2002 planeCur(I);
2003//==========================================================================
2004 ring r=0,t,Ds;
2005 ideal I=t3,t10+t14;
2006 planeCur(I);
2007//===============================
2008ideal I=t3+3t4+3t5+t6,t10+10t11+45t12+120t13+211t14+266t15+301t16+484t17+1046t18+2012t19+3004t20+
20093432t21+3003t22+2002t23+1001t24+364t25+91t26+14t27+t28
2010planeCur(I);
2011//=======================================
2012ideal I=t3+3t4+3t5+t6,t10+10t11+45t12+120t13+210t14+252t15+210t16+120t17+45t18+10t19+t20
2013 planeCur(I);
2014//===============================
2015ring r=0,t,Ds;
2016ideal I=t3+3t4+3t5+t6,t13+14t14+92t15+377t16+1079t17+2288t18+3718t19+4719t20+4719t21+3718t22+2288
2017t23+1079t24+377t25+92t26+14t27+t28,t20+20t21+190t22+1140t23+4845t24+15504t25+38760t26+77520t27+125970t28+16796
20180t29+184756t30+167960t31+125970t32+77520t33+38760t34+15504t35+4845t36+1140t37+19
20190t38+20t39+t40
2020spaceCur(I);
2021//=====================================================
2022ideal I=t3+3t4+3t5+t6,t13+14t14+92t15+377t16+1079t17+2288t18+3718t19+4719t20+4719t21+3718t22+2288
2023t23+1079t24+377t25+92t26+14t27+t28,t17+17t18+136t19+680t20+2380t21+6188t22+12376t23+19448t24+24310t25+24310t26
2024+19448t27+12376t28+6188t29+2380t30+680t31+136t32+17t33+t34
2025spaceCur(I);
2026//========================================================
2027ideal I=t3,t16,t14;
2028spaceCur(I);
2029//=============================================
2030ideal I=t3,t19,t14;
2031spaceCur(I);
2032//==============================================
2033ideal I=t3,t14+t16,t19;
2034spaceCur(I);
2035//===============================================
2036ideal I=t3,t14+t16,t25;
2037spaceCur(I);
2038//=======================================
2039ideal I=t3+3t4+3t5+t6,t14+14t15+91t16+364t17+1001t18+2002t19+3003t20+3432t21+3004t22+2024t23+1232
2040t24+1904t25+7406t26+26348t27+74614t28+170544t29+319770t30+497420t31+646646t32+70
20415432t33+646646t34+497420t35+319770t36+170544t37+74613t38+26334t39+7315t40+1540t4
20421+231t42+22t43+t44,t25+25t26+300t27+2300t28+12650t29+53130t30+177100t31+480700t32+1081575t33+2
2043042975t34+3268760t35+4457400t36+5200300t37+5200300t38+4457400t39+3268760t40+2042
2044975t41+1081575t42+480700t43+177100t44+53130t45+12650t46+2300t47+300t48+25t49+t50
2045spaceCur(I);
2046//=========================================================
2047ideal I=t3+3t4+3t5+t6,t14+14t15+91t16+364t17+1001t18+2003t19+3022t20+3603t21+3972t22+5878t23+1262
20489t24+27496t25+50479t26+75596t27+92379t28+92378t29+75582t30+50388t31+27132t32+116
204928t33+3876t34+969t35+171t36+19t37+t38,t25+25t26+300t27+2300t28+12650t29+53130t30+177100t31+480700t32+1081575t33+2
2050042975t34+3268760t35+4457400t36+5200300t37+5200300t38+4457400t39+3268760t40+2042
2051975t41+1081575t42+480700t43+177100t44+53130t45+12650t46+2300t47+300t48+25t49+t50
2052spaceCur(I);
2053//==============================================================
2054ideal I=t3+3t4+3t5+t6,t14+14t15+92t16+380t17+1121t18+2562t19+4823t20+7800t21+11011t22+13442t23+13
2055871t24+11804t25+8099t26+4382t27+1821t28+560t29+120t30+16t31+t32,t19+19t20+171t21+969t22+3876t23+11628t24+27132t25+50388t26+75582t27+92378t2
20568+92378t29+75582t30+50388t31+27132t32+11628t33+3876t34+969t35+171t36+19t37+t38
2057spaceCur(I);
2058//======================================================================
2059ideal I=t3+3t4+3t5+t6,t14+14t15+92t16+380t17+1121t18+2562t19+4823t20+7800t21+11011t22+13442t23+13
2060871t24+11804t25+8099t26+4382t27+1821t28+560t29+120t30+16t31+t32,t25+25t26+300t27+2300t28+12650t29+53130t30+177100t31+480700t32+1081575t33+2
2061042975t34+3268760t35+4457400t36+5200300t37+5200300t38+4457400t39+3268760t40+2042
2062975t41+1081575t42+480700t43+177100t44+53130t45+12650t46+2300t47+300t48+25t49+t50
2063spaceCur(I);
2064//================================================================
2065ideal I=t3+3t4+3t5+t6,t16+16t17+120t18+560t19+1820t20+4368t21+8008t22+11440t23+12870t24+11440t25+
20668008t26+4368t27+1820t28+560t29+120t30+16t31+t32
2067,t14+14t15+91t16+364t17+1001t18+2002t19+3003t20+3432t21+3003t22+2002t23+1001
2068t24+364t25+91t26+14t27+t28
2069spaceCur(I);
2070//===========================================================================================
2071*/
Note: See TracBrowser for help on using the repository browser.