Changeset 2203cf in git


Ignore:
Timestamp:
Oct 6, 2016, 3:17:04 PM (8 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '4fa496c5b814dbde0a905c54f0796301d03f6dc9')
Children:
b61d2e9bb22594df3a38d0289c0394fb3c068301
Parents:
17510a828cdc33aa541d40a9bbfb825b72812521
Message:
grobcov.lib: new version
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/grobcov.lib

    r17510a r2203cf  
    11//
    2 version="version grobcov.lib 4.0.2.0 Jul_2015 "; // $Id$
    3            // version L;  July_2015;
     2version="version grobcov.lib 4.0.3.3 Oct_2016 "; // $Id$
     3           // version M;  October_2016;
    44category="General purpose";
    55info="
    6 LIBRARY:  grobcov.lib   Groebner Cover for parametric ideals.
     6LIBRARY:  grobcovM.lib  Groebner Cover for parametric ideals.
    77
    88          Comprehensive Groebner Systems, Groebner Cover, Canonical Forms,
     
    1616          (http://www-ma2.upc.edu/~montes/).
    1717
    18 AUTHORS:  Antonio Montes , Hans Schoenemann.
     18IMPORTANT: The book,  not yet published:
     19
     20             A. Montes.
     21             \"Discussing Parametric Polynomial Systems: The Groebner Cover\"
     22
     23             can be used as a user manual of all the routines included in this library.
     24             It defines and proves all the theoretic results used here, and
     25             shows examples of all the routines.
     26
     27             I expect that it will be published soon.
     28
     29AUTHORS:  Antonio Montes (Universitat Politecnica de Catalunya),
     30          Hans Schoenemann (Technische Universitaet Kaiserslautern).
    1931
    2032OVERVIEW:
    2133            In 2010, the library was designed to contain Montes-Wibmer's
    22             algorithms for compute the canonical Groebner Cover of a
     34            algorithms for computing the canonical Groebner Cover of a
    2335            parametric ideal as described in the paper:
    2436
     
    2739            Journal of Symbolic Computation 45 (2010) 1391-1425.
    2840
    29             The central routine is grobcov. Given a parametric
    30             ideal, grobcov outputs its Canonical Groebner Cover, consisting
     41            The central routine is grobcov. Given a parametri ideal,
     42            grobcov outputs its Canonical Groebner Cover, consisting
    3143            of a set of pairs of (basis, segment). The basis (after
    3244            normalization) is the reduced Groebner basis for each point
     
    4557            reduced Comprehensive Groebner System with constant lpp.
    4658            For this purpose, in this library, the implemented algorithm is
    47             Kapur-Sun-Wang algorithm, because it is the most efficient
    48             algorithm known for this purpose.
     59            Kapur-Sun-Wang algorithm, because it is actually the most
     60            efficient algorithm known for this purpose.
    4961
    5062            D. Kapur, Y. Sun, and D.K. Wang.
     
    5971            union of locally closed sets. It is used in the new version for the
    6072            computation of loci and envelops. It determines the canonical locally closed
    61             level sets of a constructuble. They will be described in a forthcoming paper:
     73            level sets of a constructuble. It is described in:
    6274
    6375             J.M. Brunat, A. Montes,
    6476            \"Computing the canonical representation of constructible sets\".
    65             Submited to Mathematics in Computer Science. July 2015.
     77            Math.  Comput. Sci. (2016) 19: 165-178.
    6678
    6779            A new set of routines (locus, locusdg, locusto) has been included to
     
    7385            Computer-Aided Design 56 (2014) 22-33.
    7486
    75             Recently also routines for computing the envelop of a family
    76             of curves (enverlop, envelopdg), to be used in Dynamic Geometry,
    77             has been included and will be described in a forthcoming paper:
    78 
    79              Abanades, Botana, Montes, Recio:
    80              \''Envelops in Dynamic Geometry using the Groebner cover\''.
    81 
    82 
    83             This version was finished on 31/07/2015
     87            Recently also routines for computing the generalized envelop of a family
     88            of hyper-surfaces (envelop), to be used in Dynamic Geometry,
     89            has been included and is described in the book (not yet published
     90
     91             A. Montes.
     92             \"Discussing Parametric Polynomial Systems: The Groebner Cover\"
     93
     94            This version was finished on 31/010/2016
     95
    8496
    8597NOTATIONS: All given and determined polynomials and ideals are in the
     
    94106
    95107PROCEDURES:
     108setglobalrings();  Generates the global rings @R, @P and @RP that are
     109              respectively the rings Q[a][x], Q[a], Q[x,a].
     110              It is called inside each of the fundamental routines of the
     111              library: grobcov, cgsdr, locus, locusdg and killed before
     112              the output.
     113
     114              If the user want to call the routines Prep or Crep on a
     115              parametric ideal Q[a][x] for locally closed sets on Q[a],
     116              setglobalrings must be called previously.
     117              In the actual version, after the call of setglobalrings on
     118              Q[a][x], the public names of the defined ideals generated by
     119              setglobalrings are Grobcov::@R, Grobcov::@P,  Grobcov::@RP.
     120
    96121grobcov(F);  Is the basic routine giving the canonical
    97122             Groebner Cover of the parametric ideal F.
     
    103128             reduced Comprehensive Groebner System that
    104129             is used in grobcov, that can also be used
    105              independently if only the CGS is required.
     130             independently if only a CGS is required.
    106131             It is a more efficient routine than buildtree
    107132             (the own routine of 2010 that is no more used).
    108              Now, KSW algorithm is used.
     133             Now, Kapur-Sun-Wang (KSW) algorithm is used.
    109134
    110135pdivi(f,F); Performs a pseudodivision of a parametric polynomial
    111136             by a parametric ideal.
    112137
     138
    113139pnormalf(f,E,N);   Reduces a parametric polynomial f over V(E) \ V(N)
    114              ( E is the null ideal and N the non-null ideal )
     140             (E is the null ideal and N the non-null ideal )
    115141             over the parameters.
    116142
    117143Crep(N,M);  Computes the canonical C-representation of V(N) \ V(M).
     144              If the user want to call the routines Crep on a parametric
     145              ideal Q[a][x] for locally closed sets on Q[a], setglobalrings()
     146              must be called previously.
     147              In the actual version, after the call of setglobalrings on
     148              Q[a][x], the public names of the defined ideals generated by
     149              setglobalrings are Grobcov::@R, Grobcov::@P,  Grobcov::@RP.
    118150
    119151Prep(N,M);  Computes the canonical P-representation of V(N) \ V(M).
     152              If the user want to call the routines Prep on a parametric
     153              ideal Q[a][x] for locally closed sets on Q[a], setglobalrings()
     154              must be called previously.
     155              In the actual version, after the call of setglobalrings on
     156              Q[a][x], the public names of the defined ideals generated by
     157              setglobalrings are Grobcov::@R, Grobcov::@P,  Grobcov::@RP.
    120158
    121159PtoCrep(L)  Starting from the canonical Prep of a locally closed set
     
    128166             full representation of the bases. With the default
    129167             option ('ext',0) only the generic representation of
    130              the bases are computed, and one can obtain the full
     168             the bases is computed, and one can obtain the full
    131169             representation using extend.
    132170
    133 ConsLevels(L); Given a list L of locally closed sets, it returns the canonical levels
    134              of the constructible set of the union of them, as well as the levels
    135              of the complement. It is described in the paper
     171ConsLevels(L); Given a list L of locally closed sets, it returns the closures of the canonical
     172             levels of the constructible set and its complements of the union of them.
     173             It is described in the paper
    136174
    137175             J.M. Brunat, A. Montes,
    138176            \"Computing the canonical representation of constructible sets\".
    139             Submited to Mathematics in Computer Science. July 2015.
     177            Math.  Comput. Sci. (2016) 19: 165-178.
     178
     179ConsLevelsToLevels(L): Transforms the output of ConsLevels into the proper
     180             Levels of the constructible set.
    140181
    141182locus(G);    Special routine for determining the geometrical locus of points
     
    151192             Computer-Aided Design 56 (2014) 22-33.
    152193
    153              The components can be 'Normal', 'Special', 'Accumulation', 'Degenerate'.
    154              The output are the components is given in P-canonical form
     194             The components can be \"Normal\", \"Special\", \"Accumulation\", \"Degenerate\".
     195             The output are the components given in P-canonical form
    155196             It also detects automatically a possible point that is to be
    156              avoided by the mover, whose coordinates must be the last two
     197             avoided by the mover, whose coordinates must be the last
    157198             coordinates in the definition of the ring. If such a point is detected,
    158199             then it eliminates the segments of the grobcov depending on the
    159200             point that is to be avoided.
    160201
    161 locusdg(G);  Is a special routine that determines the  'Relevant' components
     202locusdg(G);  Is a special routine that determines the  \"Relevant\" components
    162203             of the locus in dynamic geometry. It is to be called to the output of locus
    163204             and selects from it the useful components.
    164205
    165 envelop(F,C); Special routine for determining the envelop of a family of curves
    166              F in Q[x,y][x_1,..xn] depending on a ideal of constraints C in Q[x_1,..,xn].
     206envelop(F,C); Special routine for determining the envelop of a family of hyper-surfaces
     207             F in Q[x1,..,xn][t1,..,tm] depending on a ideal of constraints C in Q[t1,..,tm].
    167208             It detemines the different components as well as its type:
    168              'Normal', 'Special', 'Accumulation', 'Degenerate'. And
    169              it also classifies the Special components, determining the
     209             \"Normal\", \"Special\", \"Accumulation\", \"Degenerate\". And
     210             it also classifies the \"Special\" components, determining the
    170211             zero dimensional antiimage of the component and verifying if
    171212             the component is a special curve of the family or not.
     
    173214             to obtain the complete result.
    174215             The taxonomy that it provides, as well as the algorithms involved
    175              will be described in a forthcoming paper:
    176 
    177              Abanades, Botana, Montes, Recio:
    178              \''Envelops in Dynamic Geometry using the Gr\"obner cover\''.
    179 
    180 envelopdg(ev); Is a special routine to determine the 'Relevant' components
    181              of the envelop of a family of curves to be used in Dynamic Geometry.
    182              It must be called to the output of envelop(F,C).
    183 
    184 locusto(L); Transforms the output of locus, locusdg, envelop and envelopdg
     216             are described in the book: (not yet published)
     217
     218             A. Montes.
     219             \"Discussing Parametric Polynomial Systems: The Groebner Cover\"
     220
     221locusto(L); Transforms the output of locus, locusdg, envelop
    185222             into a string that can be reed from different computational systems.
    186223
     224AssocTanToEnv(F,C,E); Having computed an envelop component E
     225             of a family of hyper-surfaces F, with constraints C,
     226             it returns the parameter values of the unique
     227             associated tangent hyper-surface of the family
     228             passing at one point of the envelop component E.
     229
     230FamElemsAtEnvCompPoints(F,C,E) Having computed the an envelop
     231             component E of a family of hyper-surfaces F, with constraints C,
     232             it returns the parameter values of all the hyper-surfaes of the family
     233             passing at one point of the envelop component E.
     234
     235discrim(f,x); Determines the factorized discriminant of a degree 2 polynomial
     236             in the variable x. The polynomial can be defined on any ring
     237             where x is a variable. The polynomial f can depend on parameters and variables
     238
     239WLemma(F,A); Given an ideal F in K[a][x] and a priime ideal A in K[a],
     240             it returns the list (lpp,B,S)  were B is the reduced Groebner basis
     241             of the specialized F over the segment computed in P-representation
     242             (or optionally in C-representation). The basis is given by I-regular
     243             functions.
    187244
    188245SEE ALSO: compregb_lib
     
    194251// ************ Begin of the grobcov library *********************
    195252
     253// Development of the library:
    196254// Library grobcov.lib
    197255// (Groebner Cover):
     
    200258// Uses buildtree for cgsdr
    201259// Final data: 3-7-2008
    202 // Release 2: (private)
     260// Release 2: (prived)
    203261// Initial data: 6-9-2009
    204262// Last version using buildtree for cgsdr
    205263// Final data: 25-10-2011
    206 // Release B: (private)
     264// Release B: (prived)
    207265// Initial data: 1-7-2012
    208266// Uses both buildtree and KSW for cgsdr
     
    213271// Final data: 21-11-2013
    214272// Release L: (public)
    215 // New routine ConsLevels: 10-7-2015
    216 // Updated locus: 10-7-2015 (uses Conslevels)
    217 // New routines for computing the envelop of a family of curves: 22-1-2015
    218 // Final data: 22-7-2015
     273// New routine ConsLevels: 25-1-2016
     274// Updated locus: 10-7-2015 (uses ConsLevels)
     275// Release M: (public)
     276// New routines for computing the envelop of a family of
     277//    hyper-surfaces and associated questions: 22-4-2016: 20-9-2016
     278// New routine WLemma for computing the result of
     279//    Wibmer's Lemma:  19-9-2016
     280// Final version October 2016
    219281
    220282//*************Auxiliary routines**************
     
    314376}
    315377
    316 
    317378// equalideals
    318379// Input: ideals F and G;
     
    325386  while ((i<=size(F)) and (t))
    326387  {
    327     if (F[i]!=G[i]){t=0;}
     388      if (F[i]!=G[i]){t=0;}
    328389    i++;
    329390  }
     
    372433  return(t);
    373434}
    374 
    375435
    376436// selectminideals
     
    404464  return(P);
    405465}
    406 
    407466
    408467// Auxiliary routine
     
    749808//}
    750809
    751 static proc setglobalrings()
    752 // "USAGE:   setglobalrings();
    753 //           No arguments
    754 // RETURN: After its call the rings Grobcov::@R=Q[a][x], Grobcov::@P=Q[a],
    755 //           Grobcov::@RP=Q[x,a] are defined as global variables.
    756 //           (a=parameters, x=variables)
    757 // NOTE: It is called internally by many basic routines of the
    758 //           library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
    759 //           envelop, envelopdg, and killed before the output.
    760 //           The user does not need to call it, except when it is interested
    761 //           in using some internal routine of the library that
    762 //           uses these rings.
    763 //           The basering R, must be of the form Q[a][x], (a=parameters,
    764 //           x=variables), and should be defined previously.
    765 // KEYWORDS: ring, rings
    766 // EXAMPLE:  setglobalrings; shows an example"
     810proc setglobalrings()
     811"USAGE:   setglobalrings();
     812          No arguments
     813RETURN: After its call the rings Grobcov::@R=Q[a][x], Grobcov::@P=Q[a],
     814          Grobcov::@RP=Q[x,a] are defined as global variables.
     815          (a=parameters, x=variables)
     816NOTE: It is called internally by many basic routines of the
     817          library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
     818          envelop, WLemma, and killed before the output.
     819          The user does not need to call it, except when it is interested
     820          in using some internal routine of the library that
     821          uses these rings.
     822          The basering R, must be of the form Q[a][x], (a=parameters,
     823          x=variables), and should be defined previously.
     824          When calling Prep or Crep of a locally closed set on a ring
     825          R[a][x], if setglobalrings() is called previously, then it is
     826          assumed that computing the P- or C-representation
     827          is considered for locally closed sets on Q[a].
     828KEYWORDS: ring; rings
     829EXAMPLE:  setglobalrings; shows an example"
    767830{
    768831  if (defined(@P))
     
    787850  setring(RR);
    788851};
    789 // example
    790 // {
    791 //   "EXAMPLE:"; echo = 2;
    792 //   ring R=(0,a,b),(x,y,z),dp;
    793 //   setglobalrings();
    794 //   " ";"R=";R;
    795 //   " ";"Grobcov::@R=";Grobcov::@R;
    796 //   " ";"Grobcov::@P=";Grobcov::@P;
    797 //   " ";"Grobcov::@RP=";Grobcov::@RP;
    798 //  " "; "ringlist(Grobcov::@R)=";  ringlist(Grobcov::@R);
    799 //  " "; "ringlist(Grobcov::@P)=";  ringlist(Grobcov::@P);
    800 //  " "; "ringlist(Grobcov::@RP)=";  ringlist(Grobcov::@RP);
    801 // }
     852example
     853{
     854  "EXAMPLE:"; echo = 2;
     855  ring R=(0,a,b),(x,y,z),dp;
     856  setglobalrings();
     857  R;
     858
     859  Grobcov::@R;
     860
     861  Grobcov::@P;
     862
     863  Grobcov::@RP;
     864
     865  ringlist(Grobcov::@R);
     866
     867  ringlist(Grobcov::@P);
     868
     869 ringlist(Grobcov::@RP);
     870}
    802871
    803872// cld : clears denominators of an ideal and normalizes to content 1
     
    859928          list (r,q,m) such that m*f=r+sum(q.G), and no lpp of a divisor
    860929          divides a pp of r.
    861 KEYWORDS: division, reduce
     930KEYWORDS: division; reduce
    862931EXAMPLE:  pdivi; shows an example"
    863932{
     
    919988  "EXAMPLE:"; echo = 2;
    920989  ring R=(0,a,b,c),(x,y),dp;
     990  // Divisor=";
    921991  poly f=(ab-ac)*xy+(ab)*x+(5c);
    922   // Divisor=";
    923   f;
     992  // Dividends=";
    924993  ideal F=ax+b,cy+a;
    925   // Dividends=";
    926   F;
    927   def r=pdivi(f,F);
    928994  // (Remainder, quotients, factor)=";
    929   r;
     995  def r=pdivi(f,F); r;
    930996  // Verifying the division:
    931997  r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2]+r[1]);
     
    10911157NOTE: Should be called from ring Q[a][x].
    10921158          Ideals E and N must be given by polynomials in Q[a].
    1093 KEYWORDS: division, pdivi, reduce
     1159KEYWORDS: division; pdivi; reduce
    10941160EXAMPLE: pnormalf; shows an example"
    10951161{
     
    11141180example
    11151181{
    1116   "EXAMPLE:"; echo = 2;
     1182"EXAMPLE:"; echo = 2;
    11171183  ring R=(0,a,b,c),(x,y),dp;
    11181184  short=0;
     
    11201186  ideal E=(c-1);
    11211187  ideal N=a-b;
    1122 
    11231188  pnormalf(f,E,N);
    11241189}
     
    13901455//    the list (a,b) of the canonical ideals
    13911456//    the Crep of V(N) \ V(M)
    1392 // Assumed to be called in the ring @R
    1393 // Works on the ring @P
     1457// Assumed to be called in the ring @R or the ring @P or a ring ring Q[a]
    13941458proc Crep(ideal N, ideal M)
    13951459 "USAGE:  Crep(N,M);
    1396  Input: ideal N (null ideal) (not necessarily radical nor maximal)
    1397            ideal M (hole ideal) (not necessarily containing N)
     1460         ideal N (null ideal) (not necessarily radical nor maximal)
     1461          ideal M (hole ideal) (not necessarily containing N)
    13981462 RETURN:  The canonical C-representation of the locally closed set.
    13991463           [ P,Q ], a pair of radical ideals with P included in Q,
    14001464           representing the set V(P) \ V(Q) = V(N) \ V(M)
    14011465 NOTE:   Operates in a ring R=Q[a] (a=parameters)
    1402  KEYWORDS: locally closed set, canoncial form
     1466           If the user want to call the routine Crep on a parametric
     1467           ideal Q[a][x] for locally closed sets on Q[a], setglobalrings()
     1468           must be called previously.
     1469
     1470 KEYWORDS: locally closed set; canoncial form
    14031471 EXAMPLE:  Crep; shows an example"
     1472{
     1473  int te;
     1474  def RR=basering;
     1475  if(defined(@P)){te=1; setring(@P); ideal Np=imap(RR,N); ideal Mp=imap(RR,M);}
     1476  else {te=0; def Np=N; def Mp=M;}
     1477  def La=Crep0(Np,Mp);
     1478  if(size(La)==0)
     1479  {
     1480    if(te==1) {setring(RR); list LL;}
     1481    if(te==0){list LL;}
     1482    return(LL);
     1483  }
     1484  else
     1485  {
     1486    if(te==1) {setring(RR); def LL=imap(@P,La);}
     1487    if(te==0){def LL=La;}
     1488  return(LL);
     1489  }
     1490}
     1491example
     1492{
     1493  "EXAMPLE:"; echo = 2;
     1494  short=0;
     1495  ring R=0,(x,y,z),lp;
     1496  ideal E=x*(x^2+y^2+z^2-25)*(x+y);
     1497  ideal N=x*(x-3)*(x+y),(y-4)*(x+y);
     1498  Crep(E,N);
     1499}
     1500
     1501// Crep0
     1502// Computes the C-representation of V(N) \ V(M).
     1503// input:
     1504//    ideal N (null ideal) (not necessarily radical nor maximal)
     1505//    ideal M (hole ideal) (not necessarily containing N)
     1506// output:
     1507//    the list (a,b) of the canonical ideals
     1508//    the Crep0 of V(N) \ V(M)
     1509// Assumed to be called in a ring Q[x] (example @P)
     1510static proc Crep0(ideal N, ideal M)
    14041511{
    14051512  list l;
     
    14151522  ideal P=ideal(1);
    14161523  L=minGTZ(Np);
     1524  //"T_Np="; Np;
     1525  //"T_minGTZ(Np)="; L;
    14171526  for(i=1;i<=size(L);i++)
    14181527  {
     1528    L[i]=std(L[i]);
    14191529    if(idcontains(L[i],Q)==0)
    14201530    {
     
    14241534  P=std(P);
    14251535  Q=std(radical(Q+P));
     1536  if(equalideals(P,Q)){return(l);}
    14261537  list T=P,Q;
    14271538  return(T);
    14281539}
     1540
     1541// Prep
     1542// Computes the P-representation of V(N) \ V(M).
     1543// input:
     1544//    ideal N (null ideal) (not necessarily radical nor maximal)
     1545//    ideal M (hole ideal) (not necessarily containing N)
     1546// output:
     1547//    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
     1548//    the Prep of V(N) \ V(M)
     1549// Assumed to be called in the ring @R or the ring @P or a ring ring Q[a]
     1550proc Prep(ideal N, ideal M)
     1551 "USAGE:  Prep(N,M);
     1552          ideal N (null ideal) (not necessarily radical nor maximal)
     1553          ideal M (hole ideal) (not necessarily containing N)
     1554 RETURN: The canonical P-representation of the locally closed set V(N) \ V(M)
     1555           Output: [ Comp_1, .. , Comp_s ] where
     1556              Comp_i=[p_i,[p_i1,..,p_is_i]]
     1557 NOTE:   To be called from the ring @R or @P or a ring Q[a]
     1558           If the user want to call the routine Prep on a parametric
     1559           ideal Q[a][x] for locally closed sets on Q[a], setglobalrings()
     1560           must be called previously.
     1561 KEYWORDS: Locally closed set; Canoncial form
     1562 EXAMPLE:  Prep; shows an example"
     1563{
     1564  int te;
     1565  def RR=basering;
     1566  if(defined(@P)){te=1; setring(@P); ideal Np=imap(RR,N); ideal Mp=imap(RR,M);}
     1567  else {te=0; def Np=N; def Mp=M;}
     1568  def La=Prep0(Np,Mp);
     1569  if(te==1) {setring(RR); def LL=imap(@P,La);}
     1570  if(te==0){def LL=La;}
     1571  return(LL);
     1572}
    14291573example
    14301574{
    14311575  "EXAMPLE:"; echo = 2;
    1432   if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
     1576  short=0;
     1577  ring R=0,(x,y,z),lp;
     1578  ideal E=x*(x^2+y^2+z^2-25)*(x+y);
     1579  ideal N=x*(x-3)*(x+y),(y-4)*(x+y);
     1580  Prep(E,N);
     1581}
     1582
     1583// Prep0
     1584// Computes the P-representation of V(N) \ V(M).
     1585// input:
     1586//    ideal N (null ideal) (not necessarily radical nor maximal)
     1587//    ideal M (hole ideal) (not necessarily containing N)
     1588// output:
     1589//    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
     1590//    the Prep of V(N) \ V(M)
     1591// Assumed to be called in a ring Q[x] (example @P)
     1592static proc Prep0(ideal N, ideal M)
     1593{
     1594  int te;
     1595  if (N[1]==1)
     1596  {
     1597    return(list(list(ideal(1),list(ideal(1)))));
     1598  }
     1599  int i; int j; list L0;
     1600  list Ni=minGTZ(N);
     1601  list prep;
     1602  for(j=1;j<=size(Ni);j++)
     1603  {
     1604    option(redSB);
     1605    Ni[j]=std(Ni[j]);
     1606  }
     1607  list Mij;
     1608  for (i=1;i<=size(Ni);i++)
     1609  {
     1610    Mij=minGTZ(Ni[i]+M);
     1611    for(j=1;j<=size(Mij);j++)
     1612    {
     1613      option(redSB);
     1614      Mij[j]=std(Mij[j]);
     1615    }
     1616    if ((size(Mij)==1) and (equalideals(Ni[i],Mij[1])==1)){;}
     1617    else
     1618    {
     1619        prep[size(prep)+1]=list(Ni[i],Mij);
     1620    }
     1621  }
     1622  //"T_before="; prep;
     1623  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
     1624  //"T_Prep="; prep;
     1625  //def Lout=CompleteA(prep,prep);
     1626  //"T_Lout="; Lout;
     1627  return(prep);
     1628}
     1629
     1630// PtoCrep
     1631// Computes the C-representation from the P-representation.
     1632// input:
     1633//    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
     1634//         the P-representation of V(N) \ V(M)
     1635// output:
     1636//    list (ideal ida, ideal idb)
     1637//    the C-representaion of V(N) \ V(M) = V(ida) \ V(idb)
     1638// Assumed to be called in the ring @R of the ring @P or a ring Q[a]
     1639proc PtoCrep(list L)
     1640"USAGE: PtoCrep(L)
     1641          L: list  [ Comp_1, .. , Comp_s ] where
     1642          Comp_i=[p_i,[p_i1,..,p_is_i] ], is
     1643          the P-representation of a locally closed set V(N) \ V(M)
     1644 RETURN:  The canonical C-representation of the locally closed set
     1645          [ P,Q ], a pair of radical ideals with P included in Q,
     1646          representing the set V(P) \ V(Q)
     1647 NOTE: To be called from the ring @R or @P or a ring Q[a]
     1648 KEYWORDS: locally closed set; canoncial form
     1649 EXAMPLE:  PtoCrep; shows an example"
     1650{
     1651  int te;
     1652  def RR=basering;
     1653  if(defined(@P)){te=1; setring(@P); list Lp=imap(RR,L);}
     1654  else {te=0; def Lp=L;}
     1655  def La=PtoCrep0(Lp);
     1656  if(te==1) {setring(RR); def LL=imap(@P,La);}
     1657  if(te==0){def LL=La;}
     1658  return(LL);
     1659}
     1660example
     1661{
     1662 "EXAMPLE:"; echo = 2;
    14331663  ring R=0,(x,y,z),lp;
    14341664  short=0;
     
    14431673}
    14441674
    1445 // Prep
    1446 // Computes the P-representation of V(N) \ V(M).
    1447 // input:
    1448 //    ideal N (null ideal) (not necessarily radical nor maximal)
    1449 //    ideal M (hole ideal) (not necessarily containing N)
    1450 // output:
    1451 //    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
    1452 //    the Prep of V(N) \ V(M)
    1453 // Assumed to be called in the ring @R
    1454 // Works on the ring @P
    1455 proc Prep(ideal N, ideal M)
    1456  "USAGE:  Prep(N,M);
    1457  Input: ideal N (null ideal) (not necessarily radical nor maximal)
    1458            ideal M (hole ideal) (not necessarily containing N)
    1459  RETURN: The canonical P-representation of the locally closed set V(N) \ V(M)
    1460            Output: [ Comp_1, .. , Comp_s ] where
    1461               Comp_i=[p_i,[p_i1,..,p_is_i]]
    1462  NOTE:   Operates in a ring R=Q[a]  (a=parameters)
    1463  KEYWORDS: Locally closed set, Canoncial form
    1464  EXAMPLE:  Prep; shows an example"
    1465 {
    1466   int te;
    1467   if (N[1]==1)
    1468   {
    1469     return(list(list(ideal(1),list(ideal(1)))));
    1470   }
    1471   int i; int j; list L0;
    1472   list Ni=minGTZ(N);
    1473   list prep;
    1474   for(j=1;j<=size(Ni);j++)
    1475   {
    1476     option(redSB);
    1477     Ni[j]=std(Ni[j]);
    1478   }
    1479   list Mij;
    1480   for (i=1;i<=size(Ni);i++)
    1481   {
    1482     Mij=minGTZ(Ni[i]+M);
    1483     for(j=1;j<=size(Mij);j++)
    1484     {
    1485       option(redSB);
    1486       Mij[j]=std(Mij[j]);
    1487     }
    1488     if ((size(Mij)==1) and (equalideals(Ni[i],Mij[1])==1)){;}
    1489     else
    1490     {
    1491         prep[size(prep)+1]=list(Ni[i],Mij);
    1492     }
    1493   }
    1494   //"T_abans="; prep;
    1495   if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
    1496   //"T_Prep="; prep;
    1497   //def Lout=CompleteA(prep,prep);
    1498   //"T_Lout="; Lout;
    1499   return(prep);
    1500 }
    1501 example
    1502 {
    1503   "EXAMPLE:"; echo = 2;
    1504   if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
    1505   short=0;
    1506   ring R=0,(x,y,z),lp;
    1507   ideal E=x*(x^2+y^2+z^2-25);
    1508   ideal N=x*(x-3),y-4;
    1509   Prep(E,N);
    1510 }
    1511 
    1512 // PtoCrep
    1513 // Computes the C-representation from the P-representation.
    1514 // input:
    1515 //    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
    1516 //         the P-representation of V(N) \ V(M)
    1517 // output:
    1518 //    list (ideal ida, ideal idb)
    1519 //    the C-representaion of V(N) \ V(M) = V(ida) \ V(idb)
    1520 // Assumed to be called in the ring @R
    1521 // Works on the ring @P
    1522 proc PtoCrep(list L)
    1523 "USAGE: PtoCrep(L)
    1524  Input L: list  [ Comp_1, .. , Comp_s ] where
    1525           Comp_i=[p_i,[p_i1,..,p_is_i] ], is
    1526           the P-representation of a locally closed set V(N) \ V(M)
    1527  RETURN:  The canonical C-representation of the locally closed set
    1528           [ P,Q ], a pair of radical ideals with P included in Q,
    1529           representing the set V(P) \ V(Q)
    1530  NOTE: Operates in a ring R=Q[a]  (a=parameters)
    1531  KEYWORDS: locally closed set, canoncial form
    1532  EXAMPLE:  PtoCrep; shows an example"
    1533 {
    1534   int te;
    1535   def RR=basering;
    1536   if(defined(@P)){te=1; setring(@P); list Lp=imap(RR,L);}
    1537   else {  te=0; def Lp=L;}
    1538   def La=PtoCrep0(Lp);
    1539   if(te==1) {setring(RR); def LL=imap(@P,La);}
    1540   if(te==0){def LL=La;}
    1541   return(LL);
    1542 }
    1543 example
    1544 {
    1545   "EXAMPLE:"; echo = 2;
    1546   if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
    1547   short=0;
    1548   ring R=0,(x,y,z),lp;
    1549 
    1550   // (P,Q) represents a locally closed set
    1551   ideal P=x^3+x*y^2+x*z^2-25*x;
    1552   ideal Q=y-4,x*z,x^2-3*x;
    1553 
    1554   // Now compute the P-representation=
    1555   def L=Prep(P,Q);
    1556   L;
    1557   // Now compute the C-representation=
    1558   def J=PtoCrep(L);
    1559   J;
    1560   // Now come back recomputing the P-represetation of the C-representation=
    1561   Prep(J[1],J[2]);
    1562 }
    1563 
    15641675// PtoCrep0
    15651676// Computes the C-representation from the P-representation.
     
    15701681//    list (ideal ida, ideal idb)
    15711682//    the C-representation of V(N) \ V(M) = V(ida) \ V(idb)
    1572 // Works in a ring Q[u_1,..,u_m] and is called on it.
     1683// Assumed to be called in a ring Q[x] (example @P)
    15731684static proc PtoCrep0(list L)
    15741685{
     
    15801691  {
    15811692    option(returnSB);
     1693    //"T_Lp[i]="; Lp[i];
    15821694    N=Lp[i][1];
     1695    Lb=Lp[i][2];
     1696    //"T_Lb="; Lb;
    15831697    ida=intersect(ida,N);
    1584     Lb=Lp[i][2];
    15851698    for(j=1;j<=size(Lb);j++)
    15861699    {
     
    15941707
    15951708// input: F a parametric ideal in Q[a][x]
    1596 // output: a rComprehensive Groebner System disjoint and reduced.
     1709// output: a disjoint and reduced Groebner System.
    15971710//      It uses Kapur-Sun-Wang algorithm, and with the options
    15981711//      can compute the homogenization before  (('can',0) or ( 'can',1))
    15991712//      and dehomogenize the result.
    16001713proc cgsdr(ideal F, list #)
    1601 "USAGE: cgsdr(F);
     1714"USAGE: cgsdr(F,options]);
    16021715          F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed.
    16031716          To compute a disjoint, reduced Comprehensive Groebner System (CGS).
     
    16101723          -option name, value- of valid options must be added to the call.
    16111724
    1612           Options:
     1725OPTIONS:
    16131726            \"can\",0-1-2: The default value is \"can\",2. In this case no
    16141727                homogenization is done. With option (\"can\",0) the given
     
    16741787          x=variables), and should be defined previously, and the ideal
    16751788          defined on R.
    1676 KEYWORDS: CGS, disjoint, reduced, Comprehensive Groebner System
     1789KEYWORDS: CGS; disjoint; reduced; Comprehensive Groebner System
    16771790EXAMPLE:  cgsdr; shows an example"
    16781791{
     
    17251838  if(comment>=1)
    17261839  {
    1727     string("Begin cgsdr with options: ",LL);
     1840    " "; string("Begin cgsdr with options: ",LL);
    17281841  }
    17291842  int ish;
     
    17311844  if (ish)
    17321845  {
    1733     if(comment>0){string("The given system is homogneous");}
     1846    if(comment>0){" "; string("The given system is homogneous");}
    17341847    def GS=KSW(B,LL);
    17351848    //can=0;
     
    17411854  {
    17421855    // WITHOUT HOMOHGENIZING
    1743     if(comment>0){string("Option of cgsdr: do not homogenize");}
     1856    if(comment>0){" "; string("Option of cgsdr: do not homogenize");}
    17441857    def GS=KSW(B,LL);
    17451858    setglobalrings();
     
    17501863    {
    17511864      // COMPUTING THE HOMOGOENIZED IDEAL
    1752       if(comment>=1){string("Homogenizing the whole ideal: option can=1");}
     1865      if(comment>0){" "; string("Homogenizing the whole ideal: option can=1");}
    17531866      list RRL=ringlist(RR);
    17541867      RRL[3][1][1]="dp";
     
    17651878      def B1=imap(RR,B);
    17661879      option(redSB);
     1880      if(comment>0){" ";string("Basis before computing its std basis="); B1;}
    17671881      B1=std(B1);
     1882      if(comment>0){" ";string("Basis after computing its std basis="); B1;}
    17681883      setring(RR);
    17691884      def B2=imap(RP,B1);
     
    17711886    else
    17721887    { // (can=0)
    1773        if(comment>0){string("Homogenizing the basis: option can=0");}
     1888       if(comment>0){" "; string( "Homogenizing the basis: option can=0");}
    17741889      def B2=B;
    17751890    }
     
    17861901      BH[i]=homog(BH[i],@t);
    17871902    }
    1788     if (comment>=2){string("Homogenized system = "); BH;}
     1903    if (comment>0){" "; string("Homogenized system = "); BH;}
    17891904    def GSH=KSW(BH,LH);
    17901905    setglobalrings();
     
    18121927    def GS=imap(RH,GSH);
    18131928    }
    1814 
    1815 
    18161929    setglobalrings();
    18171930    if(out==0)
     
    20312144{
    20322145  //"T_H="; H;
    2033   int i; int j; int k; int te; intvec notQ; int l; list sel; int used;
     2146  int i; int j; int k; int te; intvec notQ; int l; list sel;
    20342147  intvec jtesC;
    20352148  if ((size(H)==1) and (equalideals(H[1],ideal(1)))){return(H);}
     
    20622175          if (te)
    20632176          {
    2064             used++;
    20652177            if ((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;}
    20662178            else{notQ[size(notQ)+1]=i;}
     
    21002212      }
    21012213    }
    2102     //"T_Q1_0="; Q1;
    21032214    sel=selectminideals(Q1);
    21042215    kill nQ1; list nQ1;
     
    21102221  }
    21112222  if(size(Q1)==0){Q1=ideal(1),ideal(1);}
    2112   //"T_Q1_1="; Q1;
    2113   //if(used>0){string("addpartfine was ", used, " times used");}
    21142223  return(Q1);
    21152224}
    21162225
    2117 
    2118 // Auxiliary routine for grobcov: ideal F is assumed to be homogeneous
     2226// Auxiliary rutine for gcover
     2227// Deciding if combine is needed
     2228// input: list LCU=( (basis1, p_1, (p11,..p1s1)), .. (basisr, p_r, (pr1,..prsr))
     2229// output: (tes); if tes==1 then combine is needed, else not.
     2230static proc needcombine(list LCU,ideal N)
     2231{
     2232  //"Deciding if combine is needed";;
     2233  ideal BB;
     2234  int tes=0; int m=1; int j; int k; poly sp;
     2235  while((tes==0) and (m<=size(LCU[1][1])))
     2236  {
     2237    j=1;
     2238    while((tes==0) and (j<=size(LCU)))
     2239    {
     2240      k=1;
     2241      while((tes==0) and (k<=size(LCU)))
     2242      {
     2243        if(j!=k)
     2244        {
     2245          sp=pnormalf(pspol(LCU[j][1][m],LCU[k][1][m]),LCU[k][2],N);
     2246          if(sp!=0){tes=1;}
     2247        }
     2248        k++;
     2249      }
     2250      j++;
     2251    }
     2252    if(tes){break;}
     2253    m++;
     2254  }
     2255  return(tes);
     2256}
     2257
     2258// Auxiliary routine
     2259// precombine
     2260// input:  L: list of ideals (works in @P)
     2261// output: F0: ideal of polys. F0[i] is a poly in the intersection of
     2262//             all ideals in L except in the ith one, where it is not.
     2263//             L=(p1,..,ps);  F0=(f1,..,fs);
     2264//             F0[i] \in intersect_{j#i} p_i
     2265static proc precombine(list L)
     2266{
     2267  int i; int j; int tes;
     2268  def RR=basering;
     2269  setring(@P);
     2270  list L0; list L1; list L2; list L3; ideal F;
     2271  L0=imap(RR,L);
     2272  L1[1]=L0[1]; L2[1]=L0[size(L0)];
     2273  for (i=2;i<=size(L0)-1;i++)
     2274  {
     2275    L1[i]=intersect(L1[i-1],L0[i]);
     2276    L2[i]=intersect(L2[i-1],L0[size(L0)-i+1]);
     2277  }
     2278  L3[1]=L2[size(L2)];
     2279  for (i=2;i<=size(L0)-1;i++)
     2280  {
     2281    L3[i]=intersect(L1[i-1],L2[size(L0)-i]);
     2282  }
     2283  L3[size(L0)]=L1[size(L1)];
     2284  for (i=1;i<=size(L3);i++)
     2285  {
     2286    option(redSB); L3[i]=std(L3[i]);
     2287  }
     2288  for (i=1;i<=size(L3);i++)
     2289  {
     2290    tes=1; j=0;
     2291    while((tes) and (j<size(L3[i])))
     2292    {
     2293      j++;
     2294      option(redSB);
     2295      L0[i]=std(L0[i]);
     2296      if(reduce(L3[i][j],L0[i])!=0){tes=0; F[i]=L3[i][j];}
     2297    }
     2298    if (tes){"ERROR a polynomial in all p_j except p_i was not found";}
     2299  }
     2300  setring(RR);
     2301  def F0=imap(@P,F);
     2302  return(F0);
     2303}
     2304
     2305// Auxiliary routine
     2306// combine
     2307// input: a list of pairs ((p1,P1),..,(pr,Pr)) where
     2308//    ideal pi is a prime component
     2309//    poly Pi is the polynomial in Q[a][x] on V(pi)\ V(Mi)
     2310//    (p1,..,pr) are the prime decomposition of the lpp-segment
     2311//    list crep =(ideal ida,ideal idb): the Crep of the segment.
     2312//    list Pci of the intersecctions of all pj except the ith one
     2313// output:
     2314//    poly P on an open and dense set of V(p_1 int ... p_r)
     2315static proc combine(list L, ideal F)
     2316{
     2317  // ATTENTION REVISE AND USE Pci and F
     2318  int i; poly f;
     2319  f=0;
     2320  for(i=1;i<=size(L);i++)
     2321  {
     2322    f=f+F[i]*L[i][2];
     2323  }
     2324//   f=elimconstfac(f);
     2325  f=primepartZ(f);
     2326  return(f);
     2327}
     2328
     2329// Central routine for grobcov: ideal F is assumed to be homogeneous
    21192330// gcover
    21202331// input: ideal F: a generating set of a homogeneous ideal in Q[a][x]
     
    21272338{
    21282339  int i; int j; int k; ideal lpp; list GPi2; list pairspP; ideal B; int ti;
    2129   int i1; int tes; int j1; int selind; int i2; int m;
     2340  int i1; int tes; int j1; int selind; int i2; //int m;
    21302341  list prep; list crep; list LCU; poly p; poly lcp; ideal FF;
    21312342  list lpi;
     
    21722383  list S;
    21732384  poly sp;
    2174   ideal BB;
    21752385  for (i=1;i<=size(GP);i++)
    21762386  {
     
    21952405    }
    21962406    //"Deciding if combine is needed";
    2197     kill BB;
    2198     ideal BB;
    2199     tes=1; m=1;
    2200     while((tes) and (m<=size(LCU[1][1])))
    2201     {
    2202       j=1;
    2203       while((tes) and (j<=size(LCU)))
    2204       {
    2205         k=1;
    2206         while((tes) and (k<=size(LCU)))
    2207         {
    2208           if(j!=k)
    2209           {
    2210             sp=pnormalf(pspol(LCU[j][1][m],LCU[k][1][m]),LCU[k][2],N);
    2211             if(sp!=0){tes=0;}
    2212           }
    2213           k++;
    2214         }        //setglobalrings();
    2215         if(tes)
    2216         {
    2217           BB[m]=LCU[j][1][m];
    2218         }
    2219         j++;
    2220       }
    2221       if(tes==0){break;}
    2222       m++;
    2223     }
    22242407    crep=PtoCrep(prep);
    2225     if(tes==0)
     2408    if(size(LCU)>1){tes=1;}
     2409    else
     2410    {
     2411      tes=0;
     2412      for(k=1;k<=size(B);k++){B[k]=pnormalf(B[k],crep[1],crep[2]);}
     2413    }
     2414    // tes=needcombine(LCU,N);
     2415    // if(tes==1){" ";string("combine is needed for segment ",i);" ";}
     2416    //crep=PtoCrep(prep);
     2417    if(tes)
    22262418    {
    22272419      // combine is needed
     
    22312423        LL[j]=LCU[j][2];
    22322424      }
    2233       if (size(LCU)>1)
    2234       {
    2235         FF=precombint(LL);
    2236       }
     2425      FF=precombine(LL);
    22372426      for (k=1;k<=size(lpp);k++)
    22382427      {
     
    22422431          L[j]=list(LCU[j][2],LCU[j][1][k]);
    22432432        }
    2244         if (size(LCU)>1)
    2245         {
    2246           B[k]=combine(L,FF);
    2247         }
    2248         else{B[k]=L[1][2];}
    2249       }
    2250     }
    2251     else{B=BB;}
     2433        B[k]=combine(L,FF);
     2434      }
     2435    }
    22522436    for(j=1;j<=size(B);j++)
    22532437    {
     
    22772461//            N is the null conditions ideal (if desired)
    22782462//            W is the ideal of non-null conditions (if desired)
    2279 //            The value of \"can\"i s 1 by default and can be set to 0 if we do not
     2463//            The value of \"can\" is 1 by default and can be set to 0 if we do not
    22802464//            need to obtain the canonical GC, but only a GC.
    22812465//            The value of \"ext\" is 0 by default and so the generic representation
     
    22922476//            given by the lpp, the basis, and the P-representation of the segment
    22932477proc grobcov(ideal F,list #)
    2294 "USAGE:   grobcov(F); This is the fundamental routine of the
     2478"USAGE:   grobcov(F,options]);
     2479          This is the fundamental routine of the
    22952480          library. It computes the Groebner cover of a parametric ideal. See
    22962481                Montes A., Wibmer M.,
     
    23072492          -option name, value- of valid options must be added to the call.
    23082493
    2309           Options:
     2494OPTIONS:
    23102495            \"null\",ideal E: The default is (\"null\",ideal(0)).
    23112496            \"nonnull\",ideal N: The default is (\"nonnull\",ideal(1)).
     
    23322517                \"comment\" higher will provide information about the
    23332518                development of the computation.
     2519             \"showhom\",0-1: The default is (\"showhom\",0). Setting
     2520                \"showhom\",1 will output the set of homogenized lpp of
     2521                each segment as last element.
    23342522          One can give none or whatever of these options.
    23352523RETURN:   The list
     2524          (
     2525           (lpp_1,basis_1,segment_1),
     2526           ...
     2527           (lpp_s,basis_s,segment_s)
     2528          )
     2529          optionally
    23362530          (
    23372531           (lpp_1,basis_1,segment_1,lpph_1),
     
    23432537          set of lpp of the reduced Groebner basis for each point
    23442538          of the segment.
     2539          With option (\"showhom\",1) the lpph will be shown:
    23452540          The lpph corresponds to the lpp of the homogenized ideal
    23462541          and is different for each segment. It is given as a string,
     
    23792574          x=variables), and should be defined previously. The ideal must
    23802575          be defined on R.
    2381 KEYWORDS: Groebner cover, parametric ideal, canonical, discussion of
    2382           parametric ideal.
     2576KEYWORDS: Groebner cover; parametric ideal; canonical; discussion of
     2577          parametric ideal
    23832578EXAMPLE:  grobcov; shows an example"
    23842579{
     
    23902585  setglobalrings();
    23912586  list L0=#;
     2587  list Se;
    23922588  int out=0;
     2589  int showhom=0;
     2590  int hom;
    23932591  L0[size(L0)+1]="res"; L0[size(L0)+1]=ideal(1);
    23942592  // default options
     
    24212619            {
    24222620              if (L0[2*i-1]=="comment"){comment=L0[2*i];}
     2621              else
     2622              {
     2623                if (L0[2*i-1]=="showhom"){showhom=L0[2*i];}
     2624              }
    24232625            }
    24242626          }
     
    24452647  LL[11]="ext";    LL[12]=extop;
    24462648  LL[13]="rep";    LL[14]=repop;
     2649  LL[15]="showhom";    LL[16]=showhom;
    24472650  if (comment>=1)
    24482651  {
     
    24522655  def S=gcover(F,LL);
    24532656  // NOW extend
     2657  if (size(S[1])==5){hom=1;}
    24542658  if(extop)
    24552659  {
    24562660    S=extend(S,LL);
    24572661  }
     2662  // NOW representation of the segments by option repop
     2663  list Si; list nS;
     2664  if (repop==0)
     2665  {
     2666    for(i=1;i<=size(S);i++)
     2667    {
     2668      Si=list(S[i][1],S[i][2],S[i][3]);
     2669      if(hom){Si[size(Si)+1]=S[i][4];}
     2670      nS[size(nS)+1]=Si;
     2671    }
     2672    S=nS;
     2673  }
    24582674  else
    24592675  {
    2460     // NOW representation of the segments by option repop
    2461     list Si; list nS;
    2462     if(repop==0)
     2676    if (repop==1)
    24632677    {
    24642678      for(i=1;i<=size(S);i++)
    24652679      {
    2466         Si=list(S[i][1],S[i][2],S[i][3],S[i][5]);
     2680        Si=list(S[i][1],S[i][2],S[i][4]);
     2681        if(hom){Si[size(Si)+1]=S[i][5];}
    24672682        nS[size(nS)+1]=Si;
    24682683      }
    2469       kill S;
    2470       def S=nS;
     2684      S=nS;
    24712685    }
    24722686    else
    24732687    {
    2474       if(repop==1)
    2475       {
    2476         for(i=1;i<=size(S);i++)
    2477         {
    2478           Si=list(S[i][1],S[i][2],S[i][4],S[i][5]);
    2479           nS[size(nS)+1]=Si;
    2480         }
    2481         kill S;
    2482         def S=nS;
    2483       }
    2484       else
    2485       {
    2486         for(i=1;i<=size(S);i++)
    2487         {
    2488           Si=list(S[i][1],S[i][2],S[i][3],S[i][4],S[i][5]);
    2489           nS[size(nS)+1]=Si;
    2490         }
    2491         kill S;
    2492         def S=nS;
     2688      for(i=1;i<=size(S);i++)
     2689      {
     2690        Si=list(S[i][1],S[i][2],S[i][3],S[i][4]);
     2691        if(hom){Si[size(Si)+1]=S[i][5];}
     2692        nS[size(nS)+1]=Si;
    24932693      }
    24942694    }
     
    25002700  }
    25012701  if(defined(@P)==1){kill @R; kill @P; kill @RP;}
     2702  if (hom and showhom==0)
     2703  {
     2704    for(i=1;i<=size(S);i++)
     2705    {
     2706      Se=S[i];
     2707      Se=delete(Se,size(Se));
     2708      S[i]=Se;
     2709    }
     2710  }
    25022711  return(S);
    25032712}
     
    25232732// Can be called from the top
    25242733proc extend(list GC, list #);
    2525 "USAGE: extend(GC); The default option of grobcov provides
     2734"USAGE: extend(GC[,options]);
     2735          The default option of grobcov provides
    25262736          the bases in generic representation (the I-regular functions
    2527           of the bases ara given by a single polynomial. It can specialize
     2737          of the bases are given by a single polynomial. It can specialize
    25282738          to zero for some points of the segments, but in general, it
    25292739          is sufficient for many pouposes. Nevertheless the I-regular
    2530           functions allow a full representation given bey a set of
     2740          functions allow a full representation given by a set of
    25312741          polynomials specializing to the value of the function (after normalization)
    25322742          or to zero, but at least one of the polynomials specializes to non-zero.
     
    25522762          x=variables), and should be defined previously. The ideal must
    25532763          be defined on R.
    2554 KEYWORDS: Groebner cover, parametric ideal, canonical, discussion of
    2555           parametric ideal, full representation.
     2764KEYWORDS: Groebner cover; parametric ideal; canonical, discussion of
     2765          parametric ideal; full representation
    25562766EXAMPLE:  extend; shows an example"
    25572767{
     
    25672777    i++;
    25682778  }
     2779  int hom;
    25692780  if(m==1){"Warning! grobcov has already extended bases"; return(S);}
    2570   if(size(GC[1])!=5){"Warning! extend make sense only when grobcov has been called with options  'rep',2,'ext',0"; " "; return();}
     2781  if(not(size(GC[1])==5 or (size(GC[1])==4 and typeof(GC[1][4])=="list")))
     2782  {"Warning! extend make sense only when grobcov has been called with options  'rep',2,'ext',0"; " "; return();}
     2783  if(size(GC[1])==5){hom=1;}
    25712784  int repop=0;
    25722785  int start3=timer;
     
    26342847    for(i=1;i<=size(S);i++)
    26352848    {
    2636       Si=list(S[i][1],S[i][2],S[i][3],S[i][5]);
     2849      Si=list(S[i][1],S[i][2],S[i][3]);
     2850      if(hom){Si[size(Si)+1]=S[i][5];}
    26372851      nS[size(nS)+1]=Si;
    26382852    }
     
    26452859      for(i=1;i<=size(S);i++)
    26462860      {
    2647         Si=list(S[i][1],S[i][2],S[i][4],S[i][5]);
     2861        Si=list(S[i][1],S[i][2],S[i][4]);
     2862        if(hom){Si[size(Si)+1]=S[i][5];}
    26482863        nS[size(nS)+1]=Si;
    26492864      }
     
    26542869      for(i=1;i<=size(S);i++)
    26552870      {
    2656         Si=list(S[i][1],S[i][2],S[i][3],S[i][4],S[i][5]);
     2871        Si=list(S[i][1],S[i][2],S[i][3],S[i][4]);
     2872        if(hom){Si[size(Si)+1]=S[i][5];}
    26572873        nS[size(nS)+1]=Si;
    26582874      }
    2659 
    26602875    }
    26612876  }
     
    26702885  short=0;
    26712886  ideal S=a0*x^2+b0*x+c0,
    2672             a1*x^2+b1*x+c1;
     2887          a1*x^2+b1*x+c1;
    26732888  def GCS=grobcov(S,"rep",2);
    26742889  GCS;
     
    27182933  return(fr);
    27192934}
    2720 
    2721 // Auxiliary routine
    2722 // deltai
    2723 // input:
    2724 //   int i:
    2725 //   list LPr: (p1,..,pr) of prime components of an ideal in Q[a]
    2726 // output:
    2727 //   list (fr,fnr) of two polynomials that are equal on V(pi)
    2728 //       and fr=0 on V(P) \ V(pi), and fnr is nonzero on V(pj) for all j.
    2729 static proc deltai(int i, list LPr)
    2730 {
    2731   def RR=basering;
    2732   setring(@P);
    2733   def LP=imap(RR,LPr);
    2734   int j; poly p;
    2735   def F=ideal(1);
    2736   poly f;
    2737   poly fn;
    2738   ideal LPi;
    2739   for (j=1;j<=size(LP);j++)
    2740   {
    2741     if (j!=i)
    2742     {
    2743       F=idint(F,LP[j]);
    2744     }
    2745   }
    2746   p=0; j=1;
    2747   while ((p==0) and (j<=size(F)))
    2748   {
    2749     LPi=LP[i];
    2750     attrib(LPi,"isSB",1);
    2751     p=reduce(F[j],LPi);
    2752     j++;
    2753   }
    2754   f=F[j-1];
    2755   fn=nonzerodivisor(f,LP);
    2756   setring(RR);
    2757   def fr=imap(@P,f);
    2758   def fnr=imap(@P,fn);
    2759   return(list(fr,fnr));
    2760 }
    2761 
    2762 // Auxiliary routine
    2763 // combine
    2764 // input: a list of pairs ((p1,P1),..,(pr,Pr)) where
    2765 //    ideal pi is a prime component
    2766 //    poly Pi is the polynomial in Q[a][x] on V(pi)\ V(Mi)
    2767 //    (p1,..,pr) are the prime decomposition of the lpp-segment
    2768 //    list crep =(ideal ida,ideal idb): the Crep of the segment.
    2769 //    list Pci of the intersecctions of all pj except the ith one
    2770 // output:
    2771 //    poly P on an open and dense set of V(p_1 int ... p_r)
    2772 static proc combine(list L, ideal F)
    2773 {
    2774   // ATTENTION REVISE AND USE Pci and F
    2775   int i; poly f;
    2776   f=0;
    2777   for(i=1;i<=size(L);i++)
    2778   {
    2779     f=f+F[i]*L[i][2];
    2780   }
    2781 //   f=elimconstfac(f);
    2782   f=primepartZ(f);
    2783   return(f);
    2784 }
    2785 
    27862935
    27872936//Auxiliary routine
     
    30133162        M1=intersect(T0[2],T1[2]);
    30143163      }
    3015       T=list(ci,PtoCrep(Prep(N1,M1)));
     3164      T=list(ci,PtoCrep0(Prep0(N1,M1)));
    30163165      LL[size(LL)+1]=T;
    30173166      if(equalideals(T[2][1],ideal(1))){te=1; break;}
     
    30463195  return(T);
    30473196}
    3048 
    30493197
    30503198// Auxiliary routine
     
    32373385}
    32383386
    3239 // // Auxiliary routine
    3240 // // NonNull: returns 1 if the poly f is nonnull on V(E) \ V(N), 0 otherwise.
    3241 // // Input:
    3242 // //   f: polynomial
    3243 // //   E: null ideal
    3244 // //   N: nonnull ideal
    3245 // // Output:
    3246 // //   1 if f is nonnul in V(P) \ V(Q),
    3247 // //   0 if it has zeroes inside.
    3248 // //  Works in @P
    3249 // proc NonNull(poly f, ideal E, ideal N)
    3250 // {
    3251 //   int te=1; int i;
    3252 //   def RR=basering;
    3253 //   setring(@P);
    3254 //   def fp=imap(RR,f);
    3255 //   def Ep=imap(RR,E);
    3256 //   def Np=imap(RR,N);
    3257 //   ideal H;
    3258 //   ideal Ef=Ep+fp;
    3259 //   for (i=1;i<=size(Np);i++)
    3260 //   {
    3261 //     te=radicalmember(Np[i],Ef);
    3262 //     if (te==0){break;}
    3263 //   }
    3264 //   setring(RR);
    3265 //   return(te);
    3266 // }
    3267 
    32683387// Auxiliary routine
    32693388// selectextendcoef
     
    34273546}
    34283547
    3429 
    34303548// Auxiliary routine
    34313549// precombint
     
    34743592  return(F0);
    34753593}
    3476 
    34773594
    34783595// Auxiliary routine
     
    36163733  int i;
    36173734  def L=#;
    3618   if (size(L)>0)
     3735  if(size(L)>0)
    36193736  {
    36203737    for (i=1;i<=size(L) div 2;i++)
     
    36583775// Auxiliary routine
    36593776// sqf
    3660 // This is for releases of Singular before 3-5-1
    3661 // proc sqf(poly f)
    3662 // {
    3663 //  def RR=basering;
    3664 //  setring(@P);
    3665 //  def ff=imap(RR,f);
    3666 //  def G=sqrfree(ff);
    3667 //  poly fff=1;
    3668 //  int i;
    3669 //  for (i=1;i<=size(G);i++)
    3670 //  {
    3671 //    fff=fff*G[i];
    3672 //  }
    3673 //  setring(RR);
    3674 //   def ffff=imap(@P,fff);
    3675 //   return(ffff);
    3676 // }
    3677 
    3678 // Auxiliary routine
    3679 // sqf
    36803777static proc sqf(poly f)
    36813778{
     
    36883785  return(ffff);
    36893786}
    3690 
    36913787
    36923788// Auxiliary routine
     
    38033899    if (equalideals(Gr,ideal(0))==1){GrHi=H[i];}
    38043900    else {GrHi=Gr,H[i];}
    3805 //     else {for (j=1;j<=size(Gr);j++){GrHi[size(GrHi)+1]=Gr[j]*H[i];}}
    38063901    if(comment>1){"Call to KSW with arguments "; GM; GrHi;  Nh;}
    38073902    KS=KSW0(GM,GrHi,Nh,comment);
     
    39814076  int n=size(L);
    39824077  if(n<2){return(A);}
     4078  intvec w;
    39834079  for(i=1;i<=size(L);i++)
    39844080  {
     
    39904086        {
    39914087          L[i][2]=L[j][2];
     4088          w[size(w)+1]=j;
    39924089        }
    39934090      }
     4091    }
     4092  }
     4093  if(size(w)>0)
     4094  {
     4095    for(i=1; i<=size(w);i++)
     4096    {
     4097      j=w[size(w)+1-i];
     4098      L=elimfromlist(L, j);
    39944099    }
    39954100  }
     
    40134118  }
    40144119  if(equalideals(T,ideal(1))==0){L[size(L)+1]=list(std(T),ideal(1))};
    4015   //string("T_n=",n," new n0",size(L));
    40164120  return(L);
    40174121}
    40184122
    4019 // Input: list(A)
    4020 //          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
    4021 //          whose union is a constructible set from
    4022 // Output list [Lev,C]:
    4023 //          where Lev is the Crep of the canonical first level of A, and
    4024 //          C is the complement of the first level Lev wrt to the closure of A. The elements are given in Crep,
     4123// input list A=[[p1,q1],...,[pn,qn]] :
     4124//                    the list of segments of a constructible set S, where each [pi,qi] is given in C-representation
     4125// output list [topA,C]
     4126//       where topA is the closure of A
     4127//                 C is the list of segments of the complement of A given in C-representation
    40254128static proc FirstLevel(list A)
    40264129{
     
    40284131  list T=zeroone(n);
    40294132  ideal P; ideal Q;
    4030   list Cb;  ideal Cc=ideal(1);
     4133  list Cb;  ideal Cc=1;
    40314134  int i; int j;
    40324135  intvec t;
    4033   ideal Z=ideal(1);
     4136  ideal topA=1;
    40344137  list C;
    4035   for(i=1;i<=size(A);i++)
    4036   {
    4037     Z=intersect(Z,A[i][1]);
    4038   }
     4138  for(i=1;i<=n;i++)
     4139  {
     4140    topA=intersect(topA,A[i][1]);
     4141  }
     4142  //topA=std(topA);
    40394143  for(i=2; i<=size(T);i++)
    40404144  {
    40414145    t=T[i];
    4042     P=ideal(1); Q=ideal(0);
     4146    //"T_t="; t;
     4147    P=0; Q=1;
    40434148    for(j=1;j<=n;j++)
    40444149    {
    40454150      if(t[n+1-j]==1)
    40464151      {
    4047         Q=Q+A[j][2];
     4152        P=P+A[j][2];
    40484153      }
    40494154      else
    40504155      {
    4051         P=intersect(P,A[j][1]);
    4052       }
    4053     }
    4054     //"T_Q="; Q; "T_P="; P;
    4055     Cb=Crep(Q,P);
     4156        Q=intersect(Q,A[j][1]);
     4157      }
     4158    }
     4159    Cb=Crep0(P,Q);
    40564160    //"T_Cb="; Cb;
    4057     if(Cb[1][1]<>1)
    4058     {
    4059       C[size(C)+1]=Cb;
    4060       Cc=intersect(Cc,Cb[1]);
    4061     }
    4062   }
    4063   list Lev=list(Z,Cc);                //Crep(Z,Cc);
     4161    if(size(Cb)!=0)
     4162    {
     4163      if( Cb[1][1]<>1)
     4164      {
     4165        C[size(C)+1]=Cb;
     4166      }
     4167    }
     4168  }
    40644169  if(size(C)>1){C=SimplifyUnion(C);}
    4065   return(list(Lev,C));
    4066 }
    4067 
    4068 // Input: list(A)
    4069 //          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
    4070 //          whose union is a constructible set from which the algorithm computes its canonical form.
     4170  return(list(topA,C));
     4171}
     4172
     4173// Input:
    40714174// Output:
    4072 //     list [L,C] where
    4073 //          where L is the list of canonical levels of A,
    4074 //          and C is the list of canonical levels of the complement of A wrt to the closure of A.
    4075 //          All locally closed sets are given in Crep.
    4076 //          L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]] is the list of levels of A in Crep.
    4077 //          C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]  is the list of levels of C,
    4078 //                                              the complement of S wrt the closure of A, in Crep.
    4079 proc ConsLevels(list A)
    4080 "USAGE:   ConsLevels(A);
    4081           A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
    4082           whose union is a constructible set from which the algorithm computes its
    4083           canonical form.
    4084 RETURN: A list with two elements: [L,C]
    4085           where L is the list of canonical levels of A, and C is the list of canonical
    4086           levels of the
    4087           complement of A wrt to the closure of A.
    4088           All locally closed sets are given in Crep.
    4089           L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]]
    4090           C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]
    4091 NOTE: The basering R, must be of the form Q[a]
    4092 KEYWORDS: locally closed set, constructible set
    4093 EXAMPLE:  ConsLevels; shows an example"
    4094 {
    4095   list L; list C; int i;
    4096   list B; list T;
     4175static proc ConstoPrep(list L)
     4176{
     4177  list L1;
     4178  int i; int j;
     4179  list aux;
     4180  for(i=1;i<=size(L);i++)
     4181  {
     4182    aux=Prep(L[i][2][1],L[i][2][2]);
     4183    L1[size(L1)+1]=list(L[i][1],aux);
     4184  }
     4185  return(L1);
     4186}
     4187
     4188// Input:
     4189//     list A =  [[P1,Q1], .. [Pn,Qn]]
     4190//                  A constructible set as union of locally closed sets represented by pairs of ideals
     4191// Output:
     4192//     list L =[p1,p2,p3,...,pk]
     4193//        where pi is the ideal of the closure of level i alternatively of  A or its complement
     4194//        Note: the levels of A are [p1,p2], [p3,p4], [p5,p6],...
     4195//                 the levels of C are [p2,p3],[p4,p5], ...
     4196//                 expressed in C-representation
     4197//    Assumed to be called in the ring @R or the ring @P or a ring ring Q[a]
     4198proc ConsLevels(list A0)
     4199"USAGE:   ConsLevels(L);
     4200         where L=[[P1,Q1],...,[Ps,Qs]]
     4201         is a list of lists of of pairs of ideals defined over Grobcov::@P,
     4202         or over a ring without parameters and represents the constructible set
     4203         S=V(P1)\ V(Q1) u ... u V(Ps)\ V(Qs)
     4204RETURN: The list of ideals [a1,a2,...,at] representing the closures of the
     4205         canonical levels of S and its complement C wrt to the closure of S.
     4206         The canonical levels of S are represented by theirs Crep. So we have:
     4207         Levels of S:  [a1,a2],[a3,a4],...
     4208         Levels of C:  [a2,a3],[a4,a5],...
     4209         S=V(a1)\ V(a2) u V(a3)\ V(a4) u ...
     4210         C=V(a2\ V(a3) u V(a4)\ V(a5) u ...
     4211         It is used internally by locus procedure.
     4212NOTE: It is described in the paper
     4213          J.M. Brunat, A. Montes
     4214          Computing the canonical representation of constructible sets.
     4215          Math.  Comput. Sci. (2016) 19: 165-178.
     4216KEYWORDS: constructible set; locally closed set; canonical form
     4217EXAMPLE:  example ConsLevels; shows an example"
     4218{
     4219  def RR=basering;
     4220  int te;
     4221  if(defined(@P)){te=1; setring(@P); list A=imap(RR,A0);}
     4222  else {te=0; def A=A0;}
     4223
     4224  list L; list C;
     4225  list B; list T; int i;
    40974226  for(i=1; i<=size(A);i++)
    40984227  {
    4099     T=Crep(A[i][1],A[i][2]);
     4228    T=Crep0(A[i][1],A[i][2]);
    41004229    B[size(B)+1]=T;
    41014230  }
    4102   int level=0;
    41034231  list K;
    41044232  while(size(B)>0)
    41054233  {
    4106     level++;
    41074234    K=FirstLevel(B);
    4108     if(level mod 2==1){L[size(L)+1]=list(level,K[1]);}
    4109     else{C[size(C)+1]=list(level,K[1]);}
    4110     //"T_L="; L;
    4111     //"T_C="; C;
     4235    //"T_K="; K;
     4236    L[size(L)+1]=K[1];
    41124237    B=K[2];
    4113     //"T_size(B)="; size(B);
    4114   }
    4115   return(list(L,C));
     4238   }
     4239  L[size(L)+1]=ideal(1);
     4240  if(te==1) {setring(RR); def LL=imap(@P,L);}
     4241  if(te==0){def LL=L;}
     4242  return(LL);
    41164243}
    41174244example
    4118 { "EXAMPLE:"; echo = 2;
    4119 // Example of ConsLevels with nice geometrical interpretetion and 2 levels
    4120 
    4121 if (defined(R)){kill R;}
    4122 if (defined(@P)){kill @P; kill @R; kill @RP;}
    4123 
    4124   ring R=0,(x,y,z),lp;
    4125   short=0;
    4126   ideal P1=x*(x^2+y^2+z^2-1);
    4127   ideal Q1=z,x^2+y^2-1;
    4128   ideal P2=y,x^2+z^2-1;
    4129   ideal Q2=z*(z+1),y,x*(x+1);
    4130 
    4131   list Cr1=Crep(P1,Q1);
    4132   list Cr2=Crep(P2,Q2);
    4133 
    4134   list L=list(Cr1,Cr2);
    4135   L;
    4136   // ConsLevels(L)=
    4137   ConsLevels(L);
    4138 
    4139 //----------------------------------------------------------------------
    4140 // Begin Problem S93
    4141 // Automatic theorem proving
    4142 // Generalized Steiner-Lehmus Theorem
    4143 // A.Montes, T.Recio
    4144 
    4145 // Bisector of A(-1,0) = Bisector of B(1,0) varying C(a,b)
    4146 
    4147 if(defined(R1)){kill R1;}
    4148 ring R1=(0,a,b),(x1,y1,x2,y2,p,r),dp;
    4149 ideal S93=(a+1)^2+b^2-(p+1)^2,
    4150                     (a-1)^2+b^2-(1-r)^2,
    4151                     a*y1-b*x1-y1+b,
    4152                     a*y2-b*x2+y2-b,
    4153                     -2*y1+b*x1-(a+p)*y1+b,
    4154                     2*y2+b*x2-(a+r)*y2-b,
    4155                     (x1+1)^2+y1^2-(x2-1)^2-y2^2;
    4156 
     4245{
     4246"EXAMPLE:"; echo = 2;
     4247ring R=0,(x,y,z),lp;
    41574248short=0;
    4158 def GC93=grobcov(S93,"nonnull",ideal(b),"rep",1);
    4159 //"grobcov(S93,'nonnull',ideal(b),"rep",1)="; GC93;
    4160 
    4161 list L0;
    4162 for(int i=1;i<=size(GC93);i++)
    4163 {
    4164    L0[size(L0)+1]=GC93[i][3];
    4165 }
    4166 
    4167 def GC93ab=grobcov(S93,"nonnull",ideal(a*b),"rep",1);
    4168 
    4169 ring RR=0,(a,b),lp;
    4170 
    4171 list L1;
    4172 L1=imap(R1,L0);
    4173 // L1=Total elements of the grobcov(S93,'nonnull',ideal(b),'rep',1) to be added=
    4174 L1;
    4175 
    4176 // Adding all the elements of grobcov(S93,'nonnull',ideal(b),'rep',1)=
    4177 ConsLevels(L1);
     4249ideal P1=x*(x^2+y^2+z^2-1);
     4250ideal Q1=z,x^2+y^2-1;
     4251ideal P2=y,x^2+z^2-1;
     4252ideal Q2=z*(z+1),y,x*(x+1);
     4253
     4254list Cr1=Crep(P1,Q1);
     4255list Cr2=Crep(P2,Q2);
     4256
     4257list L=list(Cr1,Cr2);
     4258L;
     4259ConsLevels(L);
     4260}
     4261
     4262// Converts the output of ConsLevels, given by the set of closures of the Levels of the constructible S
     4263//     to an expression where the Levels are apparent.
     4264// Input: The ouput of ConsLevels of the form
     4265//    [A1,A2,..,Ak], where the Ai's are the closures of the levels.
     4266// Output: An expression of the form
     4267//      L1=[[1,[A1,A2]],[3,[A3,A4]],..,[2l-1,[A_{2l-1},A_{2l}]]] the list of Levels of S
     4268proc ConsLevelsToLevels(list L)
     4269"USAGE: ConsLevelsToLevels(list L);
     4270          The input list L must be the output of the call to
     4271          the routine 'ConsLevels' of a constructible set:
     4272           [A1,A2,..,Ak], where the Ai's are the closures of the levels.
     4273          'ConsLevels' outputs the closures of the levels
     4274          of the constructible and of its complements.
     4275          This routine selects the  levels of the
     4276          constructible set.
     4277RETURN: The levels of the constructible set:
     4278          Lc=[[1,[A1,A2]],[3,[A3,A4]],..,[2l-1,[A_{2l-1},A_{2l}]]] the list of
     4279          Levels of S
     4280NOTE: It must be called to the output of the 'ConsLevels' routine.
     4281KEYWORDS: constructible sets
     4282EXAMPLE:  ConsLevelsToLevels shows an example"
     4283{
     4284  int n=size(L) div 2;
     4285  int i;
     4286  list L1; list L2;
     4287  for(i=1; i<=n;i++)
     4288  {
     4289    L1[size(L1)+1]=list(2*i-1,list(L[2*i-1],L[2*i]));
     4290  }
     4291  return(L1);
     4292}
     4293example
     4294{
     4295"EXAMPLE:"; echo = 2;
     4296ring R=0,(x,y,z),lp;
     4297short=0;
     4298ideal P1=(x^2+y^2+z^2-1);
     4299ideal Q1=z,x^2+y^2-1;
     4300ideal P2=y,x^2+z^2-1;
     4301ideal Q2=z*(z+1),y,x*(x+1);
     4302ideal P3=x;
     4303ideal Q3=5*z-4,5*y-3,x;
     4304
     4305list Cr1=Crep(P1,Q1);
     4306list Cr2=Crep(P2,Q2);
     4307list Cr3=Crep(P3,Q3);
     4308
     4309list L=list(Cr1,Cr2,Cr3);
     4310L;
     4311
     4312def LL=ConsLevels(L);
     4313LL;
     4314ConsLevelsToLevels(LL);
    41784315}
    41794316
     
    42474384//   ideal N: the holes of the component
    42484385// Output:
    4249 //   int d: the dimension of B on the variables reduced by the holes.
     4386//   int d: the dimension of B on the variables (antiimage).
    42504387//     if d>0    then the component is 'Normal'
    42514388//     if d==0 then the component is 'Singular'
     
    42684405  }
    42694406  if(F!=0){Fenv=1;}
    4270   //"T_B="; B; "E="; E; "N="; N;
    42714407  list L0;
    42724408  if(dimP0(E)==0){L0=2,"Normal";} // 2 es fals pero ha de ser >0 encara que sigui 0
     
    42754411    if(version==0)
    42764412    {
    4277       //"T_B="; B;  // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part
     4413      // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part
    42784414      // independent of parameters giving the variables with dimension 0
    42794415     dd=indepparameters(B);
    4280       if (dd==1){d=0; L0=d,string(B),determineF(B,F,E);}
     4416      if (dd==1){d=0; L0=d,string(B);} // ,determineF(B,F,E)
    42814417      else{d=1; L0=2,"Normal";}
    42824418    }
     
    43404476      {
    43414477        //" ";string("Antiimage of Special component = ", GGG);
    4342          if(Fenv==1)
    4343         {
    4344           L0[3]=determineF(GGG,F,E);
    4345         }
    43464478      }
    43474479      else
     
    43514483    }
    43524484  }
     4485  //"T_L0="; L0;
    43534486  return(L0);
    43544487}
     
    43814514   ideal M=std(AA+FH);
    43824515   def rh=reduce(EH,M);
     4516   //"T_AA="; AA; "T_FH="; FH; "T_EH="; EH; "T_rh="; rh;
    43834517   if(rh==0){env=1;} else{env=0;}
    43844518   setring RR;
    43854519          //L0[3]=env;
     4520    //"T_env="; env;
    43864521    return(env);
    43874522}
    43884523
    4389 // DimPar
    4390 // Auxilliary routine to NorSing determining the dimension of a parametric ideal
    4391 // Does not use @P and define it directly because is assumes that
    4392 //                             variables and parameters have been inverted
    4393 static proc DimPar(ideal E)
    4394 {
    4395   def RRH=basering;
    4396   def RHx=ringlist(RRH);
    4397   def Prin=ring(RHx[1]);
    4398   setring(Prin);
    4399   def E2=std(imap(RRH,E));
    4400   def d=dim(E2);
    4401   setring RRH;
    4402   return(d);
    4403 }
     4524//  DimPar
     4525//  Auxilliary routine to NorSing determining the dimension of a parametric ideal
     4526//  Does not use @P and define it directly because is assumes that
     4527//                              variables and parameters have been inverted
     4528 static proc DimPar(ideal E)
     4529 {
     4530   def RRH=basering;
     4531   def RHx=ringlist(RRH);
     4532   def Prin=ring(RHx[1]);
     4533   setring(Prin);
     4534   def E2=std(imap(RRH,E));
     4535   def d=dim(E2);
     4536   setring RRH;
     4537   return(d);
     4538 }
    44044539
    44054540// locus0(G): Private routine used by locus (the public routine), that
    44064541//                builds the diferent components.
    44074542// input:      The output G of the grobcov (in generic representation, which is the default option)
    4408 //       Options: The algorithm allows the following options ar pair of arguments:
     4543//       Options: The algorithm allows the following options as pair of arguments:
    44094544//                "movdim", d  : by default movdim is 2 but it can be set to other values
     4545//                    when locus is called by envelop then as default is uses d=dim @P
    44104546//                "version", v   :  There are two versions of the algorithm. ('version',1) is
    44114547//                 a full algorithm that always distinguishes correctly between 'Normal'
     
    44184554//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
    44194555//                 The mover coordinates are to be placed as the last variables, and by default,
    4420 //                 its number is 2. If one consider locus problems in higer dimensions, the number of
    4421 //                 mover coordinates (placed as the last variables) is to be given as an option.
     4556//                 its number is dim @P, but as option it can be set to other values.
    44224557// output:
    44234558//         list, the canonical P-representation of the Normal and Non-Normal locus:
     
    44364571  int t1=1; int t2=1; int i;
    44374572  def R=basering;
    4438   //if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
    4439   //setglobalrings();
    4440   // Options
     4573  def RL=ringlist(R);
     4574  if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
     4575  setglobalrings();
     4576  //Options
    44414577  list DD=#;
    4442   int moverdim=nvars(R);
     4578  ideal vmov;
     4579  int moverdim=size(ringlist(R)[1][2]);
     4580  int dimpar=size(RL[1][2]);
     4581  int dimvar=size(RL[2]);
     4582  int nv=nvars(R);
     4583  if(moverdim>nv){moverdim=nv;}
     4584  for(i=1;i<=moverdim;i++)
     4585  {
     4586    //string("T_moverdim=",moverdim,"T_nv=",nv,"T_var(",i,")=",var(i));
     4587    vmov[size(vmov)+1]=var(i+nv-moverdim);
     4588  }
    44434589  int version=0;
    4444   int nv=nvars(R);
    44454590  if(nv<4){version=1;}
    44464591  int comment=0;
     
    44494594  for(i=1;i<=(size(DD) div 2);i++)
    44504595  {
    4451     if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
     4596    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    44524597    if(DD[2*i-1]=="version"){version=DD[2*i];}
    44534598    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
     
    44564601  }
    44574602  list HHH;
    4458   if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1){return(HHH);}
    4459    list G1; list G2;
     4603  if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1)
     4604    {return(HHH);}
     4605  list G1; list G2;
    44604606  def G=GG;
    44614607  list Q1; list Q2;
     
    44654611  {
    44664612    attrib(G[i][1],"IsSB",1);
    4467     d=locusdim(G[i][2],moverdim);
     4613    d=dim(std(G[i][1]));
    44684614    if(d==0){G1[size(G1)+1]=G[i];}
    44694615    else
     
    44844630  for(i=1;i<=size(G1RP);i++)
    44854631  {
    4486      kill B;
     4632    kill B;
    44874633     ideal B;
    44884634     for(k=1;k<=size(G1RP[i][3]);k++)
     
    44964642       P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B);
    44974643     }
    4498   }
     4644  }  //In P1RP the basis has been reduced wrt the top
    44994645  setring(R);
    45004646  ideal h;
     4647  list P1;
    45014648  if(t1)
    45024649  {
    4503     def P1=imap(@RP,P1RP);
     4650    P1=imap(@RP,P1RP);
    45044651    for(i=1;i<=size(P1);i++)
    45054652    {
     
    45144661      }
    45154662    }
    4516   }
    4517   else{list P1;}
     4663  } //In P1 factors in the basis independent of the parameters have been obtained
    45184664  ideal BB; int dd; list NS;
    45194665  for(i=1;i<=size(P1);i++)
    45204666  {
    4521        NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD);
    4522        dd=NS[1];
    4523        if(dd==0){P1[i][3]=NS;} //"Special";
    4524        else{P1[i][3]="Normal";}
     4667     NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD);
     4668     dd=NS[1];
     4669     if(dd==0){P1[i][3]=NS;} //"Special";
     4670     else{P1[i][3]="Normal";}
    45254671  }
    45264672  list P2;
     
    45484694    def C2=imap(R,P2);
    45494695    def L2=AddLocus(C2);
    4550  }
     4696  }
    45514697  else{list L2; list C2; kill P2; list P2;}
    4552     for(i=1;i<=size(L2);i++)
    4553     {
    4554       J=std(L2[i][2]);
    4555       d=dim(J); // AQUI
    4556       if(d==0)
    4557       {
    4558         L2[i][4]=string("Accumulation",L2[i][4]);
    4559       }
    4560       else{L2[i][4]=string("Degenerate",L2[i][4]);}
    4561     }
     4698  for(i=1;i<=size(L2);i++)
     4699  {
     4700    J=std(L2[i][2]);
     4701    d=dim(J);
     4702    if(d+1==dimpar)
     4703    {
     4704      L2[i][4]=string("Degenerate",L2[i][4]);
     4705    }
     4706    else{L2[i][4]=string("Accumulation",L2[i][4]);}
     4707  }
    45624708  list LN;
    45634709  if(t1==1)
     
    45694715    for(i=1;i<=size(L2);i++){LN[size(LN)+1]=L2[i];}
    45704716  }
     4717  int tLN=1;
     4718  if(size(LN)==0){tLN=0;}
    45714719  setring(R);
    4572   def L=imap(@P,LN);
    4573   for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
    4574   //if(te==0){kill @R; kill @RP; kill @P;}
    4575   list LL;
    4576   for(i=1;i<=size(L);i++)
    4577   {
    4578     if(typeof(L[i][4])=="list") {L[i][4][1]="Special";}
    4579     l[1]=L[i][2];
    4580     l[2]=L[i][3];
    4581     l[3]=L[i][4];
    4582     l[4]=L[i][5];
    4583     L[i]=l;
    4584  }
    4585  return(L);
     4720  if(tLN)
     4721  {
     4722   def L=imap(@P,LN);
     4723    for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
     4724    list LL;
     4725    for(i=1;i<=size(L);i++)
     4726    {
     4727      if(typeof(L[i][4])=="list") {L[i][4][1]="Special";}
     4728      l[1]=L[i][2];
     4729      l[2]=L[i][3];
     4730      l[3]=L[i][4];
     4731      l[4]=L[i][5];
     4732      L[i]=l;
     4733    }
     4734  }
     4735  else{list L;}
     4736  return(L);
    45864737}
    45874738
     
    45934744//               The Normal locus has two kind of components: Normal and Special.
    45944745//               The Non-normal locus has two kind of components: Accumulation and Degenerate.
    4595 //               Normal components: for each point in the component,
    4596 //               the number of solutions in the variables is finite, and
    4597 //               the solutions depend on the point in the component if the component is not 0-dimensional.
    4598 //               Special components:  for each point in the component,
    4599 //               the number of solutions in the variables is finite,
    4600 //               the component is not 0-dimensional, but the solutions do not depend on the
    4601 //               values of the parameters in the component.
    4602 //               Accumlation points: are 0-dimensional components for which it exist
    4603 //               an infinite number of solutions.
    4604 //               Degenerate components: are components of dimension greater than 0 for which
    4605 //               for every point in the component there exist infinite solutions.
     4746//          Normal component:
     4747//           - the component has non-zero antiimage
     4748//           - each point in the component has 0-dimensional antiimage.
     4749//          Special component:
     4750//           - Non-zero dimensional component whose antiimage is 0-dimensional.
     4751//          Accumulation points:
     4752//           - Zero-dimesnional  component whose antiimage is non-zero-dimensional.
     4753//          Degenerate components:
     4754//           - Non-zero-dimensional component
     4755//           - each point in the component has non-zero-dimensional antiimage.
    46064756//          The output components are given as
    46074757//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
     
    46104760//               gives the depth of the component of the constructible set.
    46114761proc locus(list GG, list #)
    4612 "USAGE: locus(G)
     4762"USAGE: locus(G[,options]);
     4763        Calling sequence:
     4764              locus(grobcov(S));
    46134765        The input must be the grobcov  of a parametrical ideal in Q[a][x],
    4614         (a=parameters, x=variables). In fact a must be the tracer coordinates
    4615         and x the mover coordinates and other auxiliary ones.
     4766        (a=parameters, x=variables). In practice a must be the tracer coordinates
     4767        and x the mover coordinates and remaining auxiliary variables.
     4768        For a parametric locus computation 'a' can contain some extra parameters.
    46164769        Special routine for determining the locus of points
    46174770        of  geometrical constructions. Given a parametric ideal J
    46184771        representing the system determining the locus of points (a)
    4619         who verify certain properties, the call to locus on the output of grobcov( J )
     4772        who verify certain properties, the call to locus on the output of grobcov(J)
    46204773        determines the different classes of locus components, following
    46214774        the taxonomy defined in
     
    46234776        \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\".
    46244777        Computer-Aided Design 56 (2014) 22-33.
    4625         The components can be Normal, Special, Accumulation, Degenerate.
     4778        The components can be Normal, Special, Accumulation or Degenerate.
    46264779OPTIONS: The algorithm allows the following options as pair of arguments:
    4627          \"movdim\", d  : by default movdim is 2 but it can be set to other values,
    4628          and represents the number of mever variables. they should be given as
    4629          the last variables of the ring.
     4780         \"vmov\", ideal(mover variables)  : by default vmov are the last n variables,
     4781         where n is the number of parameters of the ring (tracer variables plus extra
     4782         parameters). Thus, if the mover coordinates are not indicated, locus
     4783         algorithm will assume that they are the last n ring variables.
     4784         When locus is called internallt by envelop, by default, the mover variables
     4785         are assumed to be all the ring variables.
    46304786         \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is
    46314787         a full algorithm that always distinguishes correctly between 'Normal'
    4632          and 'Special' components, whereas \("version\",0) can decalre a component
    4633          as 'Normal' being really 'Special', but is more effective. By default (\"version\",1)
    4634          is used when the number of variables is less than 4 and 0 if not.
     4788         and 'Special' components, whereas \("version\",0) can declare a component
     4789         to be 'Normal' being in fact 'Special', but it is more effective. By default,
     4790         (version,1) is used when the number of variables is less than 4 and 0 if not.
    46354791         The user can force to use one or other version, but it is not recommended.
    4636          \"system\", ideal F: if the initial system is passed as an argument. This is actually not used.
    46374792         \"comments\", c: by default it is 0, but it can be set to 1.
    46384793         Usually locus problems have mover coordinates, variables and tracer coordinates.
    4639          The mover coordinates are to be placed as the last variables, and by default,
    4640          its number is 2. If one consider locus problems in higer dimensions, the number of
    4641          mover coordinates (placed as the last variables) is to be given as an option.
    4642 RETURN: The locus. The output is a list of the components ( C_1,.. C_n ) where
    4643         C_i =  (p_i,(p_i1,..p_is_i), type_i,l evel_i )   and type_i can be
    4644         'Normal', 'Special', Accumulation', 'Degenerate'. The 'Special' components
    4645         return more information, namely the antiimage of the component, that
    4646         is 0-dimensional for these kind of components.
    4647         Normal components: for each point in the component,
    4648         the number of solutions in the variables is finite, and
    4649         the solutions depend on the point in the component.
    4650         Special components:  for each point in the component,
    4651         the number of solutions in the variables is finite. The
    4652         antiimage of the component is 0-dimensional.
    4653         Accumlation points: are 0-dimensional components whose
    4654         antiimage is not zero-dimansional.
    4655         Degenerate components: are components of dimension greater than 0
    4656         whose antiimage is not-zero-diemansional.
    4657         The components are given in canonical P-representation.
    4658         The levels of a class of locus are 1,
    4659         because they represent locally closed. sets.
    4660 NOTE: It can only be called after computing the grobcov of the
    4661       parametrical ideal in generic representation ('ext',0),
    4662       which is the default.
    4663       The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
    4664 KEYWORDS: geometrical locus, locus, loci.
     4794         Example of option call:
     4795         > locus(S,\"version\",1,\"vmov\",ideal(x1,y1))
     4796RETURN:  The output is a list of the components (C_1, .. , C_n ) of the locus.
     4797          The locus is divided into two class of subsets: the normal and the non-normal
     4798          locus.
     4799            The Normal locus has two kind of components: Normal and Special.
     4800            The Non-normal locus has two kind of components: Accumulation and Degenerate.
     4801         Each component is given by
     4802              Ci=((pi,(pi1,..pis_i),type_i,level_i)
     4803         where
     4804           the first element is the canonical P-representation of the subset.
     4805           The type can be : \"Normal\", \"Special\", \"Accumulation\", \"Degenerate\".
     4806         Normal component:
     4807           - the component has non-zero antiimage
     4808           - each point in the component has 0-dimensional antiimage.
     4809         Special component:
     4810           - Non-zero dimensional component whose antiimage is 0-dimensional.
     4811           The Special components return more information, namely the antiimage of
     4812           the component, that is 0-dimensional for these kind of components.
     4813         Accumulation points:
     4814           - Zero-dimensional  component whose anti-image is non-zero-dimensional.
     4815         Degenerate components:
     4816           - Non-zero-dimensional component
     4817           - each point in the component has non-zero-dimensional anti-image.
     4818        The level is the depth of the segment of the constructible locus subset
     4819        (normal and non-normal subsets).
     4820        If all levels of a locus are 1, then both subsets are locally closed.
     4821NOTE: The input must be the grobcov of the locus system in generic
     4822        representation ('ext',0), which is the default.
     4823KEYWORDS: geometrical locus; locus; loci.
    46654824EXAMPLE: locus; shows an example"
    46664825{
     
    46694828  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
    46704829  setglobalrings();
    4671   // Options
     4830  //Options
    46724831  list DD=#;
    4673   int moverdim=nvars(R);
     4832  ideal vmov;
     4833  int moverdim=size(ringlist(R)[1][2]);
     4834  int nv=nvars(R);
     4835  if(moverdim>nv){moverdim=nv;}
     4836  for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);}
    46744837  int version=0;
    4675   int nv=nvars(R);
    46764838  if(nv<4){version=1;}
    46774839  int comment=0;
     
    46794841  for(i=1;i<=(size(DD) div 2);i++)
    46804842  {
    4681     if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
     4843    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    46824844    if(DD[2*i-1]=="version"){version=DD[2*i];}
    4683     if(DD[2*i-1]=="system"){Fm=DD[2*i];}
    46844845    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    4685     if(DD[2*i-1]=="family"){poly F=DD[2*i];}
    4686   }
    4687   int j; int k;
     4846  }
     4847  DD=list("vmov",vmov,"version",version,"comment",comment);
     4848  int j; int k; int te;
    46884849  def B0=GG[1][2];
    46894850  def H0=GG[1][3][1][1];
    4690   if (equalideals(B0,ideal(1)) or (equalideals(H0,ideal(0))==0))
     4851  if (equalideals(B0,ideal(1)) )
    46914852  {
    46924853    return(locus0(GG,DD));
     
    46954856  {
    46964857    int n=nvars(R);
    4697     ideal vmov=var(n-1),var(n);
     4858    ideal vB;
    46984859    ideal N;
    4699     intvec xw; intvec yw;
    4700     for(i=1;i<=n-1;i++){xw[i]=0;}
    4701     xw[n]=1;
    4702     for(i=1;i<=n;i++){yw[i]=0;}
    4703     yw[n-1]=1;
    4704     poly px; poly py;
    4705     int te=1;
    4706     i=1;
    4707     while( te and i<=size(B0))
    4708     {
    4709       if((deg(B0[i],xw)==1) and (deg(B0[i])==1)){px=B0[i]; te=0;}
    4710       i++;
    4711     }
    4712     i=1; te=1;
    4713     while( te and i<=size(B0))
    4714     {
    4715       if((deg(B0[i],yw)==1) and (deg(B0[i])==1)){py=B0[i]; te=0;}
    4716       i++;
    4717     }
    4718     N=px,py;
     4860    for(i=1;i<=size(B0);i++)
     4861    {
     4862      if(subset(variables(B0[i]),vmov)){N[size(N)+1]=B0[i];}
     4863    }
    47194864    te=indepparameters(N);
    47204865    if(te)
    47214866    {
    4722       string("locus detected that the mover must avoid point (",N,") in order to obtain the correct locus");
    4723       // eliminates segments of GG where N is contained in the basis
     4867      string("locus detected that the mover must avoid points (",N,") in order to obtain the correct locus");" ";
     4868      //eliminates segments of GG where N is contained in the basis
    47244869      list nGP;
    47254870      def GP=GG;
     
    47374882      }
    47384883     }
    4739     else{"Warning! Problem with more than one mover. Not able to solve it."; list L; return(L);}
     4884    else
     4885    {
     4886      " ";string("Warning! Problem with more than one mover.");
     4887      string("Try option 'vmov',ideal(of mover variables) to avoid some point of the mover");
     4888      " ";"Elements of the basis of the generic segment in mover variables=";  N;" ";
     4889      list L; return(L);
     4890    }
    47404891  }
    47414892  def LL=locus0(nGP,DD);
     
    47514902  // Concoid
    47524903  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
    4753   // System S96=
    47544904  S96;
     4905
    47554906  locus(grobcov(S96));
    4756   kill R;
    4757   // ********************************************
    4758   ring R=(0,a,b),(x4,x3,x2,x1),dp;
    4759   short=0;
    4760   ideal S=(x1-3)^2+(x2-1)^2-9,
    4761                (4-x2)*(x3-3)+(x1-3)*(x4-1),
    4762                (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
    4763                (4-x4)*a+(x3-3)*b+3*x4-4*x3,
    4764                (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
    4765   // System S=
    4766   S;
    4767   locus(grobcov(S));
    4768   kill R;
    4769   //********************************************
    4770 
    4771   ring R=(0,x,y),(x1,x2),dp;
    4772   short=0;
    4773   ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
    4774                 (x1 - 5)^2 + (x2 - 2)^2 - 4,
    4775                 -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
    4776   // System S=
    4777   S;
    4778   locus(grobcov(S));
    47794907}
    47804908
    47814909//  locusdg(G):  Special routine for determining the locus of points
    4782 //                 of  geometrical constructions. Given a parametric ideal J with
    4783 //                 parameters (a_1,..a_m) and variables (x_1,..,xn),
    4784 //                 representing the system determining
    4785 //                 the locus of points (a_1,..,a_m) who verify certain
    4786 //                 properties, computing the grobcov G of
    4787 //                 J and applying to it locus, determines the different
    4788 //                 classes of locus components. The components can be
    4789 //                 Normal, Special, Accumulation point, Degenerate.
    4790 //                 The output are the components given in P-canonical form
    4791 //                 of at most 4 constructible sets: Normal, Special, Accumulation,
    4792 //                 Degenerate.
    4793 //                 The description of the algorithm and definitions is
    4794 //                 given in a forthcoming paper by Abanades, Botana, Montes, Recio.
    4795 //                 Usually locus problems have mover coordinates, variables and tracer coordinates.
    4796 //                 The mover coordinates are to be placed as the last variables, and by default,
    4797 //                 its number is 2. If onw consider locus problems in higer dimensions, the number of
    4798 //                 mover coordinates (placed as the last variables) is to be given as an option.
    4799 //
     4910//                 of  geometrical constructions in Dynamic Geometry.
     4911//                 It is to be applied to the output of locus and selects
     4912//                 as 'Relevant' the 'Normal' and the 'Accumulation'
     4913//                 components.
    48004914//  input:      The output of locus(G);
    48014915//  output:
    4802 //          list, the canonical P-representation of the Normal and Non-Normal locus:
    4803 //               The Normal locus has two kind of components: Normal and Special.
    4804 //               The Non-normal locus has two kind of components: Accumulation and Degenerate.
    4805 //               Normal components: for each point in the component,
    4806 //               the number of solutions in the variables is finite, and
    4807 //               the solutions depend on the point in the component if the component is not 0-dimensional.
    4808 //               Special components:  for each point in the component,
    4809 //               the number of solutions in the variables is finite,
    4810 //               the component is not 0-dimensional, but the solutions do not depend on the
    4811 //               values of the parameters in the component.
    4812 //               Accumlation points: are 0-dimensional components for which it exist
    4813 //               an infinite number of solutions.
    4814 //               Degenerate components: are components of dimension greater than 0 for which
    4815 //               for every point in the component there exist infinite solutions.
     4916//          list, the canonical P-representation of the 'Relevant' components of the locus.
    48164917//          The output components are given as
    48174918//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
     
    48204921//               gives the depth of the component of the constructible set.
    48214922proc locusdg(list L)
    4822 "USAGE: locusdg(L)   The call must be locusdg(locus(grobcov(S))).
     4923"USAGE: locusdg(L);
     4924           Calling sequence:
     4925                locusdg(locus(grobcov(S))).
    48234926RETURN: The output is the list of the 'Relevant' components of the locus
    48244927           in Dynamic Geometry:(C1,..,C:m), where
     
    48264929           The 'Relevant' components are the 'Normal' and 'Accumulation' components
    48274930           of the locus. (See help for locus).
    4828 
    4829 NOTE: It can only be called after computing the locus.
    4830            Calling sequence:    locusdg(locus(grobcov(S)));
    4831 KEYWORDS: geometrical locus, locus, loci, dynamic geometry
    4832 EXAMPLE: locusdg; shows an example"
     4931KEYWORDS: geometrical locus; locus; loci; dynamic geometry
     4932EXAMPLE: example locusdg; shows an example"
    48334933{
    48344934  list LL;
     
    48514951  // Concoid
    48524952  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
    4853   // System S96=
    48544953  S96;
    48554954  locus(grobcov(S96));
    48564955  locusdg(locus(grobcov(S96)));
    4857   kill R;
    4858   //********************************************
    4859   ring R=(0,a,b),(x4,x3,x2,x1),dp;
    4860   ideal S=(x1-3)^2+(x2-1)^2-9,
    4861                (4-x2)*(x3-3)+(x1-3)*(x4-1),
    4862                (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
    4863                (4-x4)*a+(x3-3)*b+3*x4-4*x3,
    4864                (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
    4865   short=0;
    4866   locus(grobcov(S));
    4867   locusdg(locus(grobcov(S)));
    4868   kill R;
    4869   //********************************************
    4870 
    4871   ring R=(0,x,y),(x1,x2),dp;
    4872   short=0;
    4873   ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
    4874                 (x1 - 5)^2 + (x2 - 2)^2 - 4,
    4875                 -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
    4876   locus(grobcov(S));
    4877   locusdg(locus(grobcov(S)));
    4878 }
    4879 
    4880 // locusto: Transforms the output of locus, locusdg, envelop and envelopdg
     4956}
     4957
     4958// locusto: Transforms the output of locus, locusdg, envelop
    48814959//             into a string that can be reed from different computational systems.
    48824960// input:
    4883 //     list L: The output of locus
     4961//     list L: The output of locus or locusdg or envelop.
    48844962// output:
    4885 //     string s: The output of locus converted to a string readable by other programs
     4963//     string s: Converts the input into a string readable by other programs
    48864964proc locusto(list L)
    48874965"USAGE: locusto(L);
    48884966          The argument must be the output of locus or locusdg or
    4889           envelop or envelopdg.
     4967          envelop.
    48904968          It transforms the output into a string in standard form
    4891           readable in many languages (Geogebra).
     4969          readable in other languages, not only Singular (Geogebra).
    48924970RETURN: The locus in string standard form
    48934971NOTE: It can only be called after computing either
     
    48954973              - locusdg(locus(grobcov(F)))  -> locusto( locusdg(locus(grobcov(F))) )
    48964974              - envelop(F,C)                       -> locusto( envelop(F,C) )
    4897               -envelopdg(envelop(F,C))      -> locusto( envelopdg(envelop(F,C)) )
    4898 KEYWORDS: geometrical locus, locus, loci.
    4899 EXAMPLE:  locusto; shows an example"
     4975KEYWORDS: geometrical locus; locus; loci
     4976EXAMPLE:  example locusto; shows an example"
    49004977{
    49014978  int i; int j; int k;
     
    49615038  locusto(locus(grobcov(S)));
    49625039  locusto(locusdg(locus(grobcov(S))));
    4963   kill R;
    4964   //********************************************
    4965 
    4966   // 1. Take a fixed line l: x1-y1=0  and consider
    4967   //    the family F of a lines parallel to l passing through the mover point M
    4968   // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
    4969   // 3. Compute the envelop of the family of lines.
    4970 
    4971   ring R=(0,x,y),(x1,y1),lp;
    4972   poly F=(y-y1)-(x-x1);
    4973   ideal C=x1^2+y1^2-25;
    4974   short=0;
    4975 
    4976   // Curves Family F=
    4977   F;
    4978   // Conditions C=
    4979   C;
    4980 
    4981   locusto(envelop(F,C));
    4982   locusto(envelopdg(envelop(F,C)));
    4983   kill R;
    4984   //********************************************
    4985 
    4986   // Steiner Deltoid
    4987   // 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
    4988   // 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
    4989   // 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
    4990   // 4. Obtain the envelop of the lines above.
    4991 
    4992   ring R=(0,x,y),(x1,y1,x2,y2),lp;
    4993   ideal C=(x1)^2+(y1)^2-1,
    4994                x2+y2-1,
    4995                x2-y2-x1+y1;
    4996   matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
    4997   poly F=det(M);
    4998 
    4999   short=0;
    5000 
    5001   // Curves Family F=
    5002   F;
    5003   // Conditions C=
    5004   C;
    5005 
    5006   locusto(envelop(F,C));
    5007   locusto(envelopdg(envelop(F,C)));
    5008 }
    5009 
    5010 // Auxiliary routine
    5011 // locusdim
    5012 // input:
    5013 //    B:         ideal, a basis of a segment of the grobcov
    5014 //    dgdim: int, the dimension of the mover (for locus)
    5015 //                by default dgdim is equal to the number of variables
    5016 static proc locusdim(ideal B, list #)
    5017 {
    5018   def R=basering;
    5019   int dgdim;
    5020   int nv=nvars(R);
    5021   if (size(#)>0){dgdim=#[1];}
    5022   else {dgdim=nv;}
    5023   int d;
    5024   list v;
    5025   ideal vi;
    5026   int i;
    5027   for(i=1;i<=dgdim;i++)
    5028   {
    5029     v[size(v)+1]=varstr(nv-dgdim+i);
    5030     vi[size(v)+1]=var(nv-dgdim+i);
    5031   }
    5032   ideal B0;
    5033   for(i=1;i<=size(B);i++)
    5034   {
    5035     if(subset(variables(B[i]),vi))
    5036     {
    5037       B0[size(B0)+1]=B[i];
    5038     }
    5039   }
    5040   def RR=ringlist(R);
    5041   def RR0=RR;
    5042   RR0[2]=v;
    5043   def R0=ring(RR0);
    5044   setring(R0);
    5045   def B0r=imap(R,B0);
    5046   B0r=std(B0r);
    5047   d=dim(B0r);
    5048   setring R;
    5049   return(d);
    5050 }
    5051 
    5052 static proc norspec(ideal F)
    5053 {
    5054   def RR=basering;
    5055   def Rx=ringlist(RR);
    5056 
    5057    def Rx=ringlist(RR);
    5058   def @P=ring(Rx[1]);
    5059   list Lx;
    5060   Lx[1]=0;
    5061   Lx[2]=Rx[2]+Rx[1][2];
    5062   Lx[3]=Rx[1][3];
    5063   Lx[4]=Rx[1][4];
    5064   Rx[1]=0;
    5065   def D=ring(Rx);
    5066   def @RP=D+@P;
    5067   exportto(Top,@R);      // global ring Q[a][x]
    5068   exportto(Top,@P);      // global ring Q[a]
    5069   exportto(Top,@RP);     // global ring K[x,a] with product order
    5070   setring(RR);
    50715040}
    50725041
    50735042// envelop
    5074 // Input: n arguments
    5075 //   poly F: the polynomial defining the family of curves in ring R=0,(x,y),(x1,..,xn),lp;
    5076 //   ideal C=g1,..,g_{n-1}:  the set of constraints
     5043// Input:
     5044//   poly F: the polynomial defining the family of hypersurfaces in ring R=0,(x_1,..,x_n),(u_1,..,u_m),lp;
     5045//   ideal C=g1,..,g_{n-1}:  the set of constraints;
     5046//   options.
    50775047// Output: the components of the envolvent;
    50785048proc envelop(poly F, ideal C, list #)
    50795049"USAGE: envelop(F,C);
    5080           The first argument F must be the family of curves of which
    5081           on want to compute the envelop.
    5082           The second argument C must be the ideal of conditions
    5083           over the variables, and should contain as polynomials
    5084           as the number of variables -1.
     5050          The first argument poly F must be the family of hypersurfaces for which
     5051          on want to compute its envelop.
     5052          The second argument C must be the ideal of restrictions on
     5053          the variables, and should contain less polynomials
     5054          than the number of variables
     5055          (x_1,..,x_n) are the variables of the hypersurfaces of F, that are considered
     5056          as parameters of the parametric ring.
     5057          (u_1,..,u_m) are the parameteres of the hypersurfaces, that are considered
     5058          as variables of the parametric ring.
     5059          Calling sequence:
     5060              ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
     5061              poly F=F(x_1,..,x_n,u_1,..,u_m);
     5062              ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
     5063              envelop(F,C,#options);
     5064           where s<m.
     5065OPTIONS: The algorithm allows the following options as pair of arguments:
     5066         \"vmov\", ideal(mover variables)  : by default vmov are  u_1,..,u_m.
     5067         But it can be restricted by the user to the more convenient ones.
     5068         \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is
     5069         a full algorithm that always distinguishes correctly between 'Normal'
     5070         and 'Special' components, whereas \("version\",0) can decalre a component
     5071         as 'Normal' being really 'Special', but is more effective. By default (\"version\",1)
     5072         is used when the number of variables is less than 4 and 0 if not.
     5073         The user can force to use one or other version, but it is not recommended.
     5074         \"comments\", c: by default it is 0, but it can be set to 1.
    50855075RETURN: The components of the envelop with its taxonomy:
    5086            The taxonomy distinguishes 'Normal',
    5087           'Special', 'Accumulation', 'Degenerate' components.
    5088           In the case of 'Special' components, it also
    5089           outputs the antiimage of the component
    5090           and an integer (0-1). If the integer is 0
    5091           the component is not a curve of the family and is
    5092           not considered as 'Relevant' by the envelopdg routine
    5093           applied to it, but is considered as 'Relevant' if the integer is 1.
     5076           (see locus help). envelop uses locus.
     5077           The taxonomy distinguishes \"Normal\",
     5078          \"Special\", \"Accumulation\", \"Degenerate\" components.
     5079          In the case of \"Special\" components, it also
     5080          outputs the anti-image of the component
    50945081NOTE: grobcov is called internally.
    50955082          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
    5096 KEYWORDS: geometrical locus, locus, loci, envelop
    5097 EXAMPLE:  envelop; shows an example"
    5098 {
    5099   int tes=0; int i;   int j;
     5083KEYWORDS: geometrical locus; locus; loci; envelop
     5084EXAMPLE:  example envelop; shows an example"
     5085{
    51005086  def R=basering;
     5087  int tes=0; int i;   int j;  int k; int m;
     5088  int d;
     5089  int dp;
     5090  ideal BBB;
    51015091  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
    51025092  setglobalrings();
    5103   // Options
     5093  //Options
    51045094  list DD=#;
    5105   int moverdim=nvars(R);
     5095  ideal vmov;
     5096  int nv=nvars(R);
     5097  for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);}
     5098  int numpars=size(ringlist(R)[1][2]);
    51065099  int version=0;
    5107   int nv=nvars(R);
    51085100  if(nv<4){version=1;}
    51095101  int comment=0;
     5102  int familyinfo;
    51105103  ideal Fm;
    51115104  for(i=1;i<=(size(DD) div 2);i++)
    51125105  {
    5113     if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
     5106    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    51145107    if(DD[2*i-1]=="version"){version=DD[2*i];}
    5115     if(DD[2*i-1]=="system"){Fm=DD[2*i];}
    51165108    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    5117   }
    5118   int n=nvars(R);
    5119   list v;
    5120   for(i=1;i<=n;i++){v[size(v)+1]=var(i);}
    5121   def MF=jacob(F);
    5122   def TMF=transpose(MF);
    5123   def Mg=MF;
    5124   def TMg=TMF;
    5125   for(i=1;i<=n-1;i++)
    5126   {
    5127     Mg=jacob(C[i]);
    5128     TMg=transpose(Mg);
    5129     TMF=concat(TMF,TMg);
    5130   }
    5131   poly J=det(TMF);
    5132   ideal S=ideal(F)+C+ideal(J);
    5133   DD[size(DD)+1]="family";
    5134   DD[size(DD)+1]=F;
     5109    if(DD[2*i-1]=="familyinfo"){familyinfo=DD[2*i];}
     5110  };
     5111  DD=list("vmov",vmov,"version",version,"comment",comment);
     5112  int ng=size(C);
     5113  ideal S=F;
     5114  for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];}
     5115  int s=nv-ng;
     5116  if(s>0)
     5117  {
     5118    matrix M[ng+1][ng+1];
     5119    def cc=comb(nv,ng+1);
     5120    poly J;
     5121    for(k=1;k<=size(cc);k++)
     5122    {
     5123      for(j=1;j<=ng+1;j++)
     5124      {
     5125        M[1,j]=diff(F,var(cc[k][j]));
     5126      }
     5127      for(i=1;i<=ng;i++)
     5128      {
     5129        for(j=1;j<=ng+1;j++)
     5130        {
     5131          M[i+1,j]=diff(C[i],var(cc[k][j]));
     5132        }
     5133      }
     5134      J=det(M);
     5135      S[size(S)+1]=J;
     5136    }
     5137 }
     5138 if(comment>0){"System S before grobcov ="; S;}
    51355139  def G=grobcov(S,DD);
    5136    def L=locus(G, DD);
    5137   return(L);
     5140  list HHH;
     5141  if (G[1][1][1]==1 and G[1][2][1]==1 and G[1][3][1][1][1]==0 and G[1][3][1][2][1]==1)
     5142  {return(HHH);}
     5143   //DD[size(DD)+1]="vmov";
     5144   //DD[size(DD)+1]=4;
     5145  def L=locus(G,DD);
     5146  list GL;
     5147  ideal fam; ideal env;
     5148
     5149   def Rx=ringlist(R);
     5150   def P=ring(Rx[1]);
     5151   list Lx;
     5152   Lx[1]=0;
     5153   Lx[2]=Rx[2]+Rx[1][2];
     5154   Lx[3]=Rx[1][3];
     5155   Lx[4]=Rx[1][4];
     5156   Rx[1]=0;
     5157   def D=ring(Rx);
     5158   def RP=P+D;
     5159  list LL;
     5160  list NormalComp;
     5161  ideal Gi;
     5162  ideal BBBB;
     5163  poly B0;
     5164  if(familyinfo==1)
     5165  {
     5166    for(i=1;i<=size(L);i++)
     5167    {
     5168      if(typeof(L[i][3])=="string")
     5169      {
     5170        if(L[i][3]=="Normal")
     5171        {
     5172          NormalComp[size(NormalComp)+1]=L[i][1];
     5173          Gi=S;
     5174          Gi[size(Gi)+1]=L[i][1][1];
     5175          //"T_grobcov(Gi)="; grobcov(Gi);
     5176          kill HHH; list HHH;
     5177          //"T_L[i]="; L[i];
     5178          //d=DimPar(L[i][1]);
     5179          if(defined(SL)){kill SL;}
     5180          def SL=C;
     5181          SL[size(SL)+1]=F;
     5182          for(j=1;j<=size(L[i][1]);j++)
     5183          {
     5184            SL[size(SL)+1]=L[i][1][j];
     5185          }
     5186          setring RP;
     5187          if(defined(BBBB)){kill BBBB;}
     5188          ideal BBBB;
     5189          if(defined(BB)){kill BB;}
     5190          def BB=imap(R,SL);
     5191          if(defined(B0)){kill B0;}
     5192          poly B0;
     5193          if(defined(LLL)){kill LLL;}
     5194          def LLL=imap(R,L);
     5195          //"T_BB="; BB;
     5196          BB=std(BB);
     5197           for(j=1;j<=size(BB);j++)
     5198           {
     5199             B0=reduce(BB[j],LLL[i][1]);
     5200             if(not(B0==0)){BBBB[size(BBBB)+1]=B0;}
     5201           }
     5202          setring R;
     5203          BBB=imap(RP,BBBB);
     5204           L[i][5]=BBB;
     5205        }
     5206      }
     5207    }
     5208    LL[1]=L;
     5209    LL[2]=NormalComp;
     5210    list LLL; list LLLL;
     5211    int t;
     5212    for(k=1;k<=size(LL[2]);k++)
     5213    {
     5214      for(i=1;i<=size(G);i++)
     5215      {
     5216        j=1; t=0;
     5217        while(t==0 and j<=size(G[i][3]))
     5218        {
     5219          //"T_LL[2][k]="; LL[2][k];
     5220          //"T_G[i][3][j][1]="; G[i][3][j][1];
     5221          if(equalideals(LL[2][k],G[i][3][j][1]))
     5222          {
     5223            LLL[size(LLL)+1]=list(k,i,j);
     5224            t=1;
     5225          }
     5226          j++;
     5227        }
     5228      }
     5229    }
     5230    LL[3]=LLL;
     5231
     5232    for(k=1;k<=size(LLL);k++)
     5233    {
     5234      for(m=k+1;m<=size(LLL);m++)
     5235      {
     5236        for(i=1;i<=size(G[LLL[k][2]][3][LLL[k][3]][2]);i++)
     5237        {
     5238          for(j=1;j<=size(G[LLL[m][2]][3][LLL[m][3]][2]);j++)
     5239          {
     5240            //string("T_(G[",LLL[k][2],"][3][",LLL[k][3],"][2][",i,"])=");  G[LLL[k][2]][3][LLL[k][3]][2][i];
     5241            //string("T_(G[",LLL[m][2],"][3][",LLL[m][3],"][2][",j,"])=");  G[LLL[m][2]][3][LLL[m][3]][2][j];
     5242            if(equalideals( G[LLL[k][2]][3][LLL[k][3]][2][i] ,G[LLL[m][2]][3][LLL[m][3]][2][j]))
     5243            {
     5244              //"T_GGG="; LL[2][LLL[k][1]]; LL[2][LLL[m][1]]; G[LLL[k][2]][3][LLL[k][3]][2][i];
     5245              LLLL[size(LLLL)+1]=list(LL[2][LLL[k][1]],LL[2][LLL[m][1]], G[LLL[k][2]][3][LLL[k][3]][2][i]);
     5246            }
     5247          }
     5248        }
     5249      }
     5250    }
     5251    LL[3]=LLLL;
     5252  }
     5253  else{LL=L;}
     5254  return(LL);
    51385255}
    51395256example
     
    51475264
    51485265  ring R=(0,x,y),(x1,y1,x2,y2),lp;
     5266  short=0;
    51495267  ideal C=(x1)^2+(y1)^2-1,
    51505268               x2+y2-1,
     
    51525270  matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
    51535271  poly F=det(M);
    5154 
    5155   short=0;
    5156 
    5157   // Curves Family F=
     5272  // Curves Family F
    51585273  F;
    51595274  // Conditions C=
    51605275  C;
     5276  envelop(F,C);
     5277}
     5278
     5279proc AssocTanToEnv(poly F,ideal C, poly E, list #)
     5280"USAGE: AssocTanToEnv(F,C,E);
     5281          The first argument poly F must be the family of hypersurfaces for which
     5282          on want to compute its envelop.
     5283          The second argument C must be the ideal of restrictions on
     5284          the variables, and should contain s  polynomials
     5285          being s<n,
     5286          (x_1,..,x_n) are the variables of the hypersurfaces of F, that are considered
     5287          as parameters of the parametric ring.
     5288          (u_1,..,u_m) are the parameteres of the hypersurfaces, that are considered
     5289          as variables of the parametric ring.
     5290          The third parameter poly E must be the equation of a component of the
     5291          envelop of (F,C) previously determined by envelop.
     5292          Calling sequence:
     5293              ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
     5294              poly F=F(x_1,..,x_n,u_1,..,u_m);
     5295              ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
     5296              poly E(x_1,..,x_n);
     5297              AssocTanToEnv(F,C,E,#options);
     5298
     5299OPTIONS: The algorithm allows the following options as pair of arguments:
     5300         \"vmov\", ideal(mover variables)  : by default vmov are  u_1,..,u_m.
     5301         But it can be restricted by the user to the more convenient ones.
     5302         \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is
     5303         a full algorithm that always distinguishes correctly between 'Normal'
     5304         and 'Special' components, whereas \("version\",0) can decalre a component
     5305         as 'Normal' being really 'Special', but is more effective. By default (\"version\",1)
     5306         is used when the number of variables is less than 4 and 0 if not.
     5307         The user can force to use one or other version, but it is not recommended.
     5308         \"comments\", c: by default it is 0, but it can be set to 1.
     5309RETURN: The interesting segments of the grobcov each one with (lpp,basis,segment).
     5310        Fixing the values of (x_1,..,x_n) inside E, the basis allows to detemine the values of the parameters
     5311        u_1,..u_m.
     5312NOTE: grobcov is called internally.
     5313          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
     5314KEYWORDS: geometrical locus; locus; loci; envelop, associated tangent
     5315EXAMPLE:  example AssocTanToEnv; shows an example"
     5316{
     5317  def R=basering;
     5318  int tes=0; int i;   int j;  int k; int m;
     5319  int d;
     5320  int dp;
     5321  poly EE=E;
     5322  int moreinfo=1;
     5323  ideal BBB;
     5324  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
     5325  setglobalrings();
     5326  //Options
     5327  list DD=#;
     5328  ideal vmov;
     5329  int nv=nvars(R);
     5330  for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);}
     5331  int numpars=size(ringlist(R)[1][2]);
     5332  int version=0;
     5333  if(nv<4){version=1;}
     5334  int comment=0;
     5335  int familyinfo;
     5336  ideal Fm;
     5337  for(i=1;i<=(size(DD) div 2);i++)
     5338  {
     5339    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
     5340    if(DD[2*i-1]=="version"){version=DD[2*i];}
     5341    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
     5342    if(DD[2*i-1]=="familyinfo"){familyinfo=DD[2*i];}
     5343    if(DD[2*i-1]=="moreinfo"){moreinfo=DD[2*i];}
     5344  };
     5345  DD=list("vmov",vmov,"version",version,"comment",comment);
     5346  int ng=size(C);
     5347  ideal S=F;
     5348  for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];}
     5349  int s=nv-ng;
     5350  if(s>0)
     5351  {
     5352    matrix M[ng+1][ng+1];
     5353    def cc=comb(nv,ng+1);
     5354    poly J;
     5355    for(k=1;k<=size(cc);k++)
     5356    {
     5357      for(j=1;j<=ng+1;j++)
     5358      {
     5359        M[1,j]=diff(F,var(cc[k][j]));
     5360      }
     5361      for(i=1;i<=ng;i++)
     5362      {
     5363        for(j=1;j<=ng+1;j++)
     5364        {
     5365          M[i+1,j]=diff(C[i],var(cc[k][j]));
     5366        }
     5367      }
     5368      J=det(M);
     5369      S[size(S)+1]=J;
     5370    }
     5371 }
     5372 S[size(S)+1]=EE;
     5373 if(comment>0){"System S before grobcov ="; S;}
     5374  def G=grobcov(S,DD);
     5375  //"T_G=";G;
     5376  list GG;
     5377  for(i=2;i<=size(G);i++)
     5378  {
     5379    GG[size(GG)+1]=G[i];
     5380  }
     5381  G=GG;
     5382  //"T_G=";G;
     5383  if(moreinfo>0){return(G);}
     5384  else
     5385  {
     5386    int t=0;
     5387    ideal H;
     5388    i=1;
     5389    while(t==0 and i<=size(G))
     5390    {
     5391      //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1];
     5392      //string("T_EE="); EE;
     5393      if(G[i][3][1][1][1]==EE)
     5394      {
     5395         t=1;
     5396         H=G[i][2];
     5397      }
     5398      i++;
     5399    }
     5400    return(H);
     5401  }
     5402  return(G);
     5403}
     5404example
     5405{
     5406  "EXAMPLE:"; echo = 2;
     5407  ring R=(0,y0,x,y),(t),dp;
     5408  short=0;
     5409  poly F=(x-5*t)^2+y^2-3^2*t^2;
     5410  F;
     5411  ideal C;
     5412  C;
    51615413
    51625414  def Env=envelop(F,C);
    51635415  Env;
    5164 }
    5165 
    5166 // envelopdg
    5167 // Input: list L: the output of envelop(poly F, ideal C, list #)
    5168 // Output: the relevant components of the envolvent in dynamic geometry;
    5169 proc envelopdg(list L)
    5170 "USAGE: envelopdg(L);
    5171           The input list L must be the output of the call to
    5172           the routine 'envolop' of the family of curves
    5173 RETURN: The relevant components of the envelop in Dynamic Geometry.
    5174            'Normal' and 'Accumulation' components are always considered
    5175            'Relevant'. 'Special' components of the envelop outputs
    5176            three objects in its characterization: 'Special', the antiimage ideal,
    5177            and the integer 0 or 1, that indicates that the given component is
    5178            formed (1) or is not formed (0) by curves of the family. Only if yes,
    5179            'envelopdg' considers the component as 'Relevant' .
    5180 NOTE: It must be called to the output of the 'envelop' routine.
    5181           The basering R, must be of the form Q[a,b,..][x,y,..].
    5182 KEYWORDS: geometrical locus, locus, loci, envelop.
    5183 EXAMPLE:  envelop; shows an example"
    5184 {
    5185   list LL;
    5186   list Li;
     5416  // E is a component of the envelop:
     5417  poly E=Env[1][1][1];
     5418  E;
     5419  def A=AssocTanToEnv(F,C,E);
     5420  A;
     5421  // The basis of the parameter values of the associated tangent component is
     5422  A[1][2][1];
     5423  // Thus t=(5/12)*y0 and the associated tangent component at (x0,y0) is
     5424  subst(F,t,(5/12)*y0);
     5425}
     5426
     5427proc FamElemsAtEnvCompPoints(poly F,ideal C, poly E, list #)
     5428"USAGE: FamElemsAtEnvCompPoints(F,C,E[,options]);
     5429          The first argument poly F must be the family of hypersurfaces for which
     5430          on want to compute its envelop.
     5431          The second argument C must be the ideal of restrictions on
     5432          the variables, and should contain s  polynomials
     5433          being s<n,
     5434          (x_1,..,x_n) are the variables of the hypersurfaces of F, that are considered
     5435          as parameters of the parametric ring.
     5436          (u_1,..,u_m) are the parameteres of the hypersurfaces, that are considered
     5437          as variables of the parametric ring.
     5438          The third parameter poly E must be the equation of a component of the
     5439          envelop of (F,C) previously determined by envelop.
     5440          Calling sequence:
     5441              ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
     5442              poly F=F(x_1,..,x_n,u_1,..,u_m);
     5443              ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
     5444              poly E(x_1,..,x_n);
     5445              FamElemsAtEnvCompPoints(F,C,E,#options);
     5446
     5447OPTIONS: The algorithm allows the following options as pair of arguments:
     5448         \"vmov\", ideal(mover variables)  : by default vmov are  u_1,..,u_m.
     5449         But it can be restricted by the user to the more convenient ones.
     5450         \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is
     5451         a full algorithm that always distinguishes correctly between \"Normal\"
     5452         and \"Special\" components, whereas (\"version\",0) can declare a component
     5453         as \"Normal\" being really \"Special\", but is more effective. By default (\"version\",1)
     5454         is used when the number of variables is less than 4 and 0 if not.
     5455         The user can force to use one or other version, but it is not recommended.
     5456         \"comments\", c: by default it is 0, but it can be set to 1.
     5457RETURN: The interesting segments of the grobcov each one with (lpp,basis,segment).
     5458        Fixing the values of (x_1,..,x_n) inside E, the basis allows to detemine the values of the
     5459        parameters (u_1,..u_m).
     5460NOTE: grobcov is called internally.
     5461          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
     5462KEYWORDS: geometrical locus; locus; loci; envelop; associated tangent
     5463EXAMPLE: example FamElemsAtEnvCompPoints; shows an example"
     5464{
     5465  ideal S=C;
     5466  S[size(S)+1]=F;
     5467  S[size(S)+1]=E;
     5468  def G=grobcov(S);
     5469  list GG;
    51875470  int i;
    5188   for(i=1;i<=size(L);i++)
    5189   {
    5190     if(typeof(L[i][3])=="string")
    5191     {
    5192       if((L[i][3]=="Normal") or (L[i][3]=="Accumulation")){Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;}
    5193     }
     5471  for(i=2; i<=size(G); i++)
     5472  {
     5473    GG[size(GG)+1]=G[i];
     5474  }
     5475  return(GG);
     5476}
     5477example
     5478{
     5479  "EXAMPLE:"; echo = 2;
     5480  ring R=(0,y0,x,y),(t),dp;
     5481  short=0;
     5482  poly F=(x-5*t)^2+y^2-3^2*t^2;
     5483  F;
     5484  ideal C;
     5485  C;
     5486
     5487  def Env=envelop(F,C);
     5488  Env;
     5489
     5490  // E is a component of the envelop:
     5491  poly E=Env[1][1][1];
     5492  E;
     5493  def A=AssocTanToEnv(F,C,E);
     5494  A;
     5495  // The basis of the parameter values of the associated tangent component is
     5496  A[1][2][1];
     5497  // Thus t=(5/12)*y0 the assocoated tangent family element at (x0,y0) is
     5498  subst(F,t,(5/12)*y0);
     5499
     5500  FamElemsAtEnvCompPoints(F,C,E);
     5501  // Thus (12*t^2-5*y0)^2=0 and the unique circle of the family passing at (x0,y0) in E
     5502  // is the associated   tangent circle:
     5503  subst(F,t,(5/12)*y0);
     5504}
     5505
     5506// discrim
     5507proc discrim(poly F0, poly x0)
     5508"USAGE: discrim(f,x);
     5509          poly f: the polynomial in Q[a][x] or Q[x] of degree 2 in x
     5510          poly x: a variable in the ring.
     5511RETURN: the factorized discriminant of f wrt x for discussing its sign
     5512KEYWORDS: second degree; solving
     5513EXAMPLE:  discrim; shows an example"
     5514{
     5515  def RR=basering;
     5516  int i;
     5517  int te;
     5518  int d;  int dd;
     5519  if(size(ringlist(RR)[1])>0)
     5520  {
     5521    te=1;
     5522    setglobalrings();
     5523    setring @RP;
     5524    poly F=imap(RR,F0);
     5525    poly X=imap(RR,x0);
     5526  }
     5527  else
     5528  {poly F=F0; poly X=x0;}
     5529  matrix M=coef(F,X);
     5530  d=deg(M[1,1]);
     5531  if(d>2){"Degree is higher than 2. No discriminant"; setring RR; return();}
     5532    poly dis=(M[2,2])^2-4*M[2,1]*M[2,3];
     5533    def disp=factorize(dis,0);
     5534    if(te==0){return(disp);}
    51945535    else
    51955536    {
    5196       if(typeof(L[i][3])=="list")
    5197       {
    5198         if(L[i][3][3]==1)
    5199         {
    5200           Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;
    5201         }
    5202       }
    5203     }
    5204   }
    5205   return(LL);
     5537      setring RR;
     5538      def disp0=imap(@RP,disp);
     5539      return(disp0);
     5540    }
    52065541}
    52075542example
    52085543{
    52095544  "EXAMPLE:"; echo = 2;
    5210 
    5211   // 1. Take a fixed line l: x1-y1=0  and consider
    5212   //    the family F of a lines parallel to l passing through the mover point M
    5213   // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
    5214   // 3. Compute the envelop of the family of lines.
    5215 
    5216   ring R=(0,x,y),(x1,y1),lp;
    5217   short=0;
    5218   poly F=(y-y1)-(x-x1);
    5219   ideal C=x1^2+y1^2-25;
    5220   short=0;
    5221 
    5222   // Curves Family F=
    5223   F;
    5224   // Conditions C=
    5225   C;
    5226 
    5227   envelop(F,C);
    5228   envelopdg(envelop(F,C));
     5545  ring R=(0,a,b,c),(x,y),dp;
     5546  poly f=a*x^2*y+b*x*y+c*y;
     5547  discrim(f,x);
    52295548}
    52305549
    52315550// AddLocus: auxilliary routine for locus0 that computes the components of the constructible:
    52325551// Input:  the list of locally closed sets to be added, each with its type as third argument
    5233 //     L=[ [LC[11],,,LC[1k_1],  .., [LC[r1],,,LC[rk_r] ] where
     5552//     L=[ [LC[11],..,LC[1k_1],.., [LC[r1],..,LC[rk_r] ] where
    52345553//            LC[1]=[p1,[p11,..,p1k],typ]
    5235 // Output:  the list of components of the constructible union of the L, with the type of the corresponding top
     5554// Output:  the list of components of the constructible union of L, with the type of the corresponding top
    52365555//               and the level of the constructible
    52375556//     L4= [[v1,p1,[p11,..,p1l],typ_1,level]_1 ,.. [vs,ps,[ps1,..,psl],typ_s,level_s]
    52385557static proc AddLocus(list L)
    52395558{
    5240 //  int te0=0;
    5241 //  def RR=basering;
    5242 //  if(defined(@P)){te0=1;  def Rx=@R;  kill @P; setring RR;}
    52435559  list L1; int i; int j;  list L2; list L3;
    52445560  list l1; list l2;
     
    52775593      {
    52785594        v=L2[k0][2];
     5595        l4[1]=v; l4[2]=p1; l4[3]=L3[i][2][j][2];  l4[5]=level;
     5596        if(size(L2[k0])>2){l4[4]=L2[k0][3];}
     5597        L4[size(L4)+1]=l4;
    52795598      }
    52805599      else{"ERROR p1 NOT FOUND";}
    5281       l4[1]=v; l4[2]=p1; l4[3]=L3[i][2][j][2];  l4[5]=level;
    5282       if(size(L2[k0])>2){l4[4]=L2[k0][3];}
    5283       L4[size(L4)+1]=l4;
    52845600    }
    52855601  }
     
    53055621  for(i=1;i<=size(L);i++)
    53065622  {
    5307     Sc=PtoCrep(list(L[i]));
     5623    Sc=PtoCrep0(list(L[i]));
    53085624    Lc[size(Lc)+1]=Sc;
    53095625  }
    5310   list S=ConsLevels(Lc)[1];
     5626  list S=ConsLevels(Lc);
     5627  S=ConsLevelsToLevels(S);
    53115628  list Sout;
    53125629  list Lev;
    53135630  for(i=1;i<=size(S);i++)
    53145631  {
    5315     Lev=list(i,Prep(S[i][2][1],S[i][2][2]));
     5632    Lev=list(S[i][1],Prep(S[i][2][1],S[i][2][2]));
    53165633    Sout[size(Sout)+1]=Lev;
    53175634  }
     
    53215638//********************* End locus ****************************
    53225639
     5640//********************* Begin WLemma **********************
     5641
     5642// input ideal F in @R
     5643//          ideal a in @R but only depending on parameters
     5644//          F is a generating ideal in V(a);
     5645// output:  ideal b in @R but depending only on parameters
     5646//              ideal G=GBasis(F) in V(a) \ V(b)
     5647proc WLemma(ideal F,ideal a, list #)
     5648"USAGE: WLemma(F,A,#);
     5649          The first argument ideal F in K[a][x];
     5650          The second argument ideal A in K[a]
     5651              ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
     5652              ideal  F=f1(x_1,..,x_n,u_1,..,u_m),..,fs(x_1,..,x_n,u_1,..,u_m);
     5653              ideal A=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
     5654              list # : Options
     5655           Calling sequence:
     5656              WLemma(F,A,#);
     5657OPTIONS: either (\"rep\", 0) or (\"rep\",1) the representation of
     5658              the resulting segment, by default is
     5659              0 =P-representation, (default) but can be set to
     5660              1=C-representation.
     5661RETURN: The list of (lpp,B,S) = (leading power product, basis, segment)
     5662              being B the reduced Groebner Basis given by I-regular functions
     5663              of the specialized ideal F on the
     5664              segment S given in P- or C-representation
     5665NOTE: The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
     5666KEYWORDS: Wibmer's Lemma
     5667EXAMPLE:  WLemma; shows an example"
     5668{
     5669  list L=#;
     5670  int i; int j;
     5671  def RR=basering;
     5672  setglobalrings();
     5673  setring(@RP);
     5674  ideal FF=imap(RR,F);
     5675  FF=std(FF);
     5676  ideal AA=imap(RR,a);
     5677  AA=std(AA);
     5678  ideal FFa;
     5679  poly r;
     5680  for(i=1; i<=size(FF);i++)
     5681  {
     5682    r=reduce(FF[i],AA);
     5683    if(r!=0){FFa[size(FFa)+1]=r;}
     5684  }
     5685  setring RR;
     5686  ideal Fa=imap(@RP,FFa);
     5687  ideal AAA=imap(@RP,AA);
     5688  ideal lppFa;
     5689  ideal lcFa;
     5690  for(i=1;i<=size(Fa);i++)
     5691  {
     5692    lppFa[size(lppFa)+1]=leadmonom(Fa[i]);
     5693    lcFa[size(lcFa)+1]=leadcoef(Fa[i]);
     5694  }
     5695  setring @RP;
     5696  ideal lccr=imap(RR,lppFa);
     5697  lccr=std(lccr);
     5698  setring RR;
     5699  ideal lcc=imap(@RP,lccr);
     5700  list J; list Jx;
     5701  ideal Jci;
     5702  ideal Jxi;
     5703  list B;
     5704  for(i=1;i<=size(lcc);i++)
     5705  {
     5706    kill Jci; ideal Jci; kill Jxi; ideal Jxi;
     5707    for(j=1;j<=size(Fa);j++)
     5708    {
     5709      if(lppFa[j]==lcc[i])
     5710      {
     5711        Jci[size(Jci)+1]=lcFa[j];
     5712        Jxi[size(Jxi)+1]=Fa[j];
     5713      }
     5714    }
     5715    J[size(J)+1]=Jci;
     5716    B[size(B)+1]=Jxi;
     5717  }
     5718  setring @P;
     5719  list Jp=imap(RR,J);
     5720  ideal JL=product(Jp);
     5721  setring(RR);
     5722  def JLA=imap(@P,JL);
     5723  list PR;
     5724  if (size(L)>0)
     5725  {
     5726    if((L[1]=="rep") and (L[2]==1))
     5727    {
     5728      PR=Crep(AAA, JLA);
     5729    }
     5730    else
     5731    {PR=Prep(AAA, JLA);}
     5732  }
     5733  else{PR=Prep(AAA, JLA);}
     5734//  setring(RR);
     5735 // list PRR=imap(@P,PR);
     5736  return(list(lcc,B,PR));
     5737}
     5738example
     5739{
     5740"EXAMPLE:"; echo = 2;
     5741if(defined(R)){kill R;}
     5742ring R=(0,a,b,c,d,e,f),(x,y),lp;
     5743ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2;
     5744ideal A=a*e-b*d;
     5745WLemma(F,A);
     5746WLemma(F,A,"rep",1);
     5747}
     5748
     5749//********************* End WLemma ************************
Note: See TracChangeset for help on using the changeset viewer.