# source:git/Singular/LIB/grobcov.lib@e90fd1

spielwiese
Last change on this file since e90fd1 was e90fd1, checked in by Hans Schoenemann <hannes@…>, 9 years ago
new version: grobcov.lib (grobcovK)
• Property mode set to 100644
File size: 144.7 KB
Line
1//
2version="version grobcov.lib 4.0.1.2 Jan_2015 "; // $Id$
3category="General purpose";
4info="
5LIBRARY:  grobcov.lib   Groebner Cover for parametric ideals.
6
7          Comprehensive Groebner Systems, Groebner Cover, Canonical Forms,
8          Parametric Polynomial Systems, Dynamic Geometry, Loci, Envelop,
9          Constructible sets.
10          See
11
12          A. Montes A, M. Wibmer,
13          \"Groebner Bases for Polynomial Systems with parameters\",
14          Journal of Symbolic Computation 45 (2010) 1391-1425.
15          (http://www-ma2.upc.edu/~montes/).
16
17AUTHORS:  Antonio Montes , Hans Schoenemann.
18
19OVERVIEW:
20            In 2010, the library was designed to contain Montes-Wibmer's
21            algorithms for compute the canonical Groebner Cover of a
22            parametric ideal as described in the paper:
23
24            Montes A., Wibmer M.,
25            \"Groebner Bases for Polynomial Systems with parameters\".
26            Journal of Symbolic Computation 45 (2010) 1391-1425.
27
28            The central routine is grobcov. Given a parametric
29            ideal, grobcov outputs its Canonical Groebner Cover, consisting
30            of a set of pairs of (basis, segment). The basis (after
31            normalization) is the reduced Groebner basis for each point
32            of the segment. The segments are disjoint, locally closed
33            and correspond to constant lpp (leading power product)
34            of the basis, and are represented in canonical prime
35            representation. The segments are disjoint and cover the
36            whole parameter space. The output is canonical, it only
37            depends on the given parametric ideal and the monomial order.
38            This is much more than a simple Comprehensive Groebner System.
39            The algorithm grobcov allows options to solve partially the
40            problem when the whole automatic algorithm does not finish
41            in reasonable time.
42
43            grobcov uses a first algorithm cgsdr that outputs a disjoint
44            reduced Comprehensive Groebner System with constant lpp.
45            For this purpose, in this library, the implemented algorithm is
46            Kapur-Sun-Wang algorithm, because it is the most efficient
47            algorithm known for this purpose.
48
49            D. Kapur, Y. Sun, and D.K. Wang.
50            \"A New Algorithm for Computing Comprehensive Groebner Systems\".
51            Proceedings of ISSAC'2010, ACM Press, (2010), 29-36.
52
53            The library has evolved to include new applications of the
54            Groebner Cover, and new theoretical developments have been done.
55
56            A new set of routines (locus, locusdg, locusto) has been included to
57            compute loci of points. The routines are used in the Dynamic
58            Geometry software Geogebra. They are described in:
59
61            \''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\''.
62            Computer-Aided Design 56 (2014) 22-33.
63
64            Recently also routines for computing the envelop of a family
65            of curves (enverlop, envelopdg), to be used in Dynamic Geometry,
66            has been included and will be described in a forthcoming paper:
67
69             \''Envelops in Dynamic Geometry using the Groebner cover\''.
70
72            for computing the canonical form of a constructible set, given as a
73            union of locally closed sets. They are used in the new version for the
74            computation of loci and envelops. It determines the canonical locally closed
75            level sets of a constructuble. They will be described in a forthcoming paper:
76
77            A. Montes, J.M. Brunat,
78            \"Canonical representations of constructible sets\".
79
80            This version was finished on 31/01/2015
81
82NOTATIONS: All given and determined polynomials and ideals are in the
83@*         basering Q[a][x]; (a=parameters, x=variables)
84@*         After defining the ring, the main routines
85@*         grobcov, cgsdr,
86@*         generate the global rings
87@*         @R   (Q[a][x]),
88@*         @P   (Q[a]),
89@*         @RP  (Q[x,a])
90@*         that are used inside and killed before the output.
91@*         If you want to use some internal routine you must
92@*         create before the above rings by calling setglobalrings();
93@*         because some of the internal routines use these rings.
94@*         The call to the basic routines grobcov, cgsdr will
95@*         kill these rings.
96
97PROCEDURES:
98grobcov(F);  Is the basic routine giving the canonical
99             Groebner Cover of the parametric ideal F.
100             This routine accepts many options, that
101             allow to obtain results even when the canonical
102             computation does not finish in reasonable time.
103
104cgsdr(F); Is the procedure for obtaining a first disjoint,
105             reduced Comprehensive Groebner System that
106             is used in grobcov, that can also be used
107             independently if only the CGS is required.
108             It is a more efficient routine than buildtree
109             (the own routine of 2010 that is no more used).
110             Now, KSW algorithm is used.
111
112setglobalrings();  Generates the global rings @R, @P and @PR that
113              are respectively the rings Q[a][x], Q[a], Q[x,a].  (a=parameters,
114             x=variables) It is called inside each of the fundamental
115             routines of the library: grobcov, cgsdr, locus, locusdg and
116             killed before the output.
117             If the user want to use some other internal routine,
118             then setglobalrings() is to be called before, as most of
119             them use these rings.
120
121pdivi(f,F); Performs a pseudodivision of a parametric polynomial
122             by a parametric ideal.
123
124pnormalf(f,E,N);   Reduces a parametric polynomial f over V(E) \ V(N)
125             ( E is the null ideal and N the non-null ideal )
126             over the parameters.
127
128Crep(N,M);  Computes the canonical C-representation of V(N) \ V(M).
129
130Prep(N,M);  Computes the canonical P-representation of V(N) \ V(M).
131
132PtoCrep(L)  Starting from the canonical Prep of a locally closed set
133             computes its Crep.
134
135extend(GC);  When the grobcov of an ideal has been computed
136             with the default option ('ext',0) and the explicit
137             option ('rep',2) (which is not the default), then
138             one can call extend (GC) (and options) to obtain the
139             full representation of the bases. With the default
140             option ('ext',0) only the generic representation of
141             the bases are computed, and one can obtain the full
142             representation using extend.
143
144locus(G);    Special routine for determining the geometrical locus of points
145             verifying given conditions. Given a parametric ideal J with
146             parameters (x,y) and variables (x_1,..,xn), representing the
147             system determining the locus of points (x,y) who verify certain
148             properties, one can apply locus to the output of  grobcov(J),
149             locus determines the different classes of locus components.
150              described in the paper:
151
152             \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\",
153             M. Abanades, F. Botana, A. Montes, T. Recio,
154             Computer-Aided Design 56 (2014) 22-33.
155
156             The components can be 'Normal', 'Special', 'Accumulation', 'Degenerate'.
157             The output are the components is given in P-canonical form
158             It also detects automatically a possible point that is to be
159             avoided by the mover, whose coordinates must be the last two
160             coordinates in the definition of the ring. If such a point is detected,
161             then it eliminates the segments of the grobcov depending on the
162             point that is to be avoided.
163
164locusdg(G);  Is a special routine that determines the  'Relevant' components
165             of the locus in dynamic geometry. It is to be called to the output of locus
166             and selects from it the useful components.
167
168envelop(F,C); Special routine for determining the envelop of a family of curves
169             F in Q[x,y][x_1,..xn] depending on a ideal of constraints C in Q[x_1,..,xn].
170             It detemines the different components as well as its type:
171             'Normal', 'Special', 'Accumulation', 'Degenerate'. And
172             it also classifies the Special components, determining the
173             zero dimensional antiimage of the component and verifying if
174             the component is a special curve of the family or not.
175             It calls internally first grobcov and then locus with special options
176             to obtain the complete result.
177             The taxonomy that it provides, as well as the algorithms involved
178             will be described in a forthcoming paper:
179
181             \''Envelops in Dynamic Geometry using the Gr\"obner cover\''.
182
183
184envelopdg(ev); Is a special routine to determine the 'Relevant' components
185             of the envelop of a family of curves to be used in Dynamic Geometry.
186             It must be called to the output of envelop(F,C).
187
188locusto(L); Transforms the output of locus, locusdg, envelop and  envelopdg
189             into a string that can be reed from different computational systems.
190
191AddCons(L); Uses the routine AddConsP. Given a set of locally closed sets as
192             difference of of varieties (does not need to be in C-representation)
193             it returns the canonical P-representation of the constructible set
194             formed by the union of them. The result is formed by a set of embedded,
195             disjoint, locally closed sets (levels).
196
197AddConsP(L);  Given a set of locally closed P-components, it returns the
198             canonical P-representation of the constructible set
199             formed by the union of them, consisting of a set of embedded,
200             disjoint locally closed sets (levels).
202             of the constructible sets will be described in a forthcoming paper.
203
204             A. Montes, J.M. Brunat,
205             \"Canonical representations of constructible sets\".
206
208";
209
210LIB "primdec.lib";
211LIB "qhmoduli.lib";
212
213// ************ Begin of the grobcov library *********************
214
215// Library grobcov.lib
216// (Groebner Cover):
217// Release 0: (public)
218// Initial data: 21-1-2008
219// Uses buildtree for cgsdr
220// Final data: 3-7-2008
221// Release 2: (private)
222// Initial data: 6-9-2009
223// Last version using buildtree for cgsdr
224// Final data: 25-10-2011
225// Release B: (private)
226// Initial data: 1-7-2012
227// Uses both buildtree and KSW for cgsdr
228// Final data: 4-9-2012
229// Release G: (public)
230// Initial data: 4-9-2012
231// Uses KSW algorithm for cgsdr
232// Final data: 21-11-2013
233// Release K: (public)
234// Updated locus: 8-1-2015
236// Reformed many routines, examples and helps: 8-1-2015
237// New routines for computing the envelop of a family of curves: 22-1-2015
238// Final data: 22-1-2015
239
240//*************Auxiliary routines**************
241
242// elimintfromideal: elimine the constant numbers from an ideal
243//        (designed for W, nonnull conditions)
244// Input: ideal J
245// Output:ideal K with the elements of J that are non constants, in the
246//        ring K[x1,..,xm]
247static proc elimintfromideal(ideal J)
248{
249  int i;
250  int j=0;
251  ideal K;
252  if (size(J)==0){return(ideal(0));}
253  for (i=1;i<=ncols(J);i++){if (size(variables(J[i])) !=0){j++; K[j]=J[i];}}
254  return(K);
255}
256
257// delfromideal: deletes the i-th polynomial from the ideal F
258//    Works in any kind of ideal
259static proc delfromideal(ideal F, int i)
260{
261  int j;
262  ideal G;
263  if (size(F)<i){ERROR("delfromideal was called with incorrect arguments");}
264  if (size(F)<=1){return(ideal(0));}
265  if (i==0){return(F)};
266  for (j=1;j<=ncols(F);j++)
267  {
268    if (j!=i){G[ncols(G)+1]=F[j];}
269  }
270  return(G);
271}
272
273// delidfromid: deletes the polynomials in J that are in I
274// Input: ideals I, J
275// Output: the ideal J without the polynomials in I
276//   Works in any kind of ideal
277static proc delidfromid(ideal I, ideal J)
278{
279  int i; list r;
280  ideal JJ=J;
281  for (i=1;i<=size(I);i++)
282  {
283    r=memberpos(I[i],JJ);
284    if (r[1])
285    {
286      JJ=delfromideal(JJ,r[2]);
287    }
288  }
289  return(JJ);
290}
291
292// eliminates the ith element from a list or an intvec
293static proc elimfromlist(l, int i)
294{
295  if(typeof(l)=="list"){list L;}
296  if (typeof(l)=="intvec"){intvec L;}
297  if (typeof(l)=="ideal"){ideal L;}
298  int j;
299  if((size(l)==0) or (size(l)==1 and i!=1)){return(l);}
300  if (size(l)==1 and i==1){return(L);}
301  // L=l[1];
302  if(i==1)
303  {
304    for(j=2;j<=size(l);j++)
305    {
306      L[j-1]=l[j];
307    }
308  }
309  else
310  {
311    for(j=1;j<=i-1;j++)
312    {
313      L[j]=l[j];
314    }
315    for(j=i+1;j<=size(l);j++)
316    {
317      L[j-1]=l[j];
318    }
319  }
320  return(L);
321}
322
323// eliminates repeated elements form an ideal or matrix or module or intmat or bigintmat
324static proc elimrepeated(F)
325{
326  int i;
327  int nt;
328  if (typeof(F)=="ideal"){nt=ncols(F);}
329  else{nt=size(F);}
330
331  def FF=F;
332  FF=F[1];
333  for (i=2;i<=nt;i++)
334  {
335    if (not(memberpos(F[i],FF)[1]))
336    {
337      FF[size(FF)+1]=F[i];
338    }
339  }
340  return(FF);
341}
342
343
344// equalideals
345// Input: ideals F and G;
346// Output: 1 if they are identical (the same polynomials in the same order)
347//         0 else
348static proc equalideals(ideal F, ideal G)
349{
350  int i=1; int t=1;
351  if (size(F)!=size(G)){return(0);}
352  while ((i<=size(F)) and (t))
353  {
354    if (F[i]!=G[i]){t=0;}
355    i++;
356  }
357  return(t);
358}
359
360// returns 1 if the two lists of ideals are equal and 0 if not
361static proc equallistideals(list L, list M)
362{
363  int t; int i;
364  if (size(L)!=size(M)){return(0);}
365  else
366  {
367    t=1;
368    if (size(L)>0)
369    {
370      i=1;
371      while ((t) and (i<=size(L)))
372      {
373        if (equalideals(L[i],M[i])==0){t=0;}
374        i++;
375      }
376    }
377    return(t);
378  }
379}
380
381// idcontains
382// Input: ideal p, ideal q
383// Output: 1 if p contains q,  0 otherwise
384// If the routine is to be called from the top, a previous call to
385// setglobalrings() is needed.
386static proc idcontains(ideal p, ideal q)
387{
388  int t; int i;
389  t=1; i=1;
390  def P=p; def Q=q;
391  attrib(P,"isSB",1);
392  poly r;
393  while ((t) and (i<=size(Q)))
394  {
395    r=reduce(Q[i],P);
396    if (r!=0){t=0;}
397    i++;
398  }
399  return(t);
400}
401
402
403// selectminideals
404//   given a list of ideals returns the list of integers corresponding
405//   to the minimal ideals in the list
406// Input: L (list of ideals)
407// Output: the list of integers corresponding to the minimal ideals in L
408//   Works in Q[u_1,..,u_m]
409static proc selectminideals(list L)
410{
411  list P; int i; int j; int t;
412  if(size(L)==0){return(L)};
413  if(size(L)==1){P[1]=1; return(P);}
414  for (i=1;i<=size(L);i++)
415  {
416    t=1;
417    j=1;
418    while ((t) and (j<=size(L)))
419    {
420      if (i!=j)
421      {
422        if(idcontains(L[i],L[j])==1)
423        {
424          t=0;
425        }
426      }
427      j++;
428    }
429    if (t){P[size(P)+1]=i;}
430  }
431  return(P);
432}
433
434
435// Auxiliary routine
436// elimconstfac: eliminate the factors in the polynom f that are in Q[a]
437// Input:
438//   poly f:
439//   list L: of components of the segment
440// Output:
441//   poly f2  where the factors of f in Q[a] that are non-null on any component
442//   have been dropped from f
443static proc elimconstfac(poly f)
444{
445  int cond; int i; int j; int t;
446  if (f==0){return(f);}
447  def RR=basering;
448  setring(@R);
449  def ff=imap(RR,f);
450  def l=factorize(ff,0);
451  poly f1=1;
452  for(i=2;i<=size(l[1]);i++)
453  {
454      f1=f1*(l[1][i])^(l[2][i]);
455  }
456  setring(RR);
457  def f2=imap(@R,f1);
458  return(f2);
459};
460
461static proc memberpos(f,J)
462//"USAGE:  memberpos(f,J);
463//         (f,J) expected (polynomial,ideal)
464//               or       (int,list(int))
465//               or       (int,intvec)
466//               or       (intvec,list(intvec))
467//               or       (list(int),list(list(int)))
468//               or       (ideal,list(ideal))
469//               or       (list(intvec),  list(list(intvec))).
470//         Works in any kind of ideals
471//RETURN:  The list (t,pos) t int; pos int;
472//         t is 1 if f belongs to J and 0 if not.
473//         pos gives the position in J (or 0 if f does not belong).
474//EXAMPLE: memberpos; shows an example"
475{
476  int pos=0;
477  int i=1;
478  int j;
479  int t=0;
480  int nt;
481  if (typeof(J)=="ideal"){nt=ncols(J);}
482  else{nt=size(J);}
483  if ((typeof(f)=="poly") or (typeof(f)=="int"))
484  { // (poly,ideal)  or
485    // (poly,list(poly))
486    // (int,list(int)) or
487    // (int,intvec)
488    i=1;
489    while(i<=nt)
490    {
491      if (f==J[i]){return(list(1,i));}
492      i++;
493    }
494    return(list(0,0));
495  }
496  else
497  {
498    if ((typeof(f)=="intvec") or ((typeof(f)=="list") and (typeof(f[1])=="int")))
499    { // (intvec,list(intvec)) or
500      // (list(int),list(list(int)))
501      i=1;
502      t=0;
503      pos=0;
504      while((i<=nt) and (t==0))
505      {
506        t=1;
507        j=1;
508        if (size(f)!=size(J[i])){t=0;}
509        else
510        {
511          while ((j<=size(f)) and t)
512          {
513            if (f[j]!=J[i][j]){t=0;}
514            j++;
515          }
516        }
517        if (t){pos=i;}
518        i++;
519      }
520      if (t){return(list(1,pos));}
521      else{return(list(0,0));}
522    }
523    else
524    {
525      if (typeof(f)=="ideal")
526      { // (ideal,list(ideal))
527        i=1;
528        t=0;
529        pos=0;
530        while((i<=nt) and (t==0))
531        {
532          t=1;
533          j=1;
534          if (ncols(f)!=ncols(J[i])){t=0;}
535          else
536          {
537            while ((j<=ncols(f)) and t)
538            {
539              if (f[j]!=J[i][j]){t=0;}
540              j++;
541            }
542          }
543          if (t){pos=i;}
544          i++;
545        }
546        if (t){return(list(1,pos));}
547        else{return(list(0,0));}
548      }
549      else
550      {
551        if ((typeof(f)=="list") and (typeof(f[1])=="intvec"))
552        { // (list(intvec),list(list(intvec)))
553          i=1;
554          t=0;
555          pos=0;
556          while((i<=nt) and (t==0))
557          {
558            t=1;
559            j=1;
560            if (size(f)!=size(J[i])){t=0;}
561            else
562            {
563              while ((j<=size(f)) and t)
564              {
565                if (f[j]!=J[i][j]){t=0;}
566                j++;
567              }
568            }
569            if (t){pos=i;}
570            i++;
571          }
572          if (t){return(list(1,pos));}
573          else{return(list(0,0));}
574        }
575      }
576    }
577  }
578}
579//example
580//{ "EXAMPLE:"; echo = 2;
581//  list L=(7,4,5,1,1,4,9);
582//  memberpos(1,L);
583//}
584
585// Auxiliary routine
586// pos
587// Input:  intvec p of zeros and ones
588// Output: intvec W of the positions where p has ones.
589static proc pos(intvec p)
590{
591  int i;
592  intvec W; int j=1;
593  for (i=1; i<=size(p); i++)
594  {
595    if (p[i]==1){W[j]=i; j++;}
596  }
597  return(W);
598}
599
600// Input:
601//  A,B: lists of ideals
602// Output:
603//   1 if both lists of ideals are equal, or 0 if not
604static proc equallistsofideals(list A, list B)
605{
606 int i;
607 int tes=0;
608 if (size(A)!=size(B)){return(tes);}
609 tes=1; i=1;
610 while(tes==1 and i<=size(A))
611 {
612   if (equalideals(A[i],B[i])==0){tes=0; return(tes);}
613   i++;
614 }
615 return(tes);
616}
617
618// Input:
619//  A,B:  lists of P-rep, i.e. of the form [p_i,[p_{i1},..,p_{ij_i}]]
620// Output:
621//   1 if both lists of P-reps are equal, or 0 if not
622static proc equallistsA(list A, list B)
623{
624  int tes=0;
625  if(equalideals(A[1],B[1])==0){return(tes);}
626  tes=equallistsofideals(A[2],B[2]);
627  return(tes);
628}
629
630// Input:
631//  A,B:  lists lists of of P-rep, i.e. of the form [[p_1,[p_{11},..,p_{1j_1}]] .. [p_i,[p_{i1},..,p_{ij_i}]]
632// Output:
633//   1 if both lists of lists of P-rep are equal, or 0 if not
634static proc equallistsAall(list A,list B)
635{
636 int i; int tes;
637 if(size(A)!=size(B)){return(tes);}
638 tes=1; i=1;
639 while(tes and i<=size(A))
640 {
641   tes=equallistsA(A[i],B[i]);
642   i++;
643 }
644 return(tes);
645}
646
647// idint: ideal intersection
648//        in the ring @P.
649//        it works in an extended ring
650// input: two ideals in the ring @P
651// output the intersection of both (is not a GB)
652static proc idint(ideal I, ideal J)
653{
654  def RR=basering;
655  ring T=0,t,lp;
656  def K=T+RR;
657  setring(K);
658  def Ia=imap(RR,I);
659  def Ja=imap(RR,J);
660  ideal IJ;
661  int i;
662  for(i=1;i<=size(Ia);i++){IJ[i]=t*Ia[i];}
663  for(i=1;i<=size(Ja);i++){IJ[size(Ia)+i]=(1-t)*Ja[i];}
664  ideal eIJ=eliminate(IJ,t);
665  setring(RR);
666  return(imap(K,eIJ));
667}
668
669//purpose ideal intersection called in @R and computed in @P
670static proc idintR(ideal N, ideal M)
671{
672  def RR=basering;
673  setring(@P);
674  def Np=imap(RR,N);
675  def Mp=imap(RR,M);
676  def Jp=idint(Np,Mp);
677  setring(RR);
678  return(imap(@P,Jp));
679}
680
681// Auxiliary routine
682// comb: the list of combinations of elements (1,..n) of order p
683static proc comb(int n, int p)
684{
685  list L; list L0;
686  intvec c; intvec d;
687  int i; int j; int last;
688  if ((n<0) or (n<p))
689  {
690    return(L);
691  }
692  if (p==1)
693  {
694    for (i=1;i<=n;i++)
695    {
696      c=i;
697      L[size(L)+1]=c;
698    }
699    return(L);
700  }
701  else
702  {
703    L0=comb(n,p-1);
704    for (i=1;i<=size(L0);i++)
705    {
706      c=L0[i]; d=c;
707      last=c[size(c)];
708      for (j=last+1;j<=n;j++)
709      {
710        d[size(c)+1]=j;
711        L[size(L)+1]=d;
712      }
713    }
714    return(L);
715  }
716}
717
718// Auxiliary routine
719// combrep
720// Input: V=(n_1,..,n_i)
721// Output: L=(v_1,..,v_p) where p=prod_j=1^i (n_j)
722//    is the list of all intvec v_j=(v_j1,..,v_ji) where 1<=v_jk<=n_i
723static proc combrep(intvec V)
724{
725  list L; list LL;
726  int i; int j; int k;  intvec W;
727  if (size(V)==1)
728  {
729    for (i=1;i<=V[1];i++)
730    {
731      L[i]=intvec(i);
732    }
733    return(L);
734  }
735  for (i=1;i<size(V);i++)
736  {
737    W[i]=V[i];
738  }
739  LL=combrep(W);
740  for (i=1;i<=size(LL);i++)
741  {
742    W=LL[i];
743    for (j=1;j<=V[size(V)];j++)
744    {
745      W[size(V)]=j;
746      L[size(L)+1]=W;
747    }
748  }
749  return(L);
750}
751
752static proc subset(J,K)
753//"USAGE:   subset(J,K);
754//          (J,K)  expected (ideal,ideal)
755//                  or     (list, list)
756//RETURN:   1 if all the elements of J are in K, 0 if not.
757//EXAMPLE:  subset; shows an example;"
758{
759  int i=1;
760  int nt;
761  if (typeof(J)=="ideal"){nt=ncols(J);}
762  else{nt=size(J);}
763  if (size(J)==0){return(1);}
764  while(i<=nt)
765  {
766    if (memberpos(J[i],K)[1]){i++;}
767    else {return(0);}
768  }
769  return(1);
770}
771//example
772//{ "EXAMPLE:"; echo = 2;
773//  list J=list(7,3,2);
774//  list K=list(1,2,3,5,7,8);
775//  subset(J,K);
776//}
777
778proc setglobalrings()
779"USAGE:   setglobalrings();
780          No arguments
781RETURN:   After its call the rings @R=Q[a][x], @P=Q[a], @RP=Q[x,a] are
782          defined as global variables.  (a=parameters, x=variables)
783NOTE: It is called internally by many basic routines of the
784          library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
785          envelop, envelopdg, and killed before the output.
786          The user does not need to call it, except when it is interested
787          in using some internal routine of the library that
788          uses these rings.
789          The basering R, must be of the form Q[a][x], (a=parameters,
790          x=variables), and should be defined previously.
791KEYWORDS: ring, rings
792EXAMPLE:  setglobalrings; shows an example"
793{
794  if (defined(@P))
795  {
796    kill @P; kill @R; kill @RP;
797  }
798  def RR=basering;
799  def @R=basering;  // must be of the form Q[a][x], (a=parameters, x=variables)
800  def Rx=ringlist(RR);
801  def @P=ring(Rx[1]);
802  list Lx;
803  Lx[1]=0;
804  Lx[2]=Rx[2]+Rx[1][2];
805  Lx[3]=Rx[1][3];
806  Lx[4]=Rx[1][4];
807  Rx[1]=0;
808  def D=ring(Rx);
809  def @RP=D+@P;
810  export(@R);      // global ring Q[a][x]
811  export(@P);      // global ring Q[a]
812  export(@RP);     // global ring K[x,a] with product order
813  setring(RR);
814};
815example
816{
817  "EXAMPLE:"; echo = 2;
818  ring R=(0,a,b),(x,y,z),dp;
819  setglobalrings();
820  R;
821  Grobcov::@R;
822  Grobcov::@P;
823  Grobcov::@RP;
824  ringlist(Grobcov::@R);
825  ringlist(Grobcov::@P);
826  ringlist(Grobcov::@RP);
827}
828
829// cld : clears denominators of an ideal and normalizes to content 1
830//        can be used in @R or @P or @RP
831// input:
832//        ideal J (J can be also poly), but the output is an ideal;
833// output:
834//        ideal Jc (the new form of ideal J without denominators and
835//        normalized to content 1)
836static proc cld(ideal J)
837{
838  if (size(J)==0){return(ideal(0));}
839  int te=0;
840  def RR=basering;
841  if(not(defined(@RP)))
842  {
843    te=1;
844    setglobalrings();
845  }
846  setring(@RP);
847  def Ja=imap(RR,J);
848  ideal Jb;
849  if (size(Ja)==0){setring(RR); return(ideal(0));}
850  int i;
851  def j=0;
852  for (i=1;i<=ncols(Ja);i++){if (size(Ja[i])!=0){j++; Jb[j]=cleardenom(Ja[i]);}}
853  setring(RR);
854  def Jc=imap(@RP,Jb);
855  if(te){kill @R; kill @RP; kill @P;}
856  return(Jc);
857};
858
859// simpqcoeffs : simplifies a quotient of two polynomials
860// input: two coeficients (or terms), that are considered as a quotient
861// output: the two coeficients reduced without common factors
862static proc simpqcoeffs(poly n,poly m)
863{
864  def nc=content(n);
865  def mc=content(m);
866  def gc=gcd(nc,mc);
867  ideal s=n/gc,m/gc;
868  return (s);
869}
870
871// pdivi : pseudodivision of a parametric polynomial f by a parametric ideal F in Q[a][x].
872// input:
873//   poly  f
874//   ideal F
875// output:
876//   list (poly r, ideal q, poly mu)
877proc pdivi(poly f,ideal F)
878"USAGE: pdivi(f,F);
879          poly f: the polynomialin Q[a][x] to be divided
880          ideal F: the divisor ideal in Q[a][x].
881RETURN: A list (poly r, ideal q, poly m). r is the remainder of the
882          pseudodivision, q is the set of quotients, and m is the
883          coefficient factor by which f is to be multiplied.
884NOTE: pseudodivision of a poly f by an ideal F in Q[a][x]. Returns a
885          list (r,q,m) such that m*f=r+sum(q.G), and no lpp of a divisor
886          divides a pp of r.
887KEYWORDS: division, reduce
888EXAMPLE:  pdivi; shows an example"
889{
890  F=simplify(F,2);
891  int i;
892  int j;
893  poly v=1;
894  for(i=1;i<=nvars(basering);i++){v=v*var(i);}
895  poly r=0;
896  poly mu=1;
897  def p=f;
898  ideal q;
899  for (i=1; i<=ncols(F); i++){q[i]=0;};
900  ideal lpf;
901  ideal lcf;
904  poly lpp;
905  poly lcp;
906  poly qlm;
907  poly nu;
908  poly rho;
909  int divoc=0;
910  ideal qlc;
911  while (p!=0)
912  {
913    i=1;
914    divoc=0;
917    while (divoc==0 and i<=size(F))
918    {
919      qlm=lpp/lpf[i];
920      if (qlm!=0)
921      {
922        qlc=simpqcoeffs(lcp,lcf[i]);
923        nu=qlc[2];
924        mu=mu*nu;
925        rho=qlc[1]*qlm;
926        p=nu*p-rho*F[i];
927        r=nu*r;
928        for (j=1;j<=size(F);j++){q[j]=nu*q[j];}
929        q[i]=q[i]+rho;
930        divoc=1;
931      }
932      else {i++;}
933    }
934    if (divoc==0)
935    {
936      r=r+lcp*lpp;
937      p=p-lcp*lpp;
938    }
939  }
940  list res=r,q,mu;
941  return(res);
942}
943example
944{
945  "EXAMPLE:"; echo = 2;
946  ring R=(0,a,b,c),(x,y),dp;
947  poly f=(ab-ac)*xy+(ab)*x+(5c);
948  // Divisor=";
949  f;
950  ideal F=ax+b,cy+a;
951  // Dividends=";
952  F;
953  def r=pdivi(f,F);
954  // (Remainder, quotients, factor)=";
955  r;
956  // Verifying the division:
957  r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2]+r[1]);
958}
959
960// pspol : S-poly of two polynomials in @R
961// @R
962// input:
963//   poly f (given in the ring @R)
964//   poly g (given in the ring @R)
965// output:
966//   list (S, red):  S is the S-poly(f,g) and red is a Boolean variable
967//                if red then S reduces by Buchberger 1st criterion
968//                (not used)
969static proc pspol(poly f,poly g)
970{
975  def v=gcd(lpf,lpg);
976  def s=simpqcoeffs(lcf,lcg);
977  def vf=lpf/v;
978  def vg=lpg/v;
979  poly S=s[2]*vg*f-s[1]*vf*g;
980  return(S);
981}
982
983// facvar: Returns all the free-square factors of the elements
984//         of ideal J (non repeated). Integer factors are ignored,
985//         even 0 is ignored. It can be called from ideal @R, but
986//         the given ideal J must only contain poynomials in the
987//         parameters.
988//         Operates in the ring @P, but can be called from ring @R,
989//         and the ideal @P must be defined calling first setglobalrings();
990// input:  ideal J
991// output: ideal Jc: Returns all the free-square factors of the elements
992//         of ideal J (non repeated). Integer factors are ignored,
993//         even 0 is ignored. It can be called from ideal @R.
994static proc facvar(ideal J)
995//"USAGE:   facvar(J);
996//          J: an ideal in the parameters
997//RETURN:   all the free-square factors of the elements
998//          of ideal J (non repeated). Integer factors are ignored,
999//          even 0 is ignored. It can be called from ideal @R, but
1000//          the given ideal J must only contain poynomials in the
1001//          parameters.
1002//NOTE:     Operates in the ring @P, and the ideal J must contain only
1003//          polynomials in the parameters, but can be called from ring @R.
1004//KEYWORDS: factor
1005//EXAMPLE:  facvar; shows an example"
1006{
1007  int i;
1008  def RR=basering;
1009  setring(@P);
1010  def Ja=imap(RR,J);
1011  if(size(Ja)==0){setring(RR); return(ideal(0));}
1012  Ja=elimintfromideal(Ja); // also in ideal @P
1013  ideal Jb;
1014  if (size(Ja)==0){Jb=ideal(0);}
1015  else
1016  {
1017    for (i=1;i<=ncols(Ja);i++){if(size(Ja[i])!=0){Jb=Jb,factorize(Ja[i],1);}}
1018    Jb=simplify(Jb,2+4+8);
1019    Jb=cld(Jb);
1020    Jb=elimintfromideal(Jb); // also in ideal @P
1021  }
1022  setring(RR);
1023  def Jc=imap(@P,Jb);
1024  return(Jc);
1025}
1026//example
1027//{ "EXAMPLE:"; echo = 2;
1028//  ring R=(0,a,b,c),(x,y,z),dp;
1029//  setglobalrings();
1030//  ideal J=a2-b2,a2-2ab+b2,abc-bc;
1031//  facvar(J);
1032//}
1033
1034// Ered: eliminates the factors in the polynom f that are non-null.
1035//       In ring @R
1036// input:
1037//   poly f:
1038//   ideal E  of null-conditions
1039//   ideal N  of non-null conditions
1040//        (E,N) represents V(E) \ V(N),
1041//        Ered eliminates the non-null factors of f in V(E) \ V(N)
1042// output:
1043//   poly f2  where the non-null conditions have been dropped from f
1044static proc Ered(poly f,ideal E, ideal N)
1045{
1046  def RR=basering;
1047  setring(@R);
1048  poly ff=imap(RR,f);
1049  ideal EE=imap(RR,E);
1050  ideal NN=imap(RR,N);
1051  if((ff==0) or (equalideals(NN,ideal(1)))){setring(RR); return(f);}
1052  def v=variables(ff);
1053  int i;
1054  poly X=1;
1055  for(i=1;i<=size(v);i++){X=X*v[i];}
1056  matrix M=coef(ff,X);
1057  setring(@P);
1058  def RPE=imap(@R,EE);
1059  def RPN=imap(@R,NN);
1060  matrix Mp=imap(@R,M);
1061  poly g=Mp[2,1];
1062  if (size(Mp)!=2)
1063  {
1064    for(i=2;i<=size(Mp) div 2;i++)
1065    {
1066      g=gcd(g,Mp[2,i]);
1067    }
1068  }
1069  if (g==1){setring(RR); return(f);}
1070  else
1071  {
1072    def wg=factorize(g,2);
1073    if (wg[1][1]==1){setring(RR); return(f);}
1074    else
1075    {
1076      poly simp=1;
1077      int te;
1078      for(i=1;i<=size(wg[1]);i++)
1079      {
1080        te=inconsistent(RPE+wg[1][i],RPN);
1081        if(te)
1082        {
1083          simp=simp*(wg[1][i])^(wg[2][i]);
1084        }
1085      }
1086    }
1087    if (simp==1){setring(RR); return(f);}
1088    else
1089    {
1090      setring(RR);
1091      def simp0=imap(@P,simp);
1092      def f2=f/simp0;
1093      return(f2);
1094    }
1095  }
1096}
1097
1098// pnormalf: reduces a polynomial f wrt a V(E) \ V(N)
1099//           dividing by E and eliminating factors in N.
1100//           called in the ring @R,
1101//           operates in the ring @RP.
1102// input:
1103//         poly  f
1104//         ideal E  (depends only on the parameters)
1105//         ideal N  (depends only on the parameters)
1106//                  (E,N) represents V(E) \ V(N)
1107//         optional:
1108// output: poly f2 reduced wrt to V(E) \ V(N)
1109proc pnormalf(poly f, ideal E, ideal N)
1110"USAGE: pnormalf(f,E,N);
1111          f: the polynomial in Q[a][x]  (a=parameters, x=variables) to be
1112          reduced modulo V(E) \ V(N) of a segment in Q[a].
1113          E: the null conditions ideal in Q[a]
1114          N: the non-null conditions in Q[a]
1115RETURN: a reduced polynomial g of f, whose coefficients are reduced
1116          modulo E and having no factor in N.
1117NOTE: Should be called from ring Q[a][x].
1118          Ideals E and N must be given by polynomials in Q[a].
1119KEYWORDS: division, pdivi, reduce
1120EXAMPLE: pnormalf; shows an example"
1121{
1122    def RR=basering;
1123    int te=0;
1124    if (defined(@P)){te=1;}
1125    else{setglobalrings();}
1126    setring(@RP);
1127    def fa=imap(RR,f);
1128    def Ea=imap(RR,E);
1129    def Na=imap(RR,N);
1130    option(redSB);
1131    Ea=std(Ea);
1132    def r=cld(reduce(fa,Ea));
1133    poly f1=r[1];
1134    f1=Ered(r[1],Ea,Na);
1135    setring(RR);
1136    def f2=imap(@RP,f1);
1137    if(te==0){kill @R; kill @RP; kill @P;}
1138    return(f2)
1139};
1140example
1141{
1142  "EXAMPLE:"; echo = 2;
1143  ring R=(0,a,b,c),(x,y),dp;
1144  short=0;
1145  poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
1146  ideal E=(c-1);
1147  ideal N=a-b;
1148
1149  pnormalf(f,E,N);
1150}
1151
1152// lesspol: compare two polynomials by its leading power products
1153// input:  two polynomials f,g in the ring @R
1154// output: 0 if f<g,  1 if f>=g
1155static proc lesspol(poly f, poly g)
1156{
1158  else
1159  {
1161    else
1162    {
1164      else {return(0);}
1165    }
1166  }
1167};
1168
1169// sortideal: sorts the polynomials in an ideal by lm in ascending order
1170static proc sortideal(ideal Fi)
1171{
1172  def RR=basering;
1173  setring(@RP);
1174  def F=imap(RR,Fi);
1175  def H=F;
1176  ideal G;
1177  int i;
1178  int j;
1179  poly p;
1180  while (size(H)!=0)
1181  {
1182    j=1;
1183    p=H[1];
1184    for (i=1;i<=ncols(H);i++)
1185    {
1186      if(lesspol(H[i],p)){j=i;p=H[j];}
1187    }
1188    G[ncols(G)+1]=p;
1189    H=delfromideal(H,j);
1190    H=simplify(H,2);
1191  }
1192  setring(RR);
1193  def GG=imap(@RP,G);
1194  GG=simplify(GG,2);
1195  return(GG);
1196}
1197
1198// mingb: given a basis (gb reducing) it
1199// order the polynomials in ascending order and
1200// eliminates the polynomials whose lpp are divisible by some
1201// smaller one
1202static proc mingb(ideal F)
1203{
1204  int t; int i; int j;
1205  def H=sortideal(F);
1206  ideal G;
1207  if (ncols(H)<=1){return(H);}
1208  G=H[1];
1209  for (i=2; i<=ncols(H); i++)
1210  {
1211    t=1;
1212    j=1;
1213    while (t and (j<i))
1214    {
1216      j++;
1217    }
1218    if (t) {G[size(G)+1]=H[i];}
1219  }
1220  return(G);
1221}
1222
1223// redgbn: given a minimal basis (gb reducing) it
1224// reduces each polynomial wrt to V(E) \ V(N)
1225static proc redgbn(ideal F, ideal E, ideal N)
1226{
1227  int te=0;
1228  if (defined(@P)==1){te=1;}
1229  ideal G=F;
1230  ideal H;
1231  int i;
1232  if (size(G)==0){return(ideal(0));}
1233  for (i=1;i<=size(G);i++)
1234  {
1235    H=delfromideal(G,i);
1236    G[i]=pnormalf(pdivi(G[i],H)[1],E,N);
1237    G[i]=primepartZ(G[i]);
1238  }
1239  if(te==1){setglobalrings();}
1240  return(G);
1241}
1242
1243//**************Begin homogenizing************************
1244
1245// ishomog:
1246// Purpose: test if a polynomial is homogeneous in the variables or not
1247// input:  poly f
1248// output  1 if f is homogeneous, 0 if not
1249static proc ishomog(f)
1250{
1251  int i; poly r; int d; int dr;
1252  if (f==0){return(1);}
1253  d=deg(f); dr=d; r=f;
1254  while ((d==dr) and (r!=0))
1255  {
1257    dr=deg(r);
1258  }
1259  if (r==0){return(1);}
1260  else{return(0);}
1261}
1262
1263// postredgb: given a minimal basis (gb reducing) it
1264// reduces each polynomial wrt to the others
1265static proc postredgb(ideal F)
1266{
1267  int te=0;
1268  if(defined(@P)==1){te=1;}
1269  ideal G;
1270  ideal H;
1271  int i;
1272  if (size(F)==0){return(ideal(0));}
1273  for (i=1;i<=size(F);i++)
1274  {
1275    H=delfromideal(F,i);
1276    G[i]=pdivi(F[i],H)[1];
1277  }
1278  if(te==1){setglobalrings();}
1279  return(G);
1280}
1281
1282
1283//purpose reduced Groebner basis called in @R and computed in @P
1284static proc gbR(ideal N)
1285{
1286  def RR=basering;
1287  setring(@P);
1288  def Np=imap(RR,N);
1289  option(redSB);
1290  Np=std(Np);
1291  setring(RR);
1292  return(imap(@P,Np));
1293}
1294
1295//**************End homogenizing************************
1296
1297//**************Begin of Groebner Cover*****************
1298
1299// incquotient
1300// incremental quotient
1301// Input: ideal N: a Groebner basis of an ideal
1302//        poly f:
1303// Output: Na = N:<f>
1304static proc incquotient(ideal N, poly f)
1305{
1306  poly g; int i;
1307  ideal Nb; ideal Na=N;
1308  if (size(Na)==1)
1309  {
1310    g=gcd(Na[1],f);
1311    if (g!=1)
1312    {
1313      Na[1]=Na[1]/g;
1314    }
1315    attrib(Na,"IsSB",1);
1316    return(Na);
1317  }
1318  def P=basering;
1319  poly @t;
1320  ring H=0,@t,lp;
1321  def HP=H+P;
1322  setring(HP);
1323  def fh=imap(P,f);
1324  def Nh=imap(P,N);
1325  ideal Nht;
1326  for (i=1;i<=size(Nh);i++)
1327  {
1328    Nht[i]=Nh[i]*@t;
1329  }
1330  attrib(Nht,"isSB",1);
1331  def fht=(1-@t)*fh;
1332  option(redSB);
1333  Nht=std(Nht,fht);
1334  ideal Nc; ideal v;
1335  for (i=1;i<=size(Nht);i++)
1336  {
1337    v=variables(Nht[i]);
1338    if(memberpos(@t,v)[1]==0)
1339    {
1340      Nc[size(Nc)+1]=Nht[i]/fh;
1341    }
1342  }
1343  setring(P);
1344  ideal HH;
1345  def Nd=imap(HP,Nc); Nb=Nd;
1346  option(redSB);
1347  Nb=std(Nd);
1348  return(Nb);
1349}
1350
1351// Auxiliary routine to define an order for ideals
1352// Returns 1 if the ideal a is shoud precede ideal b by sorting them in idbefid order
1353//             2 if the the contrary happen
1354//             0 if the order is not relevant
1355static proc idbefid(ideal a, ideal b)
1356{
1357  poly fa; poly fb; poly la; poly lb;
1358  int te=1; int i; int j;
1359  int na=size(a);
1360  int nb=size(b);
1361  int nm;
1362  if (na<=nb){nm=na;} else{nm=nb;}
1363  for (i=1;i<=nm; i++)
1364  {
1365    fa=a[i]; fb=b[i];
1366    while((fa!=0) or (fb!=0))
1367    {
1370      fa=fa-la;
1371      fb=fb-lb;
1376    }
1377  }
1378  if(na<nb){return(1);}
1379  else
1380  {
1381    if(na>nb){return(2);}
1382    else{return(0);}
1383  }
1384}
1385
1386// sort a list of ideals using idbefid
1387static proc sortlistideals(list L)
1388{
1389  int i; int j; int n;
1390  ideal a; ideal b;
1391  list LL=L;
1392  list NL;
1393  int k; int te;
1394  i=1;
1395  while(size(LL)>0)
1396  {
1397    k=1;
1398    for(j=2;j<=size(LL);j++)
1399    {
1400      te=idbefid(LL[k],LL[j]);
1401      if (te==2){k=j;}
1402    }
1403    NL[size(NL)+1]=LL[k];
1404    n=size(LL);
1405    if (n>1){LL=elimfromlist(LL,k);} else{LL=list();}
1406  }
1407  return(NL);
1408}
1409
1410// Crep
1411// Computes the C-representation of V(N) \ V(M).
1412// input:
1413//    ideal N (null ideal) (not necessarily radical nor maximal)
1414//    ideal M (hole ideal) (not necessarily containing N)
1415// output:
1416//    the list (a,b) of the canonical ideals
1417//    the Crep of V(N) \ V(M)
1418// Assumed to be called in the ring @R
1419// Works on the ring @P
1420proc Crep(ideal N, ideal M)
1421 "USAGE:  Crep(N,M);
1422 Input: ideal N (null ideal) (not necessarily radical nor maximal)
1423           ideal M (hole ideal) (not necessarily containing N)
1424 RETURN:  The canonical C-representation of the locally closed set.
1425           [ P,Q ], a pair of radical ideals with P included in Q,
1426           representing the set V(P) \ V(Q) = V(N) \ V(M)
1427 NOTE:   Operates in a ring R=Q[a] (a=parameters)
1428 KEYWORDS: locally closed set, canoncial form
1429 EXAMPLE:  Crep; shows an example"
1430{
1431  list l;
1432  ideal Np=std(N);
1433  if (equalideals(Np,ideal(1)))
1434  {
1435    l=ideal(1),ideal(1);
1436    return(l);
1437  }
1438  int i;
1439  list L;
1440  ideal Q=Np+M;
1441  ideal P=ideal(1);
1442  L=minGTZ(Np);
1443  for(i=1;i<=size(L);i++)
1444  {
1445    if(idcontains(L[i],Q)==0)
1446    {
1447      P=intersect(P,L[i]);
1448    }
1449  }
1450  P=std(P);
1452  list T=P,Q;
1453  return(T);
1454}
1455example
1456{
1457  "EXAMPLE:"; echo = 2;
1458  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
1459  ring R=0,(x,y,z),lp;
1460  short=0;
1461  ideal E=x*(x^2+y^2+z^2-25);
1462  ideal N=x*(x-3),y-4;
1463  def Cr=Crep(E,N);
1464  Cr;
1465  def L=Prep(E,N);
1466  L;
1467  def Cr1=PtoCrep(L);
1468  Cr1;
1469}
1470
1471// Prep
1472// Computes the P-representation of V(N) \ V(M).
1473// input:
1474//    ideal N (null ideal) (not necessarily radical nor maximal)
1475//    ideal M (hole ideal) (not necessarily containing N)
1476// output:
1477//    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
1478//    the Prep of V(N) \ V(M)
1479// Assumed to be called in the ring @R
1480// Works on the ring @P
1481proc Prep(ideal N, ideal M)
1482 "USAGE:  Prep(N,M);
1483 Input: ideal N (null ideal) (not necessarily radical nor maximal)
1484           ideal M (hole ideal) (not necessarily containing N)
1485 RETURN: The canonical P-representation of the locally closed set V(N) \ V(M)
1486           Output: [ Comp_1, .. , Comp_s ] where
1487              Comp_i=[p_i,[p_i1,..,p_is_i]]
1488 NOTE:   Operates in a ring R=Q[a]  (a=parameters)
1489 KEYWORDS: Locally closed set, Canoncial form
1490 EXAMPLE:  Prep; shows an example"
1491{
1492  int te;
1493  if (N[1]==1)
1494  {
1495    return(list(list(ideal(1),list(ideal(1)))));
1496  }
1497  int i; int j; list L0;
1498  list Ni=minGTZ(N);
1499  list prep;
1500  for(j=1;j<=size(Ni);j++)
1501  {
1502    option(redSB);
1503    Ni[j]=std(Ni[j]);
1504  }
1505  list Mij;
1506  for (i=1;i<=size(Ni);i++)
1507  {
1508    Mij=minGTZ(Ni[i]+M);
1509    for(j=1;j<=size(Mij);j++)
1510    {
1511      option(redSB);
1512      Mij[j]=std(Mij[j]);
1513    }
1514    if ((size(Mij)==1) and (equalideals(Ni[i],Mij[1])==1)){;}
1515    else
1516    {
1517        prep[size(prep)+1]=list(Ni[i],Mij);
1518    }
1519  }
1520  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
1521  def Lout=CompleteA(prep,prep);
1522  return(Lout);
1523}
1524example
1525{
1526  "EXAMPLE:"; echo = 2;
1527  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
1528  short=0;
1529  ring R=0,(x,y,z),lp;
1530  ideal E=x*(x^2+y^2+z^2-25);
1531  ideal N=x*(x-3),y-4;
1532  Prep(E,N);
1533}
1534
1535// PtoCrep
1536// Computes the C-representation from the P-representation.
1537// input:
1538//    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
1539//         the P-representation of V(N) \ V(M)
1540// output:
1541//    list (ideal ida, ideal idb)
1542//    the C-representaion of V(N) \ V(M) = V(ida) \ V(idb)
1543// Assumed to be called in the ring @R
1544// Works on the ring @P
1545proc PtoCrep(list L)
1546"USAGE: PtoCrep(L)
1547 Input L: list  [ Comp_1, .. , Comp_s ] where
1548          Comp_i=[p_i,[p_i1,..,p_is_i] ], is
1549          the P-representation of a locally closed set V(N) \ V(M)
1550 RETURN:  The canonical C-representation of the locally closed set
1551          [ P,Q ], a pair of radical ideals with P included in Q,
1552          representing the set V(P) \ V(Q)
1553 NOTE: Operates in a ring R=Q[a]  (a=parameters)
1554 KEYWORDS: locally closed set, canoncial form
1555 EXAMPLE:  PtoCrep; shows an example"
1556{
1557  int te;
1558  def RR=basering;
1559  if(defined(@P)){te=1; setring(@P); list Lp=imap(RR,L);}
1560  else {  te=0; def Lp=L;}
1561  def La=PtoCrep0(Lp);
1562  if(te==1) {setring(RR); def LL=imap(@P,La);}
1563  if(te==0){def LL=La;}
1564  return(LL);
1565}
1566example
1567{
1568  "EXAMPLE:"; echo = 2;
1569  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
1570  short=0;
1571  ring R=0,(x,y,z),lp;
1572
1573  // (P,Q) represents a locally closed set
1574  ideal P=x^3+x*y^2+x*z^2-25*x;
1575  ideal Q=y-4,x*z,x^2-3*x;
1576
1577  // Now compute the P-representation=
1578  def L=Prep(P,Q);
1579  L;
1580  // Now compute the C-representation=
1581  def J=PtoCrep(L);
1582  J;
1583  // Now come back recomputing the P-represetation of the C-representation=
1584  Prep(J[1],J[2]);
1585}
1586
1587// PtoCrep0
1588// Computes the C-representation from the P-representation.
1589// input:
1590//    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
1591//         the P-representation of V(N) \ V(M)
1592// output:
1593//    list (ideal ida, ideal idb)
1594//    the C-representation of V(N) \ V(M) = V(ida) \ V(idb)
1595// Works in a ring Q[u_1,..,u_m] and is called on it.
1596static proc PtoCrep0(list L)
1597{
1598  int te=0;
1599  def Lp=L;
1600  int i; int j;
1601  ideal ida=ideal(1); ideal idb=ideal(1); list Lb; ideal N;
1602  for (i=1;i<=size(Lp);i++)
1603  {
1604    option(returnSB);
1605    N=Lp[i][1];
1606    ida=intersect(ida,N);
1607    Lb=Lp[i][2];
1608    for(j=1;j<=size(Lb);j++)
1609    {
1610      idb=intersect(idb,Lb[j]);
1611    }
1612  }
1614  def La=list(ida,idb);
1615  return(La);
1616}
1617
1618// input: F a parametric ideal in Q[a][x]
1619// output: a rComprehensive Groebner System disjoint and reduced.
1620//      It uses Kapur-Sun-Wang algorithm, and with the options
1621//      can compute the homogenization before  (('can',0) or ( 'can',1))
1622//      and dehomogenize the result.
1623proc cgsdr(ideal F, list #)
1624"USAGE: cgsdr(F);
1625          F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed.
1626          To compute a disjoint, reduced Comprehensive Groebner System (CGS).
1627          cgsdr is the starting point of the fundamental routine grobcov.
1628          Inside grobcov it is used with options 'can' set to 0,1 and
1629          not with options ('can',2).
1630          It is to be used if only a disjoint reduced CGS is required.
1631
1632          Options: To modify the default options, pairs of arguments
1633          -option name, value- of valid options must be added to the call.
1634
1635          Options:
1636            \"can\",0-1-2: The default value is \"can\",2. In this case no
1637                homogenization is done. With option (\"can\",0) the given
1638                basis is homogenized, and with option (\"can\",1) the
1639                whole given ideal is homogenized before computing the
1640                cgs and dehomogenized after.
1641                with option (\"can\",0) the homogenized basis is used
1642                with option (\"can\",1) the homogenized ideal is used
1643                with option (\"can\",2) the given basis is used
1644            \"null\",ideal E: The default is (\"null\",ideal(0)).
1645            \"nonnull\",ideal N: The default (\"nonnull\",ideal(1)).
1646                When options \"null\" and/or \"nonnull\" are given, then
1647                the parameter space is restricted to V(E) \ V(N).
1648            \"comment\",0-1: The default is (\"comment\",0). Setting (\"comment\",1)
1649                will provide information about the development of the
1650                computation.
1651            \"out\",0-1: 1 (default) the output segments are given as
1652                as difference of varieties.
1653                0: the output segments are given in P-representation
1654                and the segments grouped by lpp
1655                With options (\"can\",0) and (\"can\",1) the option (\"out\",1)
1656                is set to (\"out\",0) because it is not compatible.
1657          One can give none or whatever of these options.
1658          With the default options (\"can\",2,\"out\",1), only the
1659          Kapur-Sun-Wang algorithm is computed. This is very efficient
1660          but is only the starting point for the computation of grobcov.
1661          When grobcov is computed, the call to cgsdr inside uses
1662          specific options that are more expensive ("can",0-1,"out",0).
1663RETURN: Returns a list T describing a reduced and disjoint
1664          Comprehensive Groebner System (CGS),
1665          With option (\"out\",0)
1666          the segments are grouped by
1667          leading power products (lpp) of the reduced Groebner
1668          basis and given in P-representation.
1669          The returned list is of the form:
1670           (
1671             (lpp, (num,basis,segment),...,(num,basis,segment),lpp),
1672             ..,,
1673             (lpp, (num,basis,segment),...,(num,basis,segment),lpp)
1674           )
1675          The bases are the reduced Groebner bases (after normalization)
1676          for each point of the corresponding segment.
1677
1678          The third element of each lpp segment is the lpp of the
1679          used ideal in the CGS as a string:
1680            with option (\"can\",0) the homogenized basis is used
1681            with option (\"can\",1) the homogenized ideal is used
1682            with option (\"can\",2) the given basis is used
1683
1684          With option (\"out\",1) (default)
1685          only KSW is applied and segments are given as
1686          difference of varieties and are not grouped
1687          The returned list is of the form:
1688           (
1689             (E,N,B),..(E,N,B)
1690           )
1691          E is the null variety
1692          N is the nonnull variety
1693          segment = V(E) \ V(N)
1694          B is the reduced Groebner basis
1695
1696NOTE:  The basering R, must be of the form Q[a][x], (a=parameters,
1697          x=variables), and should be defined previously, and the ideal
1698          defined on R.
1699KEYWORDS: CGS, disjoint, reduced, Comprehensive Groebner System
1700EXAMPLE:  cgsdr; shows an example"
1701{
1702  int te;
1703  def RR=basering;
1704  if(defined(@P)){te=1;}
1705  else{te=0; setglobalrings();}
1706  // INITIALIZING OPTIONS
1707  int i; int j;
1708  def E=ideal(0);
1709  def N=ideal(1);
1710  int comment=0;
1711  int can=2;
1712  int out=1;
1713  poly f;
1714  ideal B;
1715  int start=timer;
1716  list L=#;
1717  for(i=1;i<=size(L) div 2;i++)
1718  {
1719    if(L[2*i-1]=="null"){E=L[2*i];}
1720    else
1721    {
1722      if(L[2*i-1]=="nonnull"){N=L[2*i];}
1723      else
1724      {
1725        if(L[2*i-1]=="comment"){comment=L[2*i];}
1726        else
1727        {
1728          if(L[2*i-1]=="can"){can=L[2*i];}
1729          else
1730          {
1731            if(L[2*i-1]=="out"){out=L[2*i];}
1732          }
1733        }
1734      }
1735    }
1736  }
1737  //if(can==2){out=1;}
1738  B=F;
1739  if ((printlevel) and (comment==0)){comment=printlevel;}
1740  if((can<2) and (out>0)){"Option out,1 is not compatible with can,0,1"; out=0;}
1741  // DEFINING OPTIONS
1742  list LL;
1743  LL[1]="can";     LL[2]=can;
1744  LL[3]="comment"; LL[4]=comment;
1745  LL[5]="out";     LL[6]=out;
1746  LL[7]="null";    LL[8]=E;
1747  LL[9]="nonnull"; LL[10]=N;
1748  if(comment>=1)
1749  {
1750    string("Begin cgsdr with options: ",LL);
1751  }
1752  int ish;
1753  for (i=1;i<=size(B);i++){ish=ishomog(B[i]); if(ish==0){break;};}
1754  if (ish)
1755  {
1756    if(comment>0){string("The given system is homogneous");}
1757    def GS=KSW(B,LL);
1758    //can=0;
1759  }
1760  else
1761  {
1762  // ACTING DEPENDING ON OPTIONS
1763  if(can==2)
1764  {
1765    // WITHOUT HOMOHGENIZING
1766    if(comment>0){string("Option of cgsdr: do not homogenize");}
1767    def GS=KSW(B,LL);
1768    setglobalrings();
1769  }
1770  else
1771  {
1772    if(can==1)
1773    {
1774      // COMPUTING THE HOMOGOENIZED IDEAL
1775      if(comment>=1){string("Homogenizing the whole ideal: option can=1");}
1776      list RRL=ringlist(RR);
1777      RRL[3][1][1]="dp";
1778      def Pa=ring(RRL[1]);
1779      list Lx;
1780      Lx[1]=0;
1781      Lx[2]=RRL[2]+RRL[1][2];
1782      Lx[3]=RRL[1][3];
1783      Lx[4]=RRL[1][4];
1784      RRL[1]=0;
1785      def D=ring(RRL);
1786      def RP=D+Pa;
1787      setring(RP);
1788      def B1=imap(RR,B);
1789      option(redSB);
1790      B1=std(B1);
1791      setring(RR);
1792      def B2=imap(RP,B1);
1793    }
1794    else
1795    { // (can=0)
1796       if(comment>0){string("Homogenizing the basis: option can=0");}
1797      def B2=B;
1798    }
1799    // COMPUTING HOMOGENIZED CGS
1800    poly @t;
1801    ring H=0,@t,dp;
1802    def RH=RR+H;
1803    setring(RH);
1804    setglobalrings();
1805    def BH=imap(RR,B2);
1806    def LH=imap(RR,LL);
1807    for (i=1;i<=size(BH);i++)
1808    {
1809      BH[i]=homog(BH[i],@t);
1810    }
1811    if (comment>=2){string("Homogenized system = "); BH;}
1812    def GSH=KSW(BH,LH);
1813    setglobalrings();
1814    // DEHOMOGENIZING THE RESULT
1815    if(out==0)
1816    {
1817      for (i=1;i<=size(GSH);i++)
1818      {
1819        GSH[i][1]=subst(GSH[i][1],@t,1);
1820        for(j=1;j<=size(GSH[i][2]);j++)
1821        {
1822          GSH[i][2][j][2]=subst(GSH[i][2][j][2],@t,1);
1823        }
1824      }
1825    }
1826    else
1827    {
1828      for (i=1;i<=size(GSH);i++)
1829      {
1830        GSH[i][3]=subst(GSH[i][3],@t,1);
1831        GSH[i][7]=subst(GSH[i][7],@t,1);
1832      }
1833    }
1834    setring(RR);
1835    def GS=imap(RH,GSH);
1836    }
1837
1838
1839    setglobalrings();
1840    if(out==0)
1841    {
1842      for (i=1;i<=size(GS);i++)
1843      {
1844        GS[i][1]=postredgb(mingb(GS[i][1]));
1845        for(j=1;j<=size(GS[i][2]);j++)
1846        {
1847          GS[i][2][j][2]=postredgb(mingb(GS[i][2][j][2]));
1848        }
1849      }
1850    }
1851    else
1852    {
1853      for (i=1;i<=size(GS);i++)
1854      {
1855        if(GS[i][2]==1)
1856        {
1857          GS[i][3]=postredgb(mingb(GS[i][3]));
1858          if (typeof(GS[i][7])=="ideal")
1859          { GS[i][7]=postredgb(mingb(GS[i][7]));}
1860        }
1861      }
1862    }
1863  }
1864  if(te==0){kill @P; kill @R; kill @RP;}
1865  return(GS);
1866}
1867example
1868{
1869  "EXAMPLE:"; echo = 2;
1870  // Casas conjecture for degree 4:
1871  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
1872  short=0;
1873  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
1874          x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
1875          x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
1876          x2^2+(2*a3)*x2+(a2),
1877          x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
1878          x3+(a3);
1879  cgsdr(F);
1880}
1881
1882// input:  internal routine called by cgsdr at the end to group the
1883//            lpp segments and improve the output
1884// output: grouped segments by lpp obtained in cgsdr
1885static proc grsegments(list T)
1886{
1887  int i;
1888  list L;
1889  list lpp;
1890  list lp;
1891  list ls;
1892  int n=size(T);
1893  lpp[1]=T[n][1];
1894  L[1]=list(lpp[1],list(list(T[n][2],T[n][3],T[n][4])));
1895  if (n>1)
1896  {
1897    for (i=1;i<=size(T)-1;i++)
1898    {
1899      lp=memberpos(T[n-i][1],lpp);
1900      if(lp[1]==1)
1901      {
1902        ls=L[lp[2]][2];
1903        ls[size(ls)+1]=list(T[n-i][2],T[n-i][3],T[n-i][4]);
1904        L[lp[2]][2]=ls;
1905      }
1906      else
1907      {
1908        lpp[size(lpp)+1]=T[n-i][1];
1909        L[size(L)+1]=list(T[n-i][1],list(list(T[n-i][2],T[n-i][3],T[n-i][4])));
1910      }
1911    }
1912  }
1913  return(L);
1914}
1915
1916// LCUnion
1917// Given a list of the P-representations of locally closed segments
1918// for which we know that the union is also locally closed
1919// it returns the P-representation of its union
1920// input:  L list of segments in P-representation
1921//      ((p_j^i,(p_j1^i,...,p_jk_j^i | j=1..t_i)) | i=1..s )
1922//      where i represents a segment
1923// output: P-representation of the union
1924//       ((P_j,(P_j1,...,P_jk_j | j=1..t)))
1925static proc LCUnion(list LL)
1926{
1927  def RR=basering;
1928  setring(@P);
1929  def L=imap(RR,LL);
1930  int i; int j; int k; list H; list C; list T;
1931  list L0; list P0; list P; list Q0; list Q;
1932  for (i=1;i<=size(L);i++)
1933  {
1934    for (j=1;j<=size(L[i]);j++)
1935    {
1936      P0[size(P0)+1]=L[i][j][1];
1937      L0[size(L0)+1]=intvec(i,j);
1938    }
1939  }
1940  Q0=selectminideals(P0);
1941  for (i=1;i<=size(Q0);i++)
1942  {
1943    Q[i]=L0[Q0[i]];
1944    P[i]=L[Q[i][1]][Q[i][2]];
1945  }
1946  // P is the list of the maximal components of the union
1947  //   with the corresponding initial holes.
1948  // Q is the list of intvec positions in L of the first element of the P's
1949  //   Its elements give (num of segment, num of max component (=min ideal))
1950  for (k=1;k<=size(Q);k++)
1951  {
1952    H=P[k][2]; // holes of P[k][1]
1953    for (i=1;i<=size(L);i++)
1954    {
1955      if (i!=Q[k][1])
1956      {
1957        for (j=1;j<=size(L[i]);j++)
1958        {
1959          C[size(C)+1]=L[i][j];
1960        }
1961      }
1962    }
1964  }
1965  setring(RR);
1966  def TT=imap(@P,T);
1967  return(TT);
1968}
1969
1970// Auxiliary routine
1971// called by LCUnion to modify the holes of a primepart of the union
1972// by the addition of the segments that do not correspond to that part
1973// Works on @P ring.
1974// Input:
1975//   H=(p_i1,..,p_is) the holes of a component to be transformed by the addition of
1976//        the segments C that do not correspond to that component
1977//   C=((q_1,(q_11,..,q_1l_1),pos1),..,(q_k,(q_k1,..,q_kl_k),posk))
1978// posi=(i,j) position of the component
1979//        the list of segments to be added to the holes
1980static proc addpart(list H, list C)
1981{
1982  list Q; int i; int j; int k; int l; int t; int t1;
1983  Q=H; intvec notQ; list QQ; list addq;
1984  // @Q2=list of (i,j) positions of the components that have been aded to some hole of the maximal ideals
1985  //          plus those of the components added to the holes.
1986  ideal q;
1987  i=1;
1988  while (i<=size(Q))
1989  {
1990    if (memberpos(i,notQ)[1]==0)
1991    {
1992      q=Q[i];
1993      t=1; j=1;
1994      while ((t) and (j<=size(C)))
1995      {
1996        if (equalideals(q,C[j][1]))
1997        {
1998          // \\ @Q2[size(@Q2)+1]=C[j][3];
1999          t=0;
2000          for (k=1;k<=size(C[j][2]);k++)
2001          {
2002            t1=1;
2003            l=1;
2004            while((t1) and (l<=size(Q)))
2005            {
2006              if ((l!=i) and (memberpos(l,notQ)[1]==0))
2007              {
2008                if (idcontains(C[j][2][k],Q[l]))
2009                {
2010                  t1=0;
2011                }
2012              }
2013              l++;
2014            }
2015            if (t1)
2016            {
2018              // \\ @Q2[size(@Q2)+1]=C[j][3];
2019            }
2020          }
2021          if((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;}
2022          else {notQ[size(notQ)+1]=i;}
2023        }
2024        j++;
2025      }
2027      {
2029        {
2031        }
2034      }
2035    }
2036    i++;
2037  }
2038  for (i=1;i<=size(Q);i++)
2039  {
2040    if(memberpos(i,notQ)[1]==0)
2041    {
2042      QQ[size(QQ)+1]=Q[i];
2043    }
2044  }
2045  if (size(QQ)==0){QQ[1]=ideal(1);}
2047}
2048
2049// Auxiliary routine called by addpart to finish the modification of the holes of a primepart
2050// of the union by the addition of the segments that do not correspond to
2051// that part.
2052// Works on @P ring.
2053static proc addpartfine(list H, list C0)
2054{
2055  //"T_H="; H;
2056  int i; int j; int k; int te; intvec notQ; int l; list sel; int used;
2057  intvec jtesC;
2058  if ((size(H)==1) and (equalideals(H[1],ideal(1)))){return(H);}
2059  if (size(C0)==0){return(H);}
2060  list newQ; list nQ; list Q; list nQ1; list Q0;
2061  def Q1=H;
2062  //Q1=sortlistideals(Q1,idbefid);
2063  def C=C0;
2064  while(equallistideals(Q0,Q1)==0)
2065  {
2066    Q0=Q1;
2067    i=0;
2068    Q=Q1;
2069    kill notQ; intvec notQ;
2070    while(i<size(Q))
2071    {
2072      i++;
2073      for(j=1;j<=size(C);j++)
2074      {
2075        te=idcontains(Q[i],C[j][1]);
2076        if(te)
2077        {
2078          for(k=1;k<=size(C[j][2]);k++)
2079          {
2080            if(idcontains(Q[i],C[j][2][k]))
2081            {
2082              te=0; break;
2083            }
2084          }
2085          if (te)
2086          {
2087            used++;
2088            if ((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;}
2089            else{notQ[size(notQ)+1]=i;}
2090            kill newQ; list newQ;
2091            for(k=1;k<=size(C[j][2]);k++)
2092            {
2093              nQ=minGTZ(Q[i]+C[j][2][k]);
2094              for(l=1;l<=size(nQ);l++)
2095              {
2096                option(redSB);
2097                nQ[l]=std(nQ[l]);
2098                newQ[size(newQ)+1]=nQ[l];
2099              }
2100            }
2101            sel=selectminideals(newQ);
2102            kill nQ1; list nQ1;
2103            for(l=1;l<=size(sel);l++)
2104            {
2105              nQ1[l]=newQ[sel[l]];
2106            }
2107            newQ=nQ1;
2108            for(l=1;l<=size(newQ);l++)
2109            {
2110              Q[size(Q)+1]=newQ[l];
2111            }
2112            break;
2113          }
2114        }
2115      }
2116    }
2117    kill Q1; list Q1;
2118    for(i=1;i<=size(Q);i++)
2119    {
2120      if(memberpos(i,notQ)[1]==0)
2121      {
2122        Q1[size(Q1)+1]=Q[i];
2123      }
2124    }
2125    //"T_Q1_0="; Q1;
2126    sel=selectminideals(Q1);
2127    kill nQ1; list nQ1;
2128    for(l=1;l<=size(sel);l++)
2129    {
2130      nQ1[l]=Q1[sel[l]];
2131    }
2132    Q1=nQ1;
2133  }
2134  if(size(Q1)==0){Q1=ideal(1),ideal(1);}
2135  //"T_Q1_1="; Q1;
2136  //if(used>0){string("addpartfine was ", used, " times used");}
2137  return(Q1);
2138}
2139
2140
2141// Auxiliary routine for grobcov: ideal F is assumed to be homogeneous
2142// gcover
2143// input: ideal F: a generating set of a homogeneous ideal in Q[a][x]
2144//    list #: optional
2145// output: the list
2146//   S=((lpp, generic basis, Prep, Crep),..,(lpp, generic basis, Prep, Crep))
2147//      where a Prep is ( (p1,(p11,..,p1k_1)),..,(pj,(pj1,..,p1k_j)) )
2148//            a Crep is ( ida, idb )
2149static proc gcover(ideal F,list #)
2150{
2151  int i; int j; int k; ideal lpp; list GPi2; list pairspP; ideal B; int ti;
2152  int i1; int tes; int j1; int selind; int i2; int m;
2153  list prep; list crep; list LCU; poly p; poly lcp; ideal FF;
2154  list lpi;
2155  string lpph;
2156  list L=#;
2157  int canop=1;
2158  int extop=1;
2159  int repop=0;
2160  ideal E=ideal(0);;
2161  ideal N=ideal(1);;
2162  int comment;
2163  for(i=1;i<=size(L) div 2;i++)
2164  {
2165    if(L[2*i-1]=="can"){canop=L[2*i];}
2166    else
2167    {
2168      if(L[2*i-1]=="ext"){extop=L[2*i];}
2169      else
2170      {
2171        if(L[2*i-1]=="rep"){repop=L[2*i];}
2172        else
2173        {
2174          if(L[2*i-1]=="null"){E=L[2*i];}
2175          else
2176          {
2177            if(L[2*i-1]=="nonnull"){N=L[2*i];}
2178            else
2179            {
2180              if (L[2*i-1]=="comment"){comment=L[2*i];}
2181            }
2182          }
2183        }
2184      }
2185    }
2186  }
2187  list GS; list GP;
2188  def RR=basering;
2189  GS=cgsdr(F,L); // "null",NW[1],"nonnull",NW[2],"cgs",CGS,"comment",comment);
2190  setglobalrings();
2191  int start=timer;
2192  GP=GS;
2193  ideal lppr;
2194  list LL;
2195  list S;
2196  poly sp;
2197  ideal BB;
2198  for (i=1;i<=size(GP);i++)
2199  {
2200    kill LL;
2201    list LL;
2202    lpp=GP[i][1];
2203    GPi2=GP[i][2];
2204    lpph=GP[i][3];
2205    kill pairspP; list pairspP;
2206    for(j=1;j<=size(GPi2);j++)
2207    {
2208      pairspP[size(pairspP)+1]=GPi2[j][3];
2209    }
2210    LCU=LCUnion(pairspP);
2211    kill prep; list prep;
2212    kill crep; list crep;
2213    for(k=1;k<=size(LCU);k++)
2214    {
2215      prep[k]=list(LCU[k][2],LCU[k][3]);
2216      B=GPi2[LCU[k][1][1]][2]; // ATENTION last 1 has been changed to [2]
2217      LCU[k][1]=B;
2218    }
2219    //"Deciding if combine is needed";
2220    kill BB;
2221    ideal BB;
2222    tes=1; m=1;
2223    while((tes) and (m<=size(LCU[1][1])))
2224    {
2225      j=1;
2226      while((tes) and (j<=size(LCU)))
2227      {
2228        k=1;
2229        while((tes) and (k<=size(LCU)))
2230        {
2231          if(j!=k)
2232          {
2233            sp=pnormalf(pspol(LCU[j][1][m],LCU[k][1][m]),LCU[k][2],N);
2234            if(sp!=0){tes=0;}
2235          }
2236          k++;
2237        }        //setglobalrings();
2238        if(tes)
2239        {
2240          BB[m]=LCU[j][1][m];
2241        }
2242        j++;
2243      }
2244      if(tes==0){break;}
2245      m++;
2246    }
2247    crep=PtoCrep(prep);
2248    if(tes==0)
2249    {
2250      // combine is needed
2251      kill B; ideal B;
2252      for (j=1;j<=size(LCU);j++)
2253      {
2254        LL[j]=LCU[j][2];
2255      }
2256      if (size(LCU)>1)
2257      {
2258        FF=precombint(LL);
2259      }
2260      for (k=1;k<=size(lpp);k++)
2261      {
2262        kill L; list L;
2263        for (j=1;j<=size(LCU);j++)
2264        {
2265          L[j]=list(LCU[j][2],LCU[j][1][k]);
2266        }
2267        if (size(LCU)>1)
2268        {
2269          B[k]=combine(L,FF);
2270        }
2271        else{B[k]=L[1][2];}
2272      }
2273    }
2274    else{B=BB;}
2275    for(j=1;j<=size(B);j++)
2276    {
2277      B[j]=pnormalf(B[j],crep[1],crep[2]);
2278    }
2279    S[i]=list(lpp,B,prep,crep,lpph);
2280    if(comment>=1)
2281    {
2282      lpi[size(lpi)+1]=string("[",i,"]");
2283      lpi[size(lpi)+1]=S[i][1];
2284    }
2285  }
2286  if(comment>=1)
2287  {
2288    string("Time in LCUnion + combine = ",timer-start);
2289    if(comment>=2){string("lpp=",lpi)};
2290  }
2291  if(defined(@P)==1){kill @P; kill @RP; kill @R;}
2292  return(S);
2293}
2294
2295// grobcov
2296// input:
2297//    ideal F: a parametric ideal in Q[a][x], (a=parameters, x=variables).
2298//    list #: (options) list("null",N,"nonnull",W,"can",0-1,ext",0-1, "rep",0-1-2)
2299//            where
2300//            N is the null conditions ideal (if desired)
2301//            W is the ideal of non-null conditions (if desired)
2302//            The value of \"can\"i s 1 by default and can be set to 0 if we do not
2303//            need to obtain the canonical GC, but only a GC.
2304//            The value of \"ext\" is 0 by default and so the generic representation
2305//             of the bases is given. It can be set to 1, and then the full
2306//             representation of the bases is given.
2307//            The value of \"rep\" is 0 by default, and then the segments
2308//            are given in canonical P-representation. It can be set to 1
2309//            and then they are given in canonical C-representation.
2310//            If it is set to 2, then both representations are given.
2311// output:
2312//    list S: ((lpp,basis,(idp_1,(idp_11,..,idp_1s_1))), ..
2313//             (lpp,basis,(idp_r,(idp_r1,..,idp_rs_r))) ) where
2314//            each element of S corresponds to a lpp-segment
2315//            given by the lpp, the basis, and the P-representation of the segment
2316proc grobcov(ideal F,list #)
2317"USAGE:   grobcov(F); This is the fundamental routine of the
2318          library. It computes the Groebner cover of a parametric ideal. See
2319                Montes A., Wibmer M.,
2320               \"Groebner Bases for Polynomial Systems with parameters\".
2321               JSC 45 (2010) 1391-1425.)
2322          The Groebner Cover of a parametric ideal consist of a set of
2323          pairs(S_i,B_i), where the S_i are disjoint locally closed
2324          segments of the parameter space, and the B_i are the reduced
2325          Groebner bases of the ideal on every point of S_i.
2326
2327          The ideal F must be defined on a parametric ring Q[a][x].
2328          (a=parameters, x=variables)
2329          Options: To modify the default options, pair of arguments
2330          -option name, value- of valid options must be added to the call.
2331
2332          Options:
2333            \"null\",ideal E: The default is (\"null\",ideal(0)).
2334            \"nonnull\",ideal N: The default is (\"nonnull\",ideal(1)).
2335                When options \"null\" and/or \"nonnull\" are given, then
2336                the parameter space is restricted to V(E) \ V(N).
2337            \"can\",0-1: The default is (\"can\",1). With the default option
2338                the homogenized ideal is computed before obtaining the
2339                Groebner Cover, so that the result is the canonical
2340                Groebner Cover. Setting (\"can\",0) only homogenizes the
2341                basis so the result is not exactly canonical, but the
2342                computation is shorter.
2343            \"ext\",0-1: The default is (\"ext\",0). With the default
2344                (\"ext\",0), only the generic representation of the bases is
2345                computed (single polynomials, but not specializing to non-zero
2346                for every point of the segment. With option (\"ext\",1) the
2347                full representation of the bases is computed (possible
2348                sheaves) and sometimes a simpler result is obtained,
2349                but the computation is more time consuming.
2350            \"rep\",0-1-2: The default is (\"rep\",0) and then the segments
2351                are given in canonical P-representation. Option (\"rep\",1)
2352                represents the segments in canonical C-representation,
2353                and option (\"rep\",2) gives both representations.
2354            \"comment\",0-3: The default is (\"comment\",0). Setting
2355                \"comment\" higher will provide information about the
2356                development of the computation.
2357          One can give none or whatever of these options.
2358RETURN:   The list
2359          (
2360           (lpp_1,basis_1,segment_1,lpph_1),
2361           ...
2362           (lpp_s,basis_s,segment_s,lpph_s)
2363          )
2364
2365          The lpp are constant over a segment and correspond to the
2366          set of lpp of the reduced Groebner basis for each point
2367          of the segment.
2368          The lpph corresponds to the lpp of the homogenized ideal
2369          and is different for each segment. It is given as a string,
2370          and shown only for information. With the default option
2371           \"can\",1, the segments have different lpph.
2372
2373          Basis: to each element of lpp corresponds an I-regular function
2374          given in full representation (by option (\"ext\",1)) or in
2375          generic representation (default option (\"ext\",0)). The
2376          I-regular function is the corresponding element of the reduced
2377          Groebner basis for each point of the segment with the given lpp.
2378          For each point in the segment, the polynomial or the set of
2379          polynomials representing it, if they do not specialize to 0,
2380          then after normalization, specializes to the corresponding
2381          element of the reduced Groebner basis. In the full representation
2382          at least one of the polynomials representing the I-regular
2383          function specializes to non-zero.
2384
2385          With the default option (\"rep\",0) the representation of the
2386          segment is the P-representation.
2387          With option (\"rep\",1) the representation of the segment is
2388          the C-representation.
2389          With option (\"rep\",2) both representations of the segment are
2390          given.
2391
2392          The P-representation of a segment is of the form
2393          ((p_1,(p_11,..,p_1k1)),..,(p_r,(p_r1,..,p_rkr))
2394          representing the segment Union_i ( V(p_i) \ ( Union_j V(p_ij) ) ),
2395          where the p's are prime ideals.
2396
2397          The C-representation of a segment is of the form
2398          (E,N) representing V(E) \ V(N), and the ideals E and N are
2399          radical and N contains E.
2400
2401NOTE: The basering R, must be of the form Q[a][x], (a=parameters,
2402          x=variables), and should be defined previously. The ideal must
2403          be defined on R.
2404KEYWORDS: Groebner cover, parametric ideal, canonical, discussion of
2405          parametric ideal.
2406EXAMPLE:  grobcov; shows an example"
2407{
2408  list S; int i; int ish=1; list GBR; list BR; int j; int k;
2409  ideal idp; ideal idq; int s; ideal ext; list SS;
2410  ideal E; ideal N; int canop;  int extop; int repop;
2411  int comment=0; int m;
2412  def RR=basering;
2413  setglobalrings();
2414  list L0=#;
2415  int out=0;
2416  L0[size(L0)+1]="res"; L0[size(L0)+1]=ideal(1);
2417  // default options
2418  int start=timer;
2419  E=ideal(0);
2420  N=ideal(1);
2421  canop=1; // canop=0 for homogenizing the basis but not the ideal (not canonical)
2422           // canop=1 for working with the homogenized ideal
2423  repop=0; // repop=0 for representing the segments in Prep
2424           // repop=1 for representing the segments in Crep
2425           // repop=2 for representing the segments in Prep and Crep
2426  extop=0; // extop=0 if only generic representation of the bases are to be computed
2427           // extop=1 if the full representation of the bases are to be computed
2428  for(i=1;i<=size(L0) div 2;i++)
2429  {
2430    if(L0[2*i-1]=="can"){canop=L0[2*i];}
2431    else
2432    {
2433      if(L0[2*i-1]=="ext"){extop=L0[2*i];}
2434      else
2435      {
2436        if(L0[2*i-1]=="rep"){repop=L0[2*i];}
2437        else
2438        {
2439          if(L0[2*i-1]=="null"){E=L0[2*i];}
2440          else
2441          {
2442            if(L0[2*i-1]=="nonnull"){N=L0[2*i];}
2443            else
2444            {
2445              if (L0[2*i-1]=="comment"){comment=L0[2*i];}
2446            }
2447          }
2448        }
2449      }
2450    }
2451  }
2452  if(not((canop==0) or (canop==1)))
2453  {
2454    string("Option can = ",canop," is not supported. It is changed to can = 1");
2455    canop=1;
2456  }
2457  for(i=1;i<=size(L0) div 2;i++)
2458  {
2459    if(L0[2*i-1]=="can"){L0[2*i]=canop;}
2460  }
2461  if ((printlevel) and (comment==0)){comment=printlevel;}
2462  list LL;
2463  LL[1]="can";     LL[2]=canop;
2464  LL[3]="comment"; LL[4]=comment;
2465  LL[5]="out";     LL[6]=0;
2466  LL[7]="null";    LL[8]=E;
2467  LL[9]="nonnull"; LL[10]=N;
2468  LL[11]="ext";    LL[12]=extop;
2469  LL[13]="rep";    LL[14]=repop;
2470  if (comment>=1)
2471  {
2472    string("Begin grobcov with options: ",LL);
2473  }
2474  kill S;
2475  def S=gcover(F,LL);
2476  // NOW extend
2477  if(extop)
2478  {
2479    S=extend(S,LL);
2480  }
2481  else
2482  {
2483    // NOW representation of the segments by option repop
2484    list Si; list nS;
2485    if(repop==0)
2486    {
2487      for(i=1;i<=size(S);i++)
2488      {
2489        Si=list(S[i][1],S[i][2],S[i][3],S[i][5]);
2490        nS[size(nS)+1]=Si;
2491      }
2492      kill S;
2493      def S=nS;
2494    }
2495    else
2496    {
2497      if(repop==1)
2498      {
2499        for(i=1;i<=size(S);i++)
2500        {
2501          Si=list(S[i][1],S[i][2],S[i][4],S[i][5]);
2502          nS[size(nS)+1]=Si;
2503        }
2504        kill S;
2505        def S=nS;
2506      }
2507      else
2508      {
2509        for(i=1;i<=size(S);i++)
2510        {
2511          Si=list(S[i][1],S[i][2],S[i][3],S[i][4],S[i][5]);
2512          nS[size(nS)+1]=Si;
2513        }
2514        kill S;
2515        def S=nS;
2516      }
2517    }
2518  }
2519  if (comment>=1)
2520  {
2521    string("Time in grobcov = ", timer-start);
2522    string("Number of segments of grobcov = ", size(S));
2523  }
2524  if(defined(@P)==1){kill @R; kill @P; kill @RP;}
2525  return(S);
2526}
2527example
2528{
2529  "EXAMPLE:"; echo = 2;
2530  // Casas conjecture for degree 4:
2531  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
2532  short=0;
2533  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
2534            x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
2535            x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
2536            x2^2+(2*a3)*x2+(a2),
2537            x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
2538            x3+(a3);
2539  grobcov(F);
2540}
2541
2542// Input. GC the grobcov of an ideal in generic representation of the
2543//        bases computed with option option ("rep",2).
2544// Output The grobcov in full representation.
2545// Option ("comment",1) shows the time.
2546// Can be called from the top
2547proc extend(list GC, list #);
2548"USAGE: extend(GC); The default option of grobcov provides
2549          the bases in generic representation (the I-regular functions
2550          of the bases ara given by a single polynomial. It can specialize
2551          to zero for some points of the segments, but in general, it
2552          is sufficient for many pouposes. Nevertheless the I-regular
2553          functions allow a full representation given bey a set of
2554          polynomials specializing to the value of the function (after normalization)
2555          or to zero, but at least one of the polynomials specializes to non-zero.
2556          The full representation can be obtained by computing the
2557          grobcov with option \"ext\",1. The default option is \"ext\",0.
2558          With option \"ext\",1 the computation can be much more
2559          time consuming, even if the result can be simpler.
2560          Alternatively, one can compute the full representation of the
2561          bases after computing grobcov with the defaoult option \"ext\",0
2562          and the option \"rep\",2, that outputs both the Prep and the Crep
2563          of the segments and then call \"extend\" to the output.
2564
2565RETURN: When calling extend(grobcov(S,\"rep\",2)) the result is of the form
2566              (
2567                 (lpp_1,basis_1,segment_1,lpph_1),
2568                  ...
2569                 (lpp_s,basis_s,segment_s,lpph_s)
2570              )
2571          where each function of the basis can be given by an ideal
2572          of representants.
2573
2574NOTE: The basering R, must be of the form Q[a][x], (a=parameters,
2575          x=variables), and should be defined previously. The ideal must
2576          be defined on R.
2577KEYWORDS: Groebner cover, parametric ideal, canonical, discussion of
2578          parametric ideal, full representation.
2579EXAMPLE:  extend; shows an example"
2580{
2581  list L=#;
2582  list S=GC;
2583  ideal idp;
2584  ideal idq;
2585  int i; int j; int m; int s;
2586  m=0; i=1;
2587  while((i<=size(S)) and (m==0))
2588  {
2589    if(typeof(S[i][2])=="list"){m=1;}
2590    i++;
2591  }
2592  if(m==1){"Warning! grobcov has already extended bases"; return(S);}
2593  if(size(GC[1])!=5){"Warning! extend make sense only when grobcov has been called with options  'rep',2,'ext',0"; " "; return();}
2594  int repop=0;
2595  int start3=timer;
2596  int comment;
2597  for(i=1;i<=size(L) div 2;i++)
2598  {
2599    if(L[2*i-1]=="comment"){comment=L[2*i];}
2600    else
2601    {
2602      if(L[2*i-1]=="rep"){repop=L[2*i];}
2603    }
2604  }
2606  poly ext;
2607  int te=0;
2608  list SS;
2609  def R=basering;
2610  if (defined(@R)){te=1;}
2611  else{setglobalrings();}
2612  // Now extend
2613  for (i=1;i<=size(S);i++)
2614  {
2615    m=size(S[i][2]);
2616     for (j=1;j<=m;j++)
2617    {
2618      idp=S[i][4][1];
2619      idq=S[i][4][2];
2620      if (size(idp)>0)
2621      {
2623        kill ext;
2624        def ext=extend0(S[i][2][j],idp,idq);
2625        if (typeof(ext)=="poly")
2626        {
2627          S[i][2][j]=pnormalf(ext,idp,idq);
2628        }
2629        else
2630        {
2631          if(size(ext)==1)
2632          {
2633            S[i][2][j]=ext[1];
2634          }
2635          else
2636          {
2637            kill SS; list SS;
2638            for(s=1;s<=size(ext);s++)
2639            {
2640              ext[s]=pnormalf(ext[s],idp,idq);
2641            }
2642            for(s=1;s<=size(S[i][2]);s++)
2643            {
2644              if(s!=j){SS[s]=S[i][2][s];}
2645              else{SS[s]=ext;}
2646            }
2647            S[i][2]=SS;
2648          }
2649        }
2650      }
2651    }
2652  }
2653  // NOW representation of the segments by option repop
2654  list Si; list nS;
2655  if (repop==0)
2656  {
2657    for(i=1;i<=size(S);i++)
2658    {
2659      Si=list(S[i][1],S[i][2],S[i][3],S[i][5]);
2660      nS[size(nS)+1]=Si;
2661    }
2662    S=nS;
2663  }
2664  else
2665  {
2666    if (repop==1)
2667    {
2668      for(i=1;i<=size(S);i++)
2669      {
2670        Si=list(S[i][1],S[i][2],S[i][4],S[i][5]);
2671        nS[size(nS)+1]=Si;
2672      }
2673      S=nS;
2674    }
2675    else
2676    {
2677      for(i=1;i<=size(S);i++)
2678      {
2679        Si=list(S[i][1],S[i][2],S[i][3],S[i][4],S[i][5]);
2680        nS[size(nS)+1]=Si;
2681      }
2682
2683    }
2684  }
2685  if(comment>=1){string("Time in extend = ",timer-start3);}
2686  if(te==0){kill @R; kill @RP; kill @P;}
2687  return(S);
2688}
2689example
2690{
2691  "EXAMPLE:"; echo = 2;
2692  ring R=(0,a0,b0,c0,a1,b1,c1),(x), dp;
2693  short=0;
2694  ideal S=a0*x^2+b0*x+c0,
2695            a1*x^2+b1*x+c1;
2696  def GCS=grobcov(S,"rep",2);
2697  GCS;
2698  def FGC=extend(GCS,"rep",0);
2699  // Full representation=
2700  FGC;
2701}
2702
2703// Auxiliary routine
2704// nonzerodivisor
2705// input:
2706//    poly g in Q[a],
2707//    list P=(p_1,..p_r) representing a minimal prime decomposition
2708// output
2709//    poly f such that f notin p_i for all i and
2710//           g-f in p_i for all i such that g notin p_i
2711static proc nonzerodivisor(poly gr, list Pr)
2712{
2713  def RR=basering;
2714  setring(@P);
2715  def g=imap(RR,gr);
2716  def P=imap(RR,Pr);
2717  int i; int k;  list J; ideal F;
2718  def f=g;
2719  ideal Pi;
2720  for (i=1;i<=size(P);i++)
2721  {
2722    option(redSB);
2723    Pi=std(P[i]);
2724    //attrib(Pi,"isSB",1);
2725    if (reduce(g,Pi,1)==0){J[size(J)+1]=i;}
2726  }
2727  for (i=1;i<=size(J);i++)
2728  {
2729    F=ideal(1);
2730    for (k=1;k<=size(P);k++)
2731    {
2732      if (k!=J[i])
2733      {
2734        F=idint(F,P[k]);
2735      }
2736    }
2737    f=f+F[1];
2738  }
2739  setring(RR);
2740  def fr=imap(@P,f);
2741  return(fr);
2742}
2743
2744// Auxiliary routine
2745// deltai
2746// input:
2747//   int i:
2748//   list LPr: (p1,..,pr) of prime components of an ideal in Q[a]
2749// output:
2750//   list (fr,fnr) of two polynomials that are equal on V(pi)
2751//       and fr=0 on V(P) \ V(pi), and fnr is nonzero on V(pj) for all j.
2752static proc deltai(int i, list LPr)
2753{
2754  def RR=basering;
2755  setring(@P);
2756  def LP=imap(RR,LPr);
2757  int j; poly p;
2758  def F=ideal(1);
2759  poly f;
2760  poly fn;
2761  ideal LPi;
2762  for (j=1;j<=size(LP);j++)
2763  {
2764    if (j!=i)
2765    {
2766      F=idint(F,LP[j]);
2767    }
2768  }
2769  p=0; j=1;
2770  while ((p==0) and (j<=size(F)))
2771  {
2772    LPi=LP[i];
2773    attrib(LPi,"isSB",1);
2774    p=reduce(F[j],LPi);
2775    j++;
2776  }
2777  f=F[j-1];
2778  fn=nonzerodivisor(f,LP);
2779  setring(RR);
2780  def fr=imap(@P,f);
2781  def fnr=imap(@P,fn);
2782  return(list(fr,fnr));
2783}
2784
2785// Auxiliary routine
2786// combine
2787// input: a list of pairs ((p1,P1),..,(pr,Pr)) where
2788//    ideal pi is a prime component
2789//    poly Pi is the polynomial in Q[a][x] on V(pi)\ V(Mi)
2790//    (p1,..,pr) are the prime decomposition of the lpp-segment
2791//    list crep =(ideal ida,ideal idb): the Crep of the segment.
2792//    list Pci of the intersecctions of all pj except the ith one
2793// output:
2794//    poly P on an open and dense set of V(p_1 int ... p_r)
2795static proc combine(list L, ideal F)
2796{
2797  // ATTENTION REVISE AND USE Pci and F
2798  int i; poly f;
2799  f=0;
2800  for(i=1;i<=size(L);i++)
2801  {
2802    f=f+F[i]*L[i][2];
2803  }
2804//   f=elimconstfac(f);
2805  f=primepartZ(f);
2806  return(f);
2807}
2808
2809
2810//Auxiliary routine
2811// nullin
2812// input:
2813//   poly f:  a polynomial in Q[a]
2814//   ideal P: an ideal in Q[a]
2815//   called from ring @R
2816// output:
2817//   t:  with value 1 if f reduces modulo P, 0 if not.
2818static proc nullin(poly f,ideal P)
2819{
2820  int t;
2821  def RR=basering;
2822  setring(@P);
2823  def f0=imap(RR,f);
2824  def P0=imap(RR,P);
2825  attrib(P0,"isSB",1);
2826  if (reduce(f0,P0,1)==0){t=1;}
2827  else{t=0;}
2828  setring(RR);
2829  return(t);
2830}
2831
2832// Auxiliary routine
2833// monoms
2834// Input: A polynomial f
2835// Output: The list of leading terms
2836static proc monoms(poly f)
2837{
2838  list L;
2839  poly lm; poly lc; poly lp; poly Q; poly mQ;
2840  def p=f;
2841  int i=1;
2842  while (p!=0)
2843  {
2845    p=p-lm;
2848    L[size(L)+1]=list(lc,lp);
2849    i++;
2850  }
2851  return(L);
2852}
2853
2854// Auxiliary routine called by extend
2855// extend0
2856// input:
2857//   poly f: a generic polynomial in the basis
2858//   ideal idp: such that ideal(S)=idp
2859//   ideal idq: such that S=V(idp) \ V(idq)
2860////   NW the list of ((N1,W1),..,(Ns,Ws)) of red-rep of the grouped
2861////      segments in the lpp-segment  NO MORE USED
2862// output:
2863static proc extend0(poly f, ideal idp, ideal idq)
2864{
2865  matrix CC; poly Q; list NewMonoms;
2866  int i;  int j;  poly fout; ideal idout;
2867  list L=monoms(f);
2868  int nummonoms=size(L)-1;
2869  Q=L[1][1];
2870  if (nummonoms==0){return(f);}
2871  for (i=2;i<=size(L);i++)
2872  {
2873    CC=matrix(extendcoef(L[i][1],Q,idp,idq));
2874    NewMonoms[i-1]=list(CC,L[i][2]);
2875  }
2876  if (nummonoms==1)
2877  {
2878    for(j=1;j<=ncols(NewMonoms[1][1]);j++)
2879    {
2880      fout=NewMonoms[1][1][2,j]*L[1][2]+NewMonoms[1][1][1,j]*NewMonoms[1][2];
2881      //fout=pnormalf(fout,idp,W);
2882      if(ncols(NewMonoms[1][1])>1){idout[j]=fout;}
2883    }
2884    if(ncols(NewMonoms[1][1])==1){return(fout);} else{return(idout);}
2885  }
2886  else
2887  {
2888    list cfi;
2889    list coefs;
2890    for (i=1;i<=nummonoms;i++)
2891    {
2892      kill cfi; list cfi;
2893      for(j=1;j<=ncols(NewMonoms[i][1]);j++)
2894      {
2895        cfi[size(cfi)+1]=NewMonoms[i][1][2,j];
2896      }
2897      coefs[i]=cfi;
2898    }
2899    def indexpolys=findindexpolys(coefs);
2900    for(i=1;i<=size(indexpolys);i++)
2901    {
2902      fout=L[1][2];
2903      for(j=1;j<=nummonoms;j++)
2904      {
2905        fout=fout+(NewMonoms[j][1][1,indexpolys[i][j]])/(NewMonoms[j][1][2,indexpolys[i][j]])*NewMonoms[j][2];
2906      }
2907      fout=cleardenom(fout);
2908      if(size(indexpolys)>1){idout[i]=fout;}
2909    }
2910    if (size(indexpolys)==1){return(fout);} else{return(idout);}
2911  }
2912}
2913
2914// Auxiliary routine
2915// findindexpolys
2916// input:
2917//   list coefs=( (q11,..,q1r_1),..,(qs1,..,qsr_1) )
2918//               of denominators of the monoms
2919// output:
2920//   list ind=(v_1,..,v_t) of intvec
2921//        each intvec v=(i_1,..,is) corresponds to a polynomial in the sheaf
2922//        that will be built from it in extend procedure.
2923static proc findindexpolys(list coefs)
2924{
2925  int i; int j; intvec numdens;
2926  for(i=1;i<=size(coefs);i++)
2927  {
2928    numdens[i]=size(coefs[i]);
2929  }
2930  def RR=basering;
2931  setring(@P);
2932  def coefsp=imap(RR,coefs);
2933  ideal cof; list combpolys; intvec v; int te; list mp;
2934  for(i=1;i<=size(coefsp);i++)
2935  {
2936    cof=ideal(0);
2937    for(j=1;j<=size(coefsp[i]);j++)
2938    {
2939      cof[j]=factorize(coefsp[i][j],3);
2940    }
2941    coefsp[i]=cof;
2942  }
2943  for(j=1;j<=size(coefsp[1]);j++)
2944  {
2945    v[1]=j;
2946    te=1;
2947    for (i=2;i<=size(coefsp);i++)
2948    {
2949      mp=memberpos(coefsp[1][j],coefsp[i]);
2950      if(mp[1])
2951      {
2952        v[i]=mp[2];
2953      }
2954      else{v[i]=0;}
2955    }
2956    combpolys[j]=v;
2957  }
2958  combpolys=reform(combpolys,numdens);
2959  setring(RR);
2960  return(combpolys);
2961}
2962
2963// Auxiliary routine
2964// extendcoef: given Q,P in Q[a] where P/Q specializes on an open and dense subset
2965//      of the whole V(p1 int...int pr), it returns a basis of the module
2966//      of all syzygies equivalent to P/Q,
2967static proc extendcoef(poly P, poly Q, ideal idp, ideal idq)
2968{
2969  def RR=basering;
2970  setring(@P);
2971  def PL=ringlist(@P);
2972  PL[3][1][1]="dp";
2973  def P1=ring(PL);
2974  setring(P1);
2975  ideal idp0=imap(RR,idp);
2976  option(redSB);
2977  qring q=std(idp0);
2978  poly P0=imap(RR,P);
2979  poly Q0=imap(RR,Q);
2980  ideal PQ=Q0,-P0;
2981  module C=syz(PQ);
2982  setring(@P);
2983  def idp1=imap(RR,idp);
2984  def idq1=imap(RR,idq);
2985  def C1=matrix(imap(q,C));
2986  def redC=selectregularfun(C1,idp1,idq1);
2987  setring(RR);
2988  def CC=imap(@P,redC);
2989  return(CC);
2990}
2991
2992// Auxiliary routine
2993// selectregularfun
2994// input:
2995//   list L of the polynomials matrix CC
2996//      (we assume that one of them is non-null on V(N) \ V(M))
2997//   ideal N, ideal M: ideals representing the locally closed set V(N) \ V(M)
2998// assume to work in @P
2999static proc selectregularfun(matrix CC, ideal NN, ideal MM)
3000{
3001  int numcombused;
3002  def RR=basering;
3003  setring(@P);
3004  def C=imap(RR,CC);
3005  def N=imap(RR,NN);
3006  def M=imap(RR,MM);
3007  if (ncols(C)==1){return(C);}
3008
3009  int i; int j; int k; list c; intvec ci; intvec c0; intvec c1;
3010  list T; list T0; list T1; list LL; ideal N1;ideal M1; int te=0;
3011  for(i=1;i<=ncols(C);i++)
3012  {
3013    if((C[1,i]!=0) and (C[2,i]!=0))
3014    {
3015      if(c0==intvec(0)){c0[1]=i;}
3016      else{c0[size(c0)+1]=i;}
3017    }
3018  }
3019  def C1=submat(C,1..2,c0);
3020  for (i=1;i<=ncols(C1);i++)
3021  {
3022    c=comb(ncols(C1),i);
3023    for(j=1;j<=size(c);j++)
3024    {
3025      ci=c[j];
3026      numcombused++;
3027      if(i==1){N1=N+C1[2,j]; M1=M;}
3028      if(i>1)
3029      {
3030        kill c0; intvec c0 ; kill c1; intvec c1;
3031        c1=ci[size(ci)];
3032        for(k=1;k<size(ci);k++){c0[k]=ci[k];}
3033        T0=searchinlist(c0,LL);
3034        T1=searchinlist(c1,LL);
3035        N1=T0[1]+T1[1];
3036        M1=intersect(T0[2],T1[2]);
3037      }
3038      T=list(ci,PtoCrep(Prep(N1,M1)));
3039      LL[size(LL)+1]=T;
3040      if(equalideals(T[2][1],ideal(1))){te=1; break;}
3041    }
3042    if(te){break;}
3043  }
3044  ci=T[1];
3045  def Cs=submat(C1,1..2,ci);
3046  setring(RR);
3047  return(imap(@P,Cs));
3048}
3049
3050// Auxiliary routine
3051// searchinlist
3052// input:
3053//   intvec c:
3054//   list L=( (c1,T1),..(ck,Tk) )
3055//      where the c's are assumed to be intvects
3056// output:
3057//   object T with index c
3058static proc searchinlist(intvec c,list L)
3059{
3060  int i; list T;
3061  for(i=1;i<=size(L);i++)
3062  {
3063    if (L[i][1]==c)
3064    {
3065      T=L[i][2];
3066      break;
3067    }
3068  }
3069  return(T);
3070}
3071
3072
3073// Auxiliary routine
3074// selectminsheaves
3075// Input: L=((v_11,..,v_1k_1),..,(v_s1,..,v_sk_s))
3076//    where:
3077//    The s lists correspond to the s coefficients of the polynomial f
3078//    (v_i1,..,v_ik_i) correspond to the k_i intvec v_ij of the
3079//    spezializations of the jth rekpresentant (Q,P) of the ith coefficient
3080//    v_ij is an intvec of size equal to the number of little segments
3081//    forming the lpp-segment of 0,1, where 1 represents that it specializes
3082//    to non-zedro an the whole little segment and 0 if not.
3083// Output: S=(w_1,..,w_j)
3084//    where the w_l=(n_l1,..,n_ls) are intvec of length size(L), where
3085//    n_lt fixes which element of (v_t1,..,v_tk_t) is to be
3086//    choosen to form the tth (Q,P) for the lth element of the sheaf
3087//    representing the I-regular function.
3088// The selection is done to obtian the minimal number of elements
3089//    of the sheaf that specializes to non-null everywhere.
3090static proc selectminsheaves(list L)
3091{
3092  list C=allsheaves(L);
3093  return(smsheaves(C[1],C[2]));
3094}
3095
3096// Auxiliary routine
3097// smsheaves
3098// Input:
3099//   list C of all the combrep
3100//   list L of the intvec that correesponds to each element of C
3101// Output:
3102//   list LL of the subsets of C that cover all the subsegments
3103//   (the union of the corresponding L(C) has all 1).
3104static proc smsheaves(list C, list L)
3105{
3106  int i; int i0; intvec W;
3107  int nor; int norn;
3108  intvec p;
3109  int sp=size(L[1]); int j0=1;
3110  for (i=1;i<=sp;i++){p[i]=1;}
3111  while (p!=0)
3112  {
3113    i0=0; nor=0;
3114    for (i=1; i<=size(L); i++)
3115    {
3116      norn=numones(L[i],pos(p));
3117      if (nor<norn){nor=norn; i0=i;}
3118    }
3119    W[j0]=i0;
3120    j0++;
3121    p=actualize(p,L[i0]);
3122  }
3123  list LL;
3124  for (i=1;i<=size(W);i++)
3125  {
3126    LL[size(LL)+1]=C[W[i]];
3127  }
3128  return(LL);
3129}
3130
3131// Auxiliary routine
3132// allsheaves
3133// Input: L=((v_11,..,v_1k_1),..,(v_s1,..,v_sk_s))
3134//    where:
3135//    The s lists correspond to the s coefficients of the polynomial f
3136//    (v_i1,..,v_ik_i) correspond to the k_i intvec v_ij of the
3137//    spezializations of the jth rekpresentant (Q,P) of the ith coefficient
3138//    v_ij is an intvec of size equal to the number of little segments
3139//    forming the lpp-segment of 0,1, where 1 represents that it specializes
3140//    to non-zero on the whole little segment and 1 if not.
3141// Output:
3142//    (list LL, list LLS)  where
3143//    LL is the list of all combrep
3144//    LLS is the list of intvec of the corresponding elements of LL
3145static proc allsheaves(list L)
3146{
3147  intvec V; list LL; intvec W; int r; intvec U;
3148  int i; int j; int k;
3149  int s=size(L[1][1]); // s = number of little segments of the lpp-segment
3150  list LLS;
3151  for (i=1;i<=size(L);i++)
3152  {
3153    V[i]=size(L[i]);
3154  }
3155  LL=combrep(V);
3156  for (i=1;i<=size(LL);i++)
3157  {
3158    W=LL[i];   // size(W)= number of coefficients of the polynomial
3159    kill U; intvec U;
3160    for (j=1;j<=s;j++)
3161    {
3162      k=1; r=1; U[j]=1;
3163      while((r==1) and (k<=size(W)))
3164      {
3165        if(L[k][W[k]][j]==0){r=0; U[j]=0;}
3166        k++;
3167      }
3168    }
3169    LLS[i]=U;
3170  }
3171  return(list(LL,LLS));
3172}
3173
3174// Auxiliary routine
3175// numones
3176// Input:
3177//   intvec v of (0,1) in each position
3178//   intvec pos: the positions to test
3179// Output:
3180//   int nor: the nuber of 1 of v in the positions given by pos.
3181static proc numones(intvec v, intvec pos)
3182{
3183  int i; int n;
3184  for (i=1;i<=size(pos);i++)
3185  {
3186    if (v[pos[i]]==1){n++;}
3187  }
3188  return(n);
3189}
3190
3191// Auxiliary routine
3192// actualize: actualizes zeroes of p
3193// Input:
3194//   intvec p: of zeroes and ones
3195//   intvec c: of zeroes and ones (of the same length)
3196// Output;
3197//   intvec pp: of zeroes and ones, where a 0 stays in pp[i] if either
3199static proc actualize(intvec p, intvec c)
3200{
3201  int i; intvec pp=p;
3202  for (i=1;i<=size(p);i++)
3203  {
3204    if ((pp[i]==1) and (c[i]==1)){pp[i]=0;}
3205  }
3206  return(pp);
3207}
3208
3209
3210// Auxiliary routine
3211static proc reducemodN(poly f,ideal E)
3212{
3213  def RR=basering;
3214  setring(@RPt);
3215  def fa=imap(RR,f);
3216  def Ea=imap(RR,E);
3217  attrib(Ea,"isSB",1);
3218  // option(redSB);
3219  // Ea=std(Ea);
3220  fa=reduce(fa,Ea);
3221  setring(RR);
3222  def f1=imap(@RPt,fa);
3223  return(f1);
3224}
3225
3226// Auxiliary routine
3227// intersp: computes the intersection of the ideals in S in @P
3228static proc intersp(list S)
3229{
3230  def RR=basering;
3231  setring(@P);
3232  def SP=imap(RR,S);
3233  option(returnSB);
3234  def NP=intersect(SP[1..size(SP)]);
3235  setring(RR);
3236  return(imap(@P,NP));
3237}
3238
3239// Auxiliary routine
3242{
3243  int te;
3244  def RR=basering;
3245  setring(@P);
3246  def fp=imap(RR,f);
3247  def idap=imap(RR,ida);
3248  poly @t;
3249  ring H=0,@t,dp;
3250  def PH=@P+H;
3251  setring(PH);
3252  def fH=imap(@P,fp);
3253  def idaH=imap(@P,idap);
3254  idaH[size(idaH)+1]=1-@t*fH;
3255  option(redSB);
3256  def G=std(idaH);
3257  if (G==1){te=1;} else {te=0;}
3258  setring(RR);
3259  return(te);
3260}
3261
3262// // Auxiliary routine
3263// // NonNull: returns 1 if the poly f is nonnull on V(E) \ V(N), 0 otherwise.
3264// // Input:
3265// //   f: polynomial
3266// //   E: null ideal
3267// //   N: nonnull ideal
3268// // Output:
3269// //   1 if f is nonnul in V(P) \ V(Q),
3270// //   0 if it has zeroes inside.
3271// //  Works in @P
3272// proc NonNull(poly f, ideal E, ideal N)
3273// {
3274//   int te=1; int i;
3275//   def RR=basering;
3276//   setring(@P);
3277//   def fp=imap(RR,f);
3278//   def Ep=imap(RR,E);
3279//   def Np=imap(RR,N);
3280//   ideal H;
3281//   ideal Ef=Ep+fp;
3282//   for (i=1;i<=size(Np);i++)
3283//   {
3285//     if (te==0){break;}
3286//   }
3287//   setring(RR);
3288//   return(te);
3289// }
3290
3291// Auxiliary routine
3292// selectextendcoef
3293// input:
3294//    matrix CC: CC=(p_a1 .. p_ar_a)
3295//                  (q_a1 .. q_ar_a)
3296//            the matrix of elements of a coefficient in oo[a].
3297//    (ideal ida, ideal idb): the canonical representation of the segment S.
3298// output:
3299//    list caout
3300//            the minimum set of elements of CC needed such that at least one
3301//            of the q's is non-null on S, as well as the C-rep of of the
3302//            points where the q's are null on S.
3303//            The elements of caout are of the form (p,q,prep);
3304static proc selectextendcoef(matrix CC, ideal ida, ideal idb)
3305{
3306  def RR=basering;
3307  setring(@P);
3308  def ca=imap(RR,CC);
3309  def E0=imap(RR,ida);
3310  ideal E;
3311  def N=imap(RR,idb);
3312  int r=ncols(ca);
3313  int i; int te=1; list com; int j; int k; intvec c; list prep;
3314  list cs; list caout;
3315  i=1;
3316  while ((i<=r) and (te))
3317  {
3318    com=comb(r,i);
3319    j=1;
3320    while((j<=size(com)) and (te))
3321    {
3322      E=E0;
3323      c=com[j];
3324      for (k=1;k<=i;k++)
3325      {
3326        E=E+ca[2,c[k]];
3327      }
3328      prep=Prep(E,N);
3329      if (i==1)
3330      {
3331        cs[j]=list(ca[1,j],ca[2,j],prep);
3332      }
3333      if ((size(prep)==1) and (equalideals(prep[1][1],ideal(1))))
3334      {
3335        te=0;
3336        for(k=1;k<=size(c);k++)
3337        {
3338          caout[k]=cs[c[k]];
3339        }
3340      }
3341      j++;
3342    }
3343    i++;
3344  }
3345  if (te){"error: extendcoef does not extend to the whole S";}
3346  setring(RR);
3347  return(imap(@P,caout));
3348}
3349
3350// Auxiliary routine
3351// plusP
3352// Input:
3353//   ideal E1: in some basering (depends only on the parameters)
3354//   ideal E2: in some basering (depends only on the parameters)
3355// Output:
3356//   ideal Ep=E1+E2; computed in @P
3357static proc plusP(ideal E1,ideal E2)
3358{
3359  def RR=basering;
3360  setring(@P);
3361  def E1p=imap(RR,E1);
3362  def E2p=imap(RR,E2);
3363  def Ep=E1p+E2p;
3364  setring(RR);
3365  return(imap(@P,Ep));
3366}
3367
3368// Auxiliary routine
3369// reform
3370// input:
3371//   list combpolys: (v1,..,vs)
3372//      where vi are intvec.
3373//   output outcomb: (w1,..,wt)
3374//      whre wi are intvec.
3375//      All the vi without zeroes are in outcomb, and those with zeroes are
3376//         combined to form new intvec with the rest
3377static proc reform(list combpolys, intvec numdens)
3378{
3379  list combp0; list combp1; int i; int j; int k; int l; list rest; intvec notfree;
3380  list free; intvec free1; int te; intvec v;  intvec w;
3381  int nummonoms=size(combpolys[1]);
3382  for(i=1;i<=size(combpolys);i++)
3383  {
3384    if(memberpos(0,combpolys[i])[1])
3385    {
3386      combp0[size(combp0)+1]=combpolys[i];
3387    }
3388    else {combp1[size(combp1)+1]=combpolys[i];}
3389  }
3390  for(i=1;i<=nummonoms;i++)
3391  {
3392    kill notfree; intvec notfree;
3393    for(j=1;j<=size(combpolys);j++)
3394    {
3395      if(combpolys[j][i]<>0)
3396      {
3397        if(notfree[1]==0){notfree[1]=combpolys[j][i];}
3398        else{notfree[size(notfree)+1]=combpolys[j][i];}
3399      }
3400    }
3401    kill free1; intvec free1;
3402    for(j=1;j<=numdens[i];j++)
3403    {
3404      if(memberpos(j,notfree)[1]==0)
3405      {
3406        if(free1[1]==0){free1[1]=j;}
3407        else{free1[size(free1)+1]=j;}
3408      }
3409      free[i]=free1;
3410    }
3411  }
3412  list amplcombp; list aux;
3413  for(i=1;i<=size(combp0);i++)
3414  {
3415    v=combp0[i];
3416    kill amplcombp; list amplcombp;
3417    amplcombp[1]=intvec(v[1]);
3418    for(j=2;j<=size(v);j++)
3419    {
3420      if(v[j]!=0)
3421      {
3422        for(k=1;k<=size(amplcombp);k++)
3423        {
3424          w=amplcombp[k];
3425          w[size(w)+1]=v[j];
3426          amplcombp[k]=w;
3427        }
3428      }
3429      else
3430      {
3431        kill aux; list aux;
3432        for(k=1;k<=size(amplcombp);k++)
3433        {
3434          for(l=1;l<=size(free[j]);l++)
3435          {
3436            w=amplcombp[k];
3437            w[size(w)+1]=free[j][l];
3438            aux[size(aux)+1]=w;
3439          }
3440        }
3441        amplcombp=aux;
3442      }
3443    }
3444    for(j=1;j<=size(amplcombp);j++)
3445    {
3446      combp1[size(combp1)+1]=amplcombp[j];
3447    }
3448  }
3449  return(combp1);
3450}
3451
3452
3453// Auxiliary routine
3454// precombint
3455// input:  L: list of ideals (works in @P)
3456// output: F0: ideal of polys. F0[i] is a poly in the intersection of
3457//             all ideals in L except in the ith one, where it is not.
3458//             L=(p1,..,ps);  F0=(f1,..,fs);
3459//             F0[i] \in intersect_{j#i} p_i
3460static proc precombint(list L)
3461{
3462  int i; int j; int tes;
3463  def RR=basering;
3464  setring(@P);
3465  list L0; list L1; list L2; list L3; ideal F;
3466  L0=imap(RR,L);
3467  L1[1]=L0[1]; L2[1]=L0[size(L0)];
3468  for (i=2;i<=size(L0)-1;i++)
3469  {
3470    L1[i]=intersect(L1[i-1],L0[i]);
3471    L2[i]=intersect(L2[i-1],L0[size(L0)-i+1]);
3472  }
3473  L3[1]=L2[size(L2)];
3474  for (i=2;i<=size(L0)-1;i++)
3475  {
3476    L3[i]=intersect(L1[i-1],L2[size(L0)-i]);
3477  }
3478  L3[size(L0)]=L1[size(L1)];
3479  for (i=1;i<=size(L3);i++)
3480  {
3481    option(redSB); L3[i]=std(L3[i]);
3482  }
3483  for (i=1;i<=size(L3);i++)
3484  {
3485    tes=1; j=0;
3486    while((tes) and (j<size(L3[i])))
3487    {
3488      j++;
3489      option(redSB);
3490      L0[i]=std(L0[i]);
3491      if(reduce(L3[i][j],L0[i])!=0){tes=0; F[i]=L3[i][j];}
3492    }
3493    if (tes){"ERROR a polynomial in all p_j except p_i was not found";}
3494  }
3495  setring(RR);
3496  def F0=imap(@P,F);
3497  return(F0);
3498}
3499
3500
3501// Auxiliary routine
3502// minAssGTZ eliminating denominators
3503static proc minGTZ(ideal N);
3504{
3505  int i; int j;
3506  def L=minAssGTZ(N);
3507  for(i=1;i<=size(L);i++)
3508  {
3509    for(j=1;j<=size(L[i]);j++)
3510    {
3511      L[i][j]=cleardenom(L[i][j]);
3512    }
3513  }
3514  return(L);
3515}
3516
3517//********************* Begin KapurSunWang *************************
3518
3519// Auxiliary routine
3520// inconsistent
3521// Input:
3522//   ideal E: of null conditions
3523//   ideal N: of non-null conditions representing V(E) \ V(N)
3524// Output:
3525//   1 if V(E) \ V(N) = empty
3526//   0 if not
3527static proc inconsistent(ideal E, ideal N)
3528{
3529  int j;
3530  int te=1;
3531  def R=basering;
3532  setring(@P);
3533  def EP=imap(R,E);
3534  def NP=imap(R,N);
3535  poly @t;
3536  ring H=0,@t,dp;
3537  def RH=@P+H;
3538  setring(RH);
3539  def EH=imap(@P,EP);
3540  def NH=imap(@P,NP);
3541  ideal G;
3542  j=1;
3543  while((te==1) and j<=size(NH))
3544  {
3545    G=EH+(1-@t*NH[j]);
3546    option(redSB);
3547    G=std(G);
3548    if (G[1]!=1){te=0;}
3549    j++;
3550  }
3551  setring(R);
3552  return(te);
3553}
3554
3555// Auxiliary routine
3556// MDBasis: Minimal Dickson Basis
3557static proc MDBasis(ideal G)
3558{
3559  int i; int j; int te=1;
3560  G=sortideal(G);
3561  ideal MD=G[1];
3562  poly lm;
3563  for (i=2;i<=size(G);i++)
3564  {
3565    te=1;
3567    j=1;
3568    while ((te==1) and (j<=size(MD)))
3569    {
3571      j++;
3572    }
3573    if (te==1)
3574    {
3575      MD[size(MD)+1]=(G[i]);
3576    }
3577  }
3578  return(MD);
3579}
3580
3581// Auxiliary routine
3582// primepartZ
3583static proc primepartZ(poly f);
3584{
3585  def cp=content(f);
3586  def fp=f/cp;
3587  return(fp);
3588}
3589
3590// LCMLC
3591static proc LCMLC(ideal H)
3592{
3593  int i;
3594  def R=basering;
3595  setring(@RP);
3596  def HH=imap(R,H);
3597  poly h=1;
3598  for (i=1;i<=size(HH);i++)
3599  {
3600    h=lcm(h,HH[i]);
3601  }
3602  setring(R);
3603  def hh=imap(@RP,h);
3604  return(hh);
3605}
3606
3607// KSW: Kapur-Sun-Wang algorithm for computing a CGS
3608// Input:
3609//   F:   parametric ideal to be discussed
3610//   Options:
3611//     \"out\",0 Transforms the description of the segments into
3612//     canonical P-representation form.
3613//     \"out\",1 Original KSW routine describing the segments as
3614//     difference of varieties
3615//   The ideal must be defined on C[parameters][variables]
3616// Output:
3617//   With option \"out\",0 :
3618//     ((lpp,
3619//       (1,B,((p_1,(p_11,..,p_1k_1)),..,(p_s,(p_s1,..,p_sk_s)))),
3620//       string(lpp)
3621//      )
3622//      ,..,
3623//      (lpp,
3624//       (k,B,((p_1,(p_11,..,p_1k_1)),..,(p_s,(p_s1,..,p_sk_s)))),
3625//       string(lpp))
3626//      )
3627//     )
3628//   With option \"out\",1 ((default, original KSW) (shorter to be computed,
3629//                    but without canonical description of the segments.
3630//     ((B,E,N),..,(B,E,N))
3631static proc KSW(ideal F, list #)
3632{
3633  setglobalrings();
3634  int start=timer;
3635  ideal E=ideal(0);
3636  ideal N=ideal(1);
3637  int comment=0;
3638  int out=1;
3639  int i;
3640  def L=#;
3641  if (size(L)>0)
3642  {
3643    for (i=1;i<=size(L) div 2;i++)
3644    {
3645      if (L[2*i-1]=="null"){E=L[2*i];}
3646      else
3647      {
3648        if (L[2*i-1]=="nonnull"){N=L[2*i];}
3649        else
3650        {
3651          if (L[2*i-1]=="comment"){comment=L[2*i];}
3652          else
3653          {
3654            if (L[2*i-1]=="out"){out=L[2*i];}
3655          }
3656        }
3657      }
3658    }
3659  }
3660  if (comment>0){string("Begin KSW with null = ",E," nonnull = ",N);}
3661  def CG=KSW0(F,E,N,comment);
3662  if (comment>0)
3663  {
3664    string("Number of segments in KSW (total) = ",size(CG));
3665    string("Time in KSW = ",timer-start);
3666  }
3667  if(out==0)
3668  {
3669    CG=KSWtocgsdr(CG);
3670    CG=groupKSWsegments(CG);
3671    if (comment>0)
3672    {
3673      string("Number of lpp segments = ",size(CG));
3674      string("Time in KSW + group + Prep = ",timer-start);
3675    }
3676  }
3677  if(defined(@P)){kill @P; kill @R; kill @RP;}
3678  return(CG);
3679}
3680
3681// Auxiliary routine
3682// sqf
3683// This is for releases of Singular before 3-5-1
3684// proc sqf(poly f)
3685// {
3686//  def RR=basering;
3687//  setring(@P);
3688//  def ff=imap(RR,f);
3689//  def G=sqrfree(ff);
3690//  poly fff=1;
3691//  int i;
3692//  for (i=1;i<=size(G);i++)
3693//  {
3694//    fff=fff*G[i];
3695//  }
3696//  setring(RR);
3697//   def ffff=imap(@P,fff);
3698//   return(ffff);
3699// }
3700
3701// Auxiliary routine
3702// sqf
3703static proc sqf(poly f)
3704{
3705  def RR=basering;
3706  setring(@P);
3707  def ff=imap(RR,f);
3708  poly fff=sqrfree(ff,3);
3709  setring(RR);
3710  def ffff=imap(@P,fff);
3711  return(ffff);
3712}
3713
3714
3715// Auxiliary routine
3716// KSW0: Kapur-Sun-Wang algorithm for computing a CGS, called by KSW
3717// Input:
3718//   F:   parametric ideal to be discussed
3719//   Options:
3720//   The ideal must be defined on C[parameters][variables]
3721// Output:
3722static proc KSW0(ideal F, ideal E, ideal N, int comment)
3723{
3724  def R=basering;
3725  int i; int j; list emp;
3726  list CGS;
3727  ideal N0;
3728  for (i=1;i<=size(N);i++)
3729  {
3730    N0[i]=sqf(N[i]);
3731  }
3732  ideal E0;
3733  for (i=1;i<=size(E);i++)
3734  {
3736  }
3737  setring(@P);
3738  ideal E1=imap(R,E0);
3739  E1=std(E1);
3740  ideal N1=imap(R,N0);
3741  N1=std(N1);
3742  setring(R);
3743  E0=imap(@P,E1);
3744  N0=imap(@P,N1);
3745  if (inconsistent(E0,N0)==1)
3746  {
3747    return(emp);
3748  }
3749  setring(@RP);
3750  def FRP=imap(R,F);
3751  def ERP=imap(R,E);
3752  FRP=FRP+ERP;
3753  option(redSB);
3754  def GRP=std(FRP);
3755  setring(R);
3756  def G=imap(@RP,GRP);
3757  if (memberpos(1,G)[1]==1)
3758  {
3759    if(comment>1){"Basis 1 is found"; E; N;}
3760    list KK; KK[1]=list(E0,N0,ideal(1));
3761    return(KK);
3762   }
3763  ideal Gr; ideal Gm; ideal GM;
3764  for (i=1;i<=size(G);i++)
3765  {
3766    if (variables(G[i])[1]==0){Gr[size(Gr)+1]=G[i];}
3767    else{Gm[size(Gm)+1]=G[i];}
3768  }
3769  ideal Gr0;
3770  for (i=1;i<=size(Gr);i++)
3771  {
3772    Gr0[i]=sqf(Gr[i]);
3773  }
3774
3775
3776  Gr=elimrepeated(Gr0);
3777  ideal GrN;
3778  for (i=1;i<=size(Gr);i++)
3779   {
3780    for (j=1;j<=size(N0);j++)
3781    {
3782      GrN[size(GrN)+1]=sqf(Gr[i]*N0[j]);
3783    }
3784  }
3785  if (inconsistent(E,GrN)){;}
3786  else
3787  {
3788    if(comment>1){"Basis 1 is found in a branch with arguments"; E; GrN;}
3789    CGS[size(CGS)+1]=list(E,GrN,ideal(1));
3790  }
3791  if (inconsistent(Gr,N0)){return(CGS);}
3792  GM=Gm;
3793  Gm=MDBasis(Gm);
3794  ideal H;
3795  for (i=1;i<=size(Gm);i++)
3796  {
3798  }
3799  H=facvar(H);
3800  poly h=sqf(LCMLC(H));
3801  if(comment>1){"H = "; H; "h = "; h;}
3802  ideal Nh=N0;
3803  if(size(N0)==0){Nh=h;}
3804  else
3805  {
3806    for (i=1;i<=size(N0);i++)
3807    {
3808      Nh[i]=sqf(N0[i]*h);
3809    }
3810  }
3811  if (inconsistent(Gr,Nh)){;}
3812  else
3813  {
3814    CGS[size(CGS)+1]=list(Gr,Nh,Gm);
3815  }
3816  poly hc=1;
3817  list KS;
3818  ideal GrHi;
3819  for (i=1;i<=size(H);i++)
3820  {
3821    kill GrHi;
3822    ideal GrHi;
3823    Nh=N0;
3824    if (i>1){hc=sqf(hc*H[i-1]);}
3825    for (j=1;j<=size(N0);j++){Nh[j]=sqf(N0[j]*hc);}
3826    if (equalideals(Gr,ideal(0))==1){GrHi=H[i];}
3827    else {GrHi=Gr,H[i];}
3828//     else {for (j=1;j<=size(Gr);j++){GrHi[size(GrHi)+1]=Gr[j]*H[i];}}
3829    if(comment>1){"Call to KSW with arguments "; GM; GrHi;  Nh;}
3830    KS=KSW0(GM,GrHi,Nh,comment);
3831    for (j=1;j<=size(KS);j++)
3832    {
3833      CGS[size(CGS)+1]=KS[j];
3834    }
3835    if(comment>1){"CGS after KSW = "; CGS;}
3836  }
3837  return(CGS);
3838}
3839
3840// Auxiliary routine
3841// KSWtocgsdr
3842static proc KSWtocgsdr(list L)
3843{
3844  int i; list CG; ideal B; ideal lpp; int j; list NKrep;
3845  for(i=1;i<=size(L);i++)
3846  {
3847    B=redgbn(L[i][3],L[i][1],L[i][2]);
3848    lpp=ideal(0);
3849    for(j=1;j<=size(B);j++)
3850    {
3852    }
3853    NKrep=KtoPrep(L[i][1],L[i][2]);
3854    CG[i]=list(lpp,B,NKrep);
3855  }
3856  return(CG);
3857}
3858
3859// Auxiliary routine
3860// KtoPrep
3861// Computes the P-representaion of a K-representation (N,W) of a set
3862// input:
3863//    ideal E (null conditions)
3864//    ideal N (non-null conditions ideal)
3865// output:
3866//    the ((p_1,(p_11,..,p_1k_1)),..,(p_r,(p_r1,..,p_rk_r)));
3867//    the Prep of V(N) \ V(W)
3868static proc KtoPrep(ideal N, ideal W)
3869{
3870  int i; int j;
3871  if (N[1]==1)
3872  {
3873    L0[1]=list(ideal(1),list(ideal(1)));
3874    return(L0);
3875  }
3876  def RR=basering;
3877  setring(@P);
3878  ideal B; int te; poly f;
3879  ideal Np=imap(RR,N);
3880  ideal Wp=imap(RR,W);
3881  list L;
3882  list L0; list T0;
3883  L0=minGTZ(Np);
3884  for(j=1;j<=size(L0);j++)
3885  {
3886    option(redSB);
3887    L0[j]=std(L0[j]);
3888  }
3889  for(i=1;i<=size(L0);i++)
3890  {
3891    if(inconsistent(L0[i],Wp)==0)
3892    {
3893      B=L0[i]+Wp;
3894      T0=minGTZ(B);
3895      option(redSB);
3896      for(j=1;j<=size(T0);j++)
3897      {
3898        T0[j]=std(T0[j]);
3899      }
3900      L[size(L)+1]=list(L0[i],T0);
3901    }
3902  }
3903  setring(RR);
3904  def LL=imap(@P,L);
3905  return(LL);
3906}
3907
3908// Auxiliary routine
3909// groupKSWsegments
3910// input:  the list of vertices of KSW
3911// output: the same terminal vertices grouped by lpp
3912static proc groupKSWsegments(list T)
3913{
3914  int i; int j;
3915  list L;
3916  list lpp; list lppor;
3917  list kk;
3918  lpp[1]=T[1][1]; j=1;
3919  lppor[1]=intvec(1);
3920  for(i=2;i<=size(T);i++)
3921  {
3922    kk=memberpos(T[i][1],lpp);
3923    if(kk[1]==0){j++; lpp[j]=T[i][1]; lppor[j]=intvec(i);}
3924    else{lppor[kk[2]][size(lppor[kk[2]])+1]=i;}
3925  }
3926  list ll;
3927  for (j=1;j<=size(lpp);j++)
3928  {
3929    kill ll; list ll;
3930    for(i=1;i<=size(lppor[j]);i++)
3931    {
3932      ll[size(ll)+1]=list(i,T[lppor[j][i]][2],T[lppor[j][i]][3]);
3933    }
3934    L[j]=list(lpp[j],ll,string(lpp[j]));
3935  }
3936  return(L);
3937}
3938
3939//********************* End KapurSunWang *************************
3940
3941//******************** Begin locus ******************************
3942
3943// indepparameters
3944// Auxiliary routine to detect 'Special' components of the locus
3945// Input: ideal B
3946// Output:
3947//   1 if the solutions of the ideal do not depend on the parameters
3948//   0 if they depend
3949static proc indepparameters(ideal B)
3950{
3951  def R=basering;
3952  ideal v=variables(B);
3953  setring @RP;
3954  def BP=imap(R,B);
3955  def vp=imap(R,v);
3956  ideal varpar=variables(BP);
3957  int te;
3958  te=equalideals(vp,varpar);
3959  setring(R);
3960  if(te){return(1);}
3961  else{return(0);}
3962}
3963
3964// dimP0: Auxiliary routine
3965// if the dimension in @P of an ideal in the parameters has dimension 0 then it returns 0
3966// else it retuns 1
3967static proc dimP0(ideal N)
3968{
3969  def R=basering;
3970  setring(@P);
3971  int te=1;
3972  def NP=imap(R,N);
3973  attrib(NP,"IsSB",1);
3974  int d=dim(std(NP));
3975  if(d==0){te=0;}
3976  setring(R);
3977  return(te);
3978}
3979
3980// Takes a list of intvec and sorts it and eliminates repeated elements.
3981// Auxiliary routine
3982static proc sortpairs(L)
3983{
3984  def L1=sort(L);
3985  def L2=elimrepeated(L1[1]);
3986  return(L2);
3987}
3988
3989// Eliminates the pairs of L1 that are also in L2.
3990// Auxiliary routine
3991static proc minuselements(list L1,list L2)
3992{
3993  int i;
3994  list L3;
3995  for (i=1;i<=size(L1);i++)
3996  {
3997    if(not(memberpos(L1[i],L2)[1])){L3[size(L3)+1]=L1[i];}
3998  }
3999  return(L3);
4000}
4001
4002// NorSing
4003// Input:
4004//   ideal B: the basis of a component of the grobcov
4005//   ideal E: the top of the component (assumed to be of dimension > 0 (single equation)
4006//   ideal N: the holes of the component
4007// Output:
4008//   int d: the dimension of B on the variables reduced by the holes.
4009//     if d>0    then the component is 'Normal'
4010//     if d==0 then the component is 'Singular'
4011static proc NorSing(ideal B, ideal E, ideal N, list #)
4012{
4013  int i; int j; int Fenv=0; int env; int dd;
4014  list DD=#;
4015  def RR=basering;
4016  int moverdim=2;
4017  int version=0;
4018  int nv=nvars(RR);
4019  if(nv<4){version=1;}
4020  int d;
4021  poly F;
4022  for(i=1;i<=(size(DD) div 2);i++)
4023  {
4024    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
4025    if(DD[2*i-1]=="version"){version=DD[2*i];}
4026    if(DD[2*i-1]=="family"){F=DD[2*i];}
4027  }
4028  if(F!=0){Fenv=1;}
4029  //"T_B="; B; "E="; E; "N="; N;
4030  list L0;
4031  if(dimP0(E)==0){L0=2,"Normal";} // 2 es fals pero ha de ser >0 encara que sigui 0
4032  else
4033  {
4034    if(version==0)
4035    {
4036      //"T_B="; B;  // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part
4037      // independent of parameters giving the variables with dimension 0
4038     dd=indepparameters(B);
4039      if (dd==1){d=0; L0=d,string(B),determineF(B,F,E);}
4040      else{d=1; L0=2,"Normal";}
4041    }
4042    else
4043    {
4044      def RH=ringlist(RR);
4045      //"T_RH="; RH;
4046      def H=RH;
4047      H[1]=0;
4048      H[2]=RH[1][2]+RH[2];
4049      int n=size(H[2]);
4050      intvec ll;
4051      for(i=1;i<=n;i++)
4052      {
4053        ll[i]=1;
4054      }
4055      H[3][1][1]="lp";
4056      H[3][1][2]=ll;
4057      def RRH=ring(H);
4058      setring(RRH);
4059      ideal BH=imap(RR,B);
4060      ideal EH=imap(RR,E);
4061      ideal NH=imap(RR,N);
4062      if(Fenv==1){poly FH=imap(RR,F);}
4063      for(i=1;i<=size(EH);i++){BH[size(BH)+1]=EH[i];}
4064      BH=std(BH);  // MOLT COSTOS!!!
4065      ideal G;
4066      ideal r; poly q;
4067      for(i=1;i<=size(BH);i++)
4068      {
4069        r=factorize(BH[i],1);
4070        q=1;
4071        for(j=1;j<=size(r);j++)
4072        {
4073          if((pdivi(r[j],NH)[1] != 0) or (equalideals(ideal(NH),ideal(1))))
4074          {
4075            q=q*r[j];
4076          }
4077        }
4078        if(q!=1){G[size(G)+1]=q;}
4079      }
4080      setring RR;
4081      def GG=imap(RRH,G);
4082      ideal GGG;
4083      if(defined(L0)){kill L0; list L0;}
4084      for(i=1;i<=size(GG);i++)
4085      {
4086        if(indepparameters(GG[i])){GGG[size(GGG)+1]=GG[i];}
4087      }
4088      GGG=std(GGG);
4089      ideal GLM;
4090      for(i=1;i<=size(GGG);i++)
4091      {
4093      }
4094      attrib(GLM,"IsSB",1);
4095      d=dim(std(GLM));
4096      string antiim=string(GGG);
4097      L0=d,antiim;
4098      if(d==0)
4099      {
4100        //" ";string("Antiimage of Special component = ", GGG);
4101         if(Fenv==1)
4102        {
4103          L0[3]=determineF(GGG,F,E);
4104        }
4105      }
4106      else
4107      {
4108        L0[2]="Normal";
4109      }
4110    }
4111  }
4112  return(L0);
4113}
4114
4115static proc determineF(ideal A,poly F,ideal E)
4116{
4117  int env; int i;
4118  def RR=basering;
4119  def RH=ringlist(RR);
4120  def H=RH;
4121  H[1]=0;
4122  H[2]=RH[1][2]+RH[2];
4123  int n=size(H[2]);
4124  intvec ll;
4125  for(i=1;i<=n;i++)
4126  {
4127    ll[i]=1;
4128  }
4129  H[3][1][1]="lp";
4130  H[3][1][2]=ll;
4131  def RRH=ring(H);
4132
4133        //" ";string("Antiimage of Special component = ", GGG);
4134
4135   setring(RRH);
4136   list LL;
4137   def AA=imap(RR,A);
4138   def FH=imap(RR,F);
4139   def EH=imap(RR,E);
4140   ideal M=std(AA+FH);
4141   def rh=reduce(EH,M);
4142   if(rh==0){env=1;} else{env=0;}
4143   setring RR;
4144          //L0[3]=env;
4145    return(env);
4146}
4147
4148// DimPar
4149// Auxilliary routine to NorSing determining the dimension of a parametric ideal
4150// Does not use @P and define it directly because is assumes that
4151//                             variables and parameters have been inverted
4152static proc DimPar(ideal E)
4153{
4154  def RRH=basering;
4155  def RHx=ringlist(RRH);
4156  def Prin=ring(RHx[1]);
4157  setring(Prin);
4158  def E2=std(imap(RRH,E));
4159  def d=dim(E2);
4160  setring RRH;
4161  return(d);
4162}
4163
4164// locus0(G): Private routine used by locus (the public routine), that
4165//                builds the diferent components.
4166// input:      The output G of the grobcov (in generic representation, which is the default option)
4167//       Options: The algorithm allows the following options ar pair of arguments:
4168//                "movdim", d  : by default movdim is 2 but it can be set to other values
4169//                "version", v   :  There are two versions of the algorithm. ('version',1) is
4170//                 a full algorithm that always distinguishes correctly between 'Normal'
4171//                 and 'Special' components, whereas ('version',0) can decalre a component
4172//                 as 'Normal' being really 'Special', but is more effective. By default ('version',1)
4173//                 is used when the number of variables is less than 4 and 0 if not.
4174//                 The user can force to use one or other version, but it is not recommended.
4175//                 "system", ideal F: if the initial systrem is passed as an argument. This is actually not used.
4176//                 "comments", c: by default it is 0, but it can be set to 1.
4177//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
4178//                 The mover coordinates are to be placed as the last variables, and by default,
4179//                 its number is 2. If one consider locus problems in higer dimensions, the number of
4180//                 mover coordinates (placed as the last variables) is to be given as an option.
4181// output:
4182//         list, the canonical P-representation of the Normal and Non-Normal locus:
4183//              The Normal locus has two kind of components: Normal and Special.
4184//              The Non-normal locus has two kind of components: Accumulation and Degenerate.
4185//              This routine is compemented by locus that calls it in order to eliminate problems
4186//              with degenerate points of the mover.
4187//         The output components are given as
4188//              ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
4189//         The components are given in canonical P-representation of the subset.
4190//              If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
4191//              gives the depth of the component.
4192static proc locus0(list GG, list #)
4193{
4194  int te=0;
4195  int t1=1; int t2=1; int i;
4196  def R=basering;
4197  //if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
4198  //setglobalrings();
4199  // Options
4200  list DD=#;
4201  int moverdim=nvars(R);
4202  int version=0;
4203  int nv=nvars(R);
4204  if(nv<4){version=1;}
4205  int comment=0;
4206  ideal Fm;
4207  poly F;
4208  for(i=1;i<=(size(DD) div 2);i++)
4209  {
4210    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
4211    if(DD[2*i-1]=="version"){version=DD[2*i];}
4212    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
4213    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
4214    if(DD[2*i-1]=="family"){F=DD[2*i];}
4215  }
4216  list HHH;
4217  if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1){return(HHH);}
4218   list G1; list G2;
4219  def G=GG;
4220  list Q1; list Q2;
4221  int d; int j; int k;
4222  t1=1;
4223  for(i=1;i<=size(G);i++)
4224  {
4225    attrib(G[i][1],"IsSB",1);
4226    d=locusdim(G[i][2],moverdim);
4227    if(d==0){G1[size(G1)+1]=G[i];}
4228    else
4229    {
4230      if(d>0){G2[size(G2)+1]=G[i];}
4231    }
4232  }
4233  if(size(G1)==0){t1=0;}
4234  if(size(G2)==0){t2=0;}
4235  setring(@RP);
4236  if(t1)
4237  {
4238    list G1RP=imap(R,G1);
4239  }
4240  else {list G1RP;}
4241  list P1RP;
4242  ideal B;
4243  for(i=1;i<=size(G1RP);i++)
4244  {
4245     kill B;
4246     ideal B;
4247     for(k=1;k<=size(G1RP[i][3]);k++)
4248     {
4249       attrib(G1RP[i][3][k][1],"IsSB",1);
4250       G1RP[i][3][k][1]=std(G1RP[i][3][k][1]);
4251       for(j=1;j<=size(G1RP[i][2]);j++)
4252       {
4253         B[j]=reduce(G1RP[i][2][j],G1RP[i][3][k][1]);
4254       }
4255       P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B);
4256     }
4257  }
4258  setring(R);
4259  ideal h;
4260  if(t1)
4261  {
4262    def P1=imap(@RP,P1RP);
4263    for(i=1;i<=size(P1);i++)
4264    {
4265      for(j=1;j<=size(P1[i][3]);j++)
4266      {
4267        h=factorize(P1[i][3][j],1);
4268        P1[i][3][j]=h[1];
4269        for(k=2;k<=size(h);k++)
4270        {
4271          P1[i][3][j]=P1[i][3][j]*h[k];
4272        }
4273      }
4274    }
4275  }
4276  else{list P1;}
4277  ideal BB; int dd; list NS;
4278  for(i=1;i<=size(P1);i++)
4279  {
4280       NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD);
4281       dd=NS[1];
4282       if(dd==0){P1[i][3]=NS;} //"Special";
4283       else{P1[i][3]="Normal";}
4284  }
4285  list P2;
4286  for(i=1;i<=size(G2);i++)
4287  {
4288    for(k=1;k<=size(G2[i][3]);k++)
4289    {
4290      P2[size(P2)+1]=list(G2[i][3][k][1],G2[i][3][k][2]);
4291    }
4292  }
4293  list l;
4294  for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1;
4295  for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2;
4296
4297  setring(@P);
4298  ideal J;
4299  if(t1==1)
4300  {
4301    def C1=imap(R,P1);
4303   }
4304  else{list C1; list L1; kill P1; list P1;}
4305  if(t2==1)
4306  {
4307    def C2=imap(R,P2);
4309 }
4310  else{list L2; list C2; kill P2; list P2;}
4311    for(i=1;i<=size(L2);i++)
4312    {
4313      J=std(L2[i][2]);
4314      d=dim(J); // AQUI
4315      if(d==0)
4316      {
4317        L2[i][4]=string("Accumulation",L2[i][4]);
4318      }
4319      else{L2[i][4]=string("Degenerate",L2[i][4]);}
4320    }
4321  list LN;
4322  if(t1==1)
4323  {
4324    for(i=1;i<=size(L1);i++){LN[size(LN)+1]=L1[i];}
4325  }
4326  if(t2==1)
4327  {
4328    for(i=1;i<=size(L2);i++){LN[size(LN)+1]=L2[i];}
4329  }
4330  setring(R);
4331  def L=imap(@P,LN);
4332  for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
4333  //if(te==0){kill @R; kill @RP; kill @P;}
4334  list LL;
4335  for(i=1;i<=size(L);i++)
4336  {
4337    if(typeof(L[i][4])=="list") {L[i][4][1]="Special";}
4338    l[1]=L[i][2];
4339    l[2]=L[i][3];
4340    l[3]=L[i][4];
4341    l[4]=L[i][5];
4342    L[i]=l;
4343 }
4344 return(L);
4345}
4346
4347//  locus(G):  Special routine for determining the locus of points
4348//                 of  geometrical constructions.
4349//  input:      The output G of the grobcov (in generic representation, which is the default option for grobcov)
4350//  output:
4351//          list, the canonical P-representation of the Normal and Non-Normal locus:
4352//               The Normal locus has two kind of components: Normal and Special.
4353//               The Non-normal locus has two kind of components: Accumulation and Degenerate.
4354//               Normal components: for each point in the component,
4355//               the number of solutions in the variables is finite, and
4356//               the solutions depend on the point in the component if the component is not 0-dimensional.
4357//               Special components:  for each point in the component,
4358//               the number of solutions in the variables is finite,
4359//               the component is not 0-dimensional, but the solutions do not depend on the
4360//               values of the parameters in the component.
4361//               Accumlation points: are 0-dimensional components for which it exist
4362//               an infinite number of solutions.
4363//               Degenerate components: are components of dimension greater than 0 for which
4364//               for every point in the component there exist infinite solutions.
4365//          The output components are given as
4366//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
4367//          The components are given in canonical P-representation of the subset.
4368//               If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
4369//               gives the depth of the component of the constructible set.
4370proc locus(list GG, list #)
4371"USAGE: locus(G)
4372        The input must be the grobcov  of a parametrical ideal in Q[a][x],
4373        (a=parameters, x=variables). In fact a must be the tracer coordinates
4374        and x the mover coordinates and other auxiliary ones.
4375        Special routine for determining the locus of points
4376        of  geometrical constructions. Given a parametric ideal J
4377        representing the system determining the locus of points (a)
4378        who verify certain properties, the call to locus on the output of grobcov( J )
4379        determines the different classes of locus components, following
4380        the taxonomy defined in
4382        \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\".
4383        Computer-Aided Design 56 (2014) 22-33.
4384        The components can be Normal, Special, Accumulation, Degenerate.
4385OPTIONS: The algorithm allows the following options as pair of arguments:
4386         \"movdim\", d  : by default movdim is 2 but it can be set to other values,
4387         and represents the number of mever variables. they should be given as
4388         the last variables of the ring.
4389         \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is
4390         a full algorithm that always distinguishes correctly between 'Normal'
4391         and 'Special' components, whereas \("version\",0) can decalre a component
4392         as 'Normal' being really 'Special', but is more effective. By default (\"version\",1)
4393         is used when the number of variables is less than 4 and 0 if not.
4394         The user can force to use one or other version, but it is not recommended.
4395         \"system\", ideal F: if the initial system is passed as an argument. This is actually not used.
4396         \"comments\", c: by default it is 0, but it can be set to 1.
4397         Usually locus problems have mover coordinates, variables and tracer coordinates.
4398         The mover coordinates are to be placed as the last variables, and by default,
4399         its number is 2. If one consider locus problems in higer dimensions, the number of
4400         mover coordinates (placed as the last variables) is to be given as an option.
4401RETURN: The locus. The output is a list of the components ( C_1,.. C_n ) where
4402        C_i =  (p_i,(p_i1,..p_is_i), type_i,l evel_i )   and type_i can be
4403        'Normal', 'Special', Accumulation', 'Degenerate'. The 'Special' components
4405        is 0-dimensional for these kind of components.
4406        Normal components: for each point in the component,
4407        the number of solutions in the variables is finite, and
4408        the solutions depend on the point in the component.
4409        Special components:  for each point in the component,
4410        the number of solutions in the variables is finite. The
4411        antiimage of the component is 0-dimensional.
4412        Accumlation points: are 0-dimensional components whose
4413        antiimage is not zero-dimansional.
4414        Degenerate components: are components of dimension greater than 0
4415        whose antiimage is not-zero-diemansional.
4416        The components are given in canonical P-representation.
4417        The levels of a class of locus are 1,
4418        because they represent locally closed. sets.
4419NOTE: It can only be called after computing the grobcov of the
4420      parametrical ideal in generic representation ('ext',0),
4421      which is the default.
4422      The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
4423KEYWORDS: geometrical locus, locus, loci.
4424EXAMPLE: locus; shows an example"
4425{
4426  int tes=0; int i;
4427  def R=basering;
4428  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
4429  setglobalrings();
4430  // Options
4431  list DD=#;
4432  int moverdim=nvars(R);
4433  int version=0;
4434  int nv=nvars(R);
4435  if(nv<4){version=1;}
4436  int comment=0;
4437  ideal Fm;
4438  for(i=1;i<=(size(DD) div 2);i++)
4439  {
4440    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
4441    if(DD[2*i-1]=="version"){version=DD[2*i];}
4442    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
4443    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
4444    if(DD[2*i-1]=="family"){poly F=DD[2*i];}
4445  }
4446  int j; int k;
4447  def B0=GG[1][2];
4448  def H0=GG[1][3][1][1];
4449  if (equalideals(B0,ideal(1)) or (equalideals(H0,ideal(0))==0))
4450  {
4451    return(locus0(GG,DD));
4452  }
4453  else
4454  {
4455    int n=nvars(R);
4456    ideal vmov=var(n-1),var(n);
4457    ideal N;
4458    intvec xw; intvec yw;
4459    for(i=1;i<=n-1;i++){xw[i]=0;}
4460    xw[n]=1;
4461    for(i=1;i<=n;i++){yw[i]=0;}
4462    yw[n-1]=1;
4463    poly px; poly py;
4464    int te=1;
4465    i=1;
4466    while( te and i<=size(B0))
4467    {
4468      if((deg(B0[i],xw)==1) and (deg(B0[i])==1)){px=B0[i]; te=0;}
4469      i++;
4470    }
4471    i=1; te=1;
4472    while( te and i<=size(B0))
4473    {
4474      if((deg(B0[i],yw)==1) and (deg(B0[i])==1)){py=B0[i]; te=0;}
4475      i++;
4476    }
4477    N=px,py;
4478    te=indepparameters(N);
4479    if(te)
4480    {
4481      string("locus detected that the mover must avoid point (",N,") in order to obtain the correct locus");
4482      // eliminates segments of GG where N is contained in the basis
4483      list nGP;
4484      def GP=GG;
4485      ideal BP;
4486      for(j=1;j<=size(GP);j++)
4487      {
4488        te=1; k=1;
4489        BP=GP[j][2];
4490        while((te==1) and (k<=size(N)))
4491        {
4492          if(pdivi(N[k],BP)[1]!=0){te=0;}
4493          k++;
4494        }
4495        if(te==0){nGP[size(nGP)+1]=GP[j];}
4496      }
4497     }
4498    else{"Warning! Problem with more than one mover. Not able to solve it."; list L; return(L);}
4499  }
4500  def LL=locus0(nGP,DD);
4501  kill @RP; kill @P; kill @R;
4502  return(LL);
4503}
4504example
4505{
4506  "EXAMPLE:"; echo = 2;
4507  ring R=(0,a,b),(x,y),dp;
4508  short=0;
4509
4510  // Concoid
4511  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
4512  // System S96=
4513  S96;
4514  locus(grobcov(S96));
4515  kill R;
4516  // ********************************************
4517  ring R=(0,a,b),(x4,x3,x2,x1),dp;
4518  short=0;
4519  ideal S=(x1-3)^2+(x2-1)^2-9,
4520               (4-x2)*(x3-3)+(x1-3)*(x4-1),
4521               (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
4522               (4-x4)*a+(x3-3)*b+3*x4-4*x3,
4523               (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
4524  // System S=
4525  S;
4526  locus(grobcov(S));
4527  kill R;
4528  "********************************************";
4529
4530  ring R=(0,x,y),(x1,x2),dp;
4531  short=0;
4532  ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
4533                (x1 - 5)^2 + (x2 - 2)^2 - 4,
4534                -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
4535  // System S=
4536  S;
4537  locus(grobcov(S));
4538}
4539
4540//  locusdg(G):  Special routine for determining the locus of points
4541//                 of  geometrical constructions. Given a parametric ideal J with
4542//                 parameters (a_1,..a_m) and variables (x_1,..,xn),
4543//                 representing the system determining
4544//                 the locus of points (a_1,..,a_m) who verify certain
4545//                 properties, computing the grobcov G of
4546//                 J and applying to it locus, determines the different
4547//                 classes of locus components. The components can be
4548//                 Normal, Special, Accumulation point, Degenerate.
4549//                 The output are the components given in P-canonical form
4550//                 of at most 4 constructible sets: Normal, Special, Accumulation,
4551//                 Degenerate.
4552//                 The description of the algorithm and definitions is
4553//                 given in a forthcoming paper by Abanades, Botana, Montes, Recio.
4554//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
4555//                 The mover coordinates are to be placed as the last variables, and by default,
4556//                 its number is 2. If onw consider locus problems in higer dimensions, the number of
4557//                 mover coordinates (placed as the last variables) is to be given as an option.
4558//
4559//  input:      The output of locus(G);
4560//  output:
4561//          list, the canonical P-representation of the Normal and Non-Normal locus:
4562//               The Normal locus has two kind of components: Normal and Special.
4563//               The Non-normal locus has two kind of components: Accumulation and Degenerate.
4564//               Normal components: for each point in the component,
4565//               the number of solutions in the variables is finite, and
4566//               the solutions depend on the point in the component if the component is not 0-dimensional.
4567//               Special components:  for each point in the component,
4568//               the number of solutions in the variables is finite,
4569//               the component is not 0-dimensional, but the solutions do not depend on the
4570//               values of the parameters in the component.
4571//               Accumlation points: are 0-dimensional components for which it exist
4572//               an infinite number of solutions.
4573//               Degenerate components: are components of dimension greater than 0 for which
4574//               for every point in the component there exist infinite solutions.
4575//          The output components are given as
4576//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
4577//          The components are given in canonical P-representation of the subset.
4578//               If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
4579//               gives the depth of the component of the constructible set.
4580proc locusdg(list L)
4581"USAGE: locusdg(L)   The call must be locusdg(locus(grobcov(S))).
4582RETURN: The output is the list of the 'Relevant' components of the locus
4583           in Dynamic Geometry:(C1,..,C:m), where
4584                 C_i= ( p_i,(p_i1,..p_is_i), 'Relevant', level_i )
4585           The 'Relevant' components are the 'Normal' and 'Accumulation' components
4586           of the locus. (See help for locus).
4587
4588NOTE: It can only be called after computing the locus.
4589           Calling sequence:    locusdg(locus(grobcov(S)));
4590KEYWORDS: geometrical locus, locus, loci, dynamic geometry
4591EXAMPLE: locusdg; shows an example"
4592{
4593  list LL;
4594  int i;
4595  for(i=1;i<=size(L);i++)
4596  {
4597    if(typeof(L[i][3])=="string")
4598    {
4599      if((L[i][3]=="Normal") or (L[i][3]=="Accumulation")){L[i][3]="Relevant"; LL[size(LL)+1]=L[i];}
4600    }
4601  }
4602  return(LL);
4603}
4604example
4605{
4606  "EXAMPLE:"; echo = 2;
4607  ring R=(0,a,b),(x,y),dp;
4608  short=0;
4609
4610  // Concoid
4611  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
4612  // System S96=
4613  S96;
4614  locus(grobcov(S96));
4615  locusdg(locus(grobcov(S96)));
4616  kill R;
4617  "********************************************";
4618  ring R=(0,a,b),(x4,x3,x2,x1),dp;
4619  ideal S=(x1-3)^2+(x2-1)^2-9,
4620               (4-x2)*(x3-3)+(x1-3)*(x4-1),
4621               (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
4622               (4-x4)*a+(x3-3)*b+3*x4-4*x3,
4623               (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
4624  short=0;
4625  locus(grobcov(S));
4626  locusdg(locus(grobcov(S)));
4627  kill R;
4628  "********************************************";
4629
4630  ring R=(0,x,y),(x1,x2),dp;
4631  short=0;
4632  ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
4633                (x1 - 5)^2 + (x2 - 2)^2 - 4,
4634                -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
4635  locus(grobcov(S));
4636  locusdg(locus(grobcov(S)));
4637}
4638
4639// locusto: Transforms the output of locus to a string that
4640//      can be read by different computational systems.
4641// input:
4642//     list L: The output of locus
4643// output:
4644//     string s: The output of locus converted to a string readable by other programs
4645proc locusto(list L)
4646"USAGE: locusto(L);
4647          The argument must be the output of locus or locusdg or
4648          envelop or envelopdg.
4649          It transforms the output into a string in standard form
4650          readable in many languages (Geogebra).
4651RETURN: The locus in string standard form
4652NOTE: It can only be called after computing either
4653              - locus(grobcov(F))                -> locusto( locus(grobcov(F)) )
4654              - locusdg(locus(grobcov(F)))  -> locusto( locusdg(locus(grobcov(F))) )
4655              - envelop(F,C)                       -> locusto( envelop(F,C) )
4656              -envelopdg(envelop(F,C))      -> locusto( envelopdg(envelop(F,C)) )
4657KEYWORDS: geometrical locus, locus, loci.
4658EXAMPLE:  locusto; shows an example"
4659{
4660  int i; int j; int k;
4661  string s="["; string sf="]"; string st=s+sf;
4662  if(size(L)==0){return(st);}
4663  ideal p;
4664  ideal q;
4665  for(i=1;i<=size(L);i++)
4666  {
4667    s=string(s,"[[");
4668    for (j=1;j<=size(L[i][1]);j++)
4669    {
4670      s=string(s,L[i][1][j],",");
4671    }
4672    s[size(s)]="]";
4673    s=string(s,",[");
4674    for(j=1;j<=size(L[i][2]);j++)
4675    {
4676      s=string(s,"[");
4677      for(k=1;k<=size(L[i][2][j]);k++)
4678      {
4679        s=string(s,L[i][2][j][k],",");
4680      }
4681      s[size(s)]="]";
4682      s=string(s,",");
4683    }
4684    s[size(s)]="]";
4685    s=string(s,"]");
4686    if(size(L[i])>=3)
4687    {
4688      s=string(s,",[");
4689      if(typeof(L[i][3])=="string")
4690      {
4691        s=string(s,string(L[i][3]),"]]");
4692      }
4693      else
4694      {
4695        for(k=1;k<=size(L[i][3]);k++)
4696        {
4697          s=string(s,"[",L[i][3][k],"],");
4698        }
4699        s[size(s)]="]";
4700        s=string(s,"]");
4701      }
4702    }
4703    if(size(L[i])>=4)
4704    {
4705      s[size(s)]=",";
4706      s=string(s,string(L[i][4]),"],");
4707    }
4708    s[size(s)]="]";
4709    s=string(s,",");
4710  }
4711  s[size(s)]="]";
4712  return(s);
4713}
4714example
4715{
4716  "EXAMPLE:"; echo = 2;
4717  ring R=(0,x,y),(x1,y1),dp;
4718  short=0;
4719  ideal S=x1^2+y1^2-4,(y-2)*x1-x*y1+2*x,(x-x1)^2+(y-y1)^2-1;
4720  locusto(locus(grobcov(S)));
4721  locusto(locusdg(locus(grobcov(S))));
4722  kill R;
4723  "********************************************";
4724
4725  // 1. Take a fixed line l: x1-y1=0  and consider
4726  //    the family F of a lines parallel to l passing through the mover point M
4727  // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
4728  // 3. Compute the envelop of the family of lines.
4729
4730  ring R=(0,x,y),(x1,y1),lp;
4731  poly F=(y-y1)-(x-x1);
4732  ideal C=x1^2+y1^2-25;
4733  short=0;
4734
4735  // Curves Family F=
4736  F;
4737  // Conditions C=
4738  C;
4739
4740  locusto(envelop(F,C));
4741  locusto(envelopdg(envelop(F,C)));
4742  kill R;
4743  "********************************************";
4744
4745  // Steiner Deltoid
4746  // 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
4747  // 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
4748  // 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
4749  // 4. Obtain the envelop of the lines above.
4750
4751  ring R=(0,x,y),(x1,y1,x2,y2),lp;
4752  ideal C=(x1)^2+(y1)^2-1,
4753               x2+y2-1,
4754               x2-y2-x1+y1;
4755  matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
4756  poly F=det(M);
4757
4758  short=0;
4759
4760  // Curves Family F=
4761  F;
4762  // Conditions C=
4763  C;
4764
4765  locusto(envelop(F,C));
4766  locusto(envelopdg(envelop(F,C)));
4767}
4768
4769// Auxiliary routine
4770// locusdim
4771// input:
4772//    B:         ideal, a basis of a segment of the grobcov
4773//    dgdim: int, the dimension of the mover (for locus)
4774//                by default dgdim is equal to the number of variables
4775static proc locusdim(ideal B, list #)
4776{
4777  def R=basering;
4778  int dgdim;
4779  int nv=nvars(R);
4780  if (size(#)>0){dgdim=#[1];}
4781  else {dgdim=nv;}
4782  int d;
4783  list v;
4784  ideal vi;
4785  int i;
4786  for(i=1;i<=dgdim;i++)
4787  {
4788    v[size(v)+1]=varstr(nv-dgdim+i);
4789    vi[size(v)+1]=var(nv-dgdim+i);
4790  }
4791  ideal B0;
4792  for(i=1;i<=size(B);i++)
4793  {
4794    if(subset(variables(B[i]),vi))
4795    {
4796      B0[size(B0)+1]=B[i];
4797    }
4798  }
4799  def RR=ringlist(R);
4800  def RR0=RR;
4801  RR0[2]=v;
4802  def R0=ring(RR0);
4803  setring(R0);
4804  def B0r=imap(R,B0);
4805  B0r=std(B0r);
4806  d=dim(B0r);
4807  setring R;
4808  return(d);
4809}
4810
4811// // locusdgto: Transforms the output of locusdg to a string that
4812// //      can be read by different computational systems.
4813// // input:
4814// //     list L: The output of locus
4815// // output:
4816// //     string s: The output of locus converted to a string readable by other programs
4817// //                   Outputs only the relevant dynamical geometry components.
4818// //                   Without unnecessary parenthesis
4819// proc locusdgto(list LL)
4820// "USAGE: locusdgto(L);
4821//           The argument must be the output of locusdg of a parametrical ideal
4822//           It transforms the output into a string in standard form
4823//           readable in many languages (Geogebra).
4824// RETURN: The locusdg in string standard form
4825// NOTE: It can only be called after computing the locusdg(grobcov(F)) of the
4826//           parametrical ideal.
4827//           The basering R, must be of the form Q[a,b,..][x,y,..].
4828// KEYWORDS: geometrical locus, locus, loci.
4829// EXAMPLE:  locusdgto; shows an example"
4830// {
4831//   int i; int j; int k; int short0=short; int ojo;
4832//   int te=0;
4833//   short=0;
4834//   if(size(LL)==0){ojo=1; list L;}
4835//   else
4836//   {
4837//     def L=LL;
4838//   }
4839//   string s="["; string sf="]"; string st=s+sf;
4840//   if(size(L)==0){return(st);}
4841//   ideal p;
4842//   ideal q;
4843//   for(i=1;i<=size(L);i++)
4844//   {
4845//     if(L[i][3]=="Relevant")
4846//     {
4847//       s=string(s,"[[");
4848//       for (j=1;j<=size(L[i][1]);j++)
4849//       {
4850//         s=string(s,L[i][1][j],",");
4851//       }
4852//       s[size(s)]="]";
4853//       s=string(s,",[");
4854//       for(j=1;j<=size(L[i][2]);j++)
4855//       {
4856//         s=string(s,"[");
4857//         for(k=1;k<=size(L[i][2][j]);k++)
4858//         {
4859//           s=string(s,L[i][2][j][k],",");
4860//         }
4861//         s[size(s)]="]";
4862//         s=string(s,",");
4863//       }
4864//       s[size(s)]="]";
4865//       s=string(s,"]");
4866//       s[size(s)]="]";
4867//       s=string(s,",");
4868//     }
4869//   }
4870//   if(s=="["){s="[]";}
4871//   else{s[size(s)]="]";}
4872//   short=short0;
4873//   return(s);
4874// }
4875// example
4876// {"EXAMPLE:"; echo = 2;
4877//   ring R=(0,a,b),(x,y),dp;
4878//   short=0;
4879//   ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
4880//   "System="; S96; " ";
4881//   "locusdgto(locusdg(grobcov(S96)))=";
4882//   locusdgto(locusdg(grobcov(S96)));
4883// }
4884
4885static proc norspec(ideal F)
4886{
4887  def RR=basering;
4888  def Rx=ringlist(RR);
4889
4890   def Rx=ringlist(RR);
4891  def @P=ring(Rx[1]);
4892  list Lx;
4893  Lx[1]=0;
4894  Lx[2]=Rx[2]+Rx[1][2];
4895  Lx[3]=Rx[1][3];
4896  Lx[4]=Rx[1][4];
4897  Rx[1]=0;
4898  def D=ring(Rx);
4899  def @RP=D+@P;
4900  exportto(Top,@R);      // global ring Q[a][x]
4901  exportto(Top,@P);      // global ring Q[a]
4902  exportto(Top,@RP);     // global ring K[x,a] with product order
4903  setring(RR);
4904
4905}
4906
4907// envelop
4908// Input: n arguments
4909//   poly F: the polynomial defining the family of curves in ring R=0,(x,y),(x1,..,xn),lp;
4910//   ideal C=g1,..,g_{n-1}:  the set of constraints
4911// Output: the components of the envolvent;
4912proc envelop(poly F, ideal C, list #)
4913"USAGE: envelop(F,C);
4914          The first argument F must be the family of curves of which
4915          on want to compute the envelop.
4916          The second argument C must be the ideal of conditions
4917          over the variables, and should contain as polynomials
4918          as the number of variables -1.
4919RETURN: The components of the envelop with its taxonomy:
4920           The taxonomy distinguishes 'Normal',
4921          'Special', 'Accumulation', 'Degenerate' components.
4922          In the case of 'Special' components, it also
4923          outputs the antiimage of the component
4924          and an integer (0-1). If the integer is 0
4925          the component is not a curve of the family and is
4926          not considered as 'Relevant' by the envelopdg routine
4927          applied to it, but is considered as 'Relevant' if the integer is 1.
4928NOTE: grobcov is called internally.
4929          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
4930KEYWORDS: geometrical locus, locus, loci, envelop
4931EXAMPLE:  envelop; shows an example"
4932{
4933  int tes=0; int i;   int j;
4934  def R=basering;
4935  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
4936  setglobalrings();
4937  // Options
4938  list DD=#;
4939  int moverdim=nvars(R);
4940  int version=0;
4941  int nv=nvars(R);
4942  if(nv<4){version=1;}
4943  int comment=0;
4944  ideal Fm;
4945  for(i=1;i<=(size(DD) div 2);i++)
4946  {
4947    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
4948    if(DD[2*i-1]=="version"){version=DD[2*i];}
4949    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
4950    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
4951  }
4952  int n=nvars(R);
4953  list v;
4954  for(i=1;i<=n;i++){v[size(v)+1]=var(i);}
4955  def MF=jacob(F);
4956  def TMF=transpose(MF);
4957  def Mg=MF;
4958  def TMg=TMF;
4959  for(i=1;i<=n-1;i++)
4960  {
4961    Mg=jacob(C[i]);
4962    TMg=transpose(Mg);
4963    TMF=concat(TMF,TMg);
4964  }
4965  poly J=det(TMF);
4966  ideal S=ideal(F)+C+ideal(J);
4967  DD[size(DD)+1]="family";
4968  DD[size(DD)+1]=F;
4969  def G=grobcov(S,DD);
4970   def L=locus(G, DD);
4971  return(L);
4972}
4973example
4974{
4975  "EXAMPLE:"; echo = 2;
4976  // Steiner Deltoid
4977  // 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
4978  // 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
4979  // 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
4980  // 4. Obtain the envelop of the lines above.
4981
4982  ring R=(0,x,y),(x1,y1,x2,y2),lp;
4983  ideal C=(x1)^2+(y1)^2-1,
4984               x2+y2-1,
4985               x2-y2-x1+y1;
4986  matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
4987  poly F=det(M);
4988
4989  short=0;
4990
4991  // Curves Family F=
4992  F;
4993  // Conditions C=
4994  C;
4995
4996  def Env=envelop(F,C);
4997  Env;
4998}
4999
5000// envelopdg
5001// Input: list L: the output of envelop(poly F, ideal C, list #)
5002// Output: the relevant components of the envolvent in dynamic geometry;
5003proc envelopdg(list L)
5004"USAGE: envelopdg(L);
5005          The input list L must be the output of the call to
5006          the routine 'envolop' of the family of curves
5007RETURN: The relevant components of the envelop in Dynamic Geometry.
5008           'Normal' and 'Accumulation' components are always considered
5009           'Relevant'. 'Special' components of the envelop outputs
5010           three objects in its characterization: 'Special', the antiimage ideal,
5011           and the integer 0 or 1, that indicates that the given component is
5012           formed (1) or is not formed (0) by curves of the family. Only if yes,
5013           'envelopdg' considers the component as 'Relevant' .
5014NOTE: It must be called to the output of the 'envelop' routine.
5015          The basering R, must be of the form Q[a,b,..][x,y,..].
5016KEYWORDS: geometrical locus, locus, loci, envelop.
5017EXAMPLE:  envelop; shows an example"
5018{
5019  list LL;
5020  list Li;
5021  int i;
5022  for(i=1;i<=size(L);i++)
5023  {
5024    if(typeof(L[i][3])=="string")
5025    {
5026      if((L[i][3]=="Normal") or (L[i][3]=="Accumulation")){Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;}
5027    }
5028    else
5029    {
5030      if(typeof(L[i][3])=="list")
5031      {
5032        if(L[i][3][3]==1)
5033        {
5034          Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;
5035        }
5036      }
5037    }
5038  }
5039  return(LL);
5040}
5041example
5042{
5043  "EXAMPLE:"; echo = 2;
5044
5045  // 1. Take a fixed line l: x1-y1=0  and consider
5046  //    the family F of a lines parallel to l passing through the mover point M
5047  // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
5048  // 3. Compute the envelop of the family of lines.
5049
5050  ring R=(0,x,y),(x1,y1),lp;
5051  short=0;
5052  poly F=(y-y1)-(x-x1);
5053  ideal C=x1^2+y1^2-25;
5054  short=0;
5055
5056  // Curves Family F=
5057  F;
5058  // Conditions C=
5059  C;
5060
5061  envelop(F,C);
5062  envelopdg(envelop(F,C));
5063}
5064
5065//********************* End locus ****************************
5066
5068
5069// Input: L1,L2: lists of components with common top
5070// Output L: list of the union of L1 and L2.
5071static proc Add2ComWithCommonTop(list L1, list L2)
5072{
5073  int i; int j; ideal pij; list L; list Lp; list PR; int k;
5074  for(i=1;i<=size(L1[2]);i++)
5075  {
5076    for(j=1;j<=size(L2[2]);j++)
5077    {
5078      pij=std(L1[2][i]+L2[2][j]);
5079      PR=minGTZ(pij);
5080      for(k=1;k<=size(PR);k++)
5081      {
5082        Lp[size(Lp)+1]=PR[k];
5083      }
5084    }
5085  }
5086  for(i=1; i<=size(Lp);i++)
5087  {
5088    for(j=i+1;j<=size(Lp);j++)
5089    {
5090      if(idcontains(Lp[i],Lp[j])) {Lp=delete(Lp,j);}
5091    }
5092    for(j=1;j<i;j++)
5093    {
5094      if(idcontains(Lp[i],Lp[j])){Lp=delete(Lp,j); j=j-1; i=i-1;}
5095    }
5096  }
5097  L[1]=L1[1];
5098  L[2]=Lp;
5099  return(L);
5100}
5101
5102// Input: L list od P-rep of components to be added. L[i]=[p_i,[p_{i1},...p_{ir_i}]]
5103// Output: lists A,B,L
5104//       where no top in the lists are repeated
5105//       A: list of components with higher tops
5106//       B: list of components with lower tops
5107//       L1: A,B
5108static proc SepareAB(list L)
5109{
5110  int i;  int j; list Ln=L;
5111  for(i=1;i<=size(Ln);i++)
5112  {
5113    for(j=i+1;j<=size(Ln);j++)
5114    {
5115      if (equalideals(Ln[j][1],Ln[i][1]))
5116      {
5118        Ln=delete(Ln,j);
5119        j=j-1;
5120      }
5121    }
5122  }
5123  list A; list B; int clas; list T; list L1;
5124  for(i=1;i<=size(Ln);i++)
5125  {
5126    j=1;
5127    clas=0;
5128    while((clas==0) and  (j<=size(Ln)))
5129    {
5130      if(j!=i)
5131      {
5132        if(idcontains(Ln[i][1],Ln[j][1]) ) {B[size(B)+1]=Ln[i]; clas=1;}
5133      }
5134      j++;
5135    }
5136    if(clas==0) {A[size(A)+1]=Ln[i];}
5137  }
5138  L1=A; for(j=1;j<=size(B);j++){L1[size(L1)+1]=B[j];}
5139  T[1]=A; T[2]=B; T[3]=L1;
5140  return(T);
5141}
5142
5143// Input:
5144//  A1: list of high set of P-reps to be completed by the remaining P-reps
5145//  L1: the list A1, completed with the list of lower P-reps.
5146// Output:
5147//  A: list A1 completed with all possible parts of the remaining parts of L1 with the
5148//      condition of building locally closed sets.
5149static proc CompleteA(list A1,list L1)
5150{
5151  int modif; int i; int j; int k; int l;
5152  ideal pij; ideal pk; ideal pijkl; ideal pkl;
5153  int n; list nl; int te; int clas; list vvv; list AAA;
5154  list Lp; int m; ideal Pij;
5155  list A0;
5156  modif=1;
5157  list A=A1;
5158  while(modif==1)
5159  {
5160      modif=0;
5161      A0=A;
5162      for(i=1;i<=size(A);i++)
5163      {
5164          for(j=1;j<=size(A[i][2]); j++)
5165          {
5166              pij=A[i][2][j];
5167             for(k=1;k<=size(L1);k++)
5168             {
5169                 if(k!=i)
5170                 {
5171                     pk=L1[k][1];
5172                     if(idcontains(pij,pk)==1)
5173                     {
5174                         te=0;
5175                         kill nl;
5176                         list nl; l=1;
5177                         while((te==0) and (l<=size(L1[k][2])))
5178                         {
5179                              pkl=L1[k][2][l];
5180                              if((equalideals(pij,pkl)==1) or (idcontains(pij,pkl)==1)) {te=1;}
5181                              l++;
5182                         }   // end while ((te=0) and (l>...
5183                         //"T_te="; te; pause();
5184                         if(te==0)
5185                         {
5186                           modif=1;
5187                           kill Pij; ideal Pij=1;
5188                           for(l=1; l<=size(L1[k][2]);l++)
5189                           {
5190                             pkl=L1[k][2][l];
5191                             pijkl=std(pij+pkl);
5192                             Pij=intersect(Pij,pijkl);
5193                           }
5194                           kill Lp; list Lp;
5195                           Lp=minGTZ(Pij);
5196                           for(m=1;m<=size(Lp);m++)
5197                            {
5198                               nl[size(nl)+1]=Lp[m];
5199                            }
5200                            A[i][2]=delete(A[i][2], j);
5201                            for(n=1;n<=size(nl);n++)
5202                            {
5203                              A[i][2][size(A[i][2])+1]=nl[n];
5204                            }
5205                          } // end if(te==0)
5206                     } // end if(idcontains(pij,pk)==1)
5207                 }  // end if (k!=i)
5208             } // end for k
5209         }  // end for j
5210         kill vvv; list vvv;
5211         if(modif==1)
5212         // Select the maximal ideals of the set A[I][2][j]
5213         {
5214             kill nl; list nl;
5215             nl=selectminideals(A[i][2]);
5216             kill AAA; list AAA;
5217             for(j=1;j<=size(nl);j++)
5218             {
5219               AAA[size(AAA)+1]=A[i][2][nl[j]];
5220             }
5221             A[i][2]=AAA;
5222         } // end if(modif=1)
5223      }  // end for i
5224      modif=1-equallistsAall(A,A0);
5225  } // end while(modif==1)
5226  return(A);
5227}
5228
5229// Input:
5230//   A: list of the top P-reps of one level
5231//   B: list of remaining lower P-reps that have not yeen be possible to add
5232// Output:
5233//   Bn: list B where the elements that are completely included in A are removed and the parts that are
5234//         included in A also.
5235static proc ReduceB(list A,list B)
5236{
5237     int i; int j; list nl; list Bn; int te; int k; int elim;
5238     ideal pC; ideal pD; ideal pCD; ideal pBC; list nB; int j0;
5239     for(i=1;i<=size(B);i++)
5240     {
5241         j=1; te=0;
5242         while((te==0) and (j<=size(A)))
5243         {
5244             if(idcontains(B[i][1],A[j][1])){te=1; j0=j;}
5245             else{j++;}
5246         }
5247         pD=B[i][2][1];
5248         for(k=2;k<=size(B[i][2]);k++){pD=intersect(pD,B[i][2][k]);}
5249         pC=A[j0][2][1];
5250         for(k=2;k<=size(A[j0][2]);k++) {pC=intersect(pC,A[j0][2][k]);}
5251         pCD=std(pD+pC);
5252         pBC=std(B[i][1]+pC);
5253         elim=0;
5254         if(idcontains(pBC,pCD)){elim=1;}   // B=delfromlist(B,i);i=i-1;
5255         else
5256         {
5257              nB=Prep(pBC,pCD);
5258              if(equalideals(nB[1][1],ideal(1))==0)
5259              {
5260                  for(k=1;k<=size(nB);k++){Bn[size(Bn)+1]=nB[k];}
5261              }
5262         }
5263    }   // end for i
5264    return(Bn);
5265}
5266
5267// AddConsP: given a set of components of locally closed sets in P-representation, it builds the
5268//       canonical P-representation of the corresponding constructible set, of its union,
5269//       including levels it they are.
5272      Input L: list of components in P-rep to be added
5273      [  [[p_1,[p_11,..,p_1,r1]],..[p_k,[p_k1,..,p_kr_k]]  ]
5274RETURN:
5275     list of lists of levels of the different locally closed sets of
5276     the canonical P-rep of the constructible.
5277     [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
5278        [level_s,[ [Comp_s1,..Comp_sr_1] ]
5279     ]
5280     where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
5281NOTE:     Operates in a ring R=Q[u_1,..,u_m]
5282KEYWORDS: Constructible set, Canoncial form
5284{
5285  list LL; list A; list B; list L1; list T; int level=0; list h;
5286  LL=L; int i;
5287  while(size(LL)!=0)
5288  {
5289    level++;
5290    L1=SepareAB(LL);
5291    A=L1[1]; B=L1[2]; LL=L1[3];
5292    A=CompleteA(A,LL);
5293    for(i=1;i<=size(A);i++)
5294    {
5295      LL[i]=A[i];
5296    }
5297    h[1]=level; h[2]=A;
5298    T[size(T)+1]=h;
5299    LL=ReduceB(A,B);
5300  }
5301  return(T);
5302}
5303example
5304{
5305  "EXAMPLE:"; echo = 2;
5306  if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
5307  ring R=0,(x,y,z),lp;
5308  short=0;
5309
5310  ideal P1=x;
5311  ideal Q1=x,y;
5312  ideal P2=y;
5313  ideal Q2=y,z;
5314
5315  list L=list(Prep(P1,Q1)[1],Prep(P2,Q2)[1]);
5316  L;
5318}
5319
5320// AddCons:  given a set of  locally closed sets by pairs of ideal, it builds the
5321//       canonical P-representation of the corresponding constructible set, of its union,
5322//       including levels it they are.
5323//       It makes a call to AddConsP after transforming the input.
5324// Input list of lists of pairs of ideals representing locally colsed sets:
5325//     L=  [ [E1,N1], .. , [E_s,N_s] ]
5326// Output: The canonical frepresentation of the constructible set union of the V(E_i) \ V(N_i)
5327//     T=[  [level_1,[ [p_1,[p_11,..,p_1s_1]],.., [p_k,[p_k1,..,p_ks_k]] ]],, .. , [level_r,[..       ]]  ]
5330      Calls internally AddConsP and allows a different input.
5331      Input L: list of pairs of of ideals [ [P_1,Q_1], .., [Pr,Qr] ]
5332      representing a set of locally closed setsV(P_i) \ V(Q_i)
5334RETURN:
5335      list of lists of levels of the different locally closed sets of
5336      the canonical P-rep of the constructible.
5337      [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
5338         [level_s,[ [Comp_s1,..Comp_sr_1] ]
5339      ]
5340      where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
5341NOTE: Operates in a ring R=Q[u_1,..,u_m]
5342KEYWORDS: Constructible set, Canoncial form
5344{
5345  int i; list H; list P; int j;
5346  for(i=1;i<=size(L);i++)
5347  {
5348    P=Prep(L[i][1],L[i][2]);
5349    for(j=1;j<=size(P);j++)
5350    {
5351      H[size(H)+1]=P[j];
5352    }
5353  }
5355}
5356example
5357{
5358  "EXAMPLE:"; echo = 2;
5359  if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
5360  ring R=0,(x,y,z),lp;
5361  short=0;
5362
5363  ideal P1=x;
5364  ideal Q1=x,y;
5365  ideal P2=y;
5366  ideal Q2=y,z;
5367
5368  list L=list(list(P1,Q1), list(P2,Q2) );
5369  L;
5370
5372}
5373
5374// AddLocus: auxilliary routine for locus0 that computes the components of the constructible:
5375// Input:  the list of locally closed sets to be added, each with its type as third argument
5376//     L=[ [LC[11],,,LC[1k_1],  .., [LC[r1],,,LC[rk_r] ] where
5377//            LC[1]=[p1,[p11,..,p1k],typ]
5378// Output:  the list of components of the constructible union of the L, with the type of the corresponding top
5379//               and the level of the constructible
5380//     L4= [[v1,p1,[p11,..,p1l],typ_1,level]_1 ,.. [vs,ps,[ps1,..,psl],typ_s,level_s]
5382{
5383//  int te0=0;
5384//  def RR=basering;
5385//  if(defined(@P)){te0=1;  def Rx=@R;  kill @P; setring RR;}
5386  list L1; int i; int j;  list L2; list L3;
5387  list l1; list l2;
5388  intvec v;
5389  for(i=1; i<=size(L); i++)
5390  {
5391    for(j=1;j<=size(L[i]);j++)
5392    {
5393      l1[1]=L[i][j][1];
5394      l1[2]=L[i][j][2];
5395      l2[1]=l1[1];
5396      if(size(L[i][j])>2){l2[3]=L[i][j][3];}
5397      v[1]=i; v[2]=j;
5398      l2[2]=v;
5399      L1[size(L1)+1]=l1;
5400      L2[size(L2)+1]=l2;
5401    }
5402  }
5404  list L4; int level;
5405  ideal p1; ideal pp1; int t; int k; int k0; string typ; list l4;
5406  for(i=1;i<=size(L3);i++)
5407  {
5408    level=L3[i][1];
5409    for(j=1;j<=size(L3[i][2]);j++)
5410    {
5411      p1=L3[i][2][j][1];
5412      t=1; k=1;
5413      while((t==1) and (k<=size(L2)))
5414      {
5415        pp1=L2[k][1];
5416        if(equalideals(p1,pp1)){t=0; k0=k;}
5417        k++;
5418      }
5419      if(t==0)
5420      {
5421        v=L2[k0][2];
5422      }