# source:git/Singular/LIB/grobcov.lib@1e1ec4

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