Changeset f9fc81 in git for Singular/LIB/grobcov.lib


Ignore:
Timestamp:
Dec 6, 2013, 4:57:24 PM (10 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
f1a309fae2b1e4366fb7ead1dac3f33dbfc042a5
Parents:
bdda288ddec3d88538ca426b802f82d49d384caa9e8a6c3b0f4b1814c8ea54b2bf5679f4c3e9591d
Message:
Merge pull request #444 from surface-smoothers/fix.normal.callbug.in.genus

Fix.normal.callbug.in.genus
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/grobcov.lib

    r9e8a6c3 rf9fc81  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="version grobcov.lib 4.0.0.0 Jun_2013 "; // $Id$
     2version="version grobcov.lib 4-0-0-0 Dec_2013 ";
    33category="General purpose";
    44info="
    … …  
    1515          The locus algorithm and definitions will be
    1616          published in a forthcoming paper by Abanades, Botana, Montes, Recio:
    17           \''Geometrical locus of points using the Groebner cover\''.
    18 
     17          \''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\''.
    1918
    2019          The central routine is grobcov. Given a parametric
    … …  
    6968PROCEDURES:
    7069
    71 grobcov(F);        Is the basic routine giving the canonical
    72                    Groebner cover of the parametric ideal F.
    73                    This routine accepts many options, that
    74                    allow to obtain results even when the canonical
    75                    computation does not finish in reasonable time.
    76 
    77 cgsdr(F);          Is the procedure for obtaining a first disjoint,
    78                    reduced Comprehensive Groebner System that
    79                    is used in grobcov, but that can be used
    80                    independently if only the CGS is required.
    81                    It is a more efficient routine than buildtree
    82                    (the own routine that is no more used). It uses
    83                    now KSW algorithm.
     70grobcov(F); Is the basic routine giving the canonical
     71            Groebner cover of the parametric ideal F.
     72            This routine accepts many options, that
     73            allow to obtain results even when the canonical
     74            computation does not finish in reasonable time.
     75
     76cgsdr(F);   Is the procedure for obtaining a first disjoint,
     77            reduced Comprehensive Groebner System that
     78            is used in grobcov, that can also be used
     79            independently if only the CGS is required.
     80            It is a more efficient routine than buildtree
     81            (the own routine that is no more used). It uses
     82            now KSW algorithm.
    8483
    8584setglobalrings();  Generates the global rings @R, @P and @PR that are
    86                    respectively the rings Q[a][x], Q[a], Q[x,a].
    87                    It is called inside each of the fundamental routines
    88                    of the library: grobcov, cgsdr and killed before
    89                    the output.
    90                    If the user want to use some other internal routine,
    91                    then setglobalrings() is to be called before, as
    92                    the rings @R, @P and @RP are needed in most of them.
    93                    globally, and more internal routines can be used, but
    94                    these rings are killed by the call to any of the
    95                    basic routines.
    96 
    97 pdivi(f,F);     Performs a pseudodivision of a parametric polynomial
    98                    by a parametric ideal.
     85            respectively the rings Q[a][x], Q[a], Q[x,a].
     86            It is called inside each of the fundamental routines
     87            of the library: grobcov, cgsdr, locus, locusdg and killed before
     88            the output.
     89            So, if the user want to use some other internal routine,
     90            then setglobalrings() is to be called before, as most of them use
     91            these rings.
     92
     93pdivi(f,F); Performs a pseudodivision of a parametric polynomial
     94            by a parametric ideal.
     95            Can be used without calling presiouly setglobalrings(),
    9996
    10097pnormalf(f,E,N);   Reduces a parametric polynomial f over V(E) \ V(N)
    101                    E is the null ideal and N the non-null ideal
    102                    over the parameters.
     98            E is the null ideal and N the non-null ideal
     99            over the parameters.
     100            Can be used without calling presiouly setglobalrings(),
     101
     102Prep(N,M);  Computes the P-representation of V(N) \ V(M).
     103            Can be used without calling previously setglobalrings().
    103104
    104105extend(GC); When the grobcov of an ideal has been computed
    105                    with the default option ('ext',0) and the explicit
    106                    option ('rep',2) (which is not the default), then
    107                    one can call extend (GC) (and options) to obtain the
    108                    full representation of the bases. With the default
    109                    option ('ext',0) only the generic representation of
    110                    the bases are computed, and one can obtain the full
    111                    representation using extend.
    112 
    113 locus(G):   Special routine for determining the locus of points
    114                    of  objects. Given a parametric ideal J with
    115                    parameters (a_1,..a_m) and variables (x_1,..,xn),
    116                    representing the system determining
    117                    the locus of points (a_1,..,a_m)) who verify certain
    118                    properties, computing the grobcov G of
    119                    J and applying to it locus, determines the different
    120                    classes of locus components. They can be
    121                    Normal, Special, Accumulation point, Degenerate.
    122                    The output are the components given in P-canonical form
    123                    of two constructible sets: Normal, and Non-Normal
    124                    The Normal Set has Normal and Special components
    125                    The Non-Normal set has Accumulation and Degenerate components.
    126                    The description of the algorithm and definitions will be
    127                    given in a forthcoming paper by Abanades, Botana, Montes, Recio:
    128                    ''Geometrical locus of points using the Groebner cover''
    129 
    130 locusdg(G):  Is a special routine for computing the locus in dinamical geometry.
    131                     It detects automatically a possible point that is to be avoided by the mover,
    132                     whose coordinates must be the last two coordinates in the definition of the ring.
    133                     If such a point is detected, then it eliminates the segments of the grobcov
    134                     depending on the point that is to be avoided.
    135                     Then it calls locus.
    136 
    137 locusto(L):  Transforms the output of locus to a string that
    138                   can be reed from different computational systems.
     106            with the default option ('ext',0) and the explicit
     107            option ('rep',2) (which is not the default), then
     108            one can call extend (GC) (and options) to obtain the
     109            full representation of the bases. With the default
     110            option ('ext',0) only the generic representation of
     111            the bases are computed, and one can obtain the full
     112            representation using extend.
     113            Can be used without calling presiouly setglobalrings(),
     114
     115locus(G);   Special routine for determining the locus of points
     116            of  objects. Given a parametric ideal J with
     117            parameters (a_1,..a_m) and variables (x_1,..,xn),
     118            representing the system determining
     119            the locus of points (a_1,..,a_m)) who verify certain
     120            properties, computing the grobcov G of
     121            J and applying to it locus, determines the different
     122            classes of locus components. They can be
     123            Normal, Special, Accumulation point, Degenerate.
     124            The output are the components given in P-canonical form
     125            of two constructible sets: Normal, and Non-Normal
     126            The Normal Set has Normal and Special components
     127            The Non-Normal set has Accumulation and Degenerate components.
     128            The description of the algorithm and definitions will be
     129            given in a forthcoming paper by Abanades, Botana, Montes, Recio:
     130            ''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry''.
     131            Can be used without calling presiouly setglobalrings(),
     132
     133locusdg(G); Is a special routine for computing the locus in dinamical geometry.
     134            It detects automatically a possible point that is to be avoided by the mover,
     135            whose coordinates must be the last two coordinates in the definition of the ring.
     136            If such a point is detected, then it eliminates the segments of the grobcov
     137            depending on the point that is to be avoided.
     138            Then it calls locus.
     139            Can be used without calling presiouly setglobalrings(),
     140
     141locusto(L); Transforms the output of locus to a string that
     142            can be reed from different computational systems.
     143            Can be used without calling presiouly setglobalrings(),
     144
     145addcons(L); Given a disjoint set of locally closed subsets in P-representation,
     146            it returns the canonical P-representation of the constructible
     147            set formed by the union of them.
     148            Can be used without calling presiouly setglobalrings(),
    139149
    140150SEE ALSO: compregb_lib
    … …  
    154164// Initial data: 6-9-2009
    155165// Final data: 25-10-2011
    156 // Release 3: (this release, public)
     166// Release 3:
    157167// Initial data: 1-7-2012
    158168// Final data: 4-9-2012
     169// Release 4: (this release, public)
     170// Initial data: 4-9-2012
     171// Final data: 21-11-2013
    159172// basering Q[a][x];
    160173
    … …  
    166179RETURN:   After its call the rings @R=Q[a][x], @P=Q[a], @RP=Q[x,a] are
    167180          defined as global variables.
    168 NOTE:     It is called internally by the fundamental routines of the
    169           library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusto,
     181NOTE: It is called internally by the fundamental routines of the
     182          library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg, locusto,
    170183          and killed before the output.
    171184          The user does not need to call it, except when it is interested
    … …  
    210223}
    211224
    212 //*************Auxilliary routines**************
     225//*************Auxiliary routines**************
    213226
    214227// cld : clears denominators of an ideal and normalizes to content 1
    215 //       can be used in @R or @P or @RP
     228//        can be used in @R or @P or @RP
    216229// input:
    217 //    ideal J (J can be also poly), but the output is an ideal;
     230//        ideal J (J can be also poly), but the output is an ideal;
    218231// output:
    219 //    ideal Jc (the new form of ideal J without denominators and
    220 //       normalized to content 1)
     232//        ideal Jc (the new form of ideal J without denominators and
     233//        normalized to content 1)
    221234proc cld(ideal J)
    222235{
    … …  
    369382//"USAGE:   subset(J,K);
    370383//          (J,K)  expected (ideal,ideal)
    371 //                   or     (list, list)
     384//                  or     (list, list)
    372385//RETURN:   1 if all the elements of J are in K, 0 if not.
    373386//EXAMPLE:  subset; shows an example;"
    … …  
    419432}
    420433
    421 // pdivi : pseudodivision of a poly f by a parametric ideal F in Q[a][x].
     434// pdivi : pseudodivision of a parametric polynomial f by a parametric ideal F in Q[a][x].
    422435// input:
    423 //   poly  f (in the parametric ring @R)
    424 //   ideal F (in the parametric ring @R)
     436//   poly  f
     437//   ideal F
    425438// output:
    426439//   list (poly r, ideal q, poly mu)
    … …  
    438451EXAMPLE:  pdivi; shows an example"
    439452{
     453  int i;
     454  int j;
     455//  string("T_f=",f);
     456  poly v=1;
     457  for(i=1;i<=nvars(basering);i++){v=v*var(i);}
    440458  int te=0;
     459  def R=basering;
    441460  if (defined(@P)==1){te=1;}
    442461  else{setglobalrings();}
    443   def R=basering;
    444   int i;
    445   int j;
    446462  poly r=0;
    447463  poly mu=1;
    … …  
    476492        rho=qlc[1]*qlm;
    477493        p=nu*p-rho*F[i];
     494//        string("i=",i,"  coef(p,v)=",coef(p,v));
    478495        r=nu*r;
    479496        for (j=1;j<=size(F);j++){q[j]=nu*q[j];}
    … …  
    487504      r=r+lcp*lpp;
    488505      p=p-lcp*lpp;
     506//      string("pasa al rem p=",p);
    489507    }
    490508  }
    … …  
    845863{
    846864  int i;
     865  int nt;
     866  if (typeof(F)=="ideal"){nt=ncols(F);}
     867  else{nt=size(F);}
     868
    847869  def FF=F;
    848870  FF=F[1];
    849   for (i=2;i<=ncols(F);i++)
     871  for (i=2;i<=nt;i++)
    850872  {
    851873    if (not(memberpos(F[i],FF)[1]))
    … …  
    9791001  poly g; int i;
    9801002  ideal Nb; ideal Na=N;
    981 
    982   // begins incquotient
    9831003  if (size(Na)==1)
    9841004  {
    … …  
    10551075}
    10561076
     1077// Auxiliary routine to define an order for ideals
     1078// Returns 1 if the ideal a is shoud precede ideal b by sorting them in idbefid order
     1079//             2 if the the contrary happen
     1080//             0 if the order is not relevant
    10571081proc idbefid(ideal a, ideal b)
    10581082{
    … …  
    10861110}
    10871111
     1112// sort a list of ideals using idbefid
    10881113proc sortlistideals(list L)
    10891114{
    … …  
    11391164//    the Prep of V(N)\V(M)
    11401165// Assumed to work in the ring @P of the parameters
     1166// Can be used without calling previously setglobalrings().
    11411167proc Prep(ideal N, ideal M)
    11421168{
     1169  int te;
    11431170  if (N[1]==1)
    11441171  {
    … …  
    11461173  }
    11471174  def RR=basering;
     1175  if(defined(@P)){te=1;}
     1176  else{te=0; setglobalrings();}
    11481177  setring(@P);
    11491178  ideal Np=imap(RR,N);
    … …  
    11751204  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
    11761205  setring(RR);
    1177   return(imap(@P,prep));
     1206  def pprep=imap(@P,prep);
     1207  if(te==0){kill @P; kill @R; kill @RP;}
     1208  return(pprep);
    11781209}
    11791210
    … …  
    11871218//    the C-representaion of V(N)\V(M) = V(ida)\V(idb)
    11881219// Assumed to work in the ring @P of the parameters
     1220// If the routine is to be called from the top, a previous call to
     1221// setglobalrings() is needed.
    11891222proc PtoCrep(list L)
    11901223{
    … …  
    12931326EXAMPLE:  cgsdr; shows an example"
    12941327{
     1328  int te;
    12951329  def RR=basering;
    1296   setglobalrings();
     1330  if(defined(@P)){te=1;}
     1331  else{te=0; setglobalrings();}
    12971332  // INITIALIZING OPTIONS
    12981333  int i; int j;
    … …  
    14521487    }
    14531488  }
    1454   if(defined(@P)){kill @P; kill @R; kill @RP;}
     1489  if(te==0){kill @P; kill @R; kill @RP;}
    14551490  return(GS);
    14561491}
    … …  
    14691504
    14701505// input:  internal routine called by cgsdr at the end to group the
    1471 //         lpp segments and improve the output
     1506//            lpp segments and improve the output
    14721507// output: grouped segments by lpp obtained in cgsdr
    14731508proc grsegments(list T)
    … …  
    15051540// input: ideal p, ideal q
    15061541// output: 1 if p contains q,  0 otherwise
     1542// If the routine is to be called from the top, a previous call to
     1543// setglobalrings() is needed.
    15071544proc idcontains(ideal p, ideal q)
    15081545{
    … …  
    15301567// input: L (list of ideals)
    15311568// output: the list of integers corresponding to the minimal ideals in L
     1569// If the routine is to be called from the top, a previous call to
     1570// setglobalrings() is needed.
    15321571proc selectminideals(list L)
    15331572{
    … …  
    15681607
    15691608// LCUnion
    1570 // Given a list of the P-representations of locally closed segments
     1609// Given a list of the P-representations of disjoint locally closed segments
    15711610// for which we know that the union is also locally closed
    15721611// it returns the P-representation of its union
    … …  
    15761615// output: P-representation of the union
    15771616//       ((P_j,(P_j1,...,P_jk_j | j=1..t)))
     1617// If the routine is to be called from the top, a previous call to
     1618// setglobalrings() is needed.
     1619// Auxiliary routine called by grobcov and addcons
     1620// A previous call to setglobarings() is needed if it is to be used at the top.
    15781621proc LCUnion(list LL)
    15791622{
    15801623  def RR=basering;
    15811624  setring(@P);
     1625  int care;
    15821626  def L=imap(RR,LL);
    15831627  int i; int j; int k; list H; list C; list T;
    1584   list L0; list P0; list P; list Q0; list Q;
     1628  list L0; list P0; list P; list Q0; list Q;  intvec Qi;
     1629  if(not(defined(@Q2))){list @Q2;}
     1630  else{kill @Q2; list @Q2;}
     1631  exportto(Top,@Q2);
    15851632  for (i=1;i<=size(L);i++)
    15861633  {
    … …  
    15921639  }
    15931640  Q0=selectminideals(P0);
     1641  //"T_3Q0="; Q0;
     1642  kill P; list P;
    15941643  for (i=1;i<=size(Q0);i++)
    15951644  {
    1596     Q[i]=L0[Q0[i]];
    1597     P[i]=L[Q[i][1]][Q[i][2]];
     1645    Qi=L0[Q0[i]];
     1646    Q[size(Q)+1]=Qi;
     1647      P[size(P)+1]=L[Qi[1]][Qi[2]];
     1648  }
     1649  @Q2=Q;
     1650  if(size(P)==0)
     1651  {
     1652    setring(RR);
     1653    list TT;
     1654    return(TT);
    15981655  }
    15991656  // P is the list of the maximal components of the union
    … …  
    16111668        {
    16121669          C[size(C)+1]=L[i][j];
     1670          C[size(C)][3]=intvec(i,j);
    16131671        }
    16141672      }
    16151673    }
    1616     T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C));
    1617   }
     1674    if(size(P[k])>=3){T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C),P[k][3]);}
     1675    else{T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C));}
     1676  }
     1677  @Q2=sortpairs(@Q2);
    16181678  setring(RR);
    16191679  def TT=imap(@P,T);
    … …  
    16211681}
    16221682
    1623 // Called by LCUnion to modify the holes of a primepart of the union
     1683// Auxiliary routine
     1684// called by LCUnion to modify the holes of a primepart of the union
    16241685// by the addition of the segments that do not correspond to that part
    16251686// Works on @P ring.
    … …  
    16271688//   H=(p_i1,..,p_is) the holes of a component to be transformed by the addition of
    16281689//        the segments C that do not correspond to that component
    1629 //   C=((q_1,(q_11,..,q_1l_1)),..,(q_k,(q_k1,..,q_kl_k)))
     1690//   C=((q_1,(q_11,..,q_1l_1),pos1),..,(q_k,(q_k1,..,q_kl_k),posk))
     1691// posi=(i,j) position of the component
    16301692//        the list of segments to be added to the holes
    16311693proc addpart(list H, list C)
    … …  
    16331695  list Q; int i; int j; int k; int l; int t; int t1;
    16341696  Q=H; intvec notQ; list QQ; list addq;
     1697  // @Q2=list of (i,j) positions of the components that have been aded to some hole of the maximal ideals
     1698  //          plus those of the components added to the holes.
    16351699  ideal q;
    16361700  i=1;
    … …  
    16451709        if (equalideals(q,C[j][1]))
    16461710        {
     1711          @Q2[size(@Q2)+1]=C[j][3];
    16471712          t=0;
    16481713          for (k=1;k<=size(C[j][2]);k++)
    … …  
    16661731            {
    16671732              addq[size(addq)+1]=C[j][2][k];
     1733              @Q2[size(@Q2)+1]=C[j][3];
    16681734            }
    16691735          }
    … …  
    16971763}
    16981764
    1699 // Called by addpart to finish the modification of the holes of a primepart
     1765// Auxiliary routine called by addpart to finish the modification of the holes of a primepart
    17001766// of the union by the addition of the segments that do not correspond to
    17011767// that part.
    … …  
    17871853}
    17881854
    1789 // specswellCrep
    1790 // input:
    1791 //   given two corresponding polynomials g1 and g2 with the same lpp
    1792 //   g1 belonging to the basis in the segment ida1,idb1
    1793 //   g2 belonging to the basis in the segment ida2,idb2
    1794 // output:
    1795 //   1 if g1 spezializes well to g2 on the whole (ida2,idb2) segment
    1796 //   0 if not
    1797 proc specswellCrep(poly g1, poly g2, ideal ida2)
    1798 {
    1799   poly S;
    1800   S=leadcoef(g2)*g1-leadcoef(g1)*g2;
    1801   def RR=basering;
    1802   setring(@RPt);
    1803   def SR=imap(RR,S);
    1804   def ida2R=imap(RR,ida2);
    1805   attrib(ida2R,"isSB",1);
    1806   poly S2R=reduce(SR,ida2R);
    1807   setring(RR);
    1808   def S2=imap(@RPt,S2R);
    1809   if (S2==0){return(1);}   // and (nonnullCrep(leadcoef(g1),ida2,idb2))
    1810   else {return(0);}
    1811 }
    1812 
     1855// Auxiliary routine for grobcov: ideal F is assumed to be homogeneous
    18131856// gcover
    18141857// input: ideal F: a generating set of a homogeneous ideal in Q[a][x]
    … …  
    18811924    LCU=LCUnion(pairspP);
    18821925    kill prep; list prep;
     1926    kill crep; list crep;
    18831927    for(k=1;k<=size(LCU);k++)
    18841928    {
    … …  
    19451989    for(j=1;j<=size(B);j++)
    19461990    {
    1947       B[j]=pnormalf(B[j],crep[1],N);
     1991      B[j]=pnormalf(B[j],crep[1],crep[2]);
    19481992    }
    19491993    S[i]=list(lpp,B,prep,crep,lpph);
    … …  
    20152059                each point of the segment. With option (\"ext\",1) the
    20162060                full representation of the bases is computed (possible
    2017                 shaves) and sometimes a simpler result is obtained.
     2061                sheaves) and sometimes a simpler result is obtained.
    20182062            \"rep\",0-1-2: The default is (\"rep\",0) and then the segments
    20192063                are given in canonical P-representation. Option (\"rep\",1)
    … …  
    22042248}
    22052249
    2206 // input. GC the grobcov of an ideal in generic representation of the
     2250// Input. GC the grobcov of an ideal in generic representation of the
    22072251//        bases computed with option option ("rep",2).
    2208 // output The grobcov in full representation.
    2209 // option ("comment",1) shows the time.
     2252// Output The grobcov in full representation.
     2253// Option ("comment",1) shows the time.
     2254// Can be called from the top
    22102255proc extend(list GC, list #);
    22112256"USAGE:   extend(GC); When the grobcov of an ideal has been computed
    … …  
    23942439}
    23952440
    2396 
     2441// Auxiliary routine
    23972442// nonzerodivisor
    23982443// input:
    … …  
    24352480}
    24362481
     2482// Auxiliary routine
    24372483// deltai
    24382484// input:
    … …  
    24752521}
    24762522
     2523// Auxiliary routine
    24772524// combine
    24782525// input: a list of pairs ((p1,P1),..,(pr,Pr)) where
    … …  
    24982545}
    24992546
     2547// Auxiliary routine
    25002548// elimconstfac: eliminate the factors in the polynom f that are in K[a]
    25012549// input:
    … …  
    25232571};
    25242572
     2573//Auxiliary routine
    25252574// nullin
    25262575// input:
    … …  
    25442593}
    25452594
     2595// Auxiliary routine
    25462596// monoms
     2597// Input: A polynomial f
     2598// Output: The list of leading terms
    25472599proc monoms(poly f)
    25482600{
    … …  
    25632615}
    25642616
     2617// Auxiliary routine called by extend
    25652618// extend0
    25662619// input:
    … …  
    26222675}
    26232676
     2677// Auxiliary routine
    26242678// findindexpolys
    26252679// input:
    … …  
    26702724}
    26712725
     2726// Auxiliary routine
    26722727// extendcoef: given Q,P in K[a] where P/Q specializes on an open and dense subset
    26732728//      of the whole V(p1 int...int pr), it returns a basis of the module
    … …  
    26982753}
    26992754
     2755// Auxiliary routine
    27002756// selectregularfun
    27012757// input:
    … …  
    27552811}
    27562812
     2813// Auxiliary routine
    27572814// searchinlist
    27582815// input:
    … …  
    27762833}
    27772834
     2835// Auxiliary routine
    27782836// comb: the list of combinations of elements (1,..n) of order p
    27792837proc comb(int n, int p)
    … …  
    28122870}
    28132871
     2872// Auxiliary routine
    28142873// selectminsheaves
    28152874// Input: L=((v_11,..,v_1k_1),..,(v_s1,..,v_sk_s))
    … …  
    28342893}
    28352894
     2895// Auxiliary routine
    28362896// smsheaves
    28372897// Input:
    … …  
    28682928}
    28692929
     2930// Auxiliary routine
    28702931// allsheaves
    28712932// Input: L=((v_11,..,v_1k_1),..,(v_s1,..,v_sk_s))
    … …  
    29102971}
    29112972
     2973// Auxiliary routine
    29122974// numones
    29132975// Input:
    … …  
    29262988}
    29272989
     2990// Auxiliary routine
    29282991// pos
    29292992// Input:  intvec p of zeros and ones
    … …  
    29403003}
    29413004
     3005// Auxiliary routine
    29423006// actualize: actualizes zeroes of p
    29433007// Input:
    … …  
    29573021}
    29583022
     3023// Auxiliary routine
    29593024// combrep
    29603025// Input: V=(n_1,..,n_i)
    … …  
    29903055}
    29913056
     3057// Auxiliary routine
    29923058proc reducemodN(poly f,ideal E)
    29933059{
    … …  
    30053071}
    30063072
     3073// Auxiliary routine
    30073074// intersp: computes the intersection of the ideals in S in @P
    30083075proc intersp(list S)
    … …  
    30173084}
    30183085
     3086// Auxiliary routine
    30193087// radicalmember
    30203088proc radicalmember(poly f,ideal ida)
    … …  
    30393107}
    30403108
     3109// Auxiliary routine
    30413110// NonNull: returns 1 if the poly f is nonnull on V(E)\V(N), 0 otherwise.
    30423111proc NonNull(poly f, ideal E, ideal N)
    … …  
    30593128}
    30603129
     3130// Auxiliary routine
    30613131// selectextendcoef
    30623132// input:
    … …  
    31173187}
    31183188
     3189// Auxiliary routine
    31193190// input:
    31203191//   ideal E1: in some basering (depends only on the parameters)
    … …  
    31333204}
    31343205
     3206// Auxiliary routine
    31353207// reform
    31363208// input:
    … …  
    32163288}
    32173289
     3290// Auxiliary routine
    32183291// nonnullCrep
    32193292proc nonnullCrep(poly f0,ideal ida0,ideal idb0)
    … …  
    32383311}
    32393312
     3313// Auxiliary routine
    32403314// precombint
    32413315// input:  L: list of ideals (works in @P)
    … …  
    32843358}
    32853359
    3286 // precombinediscussion
    3287 // not used, can be deleted
    3288 // input:  list L: the LCU segment with bases for each pi component
    3289 // output: intvec vv:  vv[1]=(1 if the generic polynomial of the vv[2]
    3290 //                     component already specializes well,
    3291 //                     0 if combine is to be used)
    3292 //                     vv[2]=selind, the index for which the generic basis
    3293 //                     already specializes well if combine is not to be used (vv[1]=1).
    3294 proc precombinediscussion(L,crep)
    3295 {
    3296   int tes=1; int selind; int i1; int j1; poly p; poly lcp; intvec vv;
    3297   if (size(L)==1){vv=1,1; return(vv);}
    3298   for (i1=1;i1<=size(L);i1++)
    3299   {
    3300     tes=1;
    3301     p=L[i1][2];
    3302     lcp=leadcoef(p);
    3303 
    3304 
    3305     if(nonnullCrep(lcp,crep[1],crep[2]))
    3306     {
    3307       for(j1=1;j1<=size(L);j1++)
    3308       {
    3309         if(i1!=j1)
    3310         {
    3311           if(specswellCrep(p,L[j1][2],L[j1][1])==0){tes=0; break;}
    3312         }
    3313       }
    3314     }
    3315     else{tes=0;}
    3316     if(tes){selind=i1; break;}
    3317   }
    3318   vv=tes,selind;
    3319   return(vv);
    3320 }
    3321 
     3360
     3361// Auxiliary routine
    33223362// minAssGTZ eliminating denominators
    33233363proc minGTZ(ideal N);
    … …  
    33373377//********************* Begin KapurSunWang *************************
    33383378
     3379// Auxiliary routine
    33393380// inconsistent
    33403381// Input:
    … …  
    33723413}
    33733414
     3415// Auxiliary routine
    33743416// MDBasis: Minimal Dickson Basis
    33753417proc MDBasis(ideal G)
    … …  
    33973439}
    33983440
     3441// Auxiliary routine
    33993442// primepartZ
    34003443proc primepartZ(poly f);
    … …  
    34973540}
    34983541
     3542// Auxiliary routine
    34993543// sqf
    35003544// This is for releases of Singular before 3-5-1
    … …  
    35163560// }
    35173561
     3562// Auxiliary routine
    35183563// sqf
    35193564proc sqf(poly f)
    … …  
    35293574
    35303575
     3576// Auxiliary routine
    35313577// KSW0: Kapur-Sun-Wang algorithm for computing a CGS, called by KSW
    35323578// Input:
    … …  
    35753621  {
    35763622    if(comment>1){"Basis 1 is found"; E; N;}
    3577     return(E0,N0,ideal(1));
    3578   }
     3623    list KK; KK[1]=list(E0,N0,ideal(1));
     3624    return(KK);
     3625   }
    35793626  ideal Gr; ideal Gm; ideal GM;
    35803627  for (i=1;i<=size(G);i++)
    … …  
    36543701}
    36553702
     3703// Auxiliary routine
    36563704// KSWtocgsdr
    36573705proc KSWtocgsdr(list L)
    … …  
    36723720}
    36733721
     3722// Auxiliary routine
    36743723// KtoPrep
    36753724// Computes the P-representaion of a K-representation (N,W) of a set
    … …  
    37203769}
    37213770
     3771// Auxiliary routine
    37223772// groupKSWsegments
    37233773// input:  the list of vertices of KSW
    … …  
    37543804
    37553805//******************** Begin locus ******************************
    3756 
    37573806
    37583807// indepparameters
    … …  
    37773826}
    37783827
    3779 
    3780 // dimP0: Auxiliar routine
     3828// dimP0: Auxiliary routine
    37813829// if the dimension in @P of an ideal in the parameters has dimension 0 then it returns 0
    37823830// else it retuns 1
    … …  
    37943842}
    37953843
    3796 // Auxiliar routine.
     3844// Auxiliary routine.
    37973845// input:    ideals E and F (assumed in ring @P
    37983846// returns: 1 if ideal E is contained in ideal F (assumed F is std basis)
    … …  
    38173865}
    38183866
    3819 // AddCons: given a set of locally closed components of a selection of
    3820 //       segments of the Grobner cover, it builds the canonical P-representation
    3821 //       of the corresponding constructible set, including levels it they are
    3822 // input: a list L of a selection of segments of the Groebner cover
    3823 //       given a a set of components of the form
    3824 //       ( (p_1,(p_11,..,p_1k_1).. (p_s,(p_s1,..,p_sk_s))
    3825 // output: The canonical P-representation of adding the given components.
    3826 proc AddCons(list L)
    3827 {
    3828   // First step: Selecting the top varieties
    3829   list L1; list L2; list LL; int i; int j; int t;
    3830   list Lend;
     3867// addcons: given a set of locally closed sets given in P-representation,
     3868//       (particularly a subset of components of a selection of segments
     3869//       of the Grobner cover), it builds the canonical P-representation
     3870//       of the corresponding constructible set, of its union, including levels it they are.
     3871// input: a list L of disjoint segments (for exmple a selection of segments of the Groebner cover)
     3872//       of the form
     3873//       ( ( (p1_1,(p1_11,..,p1_1k_1).. (p1_s,(p1_s1,..,p1_sk_s)),..,
     3874//         ( (pn_1,(pn_11,..,pn_1j_1).. (pn_s,(pn_s1,..,pn_sj_s))   )
     3875// output: The canonical P-representation of the  constructible set resulting by
     3876//       adding the given components.
     3877//       The first element of each component can be ignored. It is only for internal use.
     3878//       The fifth element of each component is the level of the constructible set
     3879//       of the component.
     3880//       It is no need a previous call to setglobalrings()
     3881proc addcons(list L)
     3882"USAGE:   addcons(  ( ( (p1_1,(p1_11,..,p1_1k_1).. (p1_s,(p1_s1,..,p1_sk_s)),..,
     3883  ( (pn_1,(pn_11,..,pn_1j_1).. (pn_s,(pn_s1,..,pn_sj_s))   )   )
     3884  a list L of locally closed sets in P-representation
     3885RETURN: the canonical P-representation of the constructible set of the union.
     3886NOTE:   It is called internally by the routines locus, locusdg,
     3887
     3888KEYWORDS: P-representation, constructible set
     3889EXAMPLE:  adccons; shows an example"
     3890{
     3891  int te;
     3892  if(defined(@P)){te=1;}
     3893  else{setglobalrings();}
     3894  list notadded; list P0;
     3895  int i; int j; int k; int l;
     3896  intvec v; intvec v1; intvec v0;
     3897  ideal J; list K;
    38313898  for(i=1;i<=size(L);i++)
    38323899  {
    3833     t=1;
    3834     for(j=1;j<=size(L);j++)
    3835     {
    3836       if(i!=j)
    3837       {
    3838         if(containedideal(L[j][1],L[i][1])==1)
    3839         {t=0;
    3840          // "the ideal "; L[j][1]; "is contained in ideal "; L[i][1];
    3841           j=size(L);
    3842         }
    3843       }
    3844     }
    3845     if(t==1){L1[size(L1)+1]=L[i];}
    3846     else{L2[size(L2)+1]=L[i];}
    3847   }
    3848   // Second step: Adding segments to obtain a locally closed sets for each level
    3849   int lev=1;
    3850   if(size(L2)==0)
    3851   {
    3852      for(i=1;i<=size(L1);i++)
     3900    for(j=1;j<=size(L[i]);j++)
     3901    {
     3902        v=i,j;
     3903        notadded[size(notadded)+1]=v;
     3904    }
     3905  }
     3906  //"T_notadded="; notadded;
     3907  int level=1;
     3908  list P=L;
     3909  list LL;
     3910  list Ln;
     3911  //int cont=0;
     3912  while(size(notadded)>0)
     3913  {
     3914    kill notadded; list notadded;
     3915    for(i=1;i<=size(P);i++)
     3916    {
     3917      for(j=1;j<=size(P[i]);j++)
     3918      {
     3919          v=i,j;
     3920          notadded[size(notadded)+1]=v;
     3921      }
     3922    }
     3923    //"T_notadded="; notadded;
     3924     //cont++;
     3925     P0=P;
     3926     Ln=LCUnion(P);
     3927     //"T_Ln="; Ln;
     3928     notadded=minuselements(notadded,@Q2);
     3929     //"T_@Q2="; @Q2;
     3930     //"T_notadded="; notadded;
     3931     for(i=1;i<=size(Ln);i++)
    38533932     {
    3854        if(size(L1[i])>=3)
    3855        {L1[i][3]=string(string(L1[i][3]),",",lev);}
    3856        else{L1[i][3]=string(lev);}
     3933       //Ln[i][4]=  ; JA HAURIA DE VENIR DE LCUnion
     3934       Ln[i][5]=level;
     3935       LL[size(LL)+1]=Ln[i];
    38573936     }
    3858    return(L1);
    3859   }
    3860    while(size(L2)>0)
    3861    {
    3862      LL=addtolocalclosed(L1,L2);
    3863      for(i=1;i<=size(LL[1]);i++)
     3937     i=0;j=0;
     3938     v=0,0;
     3939     v0=0,0;
     3940     kill P; list P;
     3941     for(l=1;l<=size(notadded);l++)
    38643942     {
    3865        if(size(LL[1][i])>=3)
    3866        {LL[1][i][3]=string(string(LL[1][i][3]),",",lev);}
    3867        else{LL[1][i][3]=string(lev);}
    3868        Lend[size(Lend)+1]=LL[1][i];
     3943       v1=notadded[l];
     3944       if(v1[1]<>v0[1]){i++;j=1; P[i]=K;} // list P[i];
     3945       else
     3946       {
     3947         if(v1[2]<>v0[2])
     3948         {
     3949           j=j+1;
     3950         }
     3951       }
     3952       v=i,j;
     3953       //"T_v1="; v1;
     3954       P[i][j]=K; P[i][j][1]=J; P[i][j][2]=K;
     3955       P[i][j][1]=P0[v1[1]][v1[2]][1];
     3956       //"T_P0[v1[1]][v1[2]][2]="; P0[v1[1]][v1[2]][2];
     3957       //"T_size(P0[v1[1]][v1[2]][2])="; size(P0[v1[1]][v1[2]][2]);
     3958       for(k=1;k<=size(P0[v1[1]][v1[2]][2]);k++)
     3959       {
     3960         P[i][j][2][k]=P0[v1[1]][v1[2]][2][k];
     3961         //string("T_P[",i,"][",j,"][2][",k,"]="); P[i][j][2][k];
     3962       }
     3963       v0=v1;
     3964       //"T_P_1="; P;
    38693965     }
    3870      L2=LL[2];
    3871      lev++;
    3872    }
    3873    return(Lend);
    3874 }
    3875 
    3876 // input: Two lists of P-components to be added.
    3877 //           L1 contains the top components, and L2 the remaining components
    3878 //           each of the form ( (p_1,(p_11,..,p_1k_1).. (p_s,(p_s1,..,p_sk_s))
    3879 // output: A list of two lists:
    3880 //           The first one contains the P-representation of the top components added
    3881 //           The second contains the components that have not yet been added when
    3882 //           the total subset is not locally closed. If so addtolocalclosed is to be called
    3883 //           at new after separating the new top and remaining components in order
    3884 //           to compute the next level of the constructible set.
    3885 //           The input and the output must be in the ring @P of parameters.
    3886 proc  addtolocalclosed(list L1,list L2)
    3887 {
    3888    // Second step: Adding segments to obtain a locally closed set
    3889    // L1 contains the top components (ideals and descendants)
    3890    // L2 contains the nontop components (ideals and descendants)
    3891    list LL; int i; int j; int t; intvec Added; int mesvoltes=1; intvec alreadyadded; list LN;
    3892    int k; int l; int m; ideal top; ideal hole; ideal nhole; intvec nottoadd; int t0; list h;
    3893    LL=L1;
    3894    LN=L2;
    3895    while(mesvoltes)
    3896    {
    3897      //"volta";
    3898      for(i=1;i<=size(L1);i++)
    3899      {
    3900         Added=Added,alreadyadded;
    3901         Added=sort(elimrepeatedvl(Added))[1];
    3902         kill alreadyadded; intvec alreadyadded;
    3903         top=LL[i][1][1];
    3904         j=1;
    3905         while(j<=size(LL[i][2]))
    3906         {
    3907            kill nottoadd; intvec nottoadd;
    3908            hole=LL[i][2][j];
    3909            t0=1;
    3910            k=1;
    3911            while(t0 and k<=size(LN))
    3912            {
    3913               if (equalideals(hole,LN[k][1])==1)
    3914               {
    3915                  t0=0;
    3916                   if(alreadyadded==intvec(0)){alreadyadded[1]=k;}
    3917                  else{alreadyadded[size(alreadyadded)+1]=k;}
    3918                  LL[i][2]=elimfromlist(LL[i][2],j);
    3919                  j=j-1;
    3920                  for(l=1;l<=size(LN[k][2]);l++)
    3921                  {
    3922                    nhole=LN[k][2][l],LL[i][1];
    3923                    nhole=std(nhole);
    3924                    t=1; m=1;
    3925                     while(t and m<=size(LL[i][2]))
    3926                     {
    3927                        if(containedideal(LL[i][2][m],nhole)==1)
    3928                        {
    3929                          t=0;
    3930                        }
    3931                        m++;
    3932                     }
    3933                     if(t==0){nottoadd[size(nottoadd)+1]=l;}
    3934                   }
    3935                  for(m=1;m<=size(LN[k][2]);m++)
    3936                  {
    3937                      if(memberpos(m,nottoadd)[1]==0)
    3938                     {
    3939                        LL[i][2][size(LL[i][2])+1]=LN[k][2][m];
    3940                     }
    3941                  }
    3942               }
    3943               k++;
    3944 
    3945            }
    3946            j++;
    3947         }
    3948         if(size(LL[i][2])==0 and size(LL[i][1])>0){LL[i][2][1]=ideal(1);}
    3949      }
    3950      h=1,1;
    3951      while((h[1]==1) and (alreadyadded!=intvec(0)))
    3952      {
    3953        h=memberpos(0,alreadyadded);
    3954        if(h[1]==1){alreadyadded=elimfromlist(alreadyadded,h[2]);}
    3955      }
    3956      if(alreadyadded!=intvec(0))
    3957      {alreadyadded=sort(elimrepeatedvl(alreadyadded))[1];}
    3958      if(Added==intvec(0)){Added=alreadyadded;}
    3959      else{
    3960        Added=Added,alreadyadded;
    3961        Added=sort(elimrepeatedvl(Added))[1];
    3962      }
    3963      h=1,1;
    3964      while((h[1]==1) and (Added!=intvec(0)))
    3965      {
    3966        h=memberpos(0,Added);
    3967        if(h[1]==1){Added=elimfromlist(Added,h[2]);}
    3968      }
    3969      if (alreadyadded==intvec(0))
    3970      {
    3971        mesvoltes=0;
    3972      }
    3973   }
    3974   if(Added!=intvec(0)){Added=sort(elimrepeatedvl(Added))[1]; }
    3975   if(Added!=intvec(0))
    3976   {
    3977     for(i=1;i<=size(Added);i++)
    3978     {
    3979       if(size(LN)>0){LN=elimfromlist(LN,Added[size(Added)+1-i]);}
    3980     }
    3981   }
    3982   for (i=1;i<=size(LL);i++)
    3983   {
    3984     for(j=1;j<=size(LL[i][2]);j++)
    3985     {
    3986       hole=LL[i][2][j];
    3987       for (k=1;k<=size(LL);k++)
    3988       {
    3989         if(k!=i)
    3990         {
    3991           if(containedideal(LL[k][1],hole))
    3992           {
    3993             LL[i][2]=elimfromlist(LL[i][2],j);
    3994             for(l=1;l<=size(LL[k][2]);l++)
    3995             {
    3996               nhole=hole,LL[k][2][l];
    3997               nhole=std(nhole);
    3998               if(equalideals(nhole,ideal(1))==0)
    3999               {
    4000                 m=1; t=1;
    4001                 while(t and m<size(LL[i][2]))
    4002                 {
    4003                   if(containedideal(LL[i][2][m],nhole)){t=0;}
    4004                   m++;
    4005                 }
    4006                 if(t==1){LL[i][2][size(LL[i][2])+1]=nhole;}
    4007               }
    4008             }
    4009           }
    4010         }
    4011       }
    4012     }
    4013   }
    4014   LL[1]=LL; LL[2]=LN;
     3966     //"T_P="; P;
     3967     level++;
     3968     //if(cont>3){break;}
     3969     //kill notadded; list notadded;
     3970  }
     3971  kill @Q2;
     3972  if(te==0){kill @P; kill @R; kill @RP;}
     3973  //"T_LL_sortida addcons="; LL; "Fi sortida";
    40153974  return(LL);
    40163975}
    4017 
    4018 // locus(G):  Special routine for determining the locus of points
    4019 //                of  objects. Given a parametric ideal J with
    4020 //                parameters (a_1,..a_m) and variables (x_1,..,xn),
    4021 //                representing the system determining
    4022 //                the locus of points (a_1,..,a_m)) who verify certain
    4023 //                properties, computing the grobcov G of
    4024 //                J and applying to it locus, determines the different
    4025 //                classes of locus components. They can be
    4026 //                Normal, Special, Accumulation point, Degenerate.
    4027 //                The output are the components given in P-canonical form
    4028 //                of at most 4 constructible sets: Normal, Special, Accumulation,
    4029 //                Degenerate.
    4030 //                The description of the algorithm and definitions will be
    4031 //                given in a forthcoming paper by Abanades, Botana, Montes Recio.
    4032 
    4033 // input:      The output G of the grobcov (in generic representation, which is the default)
     3976example
     3977{
     3978  "EXAMPLE:"; echo = 2;
     3979   ring R=(0,a,b),(x1,y1,x2,y2,p,r),dp;
     3980   ideal S93=(a+1)^2+b^2-(p+1)^2,
     3981                    (a-1)^2+b^2-(1-r)^2,
     3982                    a*y1-b*x1-y1+b,
     3983                    a*y2-b*x2+y2-b,
     3984                    -2*y1+b*x1-(a+p)*y1+b,
     3985                    2*y2+b*x2-(a+r)*y2-b,
     3986                    (x1+1)^2+y1^2-(x2-1)^2-y2^2;
     3987  short=0;
     3988  "System S93="; S93; " ";
     3989  def GC93=grobcov(S93);
     3990  "grobcov(S93)="; GC93; " ";
     3991  int i;
     3992  list H;
     3993  for(i=1;i<=size(GC93);i++){H[i]=GC93[i][3];}
     3994  "H="; H;
     3995  "addcons(H)="; addcons(H);
     3996}
     3997
     3998// Takes a list of intvec and sorts it and eliminates repeated elements.
     3999// Auxiliary routine used in addcons
     4000proc sortpairs(L)
     4001{
     4002  def L1=sort(L);
     4003  def L2=elimrepeated(L1[1]);
     4004  return(L2);
     4005}
     4006
     4007// Eliminates the pairs of L1 that are also in L2.
     4008// Auxiliary routine used in addcons
     4009proc minuselements(list L1,list L2)
     4010{
     4011  int i;
     4012  list L3;
     4013  for (i=1;i<=size(L1);i++)
     4014  {
     4015    if(not(memberpos(L1[i],L2)[1])){L3[size(L3)+1]=L1[i];}
     4016  }
     4017  return(L3);
     4018}
     4019
     4020// locus0(G): Private routine used by locus (the public routine), that
     4021//                builds the diferent components.
     4022// input:      The output G of the grobcov (in generic representation, which is the default option)
    40344023// output:
    40354024//         list, the canonical P-representation of the Normal and Non-Normal locus:
    40364025//              The Normal locus has two kind of components: Normal and Special.
    40374026//              The Non-normal locus has two kind of components: Accumulation and Degenerate.
    4038 //              Normal components: for each point in the component,
    4039 //              the number of solutions in the variables is finite, and
    4040 //              the solutions depend on the point in the component if the component is not 0-dimensional.
    4041 //              Special components:  for each point in the component,
    4042 //              the number of solutions in the variables is finite,
    4043 //              the component is not 0-dimensional, but the solutions do not depend on the
    4044 //              values of the parameters in the component.
    4045 //              Accumlation points: are 0-dimensional components for which it exist
    4046 //              an infinite number of solutions.
    4047 //              Degenerate components: are components of dimension greater than 0 for which
    4048 //              for every point in the component there exist infinite solutions.
     4027//              This routine is compemented by locus that calls it in order to eliminate problems
     4028//              with degenerate points of the mover.
    40494029//         The output components are given as
    40504030//              ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
    … …  
    40524032//              If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
    40534033//              gives the depth of the component.
    4054 proc locus(list GG, list #)
    4055 "USAGE:   locus(G);
    4056                The input must be the grobcov  of a parametrical ideal
    4057 RETURN:  The  locus.
    4058                The output are the components of two constructible subsets of the locus
    4059                of the parametrical system.: Normal and Non-normal.
    4060                These components are
    4061                given as a list of  (pi,(pi1,..pis_i),type_i,level_i) varying i,
    4062                where the p's are prime ideals, the type can be: Normal, Special,
    4063                Accumulation, Degenerate.
    4064 NOTE:      It can only be called after computing the grobcov of the
    4065                parametrical ideal in generic representation ('ext',0),
    4066                which is the default.
    4067                The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
    4068 KEYWORDS: geometrical locus, locus, loci.
    4069 EXAMPLE:  locus; shows an example"
     4034proc locus0(list GG, list #)
    40704035{
    40714036  int t1=1; int t2=1;
    40724037  def R=basering;
    40734038  int moverdim=nvars(R);
     4039  list HHH;
     4040  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);}
    40744041  list Lop=#;
    40754042  if(size(Lop)>0){moverdim=Lop[1];}
    … …  
    40774044  list G1; list G2;
    40784045  def G=GG;
     4046  list Q1; list Q2;
    40794047  int i; int d; int j; int k;
    40804048  t1=1;
    … …  
    40824050  {
    40834051    attrib(G[i][1],"IsSB",1);
    4084     //d=dim(std(G[i][1]));
    4085     // ULL
    40864052    d=locusdim(G[i][2],moverdim);
    40874053    if(d==0){G1[size(G1)+1]=G[i];}
    … …  
    41314097          P1[i][3][j]=P1[i][3][j]*h[k];
    41324098        }
    4133         //P1[i][3][j]=normalize(P1[i][3][j]);
    41344099      }
    41354100    }
    41364101  }
    41374102  else{list P1;}
     4103  ideal BB;
    41384104  for(i=1;i<=size(P1);i++)
    41394105  {
    … …  
    41494115    }
    41504116  }
     4117  list l;
     4118  for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1;
     4119  for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2;
     4120
    41514121  setring @P;
    4152   // CANVIAR
     4122  ideal J;
    41534123  if(t1==1)
    41544124  {
    41554125    def C1=imap(R,P1);
    4156     def L1=AddCons(C1);
     4126    def L1=addcons(C1);
    41574127   }
    41584128  else{list C1; list L1; kill P1; list P1;}
    … …  
    41604130  {
    41614131    def C2=imap(R,P2);
    4162     def L2=AddCons(C2);
     4132    def L2=addcons(C2);
    41634133  }
    41644134  else{list L2; list C2; kill P2; list P2;}
    41654135    for(i=1;i<=size(L2);i++)
    41664136    {
    4167       d=dim(std(L2[i][1]));
     4137      J=std(L2[i][2]);
     4138      d=dim(J); // AQUI
    41684139      if(d==0)
    41694140      {
    4170         L2[i][3]=string("Accumulation,",L2[i][3]);
    4171       }
    4172       else{L2[i][3]=string("Degenerate,",L2[i][3]);}
     4141        L2[i][4]=string("Accumulation",L2[i][4]);
     4142      }
     4143      else{L2[i][4]=string("Degenerate",L2[i][4]);}
    41734144    }
    41744145  list LN;
    … …  
    41854156  for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
    41864157  kill @R; kill @RP; kill @P;
    4187   return(L);
    4188 }
    4189 example
    4190 {"EXAMPLE:"; echo = 2;
    4191   ring R=(0,a,b),(x,y),dp;
    4192   short=0;
    4193   "Concoid";
    4194   ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
    4195   "System="; S96; " ";
    4196   locus(grobcov(S96));
    4197 }
    4198 
    4199 
    4200 // locusto: Transforms the output of locus to a string that
    4201 //      can be read by different computational systems.
    4202 // input:
    4203 //     list L: The output of locus
     4158  list LL;
     4159  for(i=1;i<=size(L);i++)
     4160  {
     4161    LL[i]=l;
     4162    LL[i]=elimfromlist(L[i],1);
     4163  }
     4164  return(LL);
     4165}
     4166
     4167// locus0dg(G): Private routine used by locusdg (the public routine), that
     4168//                builds the diferent components used in Dynamical Geometry
     4169// input:      The output G of the grobcov (in generic representation, which is the default option)
    42044170// output:
    4205 //     string s: The output of locus converted to a string readable by other programs
    4206 proc locusto(list L)
    4207 "USAGE:   locusto(G);
    4208           The argument must be the output of locus of a parametrical ideal
    4209           It transforms the output into a string in standard form
    4210           readable in many languages (Geogebra).
    4211 
    4212 RETURN: The locus in string standard form
    4213 NOTE:     It can only be called after computing the locus(grobcov(F)) of the
    4214               parametrical ideal.
    4215               The basering R, must be of the form Q[a,b,..][x,y,..].
    4216 KEYWORDS: geometrical locus, locus, loci.
    4217 EXAMPLE:  locusto; shows an example"
    4218 {
    4219   int i; int j; int k;
    4220   string s;
    4221   s="[";
    4222   ideal p;
    4223   ideal q;
     4171//         list, the canonical P-representation of the Relevant components of the locus in Dynamical Geometry,
     4172//              i.e. the Normal and Accumulation components.
     4173//              This routine is compemented by locusdg that calls it in order to eliminate problems
     4174//              with degenerate points of the mover.
     4175//         The output components are given as
     4176//              ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
     4177//              whwre type_i is always "Relevant".
     4178//         The components are given in canonical P-representation of the subset.
     4179//              If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
     4180//              gives the depth of the component.
     4181proc locus0dg(list GG, list #)
     4182{
     4183  int t1=1; int t2=1;
     4184  def R=basering;
     4185  int moverdim=nvars(R);
     4186  list HHH;
     4187  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);}
     4188  list Lop=#;
     4189  if(size(Lop)>0){moverdim=Lop[1];}
     4190  setglobalrings();
     4191  list G1; list G2;
     4192  def G=GG;
     4193  list Q1; list Q2;
     4194  int i; int d; int j; int k;
     4195  t1=1;
     4196  for(i=1;i<=size(G);i++)
     4197  {
     4198    attrib(G[i][1],"IsSB",1);
     4199    d=locusdim(G[i][2],moverdim);
     4200    if(d==0){G1[size(G1)+1]=G[i];}
     4201    else
     4202    {
     4203      if(d>0){G2[size(G2)+1]=G[i];}
     4204    }
     4205  }
     4206  if(size(G1)==0){t1=0;}
     4207  if(size(G2)==0){t2=0;}
     4208  setring(@RP);
     4209  if(t1)
     4210  {
     4211    list G1RP=imap(R,G1);
     4212  }
     4213  else {list G1RP;}
     4214  list P1RP;
     4215  ideal B;
     4216  for(i=1;i<=size(G1RP);i++)
     4217  {
     4218     kill B;
     4219     ideal B;
     4220     for(k=1;k<=size(G1RP[i][3]);k++)
     4221     {
     4222       attrib(G1RP[i][3][k][1],"IsSB",1);
     4223       G1RP[i][3][k][1]=std(G1RP[i][3][k][1]);
     4224       for(j=1;j<=size(G1RP[i][2]);j++)
     4225       {
     4226         B[j]=reduce(G1RP[i][2][j],G1RP[i][3][k][1]);
     4227       }
     4228       P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B);
     4229     }
     4230  }
     4231  setring(R);
     4232  ideal h;
     4233  if(t1)
     4234  {
     4235    def P1=imap(@RP,P1RP);
     4236    for(i=1;i<=size(P1);i++)
     4237    {
     4238      for(j=1;j<=size(P1[i][3]);j++)
     4239      {
     4240        h=factorize(P1[i][3][j],1);
     4241        P1[i][3][j]=h[1];
     4242        for(k=2;k<=size(h);k++)
     4243        {
     4244          P1[i][3][j]=P1[i][3][j]*h[k];
     4245        }
     4246      }
     4247    }
     4248  }
     4249  else{list P1;}
     4250  ideal BB;
     4251  for(i=1;i<=size(P1);i++)
     4252  {
     4253    if (indepparameters(P1[i][3])==1){P1[i][3]="Special";}
     4254    else{P1[i][3]="Relevant";}
     4255  }
     4256  list P2;
     4257  for(i=1;i<=size(G2);i++)
     4258  {
     4259    for(k=1;k<=size(G2[i][3]);k++)
     4260    {
     4261      P2[size(P2)+1]=list(G2[i][3][k][1],G2[i][3][k][2]);
     4262    }
     4263  }
     4264  list l;
     4265  for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1;
     4266  for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2;
     4267
     4268  setring @P;
     4269  ideal J;
     4270  if(t1==1)
     4271  {
     4272    def C1=imap(R,P1);
     4273    def L1=addcons(C1);
     4274   }
     4275  else{list C1; list L1; kill P1; list P1;}
     4276  if(t2==1)
     4277  {
     4278    def C2=imap(R,P2);
     4279    def L2=addcons(C2);
     4280  }
     4281  else{list L2; list C2; kill P2; list P2;}
     4282    for(i=1;i<=size(L2);i++)
     4283    {
     4284      J=std(L2[i][2]);
     4285      d=dim(J); // AQUI
     4286      if(d==0)
     4287      {
     4288        L2[i][4]=string("Relevant",L2[i][4]);
     4289      }
     4290      else{L2[i][4]=string("Degenerate",L2[i][4]);}
     4291    }
     4292  list LN;
     4293  list ll; list l0;
     4294  if(t1==1)
     4295  {
     4296    for(i=1;i<=size(L1);i++)
     4297    {
     4298      if(L1[i][4]=="Relevant")
     4299      {
     4300       LN[size(LN)+1]=l0;
     4301       LN[size(LN)][1]=L1[i];
     4302     }
     4303    }
     4304  }
     4305  if(t2==1)
     4306  {
     4307    for(i=1;i<=size(L2);i++)
     4308    {
     4309      if(L2[i][4]=="Relevant")
     4310      {
     4311        LN[size(LN)+1]=l0;
     4312        LN[size(LN)][1]=L2[i];
     4313      }
     4314    }
     4315  }
     4316  list LNN;
     4317  kill ll; list ll; list lll; list llll;
     4318  for(i=1;i<=size(LN);i++)
     4319  {
     4320      LNN[size(LNN)+1]=l0;
     4321      ll=LN[i][1]; lll[1]=ll[2]; lll[2]=ll[3]; lll[3]=ll[4]; LNN[size(LNN)][1]=lll;
     4322  }
     4323  kill LN; list LN;
     4324  LN=addcons(LNN);
     4325  setring(R);
     4326  def L=imap(@P,LN);
     4327  for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
     4328  kill @R; kill @RP; kill @P;
     4329  list LL;
    42244330  for(i=1;i<=size(L);i++)
    42254331  {
    4226     s=string(s,"[[");
    4227     for (j=1;j<=size(L[i][1]);j++)
    4228     {
    4229       s=string(s,L[i][1][j],",");
    4230     }
    4231     s[size(s)]="]";
    4232     s=string(s,",[");
    4233     for(j=1;j<=size(L[i][2]);j++)
    4234     {
    4235       s=string(s,"[");
    4236       for(k=1;k<=size(L[i][2][j]);k++)
    4237       {
    4238         s=string(s,L[i][2][j][k],",");
    4239       }
    4240       s[size(s)]="]";
    4241       s=string(s,",");
    4242     }
    4243     s[size(s)]="]";
    4244     s=string(s,"]");
    4245     if(size(L[i])>=3)
    4246     {
    4247       s[size(s)]=",";
    4248       s=string(s,string(L[i][3]),"]");
    4249     }
    4250     if(size(L[i])>=4)
    4251     {
    4252       s[size(s)]=",";
    4253       s=string(s,string(L[i][4]),"],");
    4254     }
    4255     s[size(s)]="]";
    4256     s=string(s,",");
    4257   }
    4258   s[size(s)]="]";
    4259   return(s);
    4260 }
    4261 example
    4262 {"EXAMPLE:"; echo = 2;
    4263   ring R=(0,a,b),(x,y),dp;
    4264   short=0;
    4265   ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
    4266   "System="; S96; " ";
    4267   locusto(locus(grobcov(S96)));
    4268 }
    4269 
    4270 // locusdg is the routine for computing the locus in dinamical geometry.
    4271 //    It detects automatically a possible point that is to be avoided by the mover,
    4272 //    whose coordinates must be the last coordinates in the definition of the ring (2 last by default,
    4273 //    or the dd last if the option locusdg(G,dd) is specified).
    4274 //    If such a point is detected, then it eliminates the segments of the grobcov that
    4275 //    have this point as solution.
    4276 //    Then it calls locus.
    4277 proc locusdg(list GG,list #)
    4278 "USAGE:  locusdg(GG) ;
    4279           The argument must be the output of grobcov
    4280 RETURN:  The  components of the locus of points in dinamical geometry
    4281 NOTE:  The basering R, must be of the form Q[a][x], a=parameters,
    4282           x=variables, and the mover coordinates must be the two last variables in the
    4283           definition of the ring.
    4284 KEYWORDS: ring, locus, grobcov
    4285 EXAMPLE:  locusdg(GG); shows an example"
     4332    LL[i]=l;
     4333    LL[i]=elimfromlist(L[i],1);
     4334  }
     4335  return(LL);
     4336}
     4337
     4338
     4339//  locus(G):  Special routine for determining the locus of points
     4340//                 of  geometrical constructions. Given a parametric ideal J with
     4341//                 parameters (a_1,..a_m) and variables (x_1,..,xn),
     4342//                 representing the system determining
     4343//                 the locus of points (a_1,..,a_m) who verify certain
     4344//                 properties, computing the grobcov G of
     4345//                 J and applying to it locus, determines the different
     4346//                 classes of locus components. The components can be
     4347//                 Normal, Special, Accumulation, Degenerate.
     4348//                 The output are the components given in P-canonical form
     4349//                 of at most 4 constructible sets: Normal, Special, Accumulation,
     4350//                 Degenerate.
     4351//                 The description of the algorithm and definitions is
     4352//                 given in a forthcoming paper by Abanades, Botana, Montes, Recio.
     4353//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
     4354//                 The mover coordinates are to be placed as the last variables, and by default,
     4355//                 its number is 2. If one consider locus problems in higer dimensions, the number of
     4356//                 mover coordinates (placed as the last variables) is to be given as an option.
     4357//
     4358//  input:      The output G of the grobcov (in generic representation, which is the default option for grobcov)
     4359//  output:
     4360//          list, the canonical P-representation of the Normal and Non-Normal locus:
     4361//               The Normal locus has two kind of components: Normal and Special.
     4362//               The Non-normal locus has two kind of components: Accumulation and Degenerate.
     4363//               Normal components: for each point in the component,
     4364//               the number of solutions in the variables is finite, and
     4365//               the solutions depend on the point in the component if the component is not 0-dimensional.
     4366//               Special components:  for each point in the component,
     4367//               the number of solutions in the variables is finite,
     4368//               the component is not 0-dimensional, but the solutions do not depend on the
     4369//               values of the parameters in the component.
     4370//               Accumlation points: are 0-dimensional components for which it exist
     4371//               an infinite number of solutions.
     4372//               Degenerate components: are components of dimension greater than 0 for which
     4373//               for every point in the component there exist infinite solutions.
     4374//          The output components are given as
     4375//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
     4376//          The components are given in canonical P-representation of the subset.
     4377//               If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
     4378//               gives the depth of the component of the constructible set.
     4379proc locus(list GG, list #)
     4380  "USAGE:   locus(G);
     4381                 The input must be the grobcov  of a parametrical ideal
     4382  RETURN:  The  locus.
     4383                 The output are the components of four constructible subsets of the locus
     4384                 of the parametrical system: Normal , Special, Accumulation and Degenerate.
     4385                 These components are
     4386                 given as a list of  (pi,(pi1,..pis_i),type_i,level_i) varying i,
     4387                 where the p's are prime ideals, the type can be: Normal, Special,
     4388                 Accumulation, Degenerate.
     4389  NOTE:      It can only be called after computing the grobcov of the
     4390                 parametrical ideal in generic representation ('ext',0),
     4391                 which is the default.
     4392                 The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
     4393  KEYWORDS: geometrical locus, locus, loci.
     4394  EXAMPLE:  locus; shows an example"
    42864395{
    42874396  def R=basering;
    … …  
    42944403  if (equalideals(B0,ideal(1)))
    42954404  {
    4296     return(locus(GG));
     4405    return(locus0(GG));
    42974406  }
    42984407  else
    … …  
    43254434    if(te)
    43264435    {
    4327       string("locusdg detected that the mover must avoid point (",N,") in order to obtain the locus");
    4328       // eliminate segments of GG where N is contained in the basis
     4436      string("locus detected that the mover must avoid point (",N,") in order to obtain the correct locus");
     4437      // eliminates segments of GG where N is contained in the basis
    43294438      list nGP;
    43304439      def GP=GG;
    … …  
    43434452   }
    43444453  }
    4345   //"T_nGP="; nGP;
    43464454  kill @RP; kill @P; kill @R;
    4347   return(locus(nGP,dd));
    4348 }
    4349 //      // Now eliminate components for which the basis has dim>0, and all elements non reducing to 0
    4350 //      // modulo N do not contain the mover coordinates
    4351 //      list nnGP; poly r;
    4352 //      for(j=1; j<=size(nGP);j++)
    4353 //      {
    4354 //        if(not(indepparameters(nGP[j][2])))
    4355 //        {
    4356 //          if(dim(std(nGP[j][1]))==0){nnGP[size(nnGP)+1]=nGP[j];}
    4357 //          else
    4358 //          {
    4359 //            te=1; k=1;
    4360 //            while(te and k<=size(nGP[2]))
    4361 //            {
    4362 //              r=pdivi(nGP[j][2][k],N)[1];
    4363 //              if(r==0){k++;}
    4364 //              else
    4365 //              {
    4366 //                if(not(subset(variables(nGP[j][2][k]),vmov )))
    4367 //                {
    4368 //                  te=0;
    4369 //                }
    4370 //                else{k++;}
    4371 //              }
    4372 //            }
    4373 //            if(te==1){nnGP[size(nnGP)+1]=nGP[j];}
    4374 //          }
    4375 //        }
    4376 //      }
    4377 //
    4378 //    }
    4379 //    kill @RP; kill @P; kill @R;
    4380 //    return(locus(nnGP));
    4381 //  }
    4382 //}
     4455  return(locus0(nGP,dd));
     4456}
    43834457example
    43844458{"EXAMPLE:"; echo = 2;
    … …  
    43904464             (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
    43914465  short=0;
    4392   locusdg(grobcov(S));
    4393 }
    4394 
     4466  locus(grobcov(S)); " ";
     4467  kill R;
     4468  ring R=(0,a,b),(x,y),dp;
     4469  short=0;
     4470  "Concoid";
     4471  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
     4472  "System="; S96; " ";
     4473  locus(grobcov(S96));
     4474  kill R; ring R=(0,x,y),(x1,x2),dp;
     4475  ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
     4476              (x1 - 5)^2 + (x2 - 2)^2 - 4,
     4477              -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
     4478 "System="; S;
     4479  locus(grobcov(S)); " ";
     4480  ideal S1=-(x - x1)*(x1 - 4) + (x2 - y)*(x2 - 4),
     4481                (x1 - 4)^2 + (x2 - 2)^2 - 4,
     4482                -2*(x1 - 4)*(2*x2 - y - 4) - 2*(x - 2*x1 + 4)*(x2 - 2);
     4483 "System="; S1;
     4484  locus(grobcov(S1));
     4485}
     4486
     4487//  locusdg(G):  Special routine for determining the locus of points
     4488//                 of  geometrical constructions. Given a parametric ideal J with
     4489//                 parameters (a_1,..a_m) and variables (x_1,..,xn),
     4490//                 representing the system determining
     4491//                 the locus of points (a_1,..,a_m) who verify certain
     4492//                 properties, computing the grobcov G of
     4493//                 J and applying to it locus, determines the different
     4494//                 classes of locus components. The components can be
     4495//                 Normal, Special, Accumulation point, Degenerate.
     4496//                 The output are the components given in P-canonical form
     4497//                 of at most 4 constructible sets: Normal, Special, Accumulation,
     4498//                 Degenerate.
     4499//                 The description of the algorithm and definitions is
     4500//                 given in a forthcoming paper by Abanades, Botana, Montes, Recio.
     4501//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
     4502//                 The mover coordinates are to be placed as the last variables, and by default,
     4503//                 its number is 2. If onw consider locus problems in higer dimensions, the number of
     4504//                 mover coordinates (placed as the last variables) is to be given as an option.
     4505//
     4506//  input:      The output G of the grobcov (in generic representation, which is the default option for grobcov)
     4507//  output:
     4508//          list, the canonical P-representation of the Normal and Non-Normal locus:
     4509//               The Normal locus has two kind of components: Normal and Special.
     4510//               The Non-normal locus has two kind of components: Accumulation and Degenerate.
     4511//               Normal components: for each point in the component,
     4512//               the number of solutions in the variables is finite, and
     4513//               the solutions depend on the point in the component if the component is not 0-dimensional.
     4514//               Special components:  for each point in the component,
     4515//               the number of solutions in the variables is finite,
     4516//               the component is not 0-dimensional, but the solutions do not depend on the
     4517//               values of the parameters in the component.
     4518//               Accumlation points: are 0-dimensional components for which it exist
     4519//               an infinite number of solutions.
     4520//               Degenerate components: are components of dimension greater than 0 for which
     4521//               for every point in the component there exist infinite solutions.
     4522//          The output components are given as
     4523//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
     4524//          The components are given in canonical P-representation of the subset.
     4525//               If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
     4526//               gives the depth of the component of the constructible set.
     4527proc locusdg(list GG, list #)
     4528  "USAGE:   locus(G);
     4529                 The input must be the grobcov  of a parametrical ideal
     4530  RETURN:  The  locus.
     4531                 The output are the components of two constructible subsets of the locus
     4532                 of the parametrical system.: Normal and Non-normal.
     4533                 These components are
     4534                 given as a list of  (pi,(pi1,..pis_i),type_i,level_i) varying i,
     4535                 where the p's are prime ideals, the type can be: Normal, Special,
     4536                 Accumulation, Degenerate.
     4537  NOTE:      It can only be called after computing the grobcov of the
     4538                 parametrical ideal in generic representation ('ext',0),
     4539                 which is the default.
     4540                 The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
     4541  KEYWORDS: geometrical locus, locus, loci.
     4542  EXAMPLE:  locusdg; shows an example"
     4543{
     4544  def R=basering;
     4545  int dd=2; int comment=0;
     4546  list DD=#;
     4547  if (size(DD)>1){comment=1;}
     4548  if(size(DD)>0){dd=DD[1];}
     4549  int i; int j; int k;
     4550  def B0=GG[1][2];
     4551  if (equalideals(B0,ideal(1)))
     4552  {
     4553    return(locus0dg(GG));
     4554  }
     4555  else
     4556  {
     4557    int n=nvars(R);
     4558    ideal vmov=var(n-1),var(n);
     4559    ideal N;
     4560    intvec xw; intvec yw;
     4561    for(i=1;i<=n-1;i++){xw[i]=0;}
     4562    xw[n]=1;
     4563    for(i=1;i<=n;i++){yw[i]=0;}
     4564    yw[n-1]=1;
     4565    poly px; poly py;
     4566    int te=1;
     4567    i=1;
     4568    while( te and i<=size(B0))
     4569    {
     4570      if((deg(B0[i],xw)==1) and (deg(B0[i])==1)){px=B0[i]; te=0;}
     4571      i++;
     4572    }
     4573    i=1; te=1;
     4574    while( te and i<=size(B0))
     4575    {
     4576      if((deg(B0[i],yw)==1) and (deg(B0[i])==1)){py=B0[i]; te=0;}
     4577      i++;
     4578    }
     4579    N=px,py;
     4580    setglobalrings();
     4581    te=indepparameters(N);
     4582    if(te)
     4583    {
     4584      string("locus detected that the mover must avoid point (",N,") in order to obtain the correct locus");
     4585      // eliminates segments of GG where N is contained in the basis
     4586      list nGP;
     4587      def GP=GG;
     4588      ideal BP;
     4589      for(j=1;j<=size(GP);j++)
     4590      {
     4591        te=1; k=1;
     4592        BP=GP[j][2];
     4593        while((te==1) and (k<=size(N)))
     4594        {
     4595          if(pdivi(N[k],BP)[1]!=0){te=0;}
     4596          k++;
     4597        }
     4598        if(te==0){nGP[size(nGP)+1]=GP[j];}
     4599      }
     4600   }
     4601  }
     4602  //"T_nGP="; nGP;
     4603  kill @RP; kill @P; kill @R;
     4604  return(locus0dg(nGP,dd));
     4605}
     4606example
     4607{"EXAMPLE:"; echo = 2;
     4608  ring R=(0,a,b),(x4,x3,x2,x1),dp;
     4609  ideal S=(x1-3)^2+(x2-1)^2-9,
     4610             (4-x2)*(x3-3)+(x1-3)*(x4-1),
     4611             (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
     4612             (4-x4)*a+(x3-3)*b+3*x4-4*x3,
     4613             (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
     4614  short=0;
     4615  locus(grobcov(S)); " ";
     4616  kill R;
     4617  ring R=(0,a,b),(x,y),dp;
     4618  short=0;
     4619  "Concoid";
     4620  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
     4621  "System="; S96; " ";
     4622  locusdg(grobcov(S96));
     4623  kill R; ring R=(0,x,y),(x1,x2),dp;
     4624  ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
     4625              (x1 - 5)^2 + (x2 - 2)^2 - 4,
     4626              -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
     4627 "System="; S;
     4628  locusdg(grobcov(S)); " ";
     4629  ideal S1=-(x - x1)*(x1 - 4) + (x2 - y)*(x2 - 4),
     4630                (x1 - 4)^2 + (x2 - 2)^2 - 4,
     4631                -2*(x1 - 4)*(2*x2 - y - 4) - 2*(x - 2*x1 + 4)*(x2 - 2);
     4632 "System="; S1;
     4633  locusdg(grobcov(S1));
     4634}
     4635
     4636// locusto: Transforms the output of locus to a string that
     4637//      can be read by different computational systems.
     4638// input:
     4639//     list L: The output of locus
     4640// output:
     4641//     string s: The output of locus converted to a string readable by other programs
     4642proc locusto(list L)
     4643"USAGE:   locusto(L);
     4644          The argument must be the output of locus of a parametrical ideal
     4645          It transforms the output into a string in standard form
     4646          readable in many languages (Geogebra).
     4647RETURN: The locus in string standard form
     4648NOTE:     It can only be called after computing the locus(grobcov(F)) of the
     4649              parametrical ideal.
     4650              The basering R, must be of the form Q[a,b,..][x,y,..].
     4651KEYWORDS: geometrical locus, locus, loci.
     4652EXAMPLE:  locusto; shows an example"
     4653{
     4654  int i; int j; int k;
     4655  string s="["; string sf="]"; string st=s+sf;
     4656  if(size(L)==0){return(st);}
     4657  ideal p;
     4658  ideal q;
     4659  for(i=1;i<=size(L);i++)
     4660  {
     4661    s=string(s,"[[");
     4662    for (j=1;j<=size(L[i][1]);j++)
     4663    {
     4664      s=string(s,L[i][1][j],",");
     4665    }
     4666    s[size(s)]="]";
     4667    s=string(s,",[");
     4668    for(j=1;j<=size(L[i][2]);j++)
     4669    {
     4670      s=string(s,"[");
     4671      for(k=1;k<=size(L[i][2][j]);k++)
     4672      {
     4673        s=string(s,L[i][2][j][k],",");
     4674      }
     4675      s[size(s)]="]";
     4676      s=string(s,",");
     4677    }
     4678    s[size(s)]="]";
     4679    s=string(s,"]");
     4680    if(size(L[i])>=3)
     4681    {
     4682      s[size(s)]=",";
     4683      s=string(s,string(L[i][3]),"]");
     4684    }
     4685    if(size(L[i])>=4)
     4686    {
     4687      s[size(s)]=",";
     4688      s=string(s,string(L[i][4]),"],");
     4689    }
     4690    s[size(s)]="]";
     4691    s=string(s,",");
     4692  }
     4693  s[size(s)]="]";
     4694  return(s);
     4695}
     4696example
     4697{"EXAMPLE:"; echo = 2;
     4698  ring R=(0,a,b),(x,y),dp;
     4699  short=0;
     4700  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
     4701  "System="; S96; " ";
     4702  locusto(locus(grobcov(S96)));
     4703}
     4704
     4705// Auxiliary routione
    43954706// locusdim
    43964707// input:
    … …  
    44374748}
    44384749
     4750// locusdgto: Transforms the output of locus to a string that
     4751//      can be read by different computational systems.
     4752// input:
     4753//     list L: The output of locus
     4754// output:
     4755//     string s: The output of locus converted to a string readable by other programs
     4756//                   Outputs only the relevant dynamical geometry components.
     4757//                   Without unnecessary parenthesis
     4758proc locusdgto(list LL)
     4759{
     4760  def RR=basering;
     4761  int i; int j; int k; int short0=short;
     4762  int te=0;
     4763  if(defined(@P)){te=1;}
     4764  else{setglobalrings();}
     4765  setring @P;
     4766  short=0;
     4767  def L=imap(RR,LL);
     4768  string s="["; string sf="]"; string st=s+sf;
     4769  if(size(L)==0){return(st);}
     4770  ideal p;
     4771  ideal q;
     4772  for(i=1;i<=size(L);i++)
     4773  {
     4774    if(L[i][3]=="Relevant")
     4775    {
     4776      s=string(s,"[[");
     4777      for (j=1;j<=size(L[i][1]);j++)
     4778      {
     4779        s=string(s,L[i][1][j],",");
     4780      }
     4781      s[size(s)]="]";
     4782      s=string(s,",[");
     4783      for(j=1;j<=size(L[i][2]);j++)
     4784      {
     4785        s=string(s,"[");
     4786        for(k=1;k<=size(L[i][2][j]);k++)
     4787        {
     4788          s=string(s,L[i][2][j][k],",");
     4789        }
     4790        s[size(s)]="]";
     4791        s=string(s,",");
     4792      }
     4793      s[size(s)]="]";
     4794      s=string(s,"]");
     4795      s[size(s)]="]";
     4796      s=string(s,",");
     4797    }
     4798  }
     4799  if(s=="["){s="[]";}
     4800  else{s[size(s)]="]";}
     4801  setring(RR);
     4802  short=short0;
     4803  if(te==0){kill @P; kill @R; kill @RP;}
     4804  return(s);
     4805}
    44394806
    44404807//********************* End locus ****************************
Note: See TracChangeset for help on using the changeset viewer.