# source:git/Singular/LIB/grobcov.lib@40df4d

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