Changeset f9fc81 in git for Singular/LIB


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

Fix.normal.callbug.in.genus
Location:
Singular/LIB
Files:
2 added
6 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 ****************************

  • Singular/LIB/modnormal.lib

    r9e8a6c3 rf9fc81  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="version modnormal.lib 4.0.0.0 Jun_2013 "; // $Id$
     2version="version modnormal.lib 4.0.0.0 Dec_2013 "; // $Id$
    33category = "Commutative Algebra";
    44info="
    … …  
    314314        modarguments[j-(k-1)*n2-sh] = list(I, L[j], condu, printTimings, #);
    315315      }
    316       modresults = parallelWaitAll("modpNormal", modarguments,
    317         list(list(list(ncores))));
     316      modresults = parallelWaitAll("modpNormal", modarguments, 0, ncores);
    318317      for(j = (k-1)*n2+1+sh; j <= k*n2+1; j++)
    319318      {

  • Singular/LIB/modstd.lib

    r9e8a6c3 rf9fc81  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="version modstd.lib 4.0.0.0 Jun_2013 "; // $Id$
     2version="version modstd.lib 4.0.0.0 Dec_2013 "; // $Id$
    33category = "Commutative Algebra";
    44info="
    … …  
    500500            arguments[i] = list("Id", p);
    501501         }
    502          parallelResults = parallelWaitAll("primeTest", arguments,
    503             list(list(list(ncores))));
     502         parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores);
    504503         for(i = size(arguments); i > 0; i--)
    505504         {
    … …  
    511510         neededPrimes = neededSize-size(L);
    512511      }
     512      kill Id;
    513513      if(size(L) > neededSize)
    514514      {
    … …  
    12391239            arguments_farey[i] = list(ideal(H[i]), N);
    12401240         }
    1241          results_farey = parallelWaitAll("farey", arguments_farey,
    1242                                          list(list(list(n1))));
     1241         results_farey = parallelWaitAll("farey", arguments_farey, 0, n1);
    12431242         for(i = ncols(H); i > 0; i--)
    12441243         {

  • Singular/LIB/normal.lib

    rbdda288 rf9fc81  
    24352435   if(size(#)!=0)
    24362436   {                              //uses the normalization to compute delta
    2437       list nor=normal(I);
     2437      list nor=normal(I,"withDelta");
    24382438      delt=nor[size(nor)][2];
    24392439      genus=genus-delt-deltainf;

  • Singular/LIB/parallel.lib

    r9e8a6c3 rf9fc81  
    1 ///////////////////////////////////////////////////////////////////
    2 version="version parallel.lib 4.0.0.0 Jun_2013 "; // $Id$
     1////////////////////////////////////////////////////////////////////
     2version="version parallel.lib 4.0.0.0 Dec_2013 "; // $Id$
    33category="General purpose";
    44info="
    5 LIBRARY:   parallel.lib  Tools for Parallelization
     5LIBRARY:   parallel.lib  An abstraction layer for parallel skeletons
     6
    67AUTHOR:    Andreas Steenpass, e-mail: steenpass@mathematik.uni-kl.de
    78
    89OVERVIEW:
    9 This library provides tools to do several computations in parallel. They
    10 are aimed at ordinary Singular users as well as authors of Singular
    11 libraries.
    12 @* Even without this library, it is possible to do execute self-defined
    13 Singular commands in parallel using @ref{link}, but the handling of
    14 such links can be quite tedious. With the pocedures described below,
    15 this can be done by one-line commands.
    16 @* There are many parallel 'skeletons' (i.e. ways in which parallel
    17 tasks rely upon and interact with each other). A few of them are already
    18 implemented. Future plans include an abstraction layer for modular
    19 techniques, 'worker farms', and parallel tests.
    20 
    21 SEE ALSO:  link, modstd_lib, assprimeszerodim_lib
    22 
    23 KEYWORDS:  parallel.lib; Parallelization; Links, user interface;
    24            Skeletons for parallelization; Distributed computing
     10This library provides implementations of several parallel 'skeletons' (i.e.
     11ways in which parallel tasks rely upon and interact with each other). It is
     12based on the library tasks.lib and aims at both ordinary Singular users as well
     13as authors of Singular libraries.
     14
     15KEYWORDS:  parallelization; parallel skeletons; distributed computing
     16
     17SEE ALSO:  resources_lib, tasks_lib, modstd_lib, modnormal_lib
    2518
    2619PROCEDURES:
    27   parallelWaitN(...)     execute several jobs in parallel,
    28                          wait for N of them to finish
    29   parallelWaitFirst(...) execute several jobs in parallel,
    30                          wait for the first to finish
    31   parallelWaitAll(...)   execute several jobs in parallel,
    32                          wait for all of them to finish
     20  parallelWaitN();      execute several jobs in parallel
     21                        and wait for N of them to finish
     22  parallelWaitFirst();  execute several jobs in parallel
     23                        and wait for the first to finish
     24  parallelWaitAll();    execute several jobs in parallel
     25                        and wait for all of them to finish
     26  parallelTestAND();    run several tests in parallel
     27                        and determine if they all succeed
     28  parallelTestOR();     run several tests in parallel
     29                        and determine if any of them succeeds
    3330";
    3431
    35 proc parallelWaitN(list commands, list args, int N, list #)
    36 "USAGE:  parallelWaitN(commands, args, N[, timeout, linktype, servers,
    37          maxmemory]); commands list, args list, N int, timeout int,
    38          linktype string, servers list, maxmemory intvec
    39 RETURN:  a list, containing the results of commands[i] applied to arg[i],
    40          i = 1, ..., size(commands).
    41          @* The procedure waits for N jobs to finish.
    42 
    43          @* OPTIONAL PARAMETERS:
    44 
    45             An optional timeout in ms can be provided. Default is 0 which
    46             disables the timeout.
    47 
    48             Supported linktypes are up to now \"ssi\" and \"mp\", see
    49             @ref{Ssi links}.
    50 
    51             Servers:
    52          @* Each server is given by a list containing the address of the server,
    53             the number of cores to use on this server and the command to start
    54             Singular.
    55          @* If the address is \"localhost\", the processes will be generated on
    56             the same machine using forks. If the command to start Singular is
    57             \"\" (the empty string), \"Singular\" will be used.
    58          @* Default is @code{list(\"localhost\", system(\"cpu\"), \"\")}.
    59          @* There are some obvious shortcuts for servers, e.g. \"myserver\" is
    60             a shortcut for
    61             @code{list(\"myserver\", [nb. of cores on myserver], \"\")}, or 3
    62             for @code{list(\"localhost\", 3, \"\")}.
    63 
    64             Memory limits:
    65          @* If an intvec maxmemory of size @code{size(commands)} is given, the
    66             i-th job will be killed if it uses more than maxmemory[i] MB of
    67             memory. If maxmemory[i] is 0, there will be no restraint for the
    68             i-th job. Default is @code{0:size(commands)}.
    69 NOTE:       The entries of the list commands must be strings.
    70          @* The entries of the list args must be lists.
    71          @* The returned list may contain more than N results if several jobs
    72             finished \"at the same time\". It may contain less than N results in
    73             the case of timeout or errors occurring.
    74          @* The check whether the jobs exceed the memory sizes given by
    75             maxmemory is only done from time to time. This feature is
    76             experimental and should be used with care.
    77 SEE ALSO: Ssi links, waitfirst, waitall
    78 KEYWORDS: parallelWaitN; Parallelization; Links, user interface;
    79           Skeletons for parallelization; Distributed computing
    80 EXAMPLE:  @code{example parallelWaitN;} shows an example."
    81 {
    82   // initialize the timer
    83   int oldtimerresolution = system("--ticks-per-sec");
    84   system("--ticks-per-sec", 1000);
    85   int t = rtimer;
    86 
    87   // auxiliary variables
    88   int i, j, m, tt;
    89 
    90   // read optional parameters
    91   list defaultserver = list("localhost", system("cpu"), "");
    92   list defaults = list(0, "ssi", list(defaultserver), 0:size(commands));
    93   for(i = 1; i <= size(defaults); i++)
    94   {
    95     if(typeof(#[i]) != typeof(defaults[i]))
    96     {
    97       # = insert(#, defaults[i], i-1);
    98     }
    99   }
    100   if(size(#) != size(defaults))
    101   {
    102     ERROR("wrong optional parameters");
    103   }
    104   for(j = size(#[3]); j > 0; j--)
    105   {
    106     if(typeof(#[3][j][1]) != typeof(defaultserver[1]))
    107     {
    108       #[3][j] = insert(#[3][j], defaultserver[1], 0);
    109     }
    110     defaultserver[3] = "";
    111     // only for ssi:tcp links, the default program is system("Singular")
    112     if(#[2] == "ssi" && #[3][j][1] != "localhost")
    113     {
    114       defaultserver[3] = system("Singular");
    115     }
    116     for(i = 2; i <= size(defaultserver); i++)
    117     {
    118       if(typeof(#[3][j][i]) != typeof(defaultserver[i]))
    119       {
    120         #[3][j] = insert(#[3][j], defaultserver[i], i-1);
    121       }
    122     }
    123     if(size(#[3][j]) != size(defaultserver))
    124     {
    125       ERROR("wrong declaration for server no. "+string(j));
    126     }
    127   }
    128   int timeout = #[1];
    129   string linktype = #[2];
    130   list servers = #[3];
    131   intvec maxmems = #[4];
    132 
    133   // error checking
    134   int njobs = size(commands);
    135   if(njobs != size(args))
    136   {
    137     ERROR("The number of commands does not match the number of lists"
    138          +newline+"of arguments.");
    139   }
    140   if(njobs == 0)
    141   {
    142     ERROR("no commands specified");
    143   }
    144   for(i = 1; i <= njobs; i++)
    145   {
    146     if(typeof(commands[i]) != "string")
    147     {
    148       ERROR("The first argument is not a list of strings.");
    149     }
    150     if(typeof(args[i]) != "list")
    151     {
    152       ERROR("The second argument is not a list of lists.");
    153     }
    154   }
    155   if(N < 0)
    156   {
    157     ERROR("The number of jobs which you want to wait for is negative.");
    158   }
    159   if(N > njobs)
    160   {
    161     ERROR("The number of jobs which you wnat to wait for is greater"
    162          +newline+"than the number of jobs itself.");
    163   }
    164   if(timeout < 0)
    165   {
    166     ERROR("The given timeout is negative.");
    167   }
    168   if(linktype != "ssi" && linktype != "mp")
    169   {
    170     ERROR("The given linktype is not recognized.");
    171   }
    172   int nservers = size(servers);
    173   if(nservers <= 0)
    174   {
    175     ERROR("no server specified");
    176   }
    177   for(i = 1; i <= nservers; i++)
    178   {
    179     if(servers[i][1] != "localhost")
    180     {
    181       if(system("sh", "ssh "+servers[i][1]+" exit"))
    182       {
    183         ERROR("Could not connect to server \""+servers[i][1]+"\"");
    184       }
    185     }
    186     if(servers[i][2] < 0)
    187     {
    188       ERROR("The number of cores to be used on server \""+servers[i][1]+"\""
    189            +newline+" is negative.");
    190     }
    191     if(servers[i][1] == "localhost")
    192     {
    193       int ncpus(i) = system("cpu");
    194     }
    195     else
    196     {
    197       //if(linktype == "ssi")
    198       //{
    199         link lcpu(i) = "ssi:tcp "+servers[i][1]+":"+servers[i][3];
    200       //}
    201       open(lcpu(i));
    202       write(lcpu(i), quote(system("cpu")));
    203       int ncpus(i) = read(lcpu(i));
    204       close(lcpu(i));
    205       kill lcpu(i);
    206     }
    207     if(servers[i][2] == 0)
    208     {
    209       servers[i][2] = ncpus(i);
    210     }
    211     else
    212     {
    213       if(servers[i][2] > ncpus(i))
    214       {
    215         ERROR("The number of cores to use on server \""+servers[i][1]+"\""
    216              +newline+"is greater than the number of available cores");
    217       }
    218     }
    219     if(servers[i][1] != "localhost")
    220     {
    221       if(system("sh", "ssh "+servers[i][1]+
    222                       " 'test -e `which "+servers[i][3]+"`'"))
    223       {
    224         ERROR("\""+servers[i][3]+"\" was not found on"
    225              +"\""+servers[i][1]+"\".");
    226       }
    227     }
    228   }
    229   if(size(maxmems) != njobs)
    230   {
    231     ERROR("The size of the intvec which specifies the maximal amount of memory"
    232          +newline+"to be used for each job does not match the number of jobs.");
    233   }
    234   int havemaxmem;
    235   for(i = 1; i <= njobs; i++)
    236   {
    237     if(maxmems[i] < 0)
    238     {
    239       ERROR("The maximal amount of memory to be used for job no. "+string(i)
    240            +"is negative.");
    241     }
    242     havemaxmem = havemaxmem+maxmems[i];
    243   }
    244 
    245   // skip those cores which won't be needed
    246   int nlinks;
    247   for(i = 1; i <= nservers; i++)
    248   {
    249     if(nlinks+servers[i][2] <= njobs)
    250     {
    251       nlinks = nlinks+servers[i][2];
    252     }
    253     else
    254     {
    255       if(nlinks == njobs)
    256       {
    257         servers = list(servers[1..(i-1)]);
    258       }
    259       else
    260       {
    261         servers = list(servers[1..i]);
    262         servers[i][2] = njobs-nlinks;
    263         nlinks = njobs;
    264       }
    265       nservers = size(servers);
    266     }
    267   }
    268 
    269   // open the links
    270   string server;
    271   int ncores;
    272   string program;
    273   int k = 1;   // the index of the link
    274   for(i = 1; i <= nservers; i++)
    275   {
    276     server = servers[i][1];
    277     ncores = servers[i][2];
    278     program = servers[i][3];
    279     for(j = 1; j <= ncores; j++)
    280     {
    281       if(server == "localhost")
    282       {
    283         //if(linktype == "ssi")
    284         //{
    285           link l(k) = "ssi:fork";
    286         //}
    287       }
    288       else
    289       {
    290         //if(linktype == "ssi")
    291         //{
    292           link l(k) = "ssi:tcp "+server+":"+program;
    293         //}
    294       }
    295       open(l(k));
    296       k++;
    297     }
    298   }
    299   list links = list(l(1..nlinks));
    300 
    301   // start a memory watchdog if needed
    302   if(havemaxmem)
    303   {
    304     //if(linktype == "ssi")
    305     //{
    306       link mempatrol = "ssi:fork";
    307     //}
    308     open(mempatrol);
    309     write(mempatrol, quote(watchlinks()));
    310     links = insert(links, mempatrol, nlinks);
    311   }
    312   int nkilled;   // the number of jobs killed by the mempatrol
    313 
    314   // distribute work to the links
    315   k = 1; // from here on, k is the index of the next job which must be
    316          // distributed to some link
    317   intvec assignment = 0:nlinks;  // link number i is currently doing
    318                                  // job number assignment[i]
    319   int pid;
    320   for(i = 1; i <= nlinks; i++)
    321   {
    322     write(l(i), quote(execute("option(noredefine);")));
    323     read(l(i));
    324     write(l(i), quote(execute("def result;")));
    325     read(l(i));
    326     write(l(i), quote(execute("list currentargs;")));
    327     read(l(i));
    328     if(status(l(i), "mode", "fork"))
    329     {
    330       write(l(i), quote(currentargs = args[eval(k)]));
    331     }
    332     else
    333     {
    334       write(l(i), quote(currentargs = eval(args[k])));
    335     }
    336     read(l(i));
    337     if(maxmems[k] > 0)
    338     {
    339       m = i;
    340       for(j = 1; j <= nservers; j++)
    341       {
    342         if(servers[j][2] > m)
    343         {
    344           server = servers[j][1];
    345           break;
    346         }
    347         else
    348         {
    349           m = m-servers[j][2];
    350         }
    351       }
    352       write(l(i), quote(system("pid")));
    353       pid = read(l(i));
    354       write(mempatrol, list(server, pid, i, maxmems[k]));
    355     }
    356     write(l(i), quote(execute("result = "+eval(commands[k])
    357       +"("+argsToString("currentargs", size(currentargs))+");")));
    358     assignment[i] = k;
    359     k++;
    360   }
    361 
    362   // distribute the rest of the work
    363   list results;
    364   for(i = njobs; i > 0; i--)
    365   {
    366     results[i] = list();  // TODO: What if the result of one of the commands is
    367                           // list()?
    368   }
    369   int nfinished;  // the number of finished jobs
    370   int wait;   // the index of the link which is finished, or 0 for timeout
    371   while(k <= njobs && nfinished < N-1)
    372   {
    373     if(timeout == 0)
    374     {
    375       wait = waitfirst(links);
    376     }
    377     else
    378     {
    379       tt = timeout-(rtimer-t);
    380       if(tt < 0)
    381       {
    382         wait = waitfirst(links, 0);
    383         wait = 0;
    384       }
    385       else
    386       {
    387         wait = waitfirst(links, tt);
    388       }
    389     }
    390     if(wait == -1)
    391     {
    392       ERROR("All links crashed.");
    393     }
    394     if(wait)
    395     {
    396       if(wait == nlinks+1)
    397       {
    398         wait = read(mempatrol);
    399         close(l(wait));
    400         open(l(wait));
    401         results[assignment[wait]] = "out of memory";
    402         nkilled++;
    403       }
    404       else
    405       {
    406         read(l(wait));
    407         write(l(wait), quote(result));
    408         results[assignment[wait]] = read(l(wait));
    409         if(maxmems[assignment[wait]] > 0)
    410         {
    411           write(mempatrol, assignment[wait]);
    412         }
    413         nfinished++;
    414       }
    415       if(status(l(wait), "mode", "fork"))
    416       {
    417         write(l(wait), quote(currentargs = args[eval(k)]));
    418       }
    419       else
    420       {
    421         write(l(wait), quote(currentargs = eval(args[k])));
    422       }
    423       read(l(wait));
    424       if(maxmems[k] > 0)
    425       {
    426         m = wait;
    427         for(j = 1; j <= nservers; j++)
    428         {
    429           if(servers[j][2] > m)
    430           {
    431             server = servers[j][1];
    432             break;
    433           }
    434           else
    435           {
    436             m = m-servers[j][2];
    437           }
    438         }
    439         write(l(wait), quote(system("pid")));
    440         pid = read(l(wait));
    441         write(mempatrol, list(server, pid, wait, maxmems[k]));
    442       }
    443       write(l(wait), quote(execute("def result = "+eval(commands[k])
    444         +"("+argsToString("currentargs", size(currentargs))+");")));
    445       assignment[wait] = k;
    446       k++;
    447     }
    448     else
    449     {
    450       break;
    451     }
    452   }
    453 
    454   // collect the rest of the results
    455   while(nfinished < N && nfinished+nkilled < njobs)
    456   {
    457     if(timeout == 0)
    458     {
    459       wait = waitfirst(links);
    460     }
    461     else
    462     {
    463       tt = timeout-(rtimer-t);
    464       if(tt < 0)
    465       {
    466         wait = waitfirst(links, 0);
    467         wait = 0;
    468       }
    469       else
    470       {
    471         wait = waitfirst(links, tt);
    472       }
    473     }
    474     if(wait == -1)
    475     {
    476       ERROR("All links crashed.");
    477     }
    478     if(wait)
    479     {
    480       if(wait == nlinks+1)
    481       {
    482         wait = read(mempatrol);
    483         close(l(wait));
    484         results[assignment[wait]] = "out of memory";
    485         nkilled++;
    486       }
    487       else
    488       {
    489         read(l(wait));
    490         write(l(wait), quote(result));
    491         results[assignment[wait]] = read(l(wait));
    492         if(maxmems[assignment[wait]] > 0)
    493         {
    494           write(mempatrol, assignment[wait]);
    495         }
    496         nfinished++;
    497       }
    498     }
    499     else
    500     {
    501       break;
    502     }
    503   }
    504 
    505   //close all links
    506   for(i = 1; i <= nlinks; i++)
    507   {
    508     if(status(l(i), "read", "ready"))
    509     {
    510       read(l(i));
    511       write(l(i), quote(result));
    512       results[assignment[i]] = read(l(i));
    513     }
    514     close(l(i));
    515   }
    516   if(havemaxmem)
    517   {
    518     close(mempatrol);
    519   }
    520 
    521   system("--ticks-per-sec", oldtimerresolution);
    522   return(results);
    523 }
    524 example
    525 {
    526   "EXAMPLE:"; echo = 2;
    527   LIB "primdec.lib";
    528   ring r = 0, (x,y,z), lp;
    529   ideal i = z8+z6+4z5+4z3+4z2+4, y-z2;
    530   ideal j = 3x3y+x3+xy3+y2z2, 2x3z-xy-xz3-y4-z2, 2x2yz-2xy2+xz2-y4;
    531   list commands = list("std", "primdecGTZ", "primdecSY",
    532                        "std", "primdecGTZ", "primdecSY");
    533   list args = list(list(i), list(i), list(i), list(j), list(j), list(j));
    534   parallelWaitN(commands, args, 3);
    535 }
    536 
    537 proc parallelWaitFirst(list commands, list args, list #)
    538 "USAGE:  parallelWaitFirst(commands, args[, timeout, linktype, servers,
    539          maxmemory]); commands list, args list, timeout int, linktype string,
    540          servers list, maxmemory intvec
    541 RETURN:  a list, containing at least one (if no timeout occurs) of the results
    542          of commands[i] applied to arg[i], i = 1, ..., size(commands).
    543          @* The command
    544          @code{parallelWaitFirst(list commands, list args, list #)} is
    545          synonymous to
    546          @code{parallelWaitN(list commands, list args, 1, list #)}. See
    547          @ref{parallelWaitN} for details on optional arguments and other
    548          remarks.
    549 SEE ALSO: Ssi links, waitfirst
    550 KEYWORDS: parallelWaitFirst; Parallelization; Links, user interface;
    551           Skeletons for parallelization; Distributed computing
    552 EXAMPLE:  @code{example parallelWaitFirst;} shows an example."
    553 {
    554   return(parallelWaitN(commands, args, 1, #));
    555 }
    556 example
    557 {
    558   "EXAMPLE:"; echo = 2;
    559   LIB "primdec.lib";
    560   ring r = 0, (x,y,z), lp;
    561   ideal i = z8+z6+4z5+4z3+4z2+4, y-z2;
    562   list commands = list("primdecGTZ", "primdecSY");
    563   list args = list(list(i), list(i));
    564   parallelWaitFirst(commands, args);
    565 }
    566 
    567 proc parallelWaitAll(def commands, list args, list #)
    568 "USAGE:  parallelWaitAll(commands, args[, timeout, linktype, servers,
    569          maxmemory]); commands list or string, args list, timeout int,
    570          linktype string, servers list, maxmemory intvec
    571 RETURN:  a list, containing the results of commands[i] applied to arg[i],
    572          i = 1, ..., size(commands).
    573          @* The command
    574          @code{parallelWaitAll(list commands, list args, list #)} is
    575          synonymous to
    576          @code{parallelWaitN(list commands, list args, size(args), list #)}. See
    577          @ref{parallelWaitN} for details on optional arguments and other
    578          remarks.
    579          If commands is of type string, this is a shortcut for a list of size
    580          @code{size(args)} whose entries are just this string.
    581 SEE ALSO: Ssi links, waitall
    582 KEYWORDS: parallelWaitAll; Parallelization; Links, user interface;
    583           Skeletons for parallelization; Distributed computing
    584 EXAMPLE:  @code{example parallelWaitAll;} shows an example."
    585 {
    586   if(typeof(commands) != "list" && typeof(commands) != "string")
    587   {
    588     ERROR("invalid type of first argument");
    589   }
    590   if(typeof(commands) == "list")
    591   {
    592     return(parallelWaitN(commands, args, size(args), #));
    593   }
    594   else
    595   {
    596     list cmds;
    597     for(int i = size(args); i > 0; i--)
    598     {
    599       cmds[i] = commands;
    600     }
    601     return(parallelWaitN(cmds, args, size(args), #));
    602   }
    603 }
    604 example
    605 {
    606   "EXAMPLE:"; echo = 2;
    607   ring r = 0, (x,y,z), dp;
    608   ideal i1 = z8+z6+4z5+4z3+4z2+4, y-z2;
    609   ideal i2 = x10+x9y2, y8-x2y7;
    610   ideal i3 = x3-2xy, x2y-2y2+x;
    611   string command = "std";
    612   list args = list(list(i1), list(i2), list(i3));
    613   parallelWaitAll(command, args);
    614 }
    615 
    616 // TODO
    617 /// http://arxiv.org/abs/1005.5663v2
    618 static proc doModular(command, args, proc deleteUnlucksPrimes, proc testInChar0)
    619 {
    620 }
    621 
    622 // TODO
    623 /* worker farm */
    624 static proc Create() {}
    625 
    626 /* auxiliary procedures */
    627 static proc watchlinks()
    628 {
    629   list parent = list(mempatrol);
    630   list watchedlinks;
    631   int wait;
    632   int i, sys;
    633   while(1)
    634   {
    635     if(size(watchedlinks) == 0)
    636     {
    637       wait = waitall(parent);
    638     }
    639     else
    640     {
    641       wait = waitall(parent, 10000);
    642     }
    643     if(wait == -1)
    644     {
    645       ERROR("The main process crashed.");
    646     }
    647     if(wait)
    648     {
    649       def query = read(mempatrol);
    650       if(typeof(query) == "list")
    651       {
    652         watchedlinks = insert(watchedlinks, query);
    653       }
    654       else // in this case, typeof(query) is assumed to be "int", the
    655            // index of the link
    656       {
    657         for(i = size(watchedlinks); i > 0; i--)
    658         {
    659           if(watchedlinks[i][3] == query)
    660           {
    661             watchedlinks = delete(watchedlinks, i);
    662             break;
    663           }
    664         }
    665       }
    666     }
    667     for(i = size(watchedlinks); i > 0; i--)
    668     {
    669       if(getusedmemory(watchedlinks[i][1], watchedlinks[i][2])
    670            > watchedlinks[i][4])
    671       {
    672         if(watchedlinks[i][1] == "localhost")
    673         {
    674           sys = system("sh", "kill "+string(watchedlinks[i][2]));
    675         }
    676         else
    677         {
    678           sys = system("sh", "ssh "+watchedlinks[i][1]+" kill "
    679                              +string(watchedlinks[i][2]));
    680         }
    681         write(mempatrol, watchedlinks[i][3]);
    682         watchedlinks = delete(watchedlinks, i);
    683       }
    684     }
    685   }
    686 }
    687 
    688 static proc getusedmemory(string server, int pid)
    689 {
    690   string s;
    691   if(server == "localhost")
    692   {
    693     s = read("|: grep VmSize /proc/"+string(pid)+"/status "+
    694              "| awk '{ print $2; }'");
    695   }
    696   else
    697   {
    698     s = read("|: ssh "+server+" grep VmSize /proc/"+string(pid)+"/status "+
    699              "| awk '{ print $2; }'");
    700   }
    701   bigint b;
    702   execute("b = "+s+";");
    703   int i = int(b/1000);
    704   return(i);
    705 }
    706 
    707 static proc argsToString(string name, int length)
    708 {
    709   string arglist;
    710   if(length > 0) {
    711     arglist = name+"[1]";
    712   }
    713   int i;
    714   for(i = 2; i <= length; i++) {
    715     arglist = arglist+", "+name+"["+string(i)+"]";
    716   }
    717   return(arglist);
    718 }
     32LIB "tasks.lib";
     33
     34proc parallelWaitN(alias list commands, alias list args, int N, list #)
     35"USAGE:   parallelWaitN(commands, arguments, N[, timeout]); commands list,
     36          arguments list, N int, timeout int
     37RETURN:   a list, containing the results of commands[i] applied to
     38          arguments[i], i = 1, ..., size(arguments).
     39       @* The procedure waits for N jobs to finish.
     40       @* An optional timeout in ms can be provided. Default is 0 which
     41          disables the timeout.
     42NOTE:     The entries of the list commands must be strings. The entries of the
     43          list arguments must be lists.
     44       @* The type of any entry of the returned list whose corresponding task
     45          did not finish (due to timeout or error) is \"none\".
     46       @* The returned list may contain more than N results if several jobs
     47          finished \"at the same time\". It may contain less than N results in
     48          the case of timeout or errors occurring.
     49SEE ALSO: parallelWaitAll, parallelWaitFirst, tasks_lib
     50EXAMPLE:  example parallelWaitN; shows an example"
     51{
     52    // auxiliary variables
     53    int i;
     54
     55    // read optional parameters
     56    int timeout;
     57    int ncores;   // obsolete, but kept for compatibility with old libraries
     58    if (size(#) > 0) {
     59        if (typeof(#[1]) != "int") {
     60            ERROR("wrong optional parameters");
     61        }
     62        timeout = #[1];
     63        if (size(#) > 1) {
     64            if (size(#) > 2 || typeof(#[2]) != "int") {
     65                ERROR("wrong optional parameters");
     66            }
     67            ncores = #[2];
     68        }
     69    }
     70
     71    // apply wrapper for obsolete optional parameter ncores
     72    if (ncores) {
     73        list semaphore_save = Resources::setcores_subtree(ncores);
     74    }
     75
     76    // error checking
     77    int njobs = size(commands);
     78    if (njobs != size(args)) {
     79        ERROR("The number of commands does not match the number of lists"
     80            +newline+"of arguments.");
     81    }
     82    if (njobs == 0) {
     83        ERROR("no commands specified");
     84    }
     85    for (i = 1; i <= njobs; i++) {
     86        if (typeof(commands[i]) != "string") {
     87            ERROR("The first argument is not a list of strings.");
     88        }
     89        if (typeof(args[i]) != "list") {
     90            ERROR("The second argument is not a list of lists.");
     91        }
     92    }
     93
     94    // compute the tasks
     95    for (i = 1; i <= njobs; i++) {
     96        task t(i) = commands[i], args[i];
     97    }
     98    startTasks(t(1..njobs));
     99    list indices = waitTasks(list(t(1..njobs)), N, timeout);
     100
     101    // wrap back to saved semaphore
     102    if (ncores) {
     103        Resources::resetcores_subtree(semaphore_save);
     104    }
     105
     106    // return results
     107    list results;
     108    for (i = size(indices); i > 0; i--) {
     109        results[indices[i]] = getResult(t(indices[i]));
     110    }
     111    for (i = 1; i <= njobs; i++) {
     112        killTask(t(i));
     113    }
     114    return(results);
     115}
     116example
     117{
     118    "EXAMPLE:";
     119    echo = 2;
     120    ring R = 0, (x,y,z), lp;
     121    ideal I = 3x3y+x3+xy3+y2z2, 2x3z-xy-xz3-y4-z2, 2x2yz-2xy2+xz2-y4;
     122    ideal J = x10+x9y2, x2y7-y8;
     123    list commands = list("std", "std");
     124    list arguments = list(list(I), list(J));
     125    parallelWaitN(commands, arguments, 1);
     126}
     127
     128proc parallelWaitFirst(alias list commands, alias list args, list #)
     129"USAGE:   parallelWaitFirst(commands, args[, timeout]); commands list,
     130          arguments list, timeout int
     131RETURN:   a list, containing at least one (if no timeout occurs) of the results
     132          of commands[i] applied to arguments[i], i = 1, ..., size(arguments).
     133       @* The command @code{parallelWaitFirst(commands, arguments[, timeout])}
     134          is synonymous to
     135          @code{parallelWaitN(commands, arguments, 1[, timeout])}. See
     136          @ref{parallelWaitN} for details on optional arguments and other
     137          remarks.
     138SEE ALSO: parallelWaitN, parallelWaitAll, tasks_lib
     139EXAMPLE:  example parallelWaitFirst; shows an example"
     140{
     141    return(parallelWaitN(commands, args, 1, #));
     142}
     143example
     144{
     145    "EXAMPLE:";
     146    echo = 2;
     147    ring R = 0, (x,y,z), lp;
     148    ideal I = 3x3y+x3+xy3+y2z2, 2x3z-xy-xz3-y4-z2, 2x2yz-2xy2+xz2-y4;
     149    ideal J = x10+x9y2, x2y7-y8;
     150    list commands = list("std", "std");
     151    list arguments = list(list(I), list(J));
     152    parallelWaitFirst(commands, arguments);
     153}
     154
     155proc parallelWaitAll(def commands, alias list args, list #)
     156"USAGE:   parallelWaitAll(commands, arguments[, timeout]); commands list or
     157          string, arguments list, timeout int
     158RETURN:   a list, containing the results of commands[i] applied to
     159          arguments[i], i = 1, ..., size(arguments).
     160       @* The command @code{parallelWaitAll(commands, arguments[, timeout])} is
     161          synonymous to @code{parallelWaitN(commands, arguments,
     162          size(arguments)[, timeout])}. See @ref{parallelWaitN} for details on
     163          optional arguments and other remarks.
     164NOTE:     As a shortcut, @code{commands} can be a string. This is synonymous to
     165          providing a list of @code{size(arguments)} copies of this string.
     166SEE ALSO: parallelWaitFirst, parallelWaitN, tasks_lib
     167EXAMPLE:  example parallelWaitAll; shows an example"
     168{
     169    if (typeof(commands) != "list" && typeof(commands) != "string") {
     170        ERROR("invalid type of first argument");
     171    }
     172    if (typeof(commands) == "list") {
     173        return(parallelWaitN(commands, args, size(args), #));
     174    }
     175    else {
     176        list cmds;
     177        for (int i = size(args); i > 0; i--) {
     178            cmds[i] = commands;
     179        }
     180        return(parallelWaitN(cmds, args, size(args), #));
     181    }
     182}
     183example
     184{
     185    "EXAMPLE:";
     186    echo = 2;
     187    ring R = 0, (x,y,z), dp;
     188    ideal I1 = z8+z6+4z5+4z3+4z2+4, -z2+y;
     189    ideal I2 = x9y2+x10, x2y7-y8;
     190    ideal I3 = x3-2xy, x2y-2y2+x;
     191    string command = "std";
     192    list arguments = list(list(I1), list(I2), list(I3));
     193    parallelWaitAll(command, arguments);
     194}
     195
     196proc parallelTestAND(def commands, alias list args, list #)
     197"USAGE:   parallelTestAND(commands, arguments[, timeout]); commands list or
     198          string, arguments list, timeout int
     199RETURN:   1, if commands[i] applied to arguments[i] is not equal to zero for
     200          all i = 1, ..., size(arguments);
     201          0, otherwise.
     202       @* An optional timeout in ms can be provided. Default is 0 which
     203          disables the timeout. In case of timeout, -1 is returned.
     204NOTE:     The entries of the list commands must be strings. The entries of the
     205          list arguments must be lists.
     206       @* commands[i] applied to arguments[i] must evaluate to an integer for
     207          i = 1, ..., size(arguments).
     208       @* As a shortcut, @code{commands} can be a string. This is synonymous to
     209          providing a list of @code{size(arguments)} copies of this string.
     210SEE ALSO: parallelTestOR, tasks_lib
     211EXAMPLE:  example parallelTestAND; shows an example"
     212{
     213    // note: this can be improved
     214    list results = parallelWaitAll(commands, args, #);
     215    int i;
     216    for (i = size(args); i > 0; i--) {
     217        if (typeof(results[i]) != "int" && typeof(results[i]) != "none") {
     218            ERROR("result no. "+string(i)+" not of type int");
     219        }
     220    }
     221    for (i = size(args); i > 0; i--) {
     222        if (typeof(results[i]) == "none") {   // timeout
     223            return(-1);
     224        }
     225    }
     226    for (i = size(results); i > 0; i--) {
     227        if (!results[i]) {
     228            return(0);
     229        }
     230    }
     231    return(1);
     232}
     233example
     234{
     235    "EXAMPLE:";
     236    echo = 2;
     237    ring R = 0, (x,y,z), dp;
     238    ideal I = x, y, z;
     239    intvec v = 0:3;
     240    list l = list(I, v);
     241    module m1 = x*gen(1);
     242    module m2;
     243    string command = "size";
     244    list arguments1 = list(list(I), list(v), list(l), list(m1));
     245    list arguments2 = list(list(I), list(v), list(l), list(m2));
     246    // test if all the arguments have non-zero size
     247    parallelTestAND(command, arguments1);
     248    parallelTestAND(command, arguments2);
     249}
     250
     251proc parallelTestOR(def commands, alias list args, list #)
     252"USAGE:   parallelTestOR(commands, arguments[, timeout]); commands list or
     253          string, arguments list, timeout int
     254RETURN:   1, if commands[i] applied to arguments[i] is not equal to zero for
     255          any i = 1, ..., size(arguments);
     256          0, otherwise.
     257       @* An optional timeout in ms can be provided. Default is 0 which
     258          disables the timeout. In case of timeout, -1 is returned.
     259NOTE:     The entries of the list commands must be strings. The entries of the
     260          list arguments must be lists.
     261       @* commands[i] applied to arguments[i] must evaluate to an integer for
     262          i = 1, ..., size(arguments).
     263       @* As a shortcut, @code{commands} can be a string. This is synonymous to
     264          providing a list of @code{size(arguments)} copies of this string.
     265SEE ALSO: parallelTestAND, tasks_lib
     266EXAMPLE:  example parallelTestAND; shows an example"
     267{
     268    // note: this can be improved
     269    list results = parallelWaitAll(commands, args, #);
     270    int i;
     271    for (i = size(args); i > 0; i--) {
     272        if (typeof(results[i]) != "int" && typeof(results[i]) != "none") {
     273            ERROR("result no. "+string(i)+" not of type int");
     274        }
     275    }
     276    for (i = size(args); i > 0; i--) {
     277        if (typeof(results[i]) == "none") {   // timeout
     278            return(-1);
     279        }
     280    }
     281    for (i = size(results); i > 0; i--) {
     282        if (results[i]) {
     283            return(1);
     284        }
     285    }
     286    return(0);
     287}
     288example
     289{
     290    "EXAMPLE:";
     291    echo = 2;
     292    ring R = 0, (x,y,z), dp;
     293    ideal I;
     294    string s;
     295    list l;
     296    module m1 = x*gen(1);
     297    module m2;
     298    string command = "size";
     299    list arguments1 = list(list(I), list(s), list(l), list(m1));
     300    list arguments2 = list(list(I), list(s), list(l), list(m2));
     301    // test if any of the arguments has non-zero size
     302    parallelTestOR(command, arguments1);
     303    parallelTestOR(command, arguments2);
     304}
     305

  • Singular/LIB/symodstd.lib

    r9e8a6c3 rf9fc81  
    11////////////////////////////////////////////////////////////////////////////
     2version="version symodstd.lib 4.0.0.0 Dec_2013 "; // $Id$
    23category = "Commutative Algebra";
    34info="
    … …  
    723724            arguments[i] = list("I", p, k);
    724725         }
    725          parallelResults = parallelWaitAll("divPrimeTest", arguments,
    726             list(list(list(ncores))));
     726         parallelResults = parallelWaitAll("divPrimeTest", arguments, 0,
     727            ncores);
    727728         for(i = size(arguments); i > 0; i--)
    728729         {
    … …  
    14361437            arguments_farey[i] = list(ideal(H[i]), N);
    14371438         }
    1438          results_farey = parallelWaitAll("farey", arguments_farey,
    1439                                          list(list(list(n1))));
     1439         results_farey = parallelWaitAll("farey", arguments_farey, 0, n1);
    14401440         for(i = ncols(H); i > 0; i--)
    14411441         {
Note: See TracChangeset for help on using the changeset viewer.