Changeset 6e3023a in git
- Timestamp:
- Dec 18, 2013, 2:57:02 PM (9 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
- Children:
- 3a9e92a717b194d4a06a429b452258c0269579076235be46a7d664256bb99b1bd6340cd8ad800e2c
- Parents:
- 74c446222122656dc7859919b51f3806768e7fd331b00db0e867fa73d4415e16969423de4337f6f7
- Files:
-
- 30 edited
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
.gitignore
r31b00db r6e3023a 8 8 *.lo 9 9 *.la 10 *.gcno 11 *.gcda 10 12 *~ 11 13 .libs -
Singular/LIB/COPYING
r31b00db r6e3023a 134 134 paramet.lib Thomas Keilen keilen@mathematik.uni-kl.de 135 135 perron.lib Oleksandr Motsak motsak at mathematik dot uni minus kl dot de 136 phindex.lib 136 phindex.lib 137 137 pointid.lib Stefan Steidel steidel@mathematik.uni-kl.de 138 138 poly.lib Olaf Bachmann obachman@mathematik.uni-kl.de 139 139 Gert-Martin Greuel greuel@mathematik.uni-kl.de 140 140 Anne Fruehbis-Krueger anne@mathematik.uni-kl.de 141 oldpolymake.libThomas Markwig keilen@mathematik.uni-kl.de141 polymake.lib Thomas Markwig keilen@mathematik.uni-kl.de 142 142 presolve.lib Gert-Martin Greuel greuel@mathematik.uni-kl.de 143 143 primdec.lib Gerhard Pfister pfister@mathematik.uni-kl.de -
Singular/LIB/divisors.lib
r31b00db r6e3023a 335 335 ASSUME: A is a divisor on X. 336 336 RETURN: a list with a basis of the space of global sections of D. 337 THEORY: We assume that the qring of X satisfies the S2-condition . We compute sat((f*J) : I) /f337 THEORY: We assume that the qring of X satisfies the S2-condition and that X is smooth. We compute sat((f*J) : I) /f 338 338 where D = (I)-(J). 339 339 KEYWORDS: divisors … … 370 370 " 371 371 { 372 return( deg(std(A.num))-deg(std(A.den)));372 return( mult(std(A.num))-mult(std(A.den))); 373 373 } 374 374 example -
Singular/LIB/gitfan.lib
r31b00db r6e3023a 1 1 /////////////////////////////////////////////////////////////////// 2 version="version gitfan.lib 4.0.0.0 Jun_2013 "; // $Id$2 version="version gitfan.lib 4.0.0.0 Jun_2013 "; // $Id$ 3 3 category="Algebraic Geometry"; 4 4 info=" … … 27 27 28 28 LIB "parallel.lib"; // for parallelWaitAll 29 LIB "general.lib";30 29 31 30 //////////////////////////////////////////////////// 32 33 proc int2face(int n, int r) 31 proc mod_init() 32 { 33 LIB"gfanlib.so"; 34 } 35 36 static proc int2face(int n, int r) 34 37 { 35 38 int k = r-1; … … 40 43 while(2^k > n0){ 41 44 k--; 45 //v[size(v)+1] = 0; 42 46 } 43 47 … … 52 56 ///////////////////////////////// 53 57 54 proc isAface(ideal a, intvec gam0 , int n)58 proc isAface(ideal a, intvec gam0) 55 59 "USAGE: isAface(a,gam0); a: ideal, gam0:intvec 56 60 PURPOSE: Checks whether the face of the positive orthant, … … 62 66 " 63 67 { 68 poly pz; 69 64 70 // special case: gam0 is the zero-cone: 65 71 if (size(gam0) == 1 and gam0[1] == 0){ 66 poly pz;ideal G;72 ideal G; 67 73 68 74 // is an a-face if and only if RL0(0,...,0) = const … … 87 93 } 88 94 89 // ring is too big: Switch to KK[T_i | e_i\in gam0] instead: 90 intvec w=ringlist(basering)[3][1][2]; 91 intvec w0; 95 96 // ring is too big: Switch to KK[T_i; e_i\in gam0] instead: 92 97 string initNewRing = "ring Rgam0 = 0,("; 93 98 94 99 for (int i=1; i<size(gam0); i++){ 95 100 initNewRing = initNewRing + string(var(gam0[i])) + ","; 96 w0[i]=w[gam0[i]]; 97 } 98 99 w0 = w0,w[gam0[i]],1; 100 initNewRing = initNewRing + string(var(gam0[size(gam0)])) + ",U),Wp("+string(w0)+");"; 101 } 102 103 initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;"; 101 104 def R = basering; 102 105 execute(initNewRing); 103 106 104 107 ideal agam0 = imap(R,a); 105 106 for (i = 1; i<=size(agam0); i=i+1)107 {108 if (size(agam0[i]) == 1)109 { return(0,0); }110 }111 108 112 109 poly p = var(1); // first entry of g; p = prod T_i with i element of g … … 117 114 118 115 agam0 = agam0, p - 1; // rad-membership 119 agam0 = timeStd(agam0,5); 120 // "timestd finished after: "+string(timer-t); 121 // int timeout = 0; 122 if (attrib(agam0,"isSB") < 1) 123 { 124 kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R; 125 return(0,1); 126 // // "timestd failed in "+string(gam0)+", falling back to saturation!"; 127 // setring R; kill Rgam0; 128 129 // initNewRing = "ring Rgam0 = 0,("; 130 131 // w0 = 0; 132 // for (i=1; i<size(gam0); i++){ 133 // initNewRing = initNewRing + string(var(gam0[i])) + ","; 134 // w0[i]=w[gam0[i]]; 135 // } 136 137 // w0 = w0,w[gam0[i]]; 138 // initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),Wp("+string(w0)+");"; 139 // execute(initNewRing); 140 141 // ideal G = imap(R,a); 142 143 // timeout = 1; 144 // int t=rtimer; 145 // for(int k=nvars(basering); k>0; k--) { G=sat(G,var(k))[1]; } 146 // t = (rtimer - t) div 1000; 147 // string(n)+": saturation successful after "+string(t)+" with: "+string(G<>1); 148 } 116 ideal G = std(agam0); 149 117 150 118 // does G contain 1?, i.e. is G = 1? 151 if(agam0 <> 1) { 152 kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R; 153 return(1,0); // true 154 } 155 156 kill agam0; kill Rgam0; kill initNewRing; kill w; kill w0; setring R; kill R; 157 return(0,0); // false 119 if(G <> 1) { 120 return(1); // true 121 } 122 123 return(0); // false 158 124 } 159 125 example … … 172 138 } 173 139 174 175 proc isAfaceNonZero(ideal a, intvec gam0)176 {177 string initNewRing = "ring Rgam0 = 0,(";178 for (int i=1; i<size(gam0); i++)179 { initNewRing = initNewRing + string(var(gam0[i])) + ","; }180 initNewRing = initNewRing + string(var(gam0[size(gam0)])) + "),dp;";181 def R = basering;182 execute(initNewRing);183 184 ideal agam0 = imap(R,a);185 186 for ( i = 1; i<=size(agam0); i=i+1)187 { if (size(agam0[i]) == 1) { return(0); } }188 189 poly p = var(1);190 for ( i = 2; i <= nvars(basering); i++ )191 { p = p * var(i); }192 193 agam0 = agam0, p - 1;194 ideal G = std(agam0);195 196 if(G <> 1)197 { return(1); }198 199 return(0);200 }201 140 //////////////////////////////////////////////////// 202 141 … … 206 145 list AF; 207 146 208 for(i = start; i <= end; i =i+1){147 for(i = start; i <= end; i++){ 209 148 if(i < 2^r){ 210 string(i)+" "+string(size(AF));211 149 gam0 = int2face(i,r); 212 150 … … 227 165 228 166 proc afaces(ideal a, list #) 229 "USAGE: afaces(a, b, c); a: ideal, b: int, c: int167 "USAGE: afaces(a, b, c); a: ideal, d: int, c: int 230 168 PURPOSE: Returns a list of all a-faces (represented by intvecs). 231 169 Moreover, it is possible to specify a dimensional bound b, … … 272 210 273 211 // do remaining ones: 274 for(i = ncores * step +1; i < 2^r; i =i+1){212 for(i = ncores * step +1; i < 2^r; i++){ 275 213 "another one needed"; 276 214 gam0 = int2face(i,r); … … 285 223 } 286 224 287 int l;288 225 // read out afaces of out into AF: 289 for( l = 1; l <= size(out); l++){290 AF = AF + out[ l];226 for(i = 1; i <= size(out); i++){ 227 AF = AF + out[i]; 291 228 } 292 229 -
Singular/LIB/grobcov.lib
r31b00db r6e3023a 113 113 Can be used without calling presiouly setglobalrings(), 114 114 115 locus(G) ;Special routine for determining the locus of points115 locus(G): Special routine for determining the locus of points 116 116 of objects. Given a parametric ideal J with 117 117 parameters (a_1,..a_m) and variables (x_1,..,xn), … … 131 131 Can be used without calling presiouly setglobalrings(), 132 132 133 locusdg(G) ;Is a special routine for computing the locus in dinamical geometry.133 locusdg(G): Is a special routine for computing the locus in dinamical geometry. 134 134 It detects automatically a possible point that is to be avoided by the mover, 135 135 whose coordinates must be the last two coordinates in the definition of the ring. … … 139 139 Can be used without calling presiouly setglobalrings(), 140 140 141 locusto(L) ;Transforms the output of locus to a string that141 locusto(L): Transforms the output of locus to a string that 142 142 can be reed from different computational systems. 143 143 Can be used without calling presiouly setglobalrings(), 144 144 145 addcons(L) ;Given a disjoint set of locally closed subsets in P-representation,145 addcons(L): Given a disjoint set of locally closed subsets in P-representation, 146 146 it returns the canonical P-representation of the constructible set 147 147 formed by the union of them. … … 1222 1222 proc PtoCrep(list L) 1223 1223 { 1224 int te=0; 1224 1225 def RR=basering; 1226 if(defined(@P)){te=1;} 1227 else{te=0; setglobalrings();} 1225 1228 setring(@P); 1226 1229 def Lp=imap(RR,L); … … 1240 1243 def La=list(ida,idb); 1241 1244 setring(RR); 1242 return(imap(@P,La)); 1245 def LL=imap(@P,La); 1246 if(te==0){kill @P; kill @R; kill @RP;} 1247 return(LL); 1243 1248 } 1244 1249 … … 3895 3900 int i; int j; int k; int l; 3896 3901 intvec v; intvec v1; intvec v0; 3897 ideal J; list K; 3902 ideal J; list K; list L1; 3903 for(i=1;i<=size(L);i++) 3904 { 3905 if(equalideals(L[i][1][1],ideal(1))==0) 3906 {L1[size(L1)+1]=L[i];} 3907 } 3908 L=L1; 3898 3909 for(i=1;i<=size(L);i++) 3899 3910 { -
Singular/LIB/modstd.lib
r31b00db r6e3023a 33 33 static proc mod_init() 34 34 { 35 newstruct("ideal _primeTest", "ideal Ideal");35 newstruct("idealPrimeTest", "ideal Ideal"); 36 36 } 37 37 … … 487 487 list arguments; 488 488 int neededPrimes = neededSize-size(L); 489 ideal _primeTest Id;489 idealPrimeTest Id; 490 490 Id.Ideal = I; 491 491 export(Id); -
Singular/LIB/normal.lib
r31b00db r6e3023a 4445 4445 printlevel = printlevel + 1; 4446 4446 ideal JDefault = 0; // Now it uses the default J; 4447 list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault) [1];4448 list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault) [1];4447 list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault); 4448 list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault); 4449 4449 printlevel = printlevel - 1; 4450 4450 option(set,save_opt); 4451 return(list(nor1, nor2)); 4451 list res = nor1 + nor2; 4452 return(res); 4452 4453 } 4453 4454 } … … 4568 4569 4569 4570 ideal JDefault = 0; // Now it uses the default J; 4570 list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault) [1];4571 list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault) [1];4571 list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault); 4572 list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault); 4572 4573 printlevel = printlevel - 1; 4573 4574 option(set,save_opt); 4574 return(list(nor1, nor2)); 4575 list res = nor1 + nor2; 4576 return(res); 4575 4577 } 4576 4578 } -
Singular/LIB/polymake.lib
r31b00db r6e3023a 1 //// 2 version="version oldpolymake.lib 4.0.0.0 Jun_2013 "; // $Id$ 1 version="version polymake.lib 4.0.0.0 Jun_2013 "; 3 2 category="Tropical Geometry"; 4 3 info=" 5 LIBRARY: oldpolymake.libComputations with polytopes and fans,6 4 LIBRARY: polymake.lib Computations with polytopes and fans, 5 interface to polymake and TOPCOM 7 6 AUTHOR: Thomas Markwig, email: keilen@mathematik.uni-kl.de 8 7 … … 10 9 Most procedures will not work unless polymake or topcom is installed and 11 10 if so, they will only work with the operating system LINUX! 12 For more detailed information see the following note orconsult the11 For more detailed information see IMPORTANT NOTE respectively consult the 13 12 help string of the procedures. 14 13 15 NOTE: 14 The conventions used in this library for polytopes and fans, e.g. the 15 length and labeling of their vertices resp. rays, differs from the conventions 16 used in polymake and thus from the conventions used in the polymake 17 extension polymake.so of Singular. We recommend to use the newer polymake.so 18 whenever possible. 19 20 IMPORTANT NOTE: 16 21 Even though this is a Singular library for computing polytopes and fans 17 22 such as the Newton polytope or the Groebner fan of a polynomial, most of … … 19 24 @* - polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, TU Darmstadt 20 25 @* (see http://www.math.tu-berlin.de/polymake/), 21 @* respectively (only in the procedure triangula rions) by the program22 @* - topcom by Joerg Rambau, Universitaet Bayreuth (see http://www.uni-bayreuth.de/23 24 @* This library should rather be seen as an interface which allows to use a26 @* respectively (only in the procedure triangulations) by the program 27 @* - topcom by Joerg Rambau, Universitaet Bayreuth (see 28 @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM); 29 @* this library should rather be seen as an interface which allows to use a 25 30 (very limited) number of options which polymake respectively topcom offers 26 31 to compute with polytopes and fans and to make the results available in … … 32 37 independent of both, polymake and topcom. 33 38 34 PROCEDURES: 35 polymakePolytope() computes the vertices of a polytope using polymake 36 newtonPolytopeP() computes the Newton polytope of a polynomial 37 newtonPolytopeLP() computes the lattice points of the Newton polytope 38 normalFan() computes the normal fan of a polytope 39 groebnerFan() computes the Groebner fan of a polynomial 40 intmatToPolymake() transforms an integer matrix into polymake format 41 polymakeToIntmat() transforms polymake output into an integer matrix 42 43 triangulations() computes all triangulations of a marked polytope 44 secondaryPolytope() computes the secondary polytope of a marked polytope 45 46 secondaryFan() computes the secondary fan of a marked polytope 47 48 cycleLength() computes the cycleLength of cycle 49 splitPolygon() splits a marked polygon into vertices, facets, interior points 50 eta() computes the eta-vector of a triangulation 51 findOrientedBoundary() computes the boundary of a convex hull 52 cyclePoints() computes lattice points connected to some lattice point 53 latticeArea() computes the lattice area of a polygon 54 picksFormula() computes the ingrediants of Pick's formula for a polygon 55 ellipticNF() computes the normal form of an elliptic polygon 56 ellipticNFDB() displays the 16 normal forms of elliptic polygons 57 58 polymakeKeepTmpFiles() determines if the files created in /tmp should be kept 39 PROCEDURES USING POLYMAKE: 40 polymakePolytope() computes the vertices of a polytope using polymake 41 newtonPolytopeP() computes the Newton polytope of a polynomial 42 newtonPolytopeLP() computes the lattice points of the Newton polytope 43 normalFanL() computes the normal fan of a polytope 44 groebnerFan() computes the Groebner fan of a polynomial 45 46 PROCEDURES USING TOPCOM: 47 triangulations() computes all triangulations of a marked polytope 48 secondaryPolytope() computes the secondary polytope of a marked polytope 49 50 PROCEDURES USING POLYMAKE AND TOPCOM: 51 secondaryFan() computes the secondary fan of a marked polytope 52 53 PROCEDURES CONERNED WITH PLANAR POLYGONS: 54 cycleLength() computes the cycleLength of cycle 55 splitPolygon() splits a marked polygon into vertices, facets, interior points 56 eta() computes the eta-vector of a triangulation 57 findOrientedBoundary() computes the boundary of a convex hull 58 cyclePoints() computes lattice points connected to some lattice point 59 latticeArea() computes the lattice area of a polygon 60 picksFormula() computes the ingrediants of Pick's formula for a polygon 61 ellipticNF() computes the normal form of an elliptic polygon 62 ellipticNFDB() displays the 16 normal forms of elliptic polygons 59 63 60 64 KEYWORDS: polytope; fan; secondary fan; secondary polytope; polymake; … … 85 89 //////////////////////////////////////////////////////////////////////////////// 86 90 91 static proc mod_init () 92 { 93 LIB "gfanlib.so"; 94 LIB "polymake.so"; 95 } 87 96 88 97 ///////////////////////////////////////////////////////////////////////////// … … 90 99 ///////////////////////////////////////////////////////////////////////////// 91 100 92 proc polymakePolytope (intmat po lytop,list #)93 "USAGE: polymakePolytope(po lytope[,#]); polytope list, # string94 ASSUME: each row of po lytopegives the coordinates of a lattice point of a101 proc polymakePolytope (intmat points) 102 "USAGE: polymakePolytope(points); polytope intmat 103 ASSUME: each row of points gives the coordinates of a lattice point of a 95 104 polytope with their affine coordinates as given by the output of 96 105 secondaryPolytope 97 106 PURPOSE: the procedure calls polymake to compute the vertices of the polytope 98 107 as well as its dimension and information on its facets 99 RETURN: list L with four entries108 RETURN: list, L with four entries 100 109 @* L[1] : an integer matrix whose rows are the coordinates of vertices 101 110 of the polytope 102 111 @* L[2] : the dimension of the polytope 103 @* L[3] : a list whose i -th entry explains to which vertices the112 @* L[3] : a list whose ith entry explains to which vertices the 104 113 ith vertex of the Newton polytope is connected 105 114 -- i.e. L[3][i] is an integer vector and an entry k in 106 115 there means that the vertex L[1][i] is connected to the 107 116 vertex L[1][k] 108 @* L[4] : an integer matrixwhose rows mulitplied by117 @* L[4] : an matrix of type bigintmat whose rows mulitplied by 109 118 (1,var(1),...,var(nvar)) give a linear system of equations 110 119 describing the affine hull of the polytope, … … 118 127 convention which starts indexing its vertices by zero while we start 119 128 with one ! 120 @* - the procedure creates the file /tmp/polytope.polymake which contains121 the polytope in polymake format; if you wish to use this for further122 computations with polymake, you have to use the procedure123 polymakeKeepTmpFiles in before124 @* - moreover, the procedure creates the file /tmp/polytope.output which125 it deletes again before ending126 @* - it is possible to provide an optional second argument a string127 which then will be used instead of 'polytope' in the name of the128 polymake output file129 129 EXAMPLE: example polymakePolytope; shows an example" 130 130 { 131 // the header for the file secendarypolytope.polymake 132 string sp="_application polytope 133 _version 2.2 134 _type RationalPolytope 135 136 POINTS 137 "; 131 // add a first column to polytope as homogenising coordinate 132 points=intmatAddFirstColumn(points,"points"); 133 polytope polytop=polytopeViaPoints(points); 134 list graph=vertexAdjacencyGraph(polytop)[2]; 138 135 int i,j; 139 // set the name for the polymake output file 140 if (size(#)>0) 141 { 142 if (typeof(#[1])=="string") 143 { 144 string dateiname=#[1]; 145 } 146 else 147 { 148 string dateiname="polytope"; 149 } 150 } 151 else 152 { 153 string dateiname="polytope"; 154 } 155 // create the lattice point list for polymake 156 sp=sp+intmatToPolymake(polytop,"points"); 157 // initialise dateiname.polymake and compute the vertices 158 write(":w /tmp/"+dateiname+".polymake",sp); 159 system("sh","cd /tmp; polymake "+dateiname+".polymake VERTICES > "+dateiname+".output"); 160 string vertices=read("/tmp/"+dateiname+".output"); 161 system("sh","/bin/rm /tmp/"+dateiname+".output"); 162 intmat np=polymakeToIntmat(vertices,"affine"); 163 // compute the dimension 164 system("sh","cd /tmp; polymake "+dateiname+".polymake DIM > "+dateiname+".output"); 165 string pdim=read("/tmp/"+dateiname+".output"); 166 system("sh","/bin/rm /tmp/"+dateiname+".output"); 167 pdim=pdim[5,size(pdim)-6]; 168 execute("int nd="+pdim+";"); 169 // compute the vertex-edge graph 170 system("sh","cd /tmp; polymake "+dateiname+".polymake GRAPH > "+dateiname+".output"); 171 string vertexedgegraph=read("/tmp/"+dateiname+".output"); 172 system("sh","/bin/rm /tmp/"+dateiname+".output"); 173 vertexedgegraph=vertexedgegraph[7,size(vertexedgegraph)-8]; 174 string newveg; 175 for (i=1;i<=size(vertexedgegraph);i++) 176 { 177 if (vertexedgegraph[i]=="{") 178 { 179 newveg=newveg+"intvec("; 180 } 181 else 182 { 183 if (vertexedgegraph[i]=="}") 184 { 185 newveg=newveg+"),"; 186 } 187 else 188 { 189 if (vertexedgegraph[i]==" ") 190 { 191 newveg=newveg+","; 192 } 193 else 194 { 195 newveg=newveg+vertexedgegraph[i]; 196 } 197 } 198 } 199 } 200 newveg=newveg[1,size(newveg)-1]; 201 execute("list nveg="+newveg+";"); 202 // raise each entry in nveg by one 203 for (i=1;i<=size(nveg);i++) 204 { 205 for (j=1;j<=size(nveg[i]);j++) 206 { 207 nveg[i][j]=nveg[i][j]+1; 208 } 209 } 210 // compute the affine hull 211 system("sh","cd /tmp; polymake "+dateiname+".polymake AFFINE_HULL > "+dateiname+".output"); 212 string equations=read("/tmp/"+dateiname+".output"); 213 system("sh","/bin/rm /tmp/"+dateiname+".output"); 214 if (size(equations)>14) 215 { 216 intmat neq=polymakeToIntmat(equations,"cleardenom"); 217 } 218 else 219 { 220 intmat neq[1][ncols(polytop)+1]; 221 } 222 // delete the tmp-files, if polymakekeeptmpfiles is not set 223 if (defined(polymakekeeptmpfiles)==0) 224 { 225 system("sh","/bin/rm /tmp/"+dateiname+".polymake"); 226 } 227 // return the files 228 return(list(np,nd,nveg,neq)); 136 for (i=1;i<=size(graph);i++) 137 { 138 for (j=1;j<=size(graph[i]);j++) 139 { 140 graph[i][j]=graph[i][j]+1; 141 } 142 } 143 return(list(intmatcoldelete(vertices(polytop),1),dimension(polytop),graph,equations(polytop))); 229 144 } 230 145 example … … 247 162 ring r=0,x(1..4),dp; 248 163 matrix M[5][1]=1,x(1),x(2),x(3),x(4); 249 np[4]*M;164 intmat(np[4])*M; 250 165 } 251 166 252 167 ///////////////////////////////////////////////////////////////////////////// 253 168 254 proc newtonPolytopeP (poly f ,list #)255 "USAGE: newtonPolytopeP(f [,#]); f poly, # string256 RETURN: list L with four entries169 proc newtonPolytopeP (poly f) 170 "USAGE: newtonPolytopeP(f); f poly 171 RETURN: list, L with four entries 257 172 @* L[1] : an integer matrix whose rows are the coordinates of vertices 258 173 of the Newton polytope of f … … 263 178 there means that the vertex L[1][i] is 264 179 connected to the vertex L[1][k] 265 @* L[4] : an integer matrixwhose rows mulitplied by180 @* L[4] : an matrix of type bigintmat whose rows mulitplied by 266 181 (1,var(1),...,var(nvar)) give a linear system of equations 267 182 describing the affine hull of the Newton polytope, i.e. the … … 269 184 NOTE: - if we replace the first column of L[4] by zeros, i.e. if we move 270 185 the affine hull to the origin, then we get the equations for the 271 orthogonal compl oment of the linearity space of the normal fan dual186 orthogonal complement of the linearity space of the normal fan dual 272 187 to the Newton polytope, i.e. we get the EQUATIONS that 273 188 we need as input for polymake when computing the normal fan … … 275 190 TU Berlin and Michael Joswig, so it only works if polymake is installed; 276 191 see http://www.math.tu-berlin.de/polymake/ 277 @* - the procedure creates the file /tmp/newtonPolytope.polymake which278 contains the polytope in polymake format and which can be used for279 further computations with polymake280 @* - moreover, the procedure creates the file /tmp/newtonPolytope.output281 and deletes it again before ending282 @* - it is possible to give as an optional second argument a string283 which then will be used instead of 'newtonPolytope' in the name of284 the polymake output file285 192 EXAMPLE: example newtonPolytopeP; shows an example" 286 193 { … … 296 203 f=f-lead(f); 297 204 } 298 if (size(#)==0)299 {300 #[1]="newtonPolytope";301 }302 205 // call polymakePolytope with exponents 303 return(polymakePolytope(exponents ,#));206 return(polymakePolytope(exponents)); 304 207 } 305 208 example … … 330 233 // the Newton polytope is contained in the affine space given 331 234 // by the equations 332 np[4]*M;235 intmat(np[4])*M; 333 236 } 334 237 … … 339 242 RETURN: list, the exponent vectors of the monomials occuring in f, 340 243 i.e. the lattice points of the Newton polytope of f 341 EXAMPLE: example n ormalFan; shows an example"244 EXAMPLE: example newtonPolytopeLP; shows an example" 342 245 { 343 246 list np; … … 363 266 ///////////////////////////////////////////////////////////////////////////// 364 267 365 proc normalFan (intmat vertices,intmataffinehull,list graph,int er,list #)366 "USAGE: normalFan (vert,aff,graph,rays,[,#]); vert,aff intmat, graph list, rays int, # string268 proc normalFanL (def vertices, def affinehull,list graph,int er,list #) 269 "USAGE: normalFanL (vert,aff,graph,rays,[,#]); vert,aff intmat, graph list, rays int, # string 367 270 ASSUME: - vert is an integer matrix whose rows are the coordinate of 368 271 the vertices of a convex lattice polytope; … … 396 299 @* - in the optional argument # it is possible to hand over other names 397 300 for the variables to be used -- be careful, the format must be correct 398 which is not tested, e.g. if you want the variable names to be 399 u00,u10,u01,u11 then you must hand over the string 'u11,u10,u01,u11' 400 EXAMPLE: example normalFan; shows an example" 401 { 301 and that is not tested, e.g. if you want the variable names to be 302 u00,u10,u01,u11 then you must hand over the string u11,u10,u01,u11 303 EXAMPLE: example normalFanL; shows an example" 304 { 305 if (typeof(affinehull) != "intmat" && typeof (affinehull) != "bigintmat") 306 { 307 ERROR("normalFanL: input affinehull has to be either intmat or bigintmat"); 308 list L; 309 return (L); 310 } 402 311 list ineq; // stores the inequalities of the cones 403 312 int i,j,k; … … 431 340 list pp; // contain strings representing the inequalities 432 341 // describing the normal cone 433 intmat ie[size(graph[i])][ncols(vertices)]; // contains the inequalities 434 // as rows 342 if (typeof (vertices) == "intmat") 343 { 344 intmat ie[size(graph[i])][ncols(vertices)]; // contains the inequalities 345 } // as rows 346 if (typeof (vertices) == "bigintmat") 347 { 348 bigintmat ie[size(graph[i])][ncols(vertices)]; // contains the inequalities 349 } // as rows 435 350 // consider all the vertices to which the ith vertex in the 436 351 // polytope is connected by an edge … … 469 384 } 470 385 // remove the first column of affine hull to compute the linearity space 471 intmat linearity=intmatcoldelete(affinehull,1); 386 bigintmat linearity[1][ncols(vertices)]; 387 if (nrows(affinehull)>0) 388 { 389 linearity=intmatcoldelete(affinehull,1); 390 } 472 391 ////////////////////////////////////////////////////////////////// 473 392 // Compute next the extreme rays of the cones … … 476 395 { 477 396 list extremerays; // keeps the result 478 string polymake; // keeps polymake output 479 // the header for ineq.polymake 480 string head="_application polytope 481 _version 2.2 482 _type RationalPolytope 483 484 INEQUALITIES 485 "; 486 // the tail for both polymake files 487 string tail=" 488 EQUATIONS 489 "; 490 tail=tail+intmatToPolymake(linearity,"rays"); 491 string ungleichungen; // keeps the inequalities for the polymake code 397 cone kegel; 398 bigintmat linearspan=intmatAddFirstColumn(linearity,"rays"); 492 399 intmat M; // the matrix keeping the inequalities 493 // create the file ineq.output494 write(":w /tmp/ineq.output","");495 int dimension; // keeps the dimension of the intersection the496 // bad cones with the u11tobeseencone497 400 for (i=1;i<=size(ineq);i++) 498 401 { 499 i,". Cone of ",nrows(vertices); // indicate how many 500 // vertices have been dealt with 501 ungleichungen=intmatToPolymake(ineq[i][1],"rays"); 502 // write the inequalities to ineq.polymake and call polymake 503 write(":w /tmp/ineq.polymake",head+ungleichungen+tail); 504 ungleichungen=""; // clear ungleichungen 505 system("sh","cd /tmp; /bin/rm ineq.output; polymake ineq.polymake VERTICES > ineq.output"); 506 // read the result of polymake 507 polymake=read("/tmp/ineq.output"); 508 intmat VERT=polymakeToIntmat(polymake,"affine"); 509 extremerays[i]=VERT; 510 kill VERT; 402 kegel=coneViaInequalities(intmatAddFirstColumn(ineq[i][1],"rays"),linearspan); 403 extremerays[i]=intmatcoldelete(rays(kegel),1); 511 404 } 512 405 for (i=1;i<=size(ineq);i++) … … 514 407 ineq[i]=ineq[i]+list(extremerays[i]); 515 408 } 516 }517 // delete the tmp-files, if polymakekeeptmpfiles is not set518 if (defined(polymakekeeptmpfiles)==0)519 {520 system("sh","/bin/rm /tmp/ineq.polymake");521 system("sh","/bin/rm /tmp/ineq.output");522 409 } 523 410 // get the linearity space … … 534 421 list np=newtonPolytopeP(f); 535 422 // the Groebner fan of f, i.e. the normal fan of the Newton polytope 536 list gf=normalFan (np[1],np[4],np[3],1,"x,y,z");423 list gf=normalFanL(np[1],np[4],np[3],1,"x,y,z"); 537 424 // the number of cones in the Groebner fan of f is: 538 425 size(gf[1]); … … 549 436 ///////////////////////////////////////////////////////////////////////////// 550 437 551 proc groebnerFan (poly f ,list #)552 "USAGE: groebnerFan(f [,#]); f poly, # string438 proc groebnerFan (poly f) 439 "USAGE: groebnerFan(f); f poly 553 440 RETURN: list, the ith entry of L[1] contains information about the ith cone 554 441 in the Groebner fan dual to the ith vertex in the Newton … … 564 451 in each cone 565 452 @* L[3] = the Newton polytope of f in the format of the procedure 566 newtonPolytope 567 @* L[4] = integer matrix where each row represents the expone t453 newtonPolytopeP 454 @* L[4] = integer matrix where each row represents the exponent 568 455 vector of one monomial occuring in the input polynomial 569 456 NOTE: - if you have already computed the Newton polytope of f then you might want 570 to use the procedure normalFan instead in order to avoid doing costly457 to use the procedure normalFanL instead in order to avoid doing costly 571 458 computation twice 572 459 @* - the procedure calls for its computation polymake by Ewgenij Gawrilow, 573 460 TU Berlin and Michael Joswig, so it only works if polymake is installed; 574 461 see http://www.math.tu-berlin.de/polymake/ 575 @* - the procedure creates the file /tmp/newtonPolytope.polymake which576 contains the Newton polytope of f in polymake format and which can577 be used for further computations with polymake578 @* - it is possible to give as an optional second argument as string which579 then will be used instead of 'newtonPolytope' in the name of the580 polymake output file581 462 EXAMPLE: example groebnerFan; shows an example" 582 463 { … … 591 472 f=f-lead(f); 592 473 } 593 if (size(#)==0)594 {595 #[1]="newtonPolytope";596 }597 474 // call polymakePolytope with exponents 598 list newtonp=polymakePolytope(exponents ,"newtonPolytope");475 list newtonp=polymakePolytope(exponents); 599 476 // get the variables as string 600 477 string variablen; … … 604 481 } 605 482 variablen=variablen[1,size(variablen)-1]; 606 // call normalFan in order to compute the Groebner fan607 list gf=normalFan (newtonp[1],newtonp[4],newtonp[3],1,variablen);483 // call normalFanL in order to compute the Groebner fan 484 list gf=normalFanL(newtonp[1],newtonp[4],newtonp[3],1,variablen); 608 485 // append newtonp to gf 609 486 gf[3]=newtonp; … … 642 519 643 520 644 /////////////////////////////////////////////////////////////////////////////645 646 proc intmatToPolymake (intmat M,string art)647 "USAGE: intmatToPolymake(M,art); M intmat, art string648 ASSUME: - M is an integer matrix which should be transformed into polymake649 format;650 @* - art is one of the following strings:651 @* + 'rays' : indicating that a first column of 0's should be added652 @* + 'points' : indicating that a first column of 1's should be added653 RETURN: string, the matrix is transformed in a string and a first column has654 been added655 EXAMPLE: example intmatToPolymake; shows an example"656 {657 if (art=="rays")658 {659 string anf="0 ";660 }661 else662 {663 string anf="1 ";664 }665 string sp;666 int i,j;667 // create the lattice point list for polymake668 for (i=1;i<=nrows(M);i++)669 {670 sp=sp+anf;671 for (j=1;j<=ncols(M);j++)672 {673 sp=sp+string(M[i,j])+" ";674 if (j==ncols(M))675 {676 sp=sp+"677 ";678 }679 }680 }681 return(sp);682 }683 example684 {685 "EXAMPLE:";686 echo=2;687 intmat M[3][4]=1,2,3,4,5,6,7,8,9,10,11,12;688 intmatToPolymake(M,"rays");689 intmatToPolymake(M,"points");690 }691 692 /////////////////////////////////////////////////////////////////////////////693 694 proc polymakeToIntmat (string pm,string art)695 "USAGE: polymakeToIntmat(pm,art); pm, art string696 ASSUME: pm is the result of calling polymake with one 'argument' like697 VERTICES, AFFINE_HULL, etc., so that the first row of the string is698 the name of the corresponding 'argument' and the further rows contain699 the result which consists of vectors either over the integers700 or over the rationals701 RETURN: intmat, the rows of the matrix are basically the vectors in pm, starting702 from the second row, where each row has been multiplied with the703 lowest common multiple of the denominators of its entries as if704 it is an integer matrix; moreover, if art=='affine', then705 the first column is omitted since we only want affine706 coordinates707 EXAMPLE: example polymakeToIntmat; shows an example"708 {709 // we need a line break710 string zeilenumbruch="711 ";712 // remove the 'argment' name, i.e. the first row of pm713 while (pm[1]!=zeilenumbruch)714 {715 pm=stringdelete(pm,1);716 }717 pm=stringdelete(pm,1);718 // find out how many entries each vector has, namely one more719 // than 'spaces' in a row720 int i=1;721 int s=1;722 int z=1;723 while (pm[i]!=zeilenumbruch)724 {725 if (pm[i]==" ")726 {727 s++;728 }729 i++;730 }731 // if we want to have affine coordinates732 if (art=="affine")733 {734 s--; // then there is one column less735 // and the entry of the first column (in the first row) has to be removed736 while (pm[1]!=" ")737 {738 pm=stringdelete(pm,1);739 }740 pm=stringdelete(pm,1);741 }742 // we add two line breaks at the end in order to have this as743 // a stopping criterion744 pm=pm+zeilenumbruch+zeilenumbruch;745 // we now have to work through each row746 for (i=1;i<=size(pm);i++)747 {748 // if there are two consecutive line breaks we are done749 if ((pm[i]==zeilenumbruch) and (pm[i+1]==zeilenumbruch))750 {751 i=size(pm)+1;752 }753 else754 {755 // a line break has to be replaced by a comma756 if (pm[i]==zeilenumbruch)757 {758 z++;759 pm[i]=",";760 // if we want to have affine coordinates,761 // then we have to delete the first entry in each row762 if (art=="affine")763 {764 while (pm[i+1]!=" ")765 {766 pm=stringdelete(pm,i+1);767 }768 pm=stringdelete(pm,i+1);769 }770 }771 // a space has to be replaced by a comma772 if (pm[i]==" ")773 {774 pm[i]=",";775 }776 }777 }778 // if we have introduced superflous commata at the end, they should be removed779 while (pm[size(pm)]==",")780 {781 pm=stringdelete(pm,size(pm));782 }783 // since the matrix could be over the rationals,784 // we need a ring with rational coefficients785 ring zwischering=0,x,lp;786 // create the matrix with the elements of pm as entries787 execute("matrix mm["+string(z)+"]["+string(s)+"]="+pm+";");788 // transform this into an integer matrix789 matrix M[1][ncols(mm)]; // takes a row of mm790 int cm; // takes a lowest common multiple791 // multiply each row by an integer such that its entries are integers792 for (int j=1;j<=nrows(mm);j++)793 {794 M=mm[j,1..ncols(mm)];795 cm=commondenominator(M);796 for (i=1;i<=ncols(mm);i++)797 {798 mm[j,i]=cm*mm[j,i];799 }800 }801 // transform the matrix mm into an integer matrix802 execute("intmat im["+string(z)+"]["+string(s)+"]="+string(mm)+";");803 return(im);804 }805 example806 {807 "EXAMPLE:";808 echo=2;809 // this is the usual output of some polymake computation810 string pm="VERTICES811 0 1 3 5/3 1/3 -1 -23/3 -1/3 5/3 1/3 1812 0 1 3 -23/3 5/3 1 5/3 1/3 1/3 -1/3 -1813 0 1 1 1/3 -1/3 -1 5/3 1/3 -23/3 5/3 3814 0 1 1 5/3 -23/3 3 1/3 5/3 -1/3 1/3 -1815 0 1 -1 1/3 5/3 3 -1/3 -23/3 1/3 5/3 1816 0 1 -1 -1/3 1/3 1 1/3 5/3 5/3 -23/3 3817 0 1 -1 1 3 -5 -1 3 -1 1 -1818 0 1 -1 -1 -1 -1 1 1 3 3 -5819 0 1 -5 3 1 -1 3 -1 1 -1 -1820 821 ";822 intmat PM=polymakeToIntmat(pm,"affine");823 // note that the first column has been removed, since we asked for824 // affine coordinates, and the denominators have been cleared825 print(PM);826 }827 828 521 /////////////////////////////////////////////////////////////////////////////// 829 522 /// PROCEDURES USING TOPCOM 830 523 /////////////////////////////////////////////////////////////////////////////// 831 524 832 proc triangulations (list polygon )833 "USAGE: triangulations(polygon ); list polygon525 proc triangulations (list polygon,list #) 526 "USAGE: triangulations(polygon[,#]); list polygon, list # 834 527 ASSUME: polygon is a list of integer vectors of the same size representing 835 528 the affine coordinates of the lattice points … … 846 539 this procedure; see 847 540 @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM 541 @* - if you only want to have the regular triangulations the procedure should 542 be called with the string 'regular' as optional argument 848 543 @* - the procedure creates the files /tmp/triangulationsinput and 849 544 /tmp/triangulationsoutput; … … 851 546 output containing the triangulations of corresponding to points in the 852 547 format of points2triangs; if you wish to use this for further 853 computations with topcom, you have to use the procedure854 polymakeKeepTmpFiles in before548 computations with topcom, you have to call the procedure with the 549 string 'keepfiles' as optional argument 855 550 @* - note that an integer i in an integer vector representing a triangle 856 551 refers to the ith lattice point, i.e. polygon[i]; this convention is … … 860 555 { 861 556 int i,j; 557 // check for optional arguments 558 int regular,keepfiles; 559 if (size(#)>0) 560 { 561 for (i=1;i<=size(#);i++) 562 { 563 if (typeof(#[i])=="string") 564 { 565 if (#[i]=="keepfiles") 566 { 567 keepfiles=1; 568 } 569 if (#[i]=="regular") 570 { 571 regular=1; 572 } 573 } 574 } 575 } 862 576 // prepare the input for points2triangs by writing the input polygon in the 863 577 // necessary format … … 875 589 write(":w /tmp/triangulationsinput",spi); 876 590 // call points2triangs 877 system("sh","cd /tmp; points2triangs < triangulationsinput > triangulationsoutput"); 591 if (regular==1) // compute only regular triangulations 592 { 593 system("sh","cd /tmp; points2triangs --regular < triangulationsinput > triangulationsoutput"); 594 } 595 else // compute all triangulations 596 { 597 system("sh","cd /tmp; points2triangs < triangulationsinput > triangulationsoutput"); 598 } 878 599 string p2t=read("/tmp/triangulationsoutput"); // takes result of points2triangs 879 // delete the tmp-files, if polymakekeeptmpfiles is not set880 if ( defined(polymakekeeptmpfiles)==0)600 // delete the tmp-files, if no second argument is given 601 if (keepfiles==0) 881 602 { 882 603 system("sh","cd /tmp; rm -f triangulationsinput; rm -f triangulationsoutput"); … … 953 674 else 954 675 { 955 np2t=np2t+p2t[i]; 676 if (p2t[i]=="[") 677 { 678 // in Topcom version 17.4 (and maybe also in earlier versions) the list 679 // of triangulations is indexed starting with index 0, in Singular 680 // we have to start with index 1 681 np2t=np2t+p2t[i]+"1+"; 682 } 683 else 684 { 685 np2t=np2t+p2t[i]; 686 } 956 687 } 957 688 } … … 962 693 list T; 963 694 execute(np2t); 695 // depending on the version of Topcom, the list T has or has not an entry T[1] 696 // if it has none, the entry should be removed 697 while (typeof(T[1])=="none") 698 { 699 T=delete(T,1); 700 } 964 701 // raise each index by one 965 702 for (i=1;i<=size(T);i++) … … 1000 737 RETURN: list, say L, such that: 1001 738 @* L[1] = intmat, each row gives the affine coordinates of a lattice 1002 point in the secondary polytope given by the marked polytope1003 corresponding to polygon739 point in the secondary polytope given by the marked 740 polytope corresponding to polygon 1004 741 @* L[2] = the list of corresponding triangulations 1005 742 NOTE: if the triangluations are not handed over as optional argument the … … 1114 851 see http://www.math.tu-berlin.de/polymake/ 1115 852 @* - in the optional argument # it is possible to hand over other names for 1116 the variables to be used -- be careful, the format must be correct 1117 whichis not tested, e.g. if you want the variable names to be853 the variables to be used -- be careful, the format must be correct and 854 that is not tested, e.g. if you want the variable names to be 1118 855 u00,u10,u01,u11 then you must hand over the string 'u11,u10,u01,u11' 1119 856 @* - if the triangluations are not handed over as optional argument the … … 1135 872 list sp=secondaryPolytope(polygon,triang); 1136 873 list spp=polymakePolytope(sp[1]); 1137 list sf=normalFan (spp[1],spp[4],spp[3],1);874 list sf=normalFanL(spp[1],spp[4],spp[3],1); 1138 875 return(list(sf[1],sf[2],spp,triang)); 1139 876 } … … 1378 1115 @* triang is a list of integer vectors all of size three describing a 1379 1116 triangulation of the polygon described by polygon; if an entry of 1380 triang is the vector (i,j,k) then the triangle is built fromthe vertices1117 triang is the vector (i,j,k) then the triangle is built by the vertices 1381 1118 with indices i, j and k 1382 1119 RETURN: intvec, the integer vector eta describing that vertex of the Newton … … 1948 1685 // points and 1949 1686 // the number b of boundary points are connected by the formula: A=b+2g-2 1950 return(list(area,bdpts,(area-bdpts+2) /2));1687 return(list(area,bdpts,(area-bdpts+2) div 2)); 1951 1688 } 1952 1689 example … … 1973 1710 "USAGE: ellipticNF(polygon); polygon list 1974 1711 ASSUME: polygon is a list of integer vectors in the plane such that their 1975 convex hull C has precisely one interior lattice point ,i.e. C is the1712 convex hull C has precisely one interior lattice point; i.e. C is the 1976 1713 Newton polygon of an elliptic curve 1977 1714 PURPOSE: compute the normal form of the polygon with respect to the unimodular … … 1991 1728 int i; // index 1992 1729 intvec edge; // stores the vector of an edge 1993 intvec boundary; // stores lattice lengths of the edges of the Newton cy lce1730 intvec boundary; // stores lattice lengths of the edges of the Newton cycle 1994 1731 // find the vertices of the Newton cycle and order it clockwise 1995 1732 list pg=findOrientedBoundary(polygon)[2]; … … 2302 2039 2303 2040 ///////////////////////////////////////////////////////////////////////////////// 2304 /// AUXILARY PROCEDURES2305 /////////////////////////////////////////////////////////////////////////////////2306 2307 proc polymakeKeepTmpFiles (int i)2308 "USAGE: polymakeKeepTmpFiles(int i); i int2309 PURPOSE: some procedures create files in the directory /tmp which are used for2310 computations with polymake respectively topcom; these will be removed2311 when the corresponding procedure is left; however, it might be2312 desireable to keep them for further computations with either polymake or2313 topcom; this can be achieved by this procedure; call the procedure as:2314 @* - polymakeKeepTmpFiles(1); - then the files will be kept2315 @* - polymakeKeepTmpFiles(0); - then files will be removed in the future2316 RETURN: none"2317 {2318 if ((i==1) and (defined(polymakekeeptmpfiles)==0))2319 {2320 int polymakekeeptmpfiles;2321 export(polymakekeeptmpfiles);2322 }2323 if (i!=1)2324 {2325 if (defined(polymakekeeptmpfiles))2326 {2327 kill polymakekeeptmpfiles;2328 }2329 }2330 }2331 2332 2333 /////////////////////////////////////////////////////////////////////////////////2334 2041 ///////////////////////////////////////////////////////////////////////////////// 2335 2042 /// AUXILARY PROCEDURES, WHICH ARE DECLARED STATIC … … 2364 2071 } 2365 2072 2366 static proc intmatcoldelete ( intmatw,int i)2073 static proc intmatcoldelete (def w,int i) 2367 2074 "USAGE: intmatcoldelete(w,i); w intmat, i int 2368 2075 RETURN: intmat, the integer matrix w with the ith comlumn deleted 2369 NOTE: the procedure is called by intmatsort and normalFan" 2370 { 2371 if ((i<1) or (i>ncols(w)) or (ncols(w)==1)) 2372 { 2373 return(w); 2374 } 2375 if (i==1) 2376 { 2377 intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),2..ncols(w)]; 2378 return(M); 2379 } 2380 if (i==ncols(w)) 2381 { 2382 intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),1..ncols(w)-1]; 2383 return(M); 2384 } 2385 else 2386 { 2387 intmat M[nrows(w)][i-1]=w[1..nrows(w),1..i-1]; 2388 intmat N[nrows(w)][ncols(w)-i]=w[1..nrows(w),i+1..ncols(w)]; 2389 return(intmatconcat(M,N)); 2076 NOTE: the procedure is called by intmatsort and normalFanL" 2077 { 2078 if (typeof(w)=="intmat") 2079 { 2080 if ((i<1) or (i>ncols(w)) or (ncols(w)==1)) 2081 { 2082 return(w); 2083 } 2084 if (i==1) 2085 { 2086 intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),2..ncols(w)]; 2087 return(M); 2088 } 2089 if (i==ncols(w)) 2090 { 2091 intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),1..ncols(w)-1]; 2092 return(M); 2093 } 2094 else 2095 { 2096 intmat M[nrows(w)][i-1]=w[1..nrows(w),1..i-1]; 2097 intmat N[nrows(w)][ncols(w)-i]=w[1..nrows(w),i+1..ncols(w)]; 2098 return(intmatconcat(M,N)); 2099 } 2100 } 2101 if (typeof(w)=="bigintmat") 2102 { 2103 if ((i<1) or (i>ncols(w)) or (ncols(w)==1)) 2104 { 2105 return(w); 2106 } 2107 if (i==1) 2108 { 2109 bigintmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),2..ncols(w)]; 2110 return(M); 2111 } 2112 if (i==ncols(w)) 2113 { 2114 bigintmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),1..ncols(w)-1]; 2115 return(M); 2116 } 2117 else 2118 { 2119 bigintmat MN[nrows(w)][ncols(w)-1]; 2120 MN[1..nrows(w),1..i-1]=w[1..nrows(w),1..i-1]; 2121 MN[1..nrows(w),i..ncols(w)-1]=w[1..nrows(w),i+1..ncols(w)]; 2122 return(MN); 2123 } 2124 } else 2125 { 2126 ERROR("intmatcoldelete: input matrix has to be of type intmat or bigintmat"); 2127 intmat M; return(M); 2390 2128 } 2391 2129 } … … 2683 2421 return(list(coord,coords,latex)); 2684 2422 } 2423 2424 static proc intmatAddFirstColumn (def M,string art) 2425 "USAGE: intmatAddFirstColumn(M,art); M intmat, art string 2426 ASSUME: - M is an integer matrix where a first column of 0's or 1's should be added 2427 @* - art is one of the following strings: 2428 @* + 'rays' : indicating that a first column of 0's should be added 2429 @* + 'points' : indicating that a first column of 1's should be added 2430 RETURN: intmat, a first column has been added to the matrix" 2431 { 2432 if (typeof (M) == "intmat") 2433 { 2434 intmat N[nrows(M)][ncols(M)+1]; 2435 int i,j; 2436 for (i=1;i<=nrows(M);i++) 2437 { 2438 if (art=="rays") 2439 { 2440 N[i,1]=0; 2441 } 2442 else 2443 { 2444 N[i,1]=1; 2445 } 2446 for (j=1;j<=ncols(M);j++) 2447 { 2448 N[i,j+1]=M[i,j]; 2449 } 2450 } 2451 return(N); 2452 } 2453 if (typeof (M) == "bigintmat") 2454 { 2455 bigintmat N[nrows(M)][ncols(M)+1]; 2456 int i,j; 2457 for (i=1;i<=nrows(M);i++) 2458 { 2459 if (art=="rays") 2460 { 2461 N[i,1]=0; 2462 } 2463 else 2464 { 2465 N[i,1]=1; 2466 } 2467 for (j=1;j<=ncols(M);j++) 2468 { 2469 N[i,j+1]=M[i,j]; 2470 } 2471 } 2472 return(N); 2473 } 2474 else 2475 { 2476 ERROR ("intmatAddFirstColumn: input matrix has to be either intmat or bigintmat"); 2477 intmat N; 2478 return (N); 2479 } 2480 } 2481 2482 2483 -
Singular/LIB/primdec.lib
r31b00db r6e3023a 709 709 return(l); 710 710 } 711 def op = option(get); 711 712 def P=basering; 712 713 int i,j,k,m,q,d,o; 713 int n =nvars(basering);714 int n = nvars(basering); 714 715 ideal s,t,u,sact; 715 716 poly ni; … … 769 770 { 770 771 setring RL; 772 771 773 option(redSB); 772 774 t=std(t); … … 785 787 rp=delete(rp,1); 786 788 keep1=keep1+rp; 789 787 790 option(noredSB); 788 791 } … … 803 806 } 804 807 l=l+keep1; 808 option(set,op); 805 809 return(l); 806 810 } … … 1898 1902 list @res,empty; 1899 1903 ideal ser; 1900 option(redSB); 1904 def op = option( get ); 1905 option( redSB ); 1901 1906 list @pr=facstd(i); 1902 1907 //if(size(@pr)==1) … … 1923 1928 // } 1924 1929 // } 1925 option(noredSB); 1930 // option( noredSB ); 1931 option( set, op ); 1926 1932 int j,k,odim,ndim,count; 1927 1933 attrib(@pr[1],"isSB",1); … … 2096 2102 } 2097 2103 2098 if(dim(i) == -1){setring P0;return(ideal(1));} 2099 if((dim(i) == 0) && (npars(P) == 0)) 2104 if( dim(i) == -1 ) 2105 { 2106 option( set,op ); 2107 setring P0; 2108 return( ideal(1) ); 2109 } 2110 if( (dim(i) == 0 ) && ( npars(P) == 0) ) 2100 2111 { 2101 2112 int di = vdim(i); … … 2546 2557 int n=nvars(R); 2547 2558 2559 def op = option(get); 2560 2548 2561 //---Anfang Provisorium 2549 2562 if((size(i)==2) && (w==2)) 2550 2563 { 2551 option( redSB);2552 ideal J =std(i);2553 option( noredSB);2554 if ((size(J)==2)&&(deg(J[1])==1))2564 option( redSB ); 2565 ideal J = std(i); 2566 option( noredSB ); 2567 if ((size(J)==2)&&(deg(J[1])==1)) 2555 2568 { 2556 2569 ideal keep; … … 2577 2590 resu[j]=L; 2578 2591 } 2592 option( set, op ); 2579 2593 return(resu); 2580 2594 } … … 2652 2666 re=convList(pr); 2653 2667 } 2654 return(re); 2668 option( set, op ); 2669 return( re ); 2655 2670 } 2656 2671 /////////////////////////////////////////////////////////////////////////////// … … 7009 7024 primary=resu; 7010 7025 } 7026 option(set,op); 7011 7027 if (intersectOption == "intersect") 7012 7028 { -
Singular/LIB/tropical.lib
r31b00db r6e3023a 1 // /////////////////////////////////////////////////////////2 version="version tropical.lib 4.0.0.0 Jun_2013 "; // $Id$1 // 2 version="version tropical.lib 4.0.0.0 Dec_2013 "; 3 3 category="Tropical Geometry"; 4 4 info=" … … 204 204 LIB "elim.lib"; 205 205 LIB "linalg.lib"; 206 LIB " oldpolymake.lib";206 LIB "polymake.lib"; 207 207 LIB "primdec.lib"; 208 208 LIB "absfact.lib"; … … 650 650 // pass first to a ring where a and @a 651 651 // are variables in order to use maps 652 string @mp= string(minpoly);652 poly mp=minpoly; 653 653 ring INTERRING=char(LIFTRING),(t,@a,a),dp; 654 execute("poly mp=" + @mp + ";");654 poly mp=imap(LIFTRING,mp); 655 655 ideal LIFT=imap(LIFTRING,LIFT); 656 656 kill LIFTRING; … … 1012 1012 { 1013 1013 def GLOBALRING=basering; 1014 string @mp= string(minpoly);1014 number mp=minpoly; 1015 1015 execute("ring LOCALRING=("+charstr(basering)+"),("+varstr(basering)+"),ds;"); 1016 execute("minpoly= " + @mp + ";");1017 1016 poly f=imap(GLOBALRING,f); 1017 minpoly=imap(GLOBALRING,mp); 1018 1018 } 1019 1019 // check if a substitution is necessary … … 1049 1049 if (minpoly!=0) 1050 1050 { 1051 string @mp= string(minpoly);1051 poly mp=minpoly; 1052 1052 def OLDRING=basering; 1053 1053 execute("ring NEWRING=0,("+varstr(basering)+","+parstr(basering)+"),ds;"); 1054 execute("ideal I = " + @mp + ";"); 1055 I = I,imap(OLDRING,f); 1054 ideal I=imap(OLDRING,mp),imap(OLDRING,f); 1056 1055 } 1057 1056 else … … 1065 1064 w=NewtP[jj]-NewtP[jj+1]; 1066 1065 ggteiler=gcd(w[1],w[2]); 1067 zw=w[1] /ggteiler;1068 w[1]=w[2] /ggteiler;1066 zw=w[1] div ggteiler; 1067 w[1]=w[2] div ggteiler; 1069 1068 w[2]=zw; 1070 1069 // if we have introduced a third variable for the parameter, then w needs a third component … … 1134 1133 else 1135 1134 { 1136 string @mp= string(minpoly);1135 poly mp=minpoly; 1137 1136 ring degreering=0,a,dp; 1138 execute("poly mp=" + @mp + ";");1137 poly mp=imap(countring,mp); 1139 1138 numberofbranchesfound=numberofbranchesfound+deg(mp); 1140 1139 setring countring; … … 1158 1157 ring r=0,(x,y),ds; 1159 1158 poly f=x2-y4+x5y7; 1160 puiseuxExpansion(puiseuxExpansion(f,3)); 1159 puiseuxExpansion(f,3,"subst"); 1160 displayPuiseuxExpansion(puiseuxExpansion(f,3)); 1161 1161 } 1162 1162 … … 1540 1540 ggt=gcd(normalvector[1],normalvector[2]); // the gcd of the entries 1541 1541 zw=normalvector[2]; // create the outward pointing normal vector 1542 normalvector[2]=-normalvector[1] /ggt;1543 normalvector[1]=zw /ggt;1542 normalvector[2]=-normalvector[1] div ggt; 1543 normalvector[1]=zw div ggt; 1544 1544 if (size(#)==0) // we are computing w.r.t. minimum 1545 1545 { … … 2718 2718 if(k<>j) // except the unit vector of the non-zero component 2719 2719 { 2720 intvec u; u[size(w)]=0; u[k]=1; 2721 O[r,1..size(w)]=u; r=r+1; 2720 intvec u; 2721 u[size(w)]=0; 2722 u[k]=1; 2723 O[r,1..size(w)]=u; 2724 r=r+1; 2725 kill u; 2722 2726 } 2723 2727 } … … 4708 4712 wneu=NewtP[1]-NewtP[2]; 4709 4713 int ggteiler=gcd(wneu[1],wneu[2]); 4710 wneu[1]=-wneu[1] /ggteiler;4711 wneu[2]=wneu[2] /ggteiler;4714 wneu[1]=-wneu[1] div ggteiler; 4715 wneu[2]=wneu[2] div ggteiler; 4712 4716 if (wneu[1]>0) 4713 4717 { … … 4807 4811 for (jj=1;jj<=anzahlvariablen+numberdeletedvariables-1;jj++) 4808 4812 { 4809 PARA[jj]=(PARA[jj]+a[jj+1])*t^(tw[jj+1]*tweight /ww[1]);4813 PARA[jj]=(PARA[jj]+a[jj+1])*t^(tw[jj+1]*tweight div ww[1]); 4810 4814 } 4811 4815 // if we have reached the stop-level, i.e. either … … 5353 5357 if (koeffizienten[j-ord(p)+1]!=0) 5354 5358 { 5355 if ((j-(N*wj) /w1)==0)5359 if ((j-(N*wj) div w1)==0) 5356 5360 { 5357 5361 dp=dp+displaycoef(koeffizienten[j-ord(p)+1]); 5358 5362 } 5359 if (((j-(N*wj) /w1)==1) and (N!=1))5363 if (((j-(N*wj) div w1)==1) and (N!=1)) 5360 5364 { 5361 5365 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t^(1/"+string(N)+")"; 5362 5366 } 5363 if (((j-(N*wj) /w1)==1) and (N==1))5367 if (((j-(N*wj) div w1)==1) and (N==1)) 5364 5368 { 5365 5369 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t"; 5366 5370 } 5367 if (((j-(N*wj) /w1)==-1) and (N==1))5371 if (((j-(N*wj) div w1)==-1) and (N==1)) 5368 5372 { 5369 5373 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*1/t"; 5370 5374 } 5371 if (((j-(N*wj) /w1)==-1) and (N!=1))5375 if (((j-(N*wj) div w1)==-1) and (N!=1)) 5372 5376 { 5373 5377 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*1/t^(1/"+string(N)+")"; 5374 5378 } 5375 if (((j-(N*wj) /w1)>1) and (N==1))5376 { 5377 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t^"+string(j-(N*wj) /w1);5378 } 5379 if (((j-(N*wj) /w1)>1) and (N!=1))5380 { 5381 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t^("+string(j-(N*wj) /w1)+"/"+string(N)+")";5382 } 5383 if (((j-(N*wj) /w1)<-1) and (N==1))5379 if (((j-(N*wj) div w1)>1) and (N==1)) 5380 { 5381 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t^"+string(j-(N*wj) div w1); 5382 } 5383 if (((j-(N*wj) div w1)>1) and (N!=1)) 5384 { 5385 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t^("+string(j-(N*wj) div w1)+"/"+string(N)+")"; 5386 } 5387 if (((j-(N*wj) div w1)<-1) and (N==1)) 5384 5388 { 5385 5389 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*1/t^"+string(wj-j); 5386 5390 } 5387 if (((j-(N*wj) /w1)<-1) and (N!=1))5388 { 5389 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*1/t^("+string((N*wj) /w1-j)+"/"+string(N)+")";5391 if (((j-(N*wj) div w1)<-1) and (N!=1)) 5392 { 5393 dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*1/t^("+string((N*wj) div w1-j)+"/"+string(N)+")"; 5390 5394 } 5391 5395 if (j<deg(p)) … … 5495 5499 { 5496 5500 setring LIFTRing; 5497 string @mp= string(minpoly);5501 poly mp=minpoly; 5498 5502 execute("ring TESTRing=("+charstr(LIFTRing)+"),("+varstr(BASERING)+"),dp;"); 5499 execute("minpoly= " + @mp + ";");5503 minpoly=number(imap(LIFTRing,mp)); 5500 5504 ideal i=imap(BASERING,i); 5501 5505 } … … 5769 5773 wneu=NewtP[jjj]-NewtP[jjj+1]; 5770 5774 int ggteiler=gcd(wneu[1],wneu[2]); 5771 wneu[1]=-wneu[1] /ggteiler;5772 wneu[2]=wneu[2] /ggteiler;5775 wneu[1]=-wneu[1] div ggteiler; 5776 wneu[2]=wneu[2] div ggteiler; 5773 5777 if (wneu[1]>0) 5774 5778 { … … 5869 5873 { 5870 5874 execute("ring PARARing=("+string(char(basering))+",@a),t,ls;"); 5871 minpoly=number(imap(PREGFANRING,m)); // ?5875 minpoly=number(imap(PREGFANRING,m)); 5872 5876 } 5873 5877 else … … 5908 5912 for (jjj=1;jjj<=anzahlvariablen+numberdeletedvariables-1;jjj++) 5909 5913 { 5910 PARA[jjj]=(PARA[jjj]+zeros[size(zeros)][jjj+1])*t^(ww[jjj+1]*tweight /ww[1]);5914 PARA[jjj]=(PARA[jjj]+zeros[size(zeros)][jjj+1])*t^(ww[jjj+1]*tweight div ww[1]); 5911 5915 } 5912 5916 // delete the last entry in zero, since that one has … … 6033 6037 { 6034 6038 I=subst(i,var(lastvar),0); 6035 while ( subst(I[1],t,0)==0)6039 while ((I[1]!=0) and (subst(I[1],t,0)==0)) 6036 6040 { 6037 6041 I[1]=I[1]/t; … … 7814 7818 poly denom=delta*hn^5; 7815 7819 poly ggt=gcd(num,denom); 7816 num=num /ggt;7817 denom=denom /ggt;7820 num=num div ggt; 7821 denom=denom div ggt; 7818 7822 setring BASERING; 7819 7823 poly num=imap(TRING,num); -
Singular/ipid.cc
r31b00db r6e3023a 433 433 else if ((IDTYP(h)==RING_CMD)||(IDTYP(h)==QRING_CMD)) 434 434 rKill(h); 435 else 435 else if (IDDATA(h)!=NULL) 436 436 s_internalDelete(IDTYP(h),IDDATA(h),r); 437 437 // general ------------------------------------------------------------- -
Singular/singular-libs
r31b00db r6e3023a 23 23 paraplanecurves.lib phindex.lib \ 24 24 pointid.lib poly.lib \ 25 oldpolymake.libpresolve.lib primdec.lib primdecint.lib \25 presolve.lib primdec.lib primdecint.lib \ 26 26 primitiv.lib qhmoduli.lib random.lib realclassify.lib \ 27 27 realrad.lib reesclos.lib resbinomial.lib \ … … 32 32 surf.lib surfacesignature.lib surfex.lib \ 33 33 teachstd.lib toric.lib triang.lib \ 34 tropical.libweierstr.lib zeroset.lib34 weierstr.lib zeroset.lib 35 35 36 36 # libs in beta testing: … … 41 41 # for new libs 42 42 43 SLIB1 = classifyci.lib classify_aeq.lib ffsolve.lib decomp.lib template.lib\ 43 SLIB1 = classifyci.lib classify_aeq.lib \ 44 divisors.lib \ 45 ffsolve.lib decomp.lib template.lib \ 44 46 findifs.lib finitediff.lib \ 47 gitfan.lib \ 45 48 locnormal.lib modnormal.lib \ 46 JMBTest.lib JMSConst.lib multigrading.lib parallel.lib \ 49 JMBTest.lib JMSConst.lib multigrading.lib\ 50 orbitparam.lib parallel.lib polymake.lib\ 47 51 realizationMatroids.lib resources.lib ringgb.lib \ 48 52 schreyer.lib symodstd.lib derham.lib polybori.lib ellipticcovers.lib \ 49 schubert.lib tasks.lib 53 schubert.lib tasks.lib tropical.lib 50 54 51 55 PLIBS = bfun.lib central.lib dmod.lib dmodapp.lib dmodvar.lib fpadim.lib \ … … 53 57 ncfactor.lib nctools.lib perron.lib qmatrix.lib \ 54 58 ncall.lib 55 56 59 60 -
Singular/subexpr.cc
r31b00db r6e3023a 455 455 void s_internalDelete(const int t, void *d, const ring r) 456 456 { 457 assume(d!=NULL); 457 458 switch (t) 458 459 { -
Tst/Short/bug_newstruct.res.gz.uu
r31b00db r6e3023a 1 1 begin 640 bug_newstruct.res.gz 2 M'XL(" )/1P5$``V)U9U]N97=S=')U8W0N<F5S`)6236O"0!"&[_D5+\%#`A+=3 M JJUMZI86+P$1P=*K^=K(0K*1_6CIO^_&F(W7GG9FWOEXAMGCYS;9`R`4N^0#4 M OE8ZJGGNQ][QICQ0V.")"ZZ#,/:Z%Y0B-^>38#]*2U/HR%J1TIEV50L*9R\C5 M S&9PR5#:5)535]$H!?ZEY-]<M=*?^C57-M<T329*-471"@:=\=H/1[C'NS%/6 M $:XENU%>T[ZLB#WKO2'MO!1<P8B255RPLA>8E*U$6Q1&2E:""U@W9U4K&?IF7 M M4W&^J5O81NF;L@SQ4"-@YM-YA2':,#?C%"$=$*WR&;`*C(!T6ID2O&SP"398 M ?[WODNTDF(?0+1K6Y,S25="_%W9=,EBN%N&_T`FQ[&ZNHR?VO!;Z&G[M.E//9 =,=_\(?/NI,2>M/L4W<F-"D@83[P_D)E)>TH"````2 M'XL("+JRJ5(``V)U9U]N97=S=')U8W0N<F5S`)62RVK#,!!%]_Z*B\G"AN!$ 3 M;=JF=:/2DHTAA$!*MXD?<A#84I#DEOY]I3B6L^U*,W/G<8;1_G.=;0$0BDWV 4 M@=!HDS2\"--@?U7N*&SPP`4W49P&[@6E*+K30;`?;517FL1:B3:Y\57W%-Y> 5 M))C-X).A35?77GU(1BD*SQ7_YEJJ<!HV7-O<KFUS4>DI2BD83,Z;,![A'F_& 6 M/"6XE&Q&>4G[LC(-K/>&H_..X!J=J%C-!:MZ@2DE%619=DJQ"ES`N@6KI6+H 7 MFS4V&<N7OH5M>/1#GBD&:NS\;#*GV"4#_FJ$(L0);I'5@%7F`D(:Y%KSD\`D 8 MVWZ];[+U))K',!(M:PMFZ6J8WS.[+!DM'A;QO]`)L>Q^KJ<G]KP6^A)^=9UI 9 >X)FO_I!Y<U)B3^H^A3MYIR,2IY/@#R.]3(I*`@`` 10 10 ` 11 11 end -
Tst/Short/polymake.tst
r31b00db r6e3023a 3 3 LIB "tst.lib"; 4 4 tst_init(); 5 LIB " oldpolymake.lib"; // TODO: update test result!5 LIB "polymake.lib"; 6 6 /////////////////////////////////////////////////////////////////////////// 7 7 // A) Test for Procedures using Polymake 8 8 /////////////////////////////////////////////////////////////////////////// 9 9 example polymakePolytope; 10 example newtonPolytope ;10 example newtonPolytopeP; 11 11 example newtonPolytopeLP; 12 example normalFan ;12 example normalFanL; 13 13 example groebnerFan; 14 example intmatToPolymake;15 example polymakeToIntmat;16 14 /////////////////////////////////////////////////////////////////////////// 17 15 // B) Test for Procedures using Topcom -
Tst/Short/sba_s.res.gz.uu
r74c4462 r6e3023a 1 begin 640 sba_s.res.gz 2 M'XL("-A/A5(``W-B85]S+G)E<P#M6UM/&T<4?N=7C%"DVO$E<Y_9.""EA#:H 3 M4:@P>4BB-"+$%&-\D6U40I7_WMF+9\:[9S:[6F@5-3P`GLLWYWQSSG?.KF!X 4 M^N+H-4*([*,G3]"CH\^/=H;9$.TC^SOKQ].K3V<?5_WU:HUZ/70]G_V)UJ/5 5 M>H4NYLMXSJ[FR>J;U6B%EN=?SJ_'YVAT>S9=7)N!B^5\BHZ&P^<'/T4:+<X6 6 MHZ7=)[P391^].OH9[9K3^M?C3[L#.Z/VD1G\.)Z-UZWV8"?^B?;W,^-FH[_Z 7 MJ_79VJ[6^PXS2C!W%_/K+]N8!+M5A/318CD_1Y.S]>IF>39KC6=K=->V\X:7 8 MO^T'0XSYBE?,T!ZZ<X@BF5B.#4DG9H91C%FW==O"_7YKUFZWNY\7;K%*43Z/ 9 MSJ[1V*S.CC8+W9HH69-\K=8)[J?1[/S2K-[][?GI\,W)\]YN)YWQ]U&2[!O= 10 M+N;+=>O$FV#^Q-B;R$P?K6^6LU9RB#=I3/UJ/WCDTBA'VV6.-X8]WA@)\,98 11 M.6_=RRWFF,@Q=SF?SO]L.?ZZGNU,U:*PL]M[Z4*$10$>.0GPR%D)CUQX/'+I 12 M>.0JXS%-G'ST<>VQR*,`BX*4L=@CN?@3+,=B>K8?1D*$N3MX>_#JZ`",/J$" 13 MK(DHP)HD):Q)YK$FN6--BFW6\L$GI4>;5`':9/0MVG+AIP@8?I8^/_H4J\/@ 14 M=O`I$:!1J0"-*BJA41./1DT=C9IE-([.YVGDS=J.-\V=/FO9SRB\,PX8;ARX 15 M[H<XO-L.NPC'*^TGND&<=J^Z$[>*]QW%;C0Y/U9RY`$F1R_,J=8<@I-#3(5J 16 M3??P8/IL[ZY'!]-.Q_E%,$UO)MZX>'S;FK;-M[NV@^#I_,3,8S<JT]$8^FIO 17 MVB&#JQC<_-@&UYMK3^!1!]VV)O$!5[E3"-XLG'0Z;C2S[:L=X)ZQJ(=:YN3V 18 MP-)(2&;6^+V9^!`O<EC:;74$T>1<BT_3:WAOKO0#,@N-C2;P$U-[OKF4;[AV 19 ME%!9X)KDN*:ZGR-CZC#9MBF,]KV5SF+&,Q.WW6/2AB=AEO5BLAT>'+M,NXO9 20 MV^SBR?F%9"*<>N,GWKA+",(3W\%\([%H6Z]X9/.-".SEVV4^X4@LX_9W%D@X 21 MDDIS(.&V!8ND>FP_16#*D52#<RE')"NF')&BF'*IOI:F7*:T^91[?&E!,FG- 22 M)=U&1+^1=$I43#JEH*3+Q-/=F29`TCV^_(-YB:=9,/&T@!)/JZUHUU%IXIG3 23 MJ-T:D6+R1>R;R1>)DN2+MLV)(C_Y[+U03*#DHYC94*6XI%<H)-]6H:-8@1E( 24 M<01F("4N0VC:BX,92(G7;E'BVBT:M]YQ`DS<4M/4Q@\TK<F>N=[)LSTYB(.C 25 M;><CUTK05#QS3FXZMXD[/U75U<*L6U_DC:->A;L86-,R-36`YM&F->Z2+O;V 26 MZ/RDTS&:ZJB;Q/[.5%7]27]G8LID?&T2OWO1/>DFIKII&7;WTO<WU=^`OZG, 27 M%OSEM,1?SDO\Y;+$7ZY+_!6XU%]!`7_M<Z'GK^`E_@H)^BL*5^A9+7&)O[)P 28 MA?Y.7N*OE*7^2EWB[]8%*USBL**@PZIPAY[9JA#MGMFJ<(?>3EV(=F^GIJ4. 29 M:PXXG#2_GJ]IKQOP-6UW"[Y&A?OS+(X*D>Y9'!7NS]]9B'1_IR[SE6%(J]+& 30 MPSG+<(E2,0PJ%<,E2L5PB5(Q4J)4C)0H%2.E2L7B-OCK3O9<Q7;>'9X<HY/# 31 M%V\.3H^.7P^?FJUXY_3X=S3L;8WRG>'1NT-T_`OZ]:FM8<P.#M^^VPSS_S7F 32 MAMC>R^_&Y.\%<_,FJGK0*@B=-[/X(3#E?XQIF:T1M3^XK889]_75(Y9!R!I" 33 MEO\6)H$P(PA30)CB03`35FM$ZP]>JV%FY0M,!1"<T<HLJ,J@%*16-0.M8:E^ 34 M"$NU)1>.6Q@?E,2&]!(P(AK26\/2ZO36L%3;*E9=QRDHY-5M?A!0D`@)LAL@ 35 MXD%`+;LU5/<[XA<$;<IO'="XHE6/7`)">Z,>-JE>TIK#,K#+%R`L!F$)!,LA 36 M6`$^/01@:<IP'>W%(#P8::2I^M:!A3D&*U%SCD&-"'*<E3BP.0$C3DBP++,( 37 M/""J_F"I09KKX"H(5RO07D&;XO)ZN!NB:ZBQD/IAJ([`&&E.-07M;4XUK')! 38 MW$WEJQ[4&M:[INK<'!94#J;!T"/@:Q@"UB@85\`"&L*U/->(:0V7@,9,-X:% 39 M&8$SL#G3DM3"C2MA]6B6\$LZD`R.J\."!8N!#Z8,;KK@_I"!9%`"`5/X21?N 40 M!^!^F8*UD*;O*`*Q#);:"$QO#K\-!>T&837\5@VD@X.-08!G\5`\1Z!N!'G. 41 M2B&8AAI\S:@"RB3!F%;5D:6`7P@0`8:UA*4#K%I$J``O`FP=!=PZPMA,<;`F 42 M4A["WI`.QS=,NU:P3-T#[8$V]1YHESB@)_=!N\8U:=_4R!IO*2EH/@.['3-: 43 M^3%<@7)5!S<"<8F&.0&!*5C&`L`4['."P);K.N^9,5PL&[,-UY[F;%,,"M8] 44 ML,WA:AD"CJMEC:B&'^\TW/'@ZLU.@&@%`FL0U[@()3N&VW?X,0FN[Q0L\(00 45 M^,$.?--DA#ZA.O`B!*S(7,#O+#!8ZPG<$</0C,,7"=H>@8$=X)LHN,%LSC<E 46 G<#*&^-[\Y86R?UW'XK\XCO]7)_Z/G)M5B[0'CW;^`87ADYAR-``` 1 begin 644 sba_s.res.gz 2 M'XL("%K7JE(``W-B85]S+G)E<P"5F&UOVS80Q]_G4Q!&@,FQK8J/HN#%0)85 3 M:+!B`Y+NQ5!TA>,XM>-'6`J6>LAW'R7*Y)&2F-4O9(K'(T\_W?U)^^[3KS>_ 4 M(X3P!+U[A\YO'L[/[NHN$B/3IG%ISN^G7_.XR`LT&J'U;OL-%?.\R-'C[E#: 5 MS&A6C7[.YSDZS+[/ULL9FK],-_NUZG@\[#;HYN[NZOJG3*+]=#\_&#\.5A0Q 6 M^GCS"^JIU>+U\KXW-I9T@E3GU^5V643]\5GYC2:3.KCM_)\X+Z:%&2TG=LZL 7 MFK.WWZV_NW/BQ([".$;[PVZ&5M,B?SY,M]%R6Z!CW]@5EW_-C0*C/N6(+;I$ 8 M1SLCKPR'I8)TJRR4)`D=1B]1$L?1MM_O#Q_V=G"J9WF83]=HJ4;72ZN!=DQ6 9 MC:D^>5'->S_?SA9J=.^WJT]W?]Y>C7H#;8%^!%=^\Y?][E!$M\!`H6$)#'7H 10 M\^+YL(VJ18!1A?IJ;@!<DGG8%AXWF@!N%'=PHS3,;;APR%'ND5OL-KMOD>4W 11 M!+'3](<0#GJC#S9%:-;!D>$.CHP&.#(..#)A.;*TYJ@+Q\\^)@%%EG50Y#A$ 12 M<82]_./4HZC7AFG$>3>[Z[^N/]Y<MV8?3SNH\:R#FL`!:H(":H)9:H*[U/SD 13 M$P)@$VD'-I&]A<U+OQ2WII_!![,OI3]"T$V^E'=@3-,.C&D6P"@QP"B)Q2AI 14 MC7$^V^G,V_8M-\FL/DL1UPB/Z@$4&SNYC+L8'MVTRY)RI+DCIQDWPZ?ARHYB 15 ML45L>ZOU2R5'8,)JZ;U:U82#DVH1M4-%F\MDO/GY\C@BX\U@8)\+)T2_F=)Q 16 M?_$2;?KJ<NS;*9BVKY0]L;U"]Y93/UUN!GC\5$ZNOMS)Y>FU5].C`7J)5N4" 17 M3]XJ.#D-7`T&MK>.[=5T,!`L&J%(K=P?&XP8UV$M/RO#EW*0G4M:5PN(5.N: 18 M^8E^#9_5*_V"U$`5HTK\*M01#)>P$VN+A(@&:^RQ)C+V8&SLG-0-A9(8C+01 19 M4U:'Z#X>%28],374F\7V_OH/6VG'DM[)BU7K-XH),P+Z;T&_+0C,JF=OK3=< 20 MBK9Y*I:9>L,\`?6V\`L.ES)NVK2CX+"6YHZ"<P4+:STV=UEKR6&MP5[)84&; 21 M)8<%;Y:<UM=@R=5*ZY?<Q<),4DNK5W0G$7VCZ%+^/XLN3=N*KA9/^\XD;BFZ 22 MB\7?%!2>I)V%)WE;X<G4R7:9!0M/K4:,:X:;Q9?1-XLOXX'BR]QPL@P6GWDO 23 M),%MQ4<2:E*5)(&S0J/XG(V.)&EK!9(D:ZU`@FV%$'T6;ZU`@L%QBV![W"+E 24 MT;LL@)4=J@ZUY0^::'4IQF@U(>HR&O6-.;,G":*UTWO&T\%M99?7HIKOU;CB 25 MT8^-@`WN<6PBJ\543:A^V43+(1TFP$?Z1BMC1,NH-1+HJ445&J$G<XW8\12^ 26 M$7IZ`270DR6^$7AJ<5TMUTIKAH_#VV&%QYI9-^(%9*P%N(,QDZV,>1)@S$F` 27 M,6<!QEP$&',98"R2`&-!`HP%"S`6(LA8R!;&YN<O8)PF`<8I:66<L@#CM)'D 28 M(.I4!AC+1I(#3TD"C&4CR:&G"#"6C20'GOHTV\E8'V\[&#N)K`^\'9#UN;<! 29 M.6OH@0V;)HTLQ\#8T`/HV<ARZ-G0`^C9R'+@B1MZ`#QQ(\NA)PM!IOKPZT&N 30 M?LNLP!RRFR_5BN[SI:2A!2!BTLAP$+$OXPY?7\8=OKZ,.WQ]&7?X^C+N\/5E 31 MW.&K9;R3+VO;[_39%0!F@=V.LM;=CK+`;D=98+>COH8[@'T-=P#[&NX`]C7< 32 M`>QKN`/8UW`'L`CN=K14\=>S^N\`81JC#V>G/ZB$;:G>\C@E]%7=U<.Y:0`_ 33 M;ENU']=7Z\=,`_@QVZK]F+Y:/VH:P(_:5NU']57=G1Z6FW,8+7>H\E_=\K_; 34 0YSQ2Q,[/_@-@&(*'G!8````` 47 35 ` 48 36 end -
Tst/Short/sba_s.stat
r74c4462 r6e3023a 1 1 >> tst_memory_0 :: 138 4468440:3.1.3.sw, 32 bit:spielwiese:x86_64-Darwin:schmetterhemd.mathematik.uni-kl.de:195212522 1 >> tst_memory_1 :: 138 4468440:3.1.3.sw, 32 bit:spielwiese:x86_64-Darwin:schmetterhemd.mathematik.uni-kl.de:209715203 1 >> tst_memory_2 :: 138 4468440:3.1.3.sw, 32 bit:spielwiese:x86_64-Darwin:schmetterhemd.mathematik.uni-kl.de:212694524 1 >> tst_timer_1 :: 138 4468440:3.1.3.sw, 32 bit:spielwiese:x86_64-Darwin:schmetterhemd.mathematik.uni-kl.de:381 1 >> tst_memory_0 :: 1385654086:3.1.3.sw, 64 1385654090:3.1.3.sw, 64 1385655331:3.1.3.sw, 64 1386104581:3.1.3.sw, 64 1386927962:3.1.3.sw, 64 bit:spielwiese:x86_64-Linux:schosch:245535672 2 1 >> tst_memory_1 :: 1385654086:3.1.3.sw, 64 1385654090:3.1.3.sw, 64 1385655331:3.1.3.sw, 64 1386104581:3.1.3.sw, 64 1386927962:3.1.3.sw, 64 bit:spielwiese:x86_64-Linux:schosch:249974784 3 1 >> tst_memory_2 :: 1385654086:3.1.3.sw, 64 1385654090:3.1.3.sw, 64 1385655331:3.1.3.sw, 64 1386104581:3.1.3.sw, 64 1386927962:3.1.3.sw, 64 bit:spielwiese:x86_64-Linux:schosch:250400768 4 1 >> tst_timer_1 :: 1385654086:3.1.3.sw, 64 1385654090:3.1.3.sw, 64 1385655331:3.1.3.sw, 64 1386104581:3.1.3.sw, 64 1386927962:3.1.3.sw, 64 bit:spielwiese:x86_64-Linux:schosch:524 -
Tst/Short/sba_s.tst
r74c4462 r6e3023a 124 124 125 125 int k; 126 for (k= 3; k<=6; k++)126 for (k=6; k>2; k--) 127 127 { 128 128 string bench = cyclicn(k); 129 129 sprintf(bench); 130 130 ideal f; 131 f = sba(i,3,0); 132 f = sba(i,3,1); 133 f = sba(i,2,0); 134 f = sba(i,2,1); 131 135 f = sba(i,1,0); 132 136 f = sba(i,1,1); … … 137 141 sprintf(bench); 138 142 ideal f; 143 f = sba(i,3,0); 144 f = sba(i,3,1); 145 f = sba(i,2,0); 146 f = sba(i,2,1); 139 147 f = sba(i,1,0); 140 148 f = sba(i,1,1); … … 145 153 sprintf(bench); 146 154 ideal f; 155 f = sba(i,3,0); 156 f = sba(i,3,1); 157 f = sba(i,2,0); 158 f = sba(i,2,1); 147 159 f = sba(i,1,0); 148 160 f = sba(i,1,1); … … 153 165 sprintf(bench); 154 166 ideal f; 167 f = sba(i,3,0); 168 f = sba(i,3,1); 169 f = sba(i,2,0); 170 f = sba(i,2,1); 155 171 f = sba(i,1,0); 156 172 f = sba(i,1,1); … … 161 177 sprintf(bench); 162 178 ideal f; 179 f = sba(i,3,0); 180 f = sba(i,3,1); 181 f = sba(i,2,0); 182 f = sba(i,2,1); 163 183 f = sba(i,1,0); 164 184 f = sba(i,1,1); … … 169 189 sprintf(bench); 170 190 ideal f; 191 f = sba(i,3,0); 192 f = sba(i,3,1); 193 f = sba(i,2,0); 194 f = sba(i,2,1); 171 195 f = sba(i,1,0); 172 196 f = sba(i,1,1); -
debian/changelog
r31b00db r6e3023a 1 singular ( 3.1.3.sw-1) unstable; urgency=low1 singular (4.0.0-1) unstable; urgency=low 2 2 3 3 * Initial release. -
doc/COPYING.texi
r31b00db r6e3023a 131 131 @copyright{} Winfried Bruns and Bogdan Ichim 132 132 @* @uref{http://www.mathematik.uni-osnabrueck.de/normaliz/} 133 @item polymake (used by oldpolymake.lib, @pxref{oldpolymake_lib})133 @item polymake (used by polymake.lib, @pxref{polymake_lib}) 134 134 @copyright{} Ewgenij Gawrilow and Michael Joswig 135 135 @* @uref{http://www.polymake.de/} … … 142 142 @copyright{} Oliver Labs (2001-2008), Stephan Holzer (2004-2005) 143 143 @* @uref{http://surfex.AlgebraicSurface.net} 144 @item TOPCOM (used by oldpolymake.lib, @pxref{oldpolymake_lib})144 @item TOPCOM (used by polymake.lib, @pxref{polymake_lib}) 145 145 @copyright{} J@"org Rambau 146 146 @* @uref{http://www.rambau.wm.uni-bayreuth.de/TOPCOM/} -
doc/Makefile.bsp
r31b00db r6e3023a 63 63 CHKSUM_DB = ${DOC_SUBDIR}/chksum 64 64 # which tags to avoid: 65 DOC2TEX_EXAMPLE_EXCLUSIONS = -exclude oldpolymake65 DOC2TEX_EXAMPLE_EXCLUSIONS = -exclude polymake 66 66 # which tags to avoid: 67 TAG = oldpolymake67 TAG = polymake 68 68 DOC2TEX = ${PERL} ./doc2tex.pl -docdir ${DOC_SUBDIR} \ 69 69 -Singular ${SINGULAR} -verbose ${VERBOSE} -make ${MAKE} \ -
doc/NEWS.texi
r31b00db r6e3023a 34 34 @item 35 35 name spaces (@nref{package}) 36 @item 37 algebraic/transcendental extension rewritten 36 38 @end itemize 37 39 Besides theses internal changes, @sc{Singular} version 4 offers many new … … 42 44 @itemize 43 45 @item @nref{sba}: an F5 like Groebner base computation 46 @item @nref{ASSUME}: support for debugging libraries 44 47 @end itemize 45 48 … … 47 50 @itemize 48 51 @item @nref{classifyci_lib}: Isolated complete intersection singularities in characteristic 0 49 @item @ xref{schubert_lib}: Proceduces for Intersection Theory52 @item @nref{derham_lib}: Computation of deRham cohomology 50 53 @end itemize 51 54 … … 54 57 @item 55 58 improved grobcov.lib (@nref{grobcov_lib}) 59 @item 60 fixed normal.lib (@nref{normal_lib}) 56 61 @end itemize 57 62 -
kernel/kInline.h
r74c4462 r6e3023a 1156 1156 // dummy function for function pointer strat->rewCrit being usable in all 1157 1157 // possible choices for criteria 1158 KINLINE BOOLEAN arriRewDummy(poly /*sig*/, unsigned long /*not_sevSig*/, kStrategy /*strat*/, int /*start=0*/)1158 KINLINE BOOLEAN arriRewDummy(poly /*sig*/, unsigned long /*not_sevSig*/, poly /*lm*/, kStrategy /*strat*/, int /*start=0*/) 1159 1159 { 1160 1160 return FALSE; -
kernel/kspoly.cc
r74c4462 r6e3023a 201 201 * TODO: 202 202 * -------------------------------------------- 203 * if strat-> incremental203 * if strat->sbaOrder == 1 204 204 * Since we are subdividing lower index and 205 205 * current index reductions it is enough to … … 207 207 * for a check. This should speed-up checking 208 208 * a lot! 209 * if !strat-> incremental209 * if !strat->sbaOrder == 0 210 210 * We are not subdividing lower and current index 211 211 * due to the fact that we are using the induced … … 223 223 if (!PW->is_sigsafe) 224 224 { 225 poly f1 = p_Copy(PR->GetLmCurrRing(),currRing); 225 poly ftmp; 226 ring rtmp; 227 PR->SetLmCurrRing(); 228 poly f1 = pCopy(PR->GetLmCurrRing()); 229 //PR->GetLm(ftmp,rtmp); 230 //poly f1 = k_LmInit_tailRing_2_currRing(ftmp,rtmp); 226 231 poly f2 = PW->GetLmCurrRing(); 227 232 poly sigMult = pCopy(PW->sig); // copy signature of reducer … … 235 240 printf("--------------\n"); 236 241 #endif 237 sigMult = p p_Mult_qq(f1,sigMult,currRing);242 sigMult = p_Mult_q(f1,sigMult,currRing); 238 243 //#if 1 239 244 #ifdef DEBUGF5 … … 253 258 254 259 #endif 255 pDelete(&f1);260 //pDelete(&f1); 256 261 pDelete(&sigMult); 257 262 // go on with the computations only if the signature of p2 is greater than the … … 262 267 return 3; 263 268 } 264 PW->is_sigsafe = TRUE;269 //PW->is_sigsafe = TRUE; 265 270 } 266 271 PR->is_redundant = FALSE; -
kernel/kstd1.cc
r74c4462 r6e3023a 1422 1422 // standard basis computations 1423 1423 strat->red = redSig; 1424 //strat-> incremental = TRUE;1424 //strat->sbaOrder = 1; 1425 1425 strat->currIdx = 1; 1426 1426 } … … 2253 2253 } 2254 2254 2255 ideal kSba(ideal F, ideal Q, tHomog h,intvec ** w, int incremental, int arri, intvec *hilb,int syzComp,2255 ideal kSba(ideal F, ideal Q, tHomog h,intvec ** w, int sbaOrder, int arri, intvec *hilb,int syzComp, 2256 2256 int newIdeal, intvec *vw) 2257 2257 { … … 2263 2263 BOOLEAN delete_w=(w==NULL); 2264 2264 kStrategy strat=new skStrategy; 2265 if (incremental!=0) 2266 { 2267 strat->incremental = TRUE; 2268 } 2269 else 2270 { 2271 strat->incremental = FALSE; 2272 } 2265 strat->sbaOrder = sbaOrder; 2273 2266 if (arri!=0) 2274 2267 { 2275 2268 strat->rewCrit1 = arriRewDummy; 2276 2269 strat->rewCrit2 = arriRewCriterion; 2270 strat->rewCrit3 = arriRewCriterionPre; 2277 2271 } 2278 2272 else … … 2280 2274 strat->rewCrit1 = faugereRewCriterion; 2281 2275 strat->rewCrit2 = faugereRewCriterion; 2276 strat->rewCrit3 = faugereRewCriterion; 2282 2277 } 2283 2278 -
kernel/kstd2.cc
r74c4462 r6e3023a 48 48 #endif 49 49 50 #define F5C 050 #define F5C 1 51 51 #if F5C 52 52 #define F5CTAILRED 1 53 53 #endif 54 54 55 #define SBA_PRODUCT_CRITERION 0 56 #define SBA_PRINT_ZERO_REDUCTIONS 1 57 #define SBA_PRINT_REDUCTION_STEPS 1 58 #define SBA_PRINT_SIZE_G 1 59 #define SBA_PRINT_SIZE_SYZ 1 60 #define SBA_PRINT_PRODUCT_CRITERION 0 61 55 #define SBA_INTERRED_START 0 56 #define SBA_TAIL_RED 1 57 #define SBA_PRODUCT_CRITERION 0 58 #define SBA_PRINT_ZERO_REDUCTIONS 0 59 #define SBA_PRINT_REDUCTION_STEPS 0 60 #define SBA_PRINT_OPERATIONS 0 61 #define SBA_PRINT_SIZE_G 0 62 #define SBA_PRINT_SIZE_SYZ 0 63 #define SBA_PRINT_PRODUCT_CRITERION 0 64 65 // counts sba's reduction steps 66 #if SBA_PRINT_REDUCTION_STEPS 62 67 long sba_reduction_steps; 68 long sba_interreduction_steps; 69 #endif 70 #if SBA_PRINT_OPERATIONS 71 long sba_operations; 72 long sba_interreduction_operations; 73 #endif 74 63 75 /*********************************************** 64 76 * SBA stuff -- done … … 466 478 467 479 ksReducePoly(h, &(strat->T[ii]), NULL, NULL, strat); 480 #if SBA_PRINT_REDUCTION_STEPS 481 sba_interreduction_steps++; 482 #endif 483 #if SBA_PRINT_OPERATIONS 484 sba_interreduction_operations += pLength(strat->T[ii].p); 485 #endif 468 486 469 487 #ifdef KDEBUG … … 513 531 } 514 532 } 533 } 534 535 KINLINE int ksReducePolyTailSig(LObject* PR, TObject* PW, LObject* Red) 536 { 537 BOOLEAN ret; 538 number coef; 539 540 assume(PR->GetLmCurrRing() != PW->GetLmCurrRing()); 541 Red->HeadNormalize(); 542 /* 543 printf("------------------------\n"); 544 pWrite(Red->GetLmCurrRing()); 545 */ 546 ret = ksReducePolySig(Red, PW, 1, NULL, &coef); 547 548 549 if (!ret) 550 { 551 if (! n_IsOne(coef, currRing->cf)) 552 { 553 PR->Mult_nn(coef); 554 // HANNES: mark for Normalize 555 } 556 n_Delete(&coef, currRing->cf); 557 } 558 return ret; 515 559 } 516 560 … … 623 667 sigSafe = ksReducePolySig(h, &(strat->T[ii]), strat->S_2_R[ii], NULL, NULL, strat); 624 668 #if SBA_PRINT_REDUCTION_STEPS 625 sba_reduction_steps++; 669 if (sigSafe != 3) 670 sba_reduction_steps++; 671 #endif 672 #if SBA_PRINT_OPERATIONS 673 if (sigSafe != 3) 674 sba_operations += pLength(strat->T[ii].p); 626 675 #endif 627 676 // if reduction has taken place, i.e. the reduction was sig-safe … … 685 734 } 686 735 } 736 } 737 738 // tail reduction for SBA 739 poly redtailSba (LObject* L, int pos, kStrategy strat, BOOLEAN withT, BOOLEAN normalize) 740 { 741 #define REDTAIL_CANONICALIZE 100 742 strat->redTailChange=FALSE; 743 if (strat->noTailReduction) return L->GetLmCurrRing(); 744 poly h, p; 745 p = h = L->GetLmTailRing(); 746 if ((h==NULL) || (pNext(h)==NULL)) 747 return L->GetLmCurrRing(); 748 749 TObject* With; 750 // placeholder in case strat->tl < 0 751 TObject With_s(strat->tailRing); 752 753 LObject Ln(pNext(h), strat->tailRing); 754 Ln.sig = L->sig; 755 Ln.sevSig = L->sevSig; 756 Ln.pLength = L->GetpLength() - 1; 757 758 pNext(h) = NULL; 759 if (L->p != NULL) pNext(L->p) = NULL; 760 L->pLength = 1; 761 762 Ln.PrepareRed(strat->use_buckets); 763 764 int cnt=REDTAIL_CANONICALIZE; 765 while(!Ln.IsNull()) 766 { 767 loop 768 { 769 Ln.SetShortExpVector(); 770 if (withT) 771 { 772 int j; 773 j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, &Ln); 774 if (j < 0) break; 775 With = &(strat->T[j]); 776 } 777 else 778 { 779 With = kFindDivisibleByInS(strat, pos, &Ln, &With_s); 780 if (With == NULL) break; 781 } 782 cnt--; 783 if (cnt==0) 784 { 785 cnt=REDTAIL_CANONICALIZE; 786 /*poly tmp=*/Ln.CanonicalizeP(); 787 if (normalize) 788 { 789 Ln.Normalize(); 790 //pNormalize(tmp); 791 //if (TEST_OPT_PROT) { PrintS("n"); mflush(); } 792 } 793 } 794 if (normalize && (!TEST_OPT_INTSTRATEGY) && (!nIsOne(pGetCoeff(With->p)))) 795 { 796 With->pNorm(); 797 } 798 strat->redTailChange=TRUE; 799 int ret = ksReducePolyTailSig(L, With, &Ln); 800 #if SBA_PRINT_REDUCTION_STEPS 801 if (ret != 3) 802 sba_reduction_steps++; 803 #endif 804 #if SBA_PRINT_OPERATIONS 805 if (ret != 3) 806 sba_operations += pLength(With->p); 807 #endif 808 if (ret) 809 { 810 // reducing the tail would violate the exp bound 811 // set a flag and hope for a retry (in bba) 812 strat->completeReduce_retry=TRUE; 813 if ((Ln.p != NULL) && (Ln.t_p != NULL)) Ln.p=NULL; 814 do 815 { 816 pNext(h) = Ln.LmExtractAndIter(); 817 pIter(h); 818 L->pLength++; 819 } while (!Ln.IsNull()); 820 goto all_done; 821 } 822 if (Ln.IsNull()) goto all_done; 823 if (! withT) With_s.Init(currRing); 824 } 825 pNext(h) = Ln.LmExtractAndIter(); 826 pIter(h); 827 pNormalize(h); 828 L->pLength++; 829 } 830 831 all_done: 832 Ln.Delete(); 833 if (L->p != NULL) pNext(L->p) = pNext(p); 834 835 if (strat->redTailChange) 836 { 837 L->length = 0; 838 } 839 840 //if (TEST_OPT_PROT) { PrintS("N"); mflush(); } 841 //L->Normalize(); // HANNES: should have a test 842 assume(kTest_L(L)); 843 return L->GetLmCurrRing(); 687 844 } 688 845 … … 772 929 773 930 ksReducePoly(h, &(strat->T[ii]), NULL, NULL, strat); 931 #if SBA_PRINT_REDUCTION_STEPS 932 sba_interreduction_steps++; 933 #endif 934 #if SBA_PRINT_OPERATIONS 935 sba_interreduction_operations += pLength(strat->T[ii].p); 936 #endif 774 937 775 938 #ifdef KDEBUG … … 942 1105 number coef; 943 1106 ksReducePoly(h,&(strat->T[ii]),strat->kNoetherTail(),&coef,strat); 1107 #if SBA_PRINT_REDUCTION_STEPS 1108 sba_interreduction_steps++; 1109 #endif 1110 #if SBA_PRINT_OPERATIONS 1111 sba_interreduction_operations += pLength(strat->T[ii].p); 1112 #endif 944 1113 #ifdef KDEBUG 945 1114 if (TEST_OPT_DEBUG) … … 1509 1678 // induced Schreyer order. 1510 1679 // The corresponding orders are computed in sbaRing(), depending 1511 // on the flag strat->incremental 1512 long zeroreductions = 0; 1513 long product_criterion = 0; 1514 long size_g = 0; 1515 long size_syz = 0; 1680 // on the flag strat->sbaOrder 1681 #if SBA_PRINT_ZERO_REDUCTIONS 1682 long zeroreductions = 0; 1683 #endif 1684 #if SBA_PRINT_PRODUCT_CRITERION 1685 long product_criterion = 0; 1686 #endif 1687 #if SBA_PRINT_SIZE_G 1688 long size_g = 0; 1689 #endif 1690 #if SBA_PRINT_SIZE_SYZ 1691 long size_syz = 0; 1692 #endif 1516 1693 // global variable 1517 sba_reduction_steps = 0; 1518 1519 ideal F = F0; 1694 #if SBA_PRINT_REDUCTION_STEPS 1695 sba_reduction_steps = 0; 1696 sba_interreduction_steps = 0; 1697 #endif 1698 #if SBA_PRINT_OPERATIONS 1699 sba_operations = 0; 1700 sba_interreduction_operations = 0; 1701 #endif 1702 1703 ideal F1 = F0; 1520 1704 ring sRing, currRingOld; 1521 1705 currRingOld = currRing; 1522 if (strat-> incremental)1706 if (strat->sbaOrder == 1 || strat->sbaOrder == 3) 1523 1707 { 1524 1708 sRing = sbaRing(strat); … … 1526 1710 { 1527 1711 rChangeCurrRing (sRing); 1528 F = idrMoveR (F0, currRingOld, currRing); 1529 } 1530 } 1531 #if 0 1712 F1 = idrMoveR (F0, currRingOld, currRing); 1713 } 1714 } 1715 // sort ideal F 1716 ideal F = idInit(IDELEMS(F1),F1->rank); 1717 intvec *sort = idSort(F1); 1718 for (int i=0; i<sort->length();++i) 1719 F->m[i] = F1->m[(*sort)[i]-1]; 1720 #if SBA_INTERRED_START 1721 F = kInterRed(F,NULL); 1722 #endif 1723 #if F5DEBUG 1532 1724 printf("SBA COMPUTATIONS DONE IN THE FOLLOWING RING:\n"); 1533 1725 rWrite (currRing); 1726 printf("ordSgn = %d\n",currRing->OrdSgn); 1534 1727 printf("\n"); 1535 1728 #endif … … 1542 1735 int hilbeledeg=1,hilbcount=0,minimcnt=0; 1543 1736 LObject L; 1544 // BOOLEAN withT = FALSE;1737 BOOLEAN withT = TRUE; 1545 1738 strat->max_lower_index = 0; 1546 1739 … … 1605 1798 #endif 1606 1799 if (strat->Ll== 0) strat->interpt=TRUE; 1800 /* 1607 1801 if (TEST_OPT_DEGBOUND 1608 1802 && ((strat->honey && (strat->L[strat->Ll].ecart+currRing->pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg)) 1609 1803 || ((!strat->honey) && (currRing->pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg)))) 1610 1804 { 1611 /* 1612 *stops computation if 1613 * 24 IN test and the degree +ecart of L[strat->Ll] is bigger then 1614 *a predefined number Kstd1_deg 1615 */ 1805 1806 //stops computation if 1807 // 24 IN test and the degree +ecart of L[strat->Ll] is bigger then 1808 //a predefined number Kstd1_deg 1616 1809 while ((strat->Ll >= 0) 1617 1810 && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL) … … 1623 1816 else strat->noClearS=TRUE; 1624 1817 } 1625 if (strat->incremental && pGetComp(strat->L[strat->Ll].sig) != strat->currIdx) 1818 */ 1819 if (strat->sbaOrder == 1 && pGetComp(strat->L[strat->Ll].sig) != strat->currIdx) 1626 1820 { 1627 1821 strat->currIdx = pGetComp(strat->L[strat->Ll].sig); … … 1630 1824 // 2. generation of new principal syzygy rules for syzCriterion 1631 1825 f5c ( strat, olddeg, minimcnt, hilbeledeg, hilbcount, srmax, 1632 1826 lrmax, reduc, Q, w, hilb ); 1633 1827 #endif 1634 1828 // initialize new syzygy rules for the next iteration step 1635 1829 initSyzRules(strat); 1830 1636 1831 } 1637 1832 /********************************************************************* 1638 * interrreduction step is done, we can go on with the next iteration1639 * step of the signature-based algorithm1640 ********************************************************************/1833 * interrreduction step is done, we can go on with the next iteration 1834 * step of the signature-based algorithm 1835 ********************************************************************/ 1641 1836 /* picks the last element from the lazyset L */ 1642 1837 strat->P = strat->L[strat->Ll]; 1643 1838 strat->Ll--; 1644 //#if 1 1839 /* reduction of the element choosen from L */ 1840 1841 if (!strat->rewCrit2(strat->P.sig, ~strat->P.sevSig, strat->P.GetLmCurrRing(), strat, strat->P.checked+1)) { 1842 //#if 1 1645 1843 #ifdef DEBUGF5 1646 Print("SIG OF NEXT PAIR TO HANDLE IN SIG-BASED ALGORITHM\n"); 1647 Print("-------------------------------------------------\n"); 1648 pWrite(strat->P.sig); 1649 pWrite(pHead(strat->P.p)); 1650 pWrite(pHead(strat->P.p1)); 1651 pWrite(pHead(strat->P.p2)); 1652 Print("-------------------------------------------------\n"); 1653 #endif 1654 if (pNext(strat->P.p) == strat->tail) 1655 { 1656 // deletes the short spoly 1844 Print("SIG OF NEXT PAIR TO HANDLE IN SIG-BASED ALGORITHM\n"); 1845 Print("-------------------------------------------------\n"); 1846 pWrite(strat->P.sig); 1847 pWrite(pHead(strat->P.p)); 1848 pWrite(pHead(strat->P.p1)); 1849 pWrite(pHead(strat->P.p2)); 1850 Print("-------------------------------------------------\n"); 1851 #endif 1852 if (pNext(strat->P.p) == strat->tail) 1853 { 1854 // deletes the short spoly 1855 /* 1657 1856 #ifdef HAVE_RINGS 1658 if (rField_is_Ring(currRing)) 1659 pLmDelete(strat->P.p); 1857 if (rField_is_Ring(currRing)) 1858 pLmDelete(strat->P.p); 1859 else 1860 #endif 1861 pLmFree(strat->P.p); 1862 */ 1863 // TODO: needs some masking 1864 // TODO: masking needs to vanish once the signature 1865 // sutff is completely implemented 1866 strat->P.p = NULL; 1867 poly m1 = NULL, m2 = NULL; 1868 1869 // check that spoly creation is ok 1870 while (strat->tailRing != currRing && 1871 !kCheckSpolyCreation(&(strat->P), strat, m1, m2)) 1872 { 1873 assume(m1 == NULL && m2 == NULL); 1874 // if not, change to a ring where exponents are at least 1875 // large enough 1876 if (!kStratChangeTailRing(strat)) 1877 { 1878 WerrorS("OVERFLOW..."); 1879 break; 1880 } 1881 } 1882 // create the real one 1883 ksCreateSpoly(&(strat->P), NULL, strat->use_buckets, 1884 strat->tailRing, m1, m2, strat->R); 1885 1886 } 1887 else if (strat->P.p1 == NULL) 1888 { 1889 if (strat->minim > 0) 1890 strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing); 1891 // for input polys, prepare reduction 1892 strat->P.PrepareRed(strat->use_buckets); 1893 } 1894 if (strat->P.p == NULL && strat->P.t_p == NULL) 1895 { 1896 red_result = 0; 1897 } 1660 1898 else 1661 #endif 1662 pLmFree(strat->P.p); 1663 1664 // TODO: needs some masking 1665 // TODO: masking needs to vanish once the signature 1666 // sutff is completely implemented 1667 strat->P.p = NULL; 1668 poly m1 = NULL, m2 = NULL; 1669 1670 // check that spoly creation is ok 1671 while (strat->tailRing != currRing && 1672 !kCheckSpolyCreation(&(strat->P), strat, m1, m2)) 1673 { 1674 assume(m1 == NULL && m2 == NULL); 1675 // if not, change to a ring where exponents are at least 1676 // large enough 1677 if (!kStratChangeTailRing(strat)) 1678 { 1679 WerrorS("OVERFLOW..."); 1680 break; 1681 } 1682 } 1683 // create the real one 1684 ksCreateSpoly(&(strat->P), NULL, strat->use_buckets, 1685 strat->tailRing, m1, m2, strat->R); 1686 1687 } 1688 else if (strat->P.p1 == NULL) 1689 { 1690 if (strat->minim > 0) 1691 strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing); 1692 // for input polys, prepare reduction 1693 strat->P.PrepareRed(strat->use_buckets); 1694 } 1695 1696 if (strat->P.p == NULL && strat->P.t_p == NULL) 1697 { 1698 red_result = 0; 1699 } 1700 else 1701 { 1702 if (TEST_OPT_PROT) 1703 message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(), 1704 &olddeg,&reduc,strat, red_result); 1705 1706 //#if 1 1899 { 1900 //#if 1 1707 1901 #ifdef DEBUGF5 1708 Print("Poly before red: "); 1709 pWrite(strat->P.p); 1710 #endif 1711 /* reduction of the element choosen from L */ 1712 if (!strat->rewCrit2(strat->P.sig, ~strat->P.sevSig, strat, strat->P.checked+1)) { 1902 Print("Poly before red: "); 1903 pWrite(pHead(strat->P.p)); 1904 pWrite(strat->P.sig); 1905 #endif 1713 1906 #if SBA_PRODUCT_CRITERION 1714 if (strat->P.checked == 3) { 1907 if (strat->P.prod_crit) { 1908 #if SBA_PRINT_PRODUCT_CRITERION 1715 1909 product_criterion++; 1716 enterSyz(strat->P, strat); 1910 #endif 1911 int pos = posInSyz(strat, strat->P.sig); 1912 enterSyz(strat->P, strat, pos); 1717 1913 if (strat->P.lcm!=NULL) 1718 1914 pLmFree(strat->P.lcm); … … 1722 1918 } 1723 1919 #else 1724 red_result = strat->red(&strat->P,strat); 1725 #endif 1726 } else { 1727 if (strat->P.lcm!=NULL) 1728 pLmFree(strat->P.lcm); 1729 red_result = 2; 1730 } 1731 if (errorreported) break; 1732 } 1920 red_result = strat->red(&strat->P,strat); 1921 #endif 1922 } 1923 } else { 1924 /* 1925 if (strat->P.lcm != NULL) 1926 pLmFree(strat->P.lcm); 1927 */ 1928 red_result = 2; 1929 } 1930 if (errorreported) break; 1931 1932 //#if 1 1933 #ifdef DEBUGF5 1934 if (red_result != 0) { 1935 Print("Poly after red: "); 1936 pWrite(pHead(strat->P.p)); 1937 pWrite(strat->P.GetLmCurrRing()); 1938 pWrite(strat->P.sig); 1939 printf("%d\n",red_result); 1940 } 1941 #endif 1733 1942 1734 1943 if (strat->overflow) 1735 1944 { 1736 1945 if (!kStratChangeTailRing(strat)) { Werror("OVERFLOW.."); break;} 1737 }1738 if (strat->incremental)1739 {1740 for (int jj = 0; jj<strat->tl+1; jj++)1741 {1742 if (pGetComp(strat->T[jj].sig) == strat->currIdx)1743 {1744 strat->T[jj].is_sigsafe = FALSE;1745 }1746 }1747 }1748 else1749 {1750 for (int jj = 0; jj<strat->tl+1; jj++)1751 {1752 strat->T[jj].is_sigsafe = FALSE;1753 }1754 1946 } 1755 1947 … … 1785 1977 // in the ring case we cannot expect LC(f) = 1, 1786 1978 // therefore we call pContent instead of pNorm 1787 /* 1788 if ((TEST_OPT_INTSTRATEGY) || (rField_is_Ring(currRing))) 1789 { 1790 strat->P.pCleardenom(); 1791 if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL)) 1792 { 1793 strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT); 1979 #if SBA_TAIL_RED 1980 if (strat->sbaOrder != 2) { 1981 if ((TEST_OPT_INTSTRATEGY) || (rField_is_Ring(currRing))) 1982 { 1794 1983 strat->P.pCleardenom(); 1795 } 1796 } 1797 else 1798 { 1799 strat->P.pNorm(); 1800 if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL)) 1801 strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT); 1802 } 1803 */ 1984 if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL)) 1985 { 1986 strat->P.p = redtailSba(&(strat->P),pos-1,strat, withT); 1987 strat->P.pCleardenom(); 1988 } 1989 } 1990 else 1991 { 1992 strat->P.pNorm(); 1993 if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL)) 1994 strat->P.p = redtailSba(&(strat->P),pos-1,strat, withT); 1995 } 1996 } 1997 #endif 1998 1999 // remove sigsafe label since it is no longer valid for the next element to 2000 // be reduced 2001 if (strat->sbaOrder == 1) 2002 { 2003 for (int jj = 0; jj<strat->tl+1; jj++) 2004 { 2005 if (pGetComp(strat->T[jj].sig) == strat->currIdx) 2006 { 2007 strat->T[jj].is_sigsafe = FALSE; 2008 } 2009 } 2010 } 2011 else 2012 { 2013 for (int jj = 0; jj<strat->tl+1; jj++) 2014 { 2015 strat->T[jj].is_sigsafe = FALSE; 2016 } 2017 } 1804 2018 #ifdef KDEBUG 1805 2019 if (TEST_OPT_DEBUG){PrintS("new s:");strat->P.wrp();PrintLn();} … … 1833 2047 // enter into S, L, and T 1834 2048 //if ((!TEST_OPT_IDLIFT) || (pGetComp(strat->P.p) <= strat->syzComp)) 1835 if(!strat->incremental) 1836 { 1837 BOOLEAN overwrite = TRUE; 2049 enterT(strat->P, strat); 2050 strat->T[strat->tl].is_sigsafe = FALSE; 2051 /* 2052 printf("hier\n"); 2053 pWrite(strat->P.GetLmCurrRing()); 2054 pWrite(strat->P.sig); 2055 */ 2056 #ifdef HAVE_RINGS 2057 if (rField_is_Ring(currRing)) 2058 superenterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl); 2059 else 2060 #endif 2061 enterpairsSig(strat->P.p,strat->P.sig,strat->sl+1,strat->sl,strat->P.ecart,pos,strat, strat->tl); 2062 // posInS only depends on the leading term 2063 strat->enterS(strat->P, pos, strat, strat->tl); 2064 if(strat->sbaOrder != 1) 2065 { 2066 BOOLEAN overwrite = FALSE; 1838 2067 for (int tk=0; tk<strat->sl+1; tk++) 1839 2068 { … … 1855 2084 1856 2085 strat->P.sevSig = pGetShortExpVector (strat->P.sig); 2086 int i; 2087 LObject Q; 1857 2088 for(int ps=0;ps<strat->sl+1;ps++) 1858 2089 { 1859 int i = strat->syzl;1860 2090 1861 2091 strat->newt = TRUE; … … 1869 2099 strat->syzmax += setmaxTinc; 1870 2100 } 1871 strat->syz[i]= pCopy(strat->P.sig);2101 Q.sig = pCopy(strat->P.sig); 1872 2102 // add LM(F->m[i]) to the signature to get a Schreyer order 1873 2103 // without changing the underlying polynomial ring at all 1874 p_ExpVectorAdd (strat->syz[i],strat->S[ps],currRing); 2104 if (strat->sbaOrder == 0) 2105 p_ExpVectorAdd (Q.sig,strat->S[ps],currRing); 1875 2106 // since p_Add_q() destroys all input 1876 2107 // data we need to recreate help … … 1881 2112 // the corresponding principal syzygy 1882 2113 // => we do not need to compute the "real" syzygy completely 1883 poly help = p Copy(strat->sig[ps]);2114 poly help = p_Copy(strat->sig[ps],currRing); 1884 2115 p_ExpVectorAdd (help,strat->P.p,currRing); 1885 strat->syz[i] = p_Add_q(strat->syz[i],help,currRing);2116 Q.sig = p_Add_q(Q.sig,help,currRing); 1886 2117 //printf("%d. SYZ ",i+1); 1887 2118 //pWrite(strat->syz[i]); 1888 strat->sevSyz[i] = p_GetShortExpVector(strat->syz[i],currRing); 1889 strat->syzl++; 2119 Q.sevSig = p_GetShortExpVector(Q.sig,currRing); 2120 i = posInSyz(strat, Q.sig); 2121 enterSyz(Q, strat, i); 1890 2122 } 1891 2123 } 1892 2124 } 1893 enterT(strat->P, strat); 1894 strat->T[strat->tl].is_sigsafe = FALSE; 1895 #ifdef HAVE_RINGS 1896 if (rField_is_Ring(currRing)) 1897 superenterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl); 1898 else 1899 #endif 1900 enterpairsSig(strat->P.p,strat->P.sig,strat->sl+1,strat->sl,strat->P.ecart,pos,strat, strat->tl); 1901 // posInS only depends on the leading term 1902 strat->enterS(strat->P, pos, strat, strat->tl); 2125 // deg - idx - lp/rp 2126 // => we need to add syzygies with indices > pGetComp(strat->P.sig) 2127 if(strat->sbaOrder == 0 || strat->sbaOrder == 3) 2128 { 2129 int cmp = pGetComp(strat->P.sig); 2130 int max_cmp = IDELEMS(F); 2131 int* vv = (int*)omAlloc((currRing->N+1)*sizeof(int)); 2132 pGetExpV (strat->P.p,vv); 2133 LObject Q; 2134 int pos; 2135 int idx = p_GetComp(strat->P.sig,currRing); 2136 //printf("++ -- adding syzygies -- ++\n"); 2137 // if new element is the first one in this index 2138 if (strat->currIdx < idx) { 2139 for (int i=0; i<strat->sl; ++i) { 2140 Q.sig = p_Copy(strat->P.sig,currRing); 2141 p_ExpVectorAdd(Q.sig,strat->S[i],currRing); 2142 poly help = p_Copy(strat->sig[i],currRing); 2143 p_ExpVectorAdd(help,strat->P.p,currRing); 2144 Q.sig = p_Add_q(Q.sig,help,currRing); 2145 //pWrite(Q.sig); 2146 pos = posInSyz(strat, Q.sig); 2147 enterSyz(Q, strat, pos); 2148 } 2149 strat->currIdx = idx; 2150 } else { 2151 // if the element is not the first one in the given index we build all 2152 // possible syzygies with elements of higher index 2153 for (int i=cmp+1; i<=max_cmp; ++i) { 2154 pos = -1; 2155 for (int j=0; j<strat->sl; ++j) { 2156 if (p_GetComp(strat->sig[j],currRing) == i) { 2157 pos = j; 2158 break; 2159 } 2160 } 2161 if (pos != -1) { 2162 Q.sig = p_One(currRing); 2163 p_SetExpV(Q.sig, vv, currRing); 2164 // F->m[i-1] corresponds to index i 2165 p_ExpVectorAdd(Q.sig,F->m[i-1],currRing); 2166 p_SetComp(Q.sig, i, currRing); 2167 poly help = p_Copy(strat->P.sig,currRing); 2168 p_ExpVectorAdd(help,strat->S[pos],currRing); 2169 Q.sig = p_Add_q(Q.sig,help,currRing); 2170 if (strat->sbaOrder == 0) { 2171 if (p_LmCmp(Q.sig,strat->syz[strat->syzl-1],currRing) == -currRing->OrdSgn) { 2172 pos = posInSyz(strat, Q.sig); 2173 enterSyz(Q, strat, pos); 2174 } 2175 } else { 2176 pos = posInSyz(strat, Q.sig); 2177 enterSyz(Q, strat, pos); 2178 } 2179 } 2180 } 2181 //printf("++ -- done adding syzygies -- ++\n"); 2182 } 2183 } 1903 2184 //#if 1 1904 2185 #if DEBUGF50 … … 1945 2226 // pair was not detected by the rewritten criterion in strat->red = redSig 1946 2227 if (red_result!=2) { 2228 #if SBA_PRINT_ZERO_REDUCTIONS 1947 2229 zeroreductions++; 1948 enterSyz(strat->P,strat); 2230 #endif 2231 int pos = posInSyz(strat, strat->P.sig); 2232 enterSyz(strat->P, strat, pos); 1949 2233 //#if 1 1950 2234 #ifdef DEBUGF5 1951 2235 Print("ADDING STUFF TO SYZ : "); 1952 pWrite(strat->P.p);2236 //pWrite(strat->P.p); 1953 2237 pWrite(strat->P.sig); 1954 2238 #endif … … 2012 2296 #endif 2013 2297 #if SBA_PRINT_SIZE_SYZ 2014 size_syz = strat->syzl+1;2015 #endif 2016 2298 // that is correct, syzl is counting one too far 2299 size_syz = strat->syzl; 2300 #endif 2017 2301 exitSba(strat); 2018 2302 // if (TEST_OPT_WEIGHTM) … … 2043 2327 } 2044 2328 #endif 2045 if ( strat->incremental&& sRing!=currRingOld)2329 if ((strat->sbaOrder == 1 || strat->sbaOrder == 3) && sRing!=currRingOld) 2046 2330 { 2047 2331 rChangeCurrRing (currRingOld); 2048 F0 = idrMoveR (F , sRing, currRing);2332 F0 = idrMoveR (F1, sRing, currRing); 2049 2333 strat->Shdl = idrMoveR_NoSort (strat->Shdl, sRing, currRing); 2050 2334 rDelete (sRing); … … 2063 2347 #endif 2064 2348 #if SBA_PRINT_ZERO_REDUCTIONS 2065 printf("ZERO REDUCTIONS: %ld\n",zeroreductions); 2349 printf("----------------------------------------------------------\n"); 2350 printf("ZERO REDUCTIONS: %ld\n",zeroreductions); 2351 zeroreductions = 0; 2066 2352 #endif 2067 2353 #if SBA_PRINT_REDUCTION_STEPS 2068 printf("TOP S-REDUCTIONS: %ld\n",sba_reduction_steps); 2354 printf("----------------------------------------------------------\n"); 2355 printf("S-REDUCTIONS: %ld\n",sba_reduction_steps); 2356 #endif 2357 #if SBA_PRINT_OPERATIONS 2358 printf("OPERATIONS: %ld\n",sba_operations); 2359 #endif 2360 #if SBA_PRINT_REDUCTION_STEPS 2361 printf("- - - - - - - - - - - - - - - - - - - - - - - - - - - - - \n"); 2362 printf("INTERREDUCTIONS: %ld\n",sba_interreduction_steps); 2363 #endif 2364 #if SBA_PRINT_OPERATIONS 2365 printf("INTERREDUCTION OPERATIONS: %ld\n",sba_interreduction_operations); 2366 #endif 2367 #if SBA_PRINT_REDUCTION_STEPS 2368 printf("- - - - - - - - - - - - - - - - - - - - - - - - - - - - - \n"); 2369 printf("ALL REDUCTIONS: %ld\n",sba_reduction_steps+sba_interreduction_steps); 2370 sba_interreduction_steps = 0; 2371 sba_reduction_steps = 0; 2372 #endif 2373 #if SBA_PRINT_OPERATIONS 2374 printf("ALL OPERATIONS: %ld\n",sba_operations+sba_interreduction_operations); 2375 sba_interreduction_operations = 0; 2376 sba_operations = 0; 2069 2377 #endif 2070 2378 #if SBA_PRINT_SIZE_G 2071 printf("SIZE OF G: %ld\n",size_g); 2379 printf("----------------------------------------------------------\n"); 2380 printf("SIZE OF G: %ld\n",size_g); 2381 size_g = 0; 2072 2382 #endif 2073 2383 #if SBA_PRINT_SIZE_SYZ 2074 printf("SIZE OF SYZ: %ld\n",size_syz); 2384 printf("SIZE OF SYZ: %ld\n",size_syz); 2385 printf("----------------------------------------------------------\n"); 2386 size_syz = 0; 2075 2387 #endif 2076 2388 #if SBA_PRINT_PRODUCT_CRITERION 2077 printf("PRODUCT CRITERIA: %ld\n",product_criterion); 2078 #endif 2079 zeroreductions = 0; 2080 size_g = 0; 2081 size_syz = 0; 2082 product_criterion = 0; 2083 sba_reduction_steps = 0; 2389 printf("PRODUCT CRITERIA: %ld\n",product_criterion); 2390 product_criterion = 0; 2391 #endif 2084 2392 return (strat->Shdl); 2085 2393 } … … 2404 2712 // therefore we call pContent instead of pNorm 2405 2713 #if F5CTAILRED 2714 BOOLEAN withT = TRUE; 2406 2715 if ((TEST_OPT_INTSTRATEGY) || (rField_is_Ring(currRing))) 2407 2716 { -
kernel/kutil.cc
r74c4462 r6e3023a 987 987 memmove(&(strat->sevSig[i]),&(strat->sevSig[i+1]),(strat->sl - i)*sizeof(unsigned long)); 988 988 memmove(&(strat->S_2_R[i]),&(strat->S_2_R[i+1]),(strat->sl - i)*sizeof(int)); 989 memmove(&(strat->fromS[i]),&(strat->fromS[i+1]),(strat->sl - i)*sizeof(int));990 989 #else 991 990 int j; … … 998 997 strat->sevSig[j] = strat->sevSig[j+1]; 999 998 strat->S_2_R[j] = strat->S_2_R[j+1]; 1000 strat->fromS[j] = strat->fromS[j+1];1001 999 } 1002 1000 #endif … … 1045 1043 #endif 1046 1044 pLmFree(set[j].lcm); 1045 } 1046 if (set[j].sig!=NULL) 1047 { 1048 #ifdef HAVE_RINGS 1049 if (pGetCoeff(set[j].sig) != NULL) 1050 pLmDelete(set[j].sig); 1051 else 1052 #endif 1053 pLmFree(set[j].sig); 1047 1054 } 1048 1055 if (set[j].p!=NULL) … … 1779 1786 // the corresponding signatures for criteria checks 1780 1787 LObject Lp; 1781 // poly last;1782 1788 poly pSigMult = p_Copy(pSig,currRing); 1783 1789 poly sSigMult = p_Copy(strat->sig[i],currRing); … … 1806 1812 Print("P1 "); 1807 1813 pWrite(pHead(p)); 1808 Print("FROM: %d\n", from);1809 1814 Print("P2 "); 1810 1815 pWrite(pHead(strat->S[i])); 1811 Print("FROM: %d\n", strat->fromS[i]);1812 1816 Print("M1 "); 1813 1817 pWrite(m1); … … 1827 1831 pWrite(sSigMult); 1828 1832 Print("----------------\n"); 1829 #endif 1830 // testing by syzCrit = F5 Criterion 1831 // testing by rewCrit1 = Rewritten Criterion 1832 if ( strat->syzCrit(pSigMult,pSigMultNegSev,strat) || 1833 strat->syzCrit(sSigMult,sSigMultNegSev,strat) 1834 || strat->rewCrit1(sSigMult,sSigMultNegSev,strat,i+1) 1835 ) 1836 { 1837 pDelete(&pSigMult); 1838 pDelete(&sSigMult); 1839 strat->cp++; 1840 pLmFree(Lp.lcm); 1841 Lp.lcm=NULL; 1842 pDelete (&m1); 1843 pDelete (&m2); 1844 return; 1845 } 1846 // in any case Lp is checked up to the next strat->P which is added 1847 // to S right after this critical pair creation. 1848 // NOTE: this even holds if the 2nd generator gives the bigger signature 1849 // moreover, this improves rewCriterion, 1850 // i.e. strat->checked > strat->from if and only if the 2nd generator 1851 // gives the bigger signature. 1852 Lp.checked = strat->sl+1; 1833 Lp.checked = 0; 1834 #endif 1853 1835 int sigCmp = p_LmCmp(pSigMult,sSigMult,currRing); 1854 1836 //#if 1 … … 1862 1844 // printf("!!!! EQUAL SIGS !!!!\n"); 1863 1845 // pSig = sSig, delete element due to Rewritten Criterion 1864 strat->cp++;1865 1846 pDelete(&pSigMult); 1866 1847 pDelete(&sSigMult); … … 1871 1852 return; 1872 1853 } 1873 // at this point it is clear that the pair will be added to L, since it has 1874 // passed all tests up to now 1875 1876 // store from which element this pair comes from for further tests 1877 Lp.from = strat->sl+1; 1878 if(sigCmp==currRing->OrdSgn) 1879 { 1880 // pSig > sSig 1881 pDelete (&sSigMult); 1882 Lp.sig = pSigMult; 1883 Lp.sevSig = ~pSigMultNegSev; 1884 } 1885 else 1886 { 1887 // pSig < sSig 1888 pDelete (&pSigMult); 1889 Lp.sig = sSigMult; 1890 Lp.sevSig = ~sSigMultNegSev; 1891 } 1892 // adds buchberger's first criterion 1893 if (pLmCmp(m2,pHead(p)) == 0) { 1894 Lp.checked = 3; // 3 == Product Criterion 1895 #if 0 1896 enterSyz(Lp, strat); 1854 // testing by syzCrit = F5 Criterion 1855 // testing by rewCrit1 = Rewritten Criterion 1856 // NOTE: Arri's Rewritten Criterion is tested below, we need Lp.p for it! 1857 if ( strat->syzCrit(pSigMult,pSigMultNegSev,strat) || 1858 strat->syzCrit(sSigMult,sSigMultNegSev,strat) 1859 || strat->rewCrit1(sSigMult,sSigMultNegSev,Lp.lcm,strat,i+1) 1860 ) 1861 { 1862 pDelete(&pSigMult); 1863 pDelete(&sSigMult); 1864 pLmFree(Lp.lcm); 1897 1865 Lp.lcm=NULL; 1898 1866 pDelete (&m1); 1899 1867 pDelete (&m2); 1900 1868 return; 1901 #endif 1902 } 1903 pDelete (&m1); 1904 pDelete (&m2); 1905 #if DEBUGF5 1906 printf("SIGNATURE OF PAIR: "); 1907 pWrite(Lp.sig); 1908 #endif 1869 } 1909 1870 /* 1910 1871 *the pair (S[i],p) enters B if the spoly != 0 … … 1985 1946 } 1986 1947 } 1948 // store from which element this pair comes from for further tests 1949 //Lp.from = strat->sl+1; 1950 if(sigCmp==currRing->OrdSgn) 1951 { 1952 // pSig > sSig 1953 pDelete (&sSigMult); 1954 Lp.sig = pSigMult; 1955 Lp.sevSig = ~pSigMultNegSev; 1956 } 1957 else 1958 { 1959 // pSig < sSig 1960 pDelete (&pSigMult); 1961 Lp.sig = sSigMult; 1962 Lp.sevSig = ~sSigMultNegSev; 1963 } 1987 1964 if (Lp.p == NULL) 1988 1965 { 1989 /*- the case that the s-poly is 0 -*/1990 if (strat->pairtest==NULL) initPairtest(strat);1991 strat->pairtest[i] = TRUE;/*- hint for spoly(S^[i],p)=0 -*/1992 strat->pairtest[strat->sl+1] = TRUE;1993 /*hint for spoly(S[i],p) == 0 for some i,0 <= i <= sl*/1994 /*1995 *suppose we have (s,r),(r,p),(s,p) and spoly(s,p) == 0 and (r,p) is1996 *still in B (i.e. lcm(r,p) == lcm(s,p) or the leading term of s does not1997 *devide lcm(r,p)). In the last case (s,r) can be canceled if the leading1998 *term of p devides the lcm(s,r)1999 *(this canceling should be done here because2000 *the case lcm(s,p) == lcm(s,r) is not covered in chainCrit)2001 *the first case is handeled in chainCrit2002 */2003 1966 if (Lp.lcm!=NULL) pLmFree(Lp.lcm); 1967 int pos = posInSyz(strat, Lp.sig); 1968 enterSyz(Lp, strat, pos); 2004 1969 } 2005 1970 else 2006 1971 { 1972 // testing by rewCrit3 = Arris Rewritten Criterion (for F5 nothing happens!) 1973 if (strat->rewCrit3(Lp.sig,~Lp.sevSig,Lp.p,strat,strat->sl+1)) { 1974 pLmFree(Lp.lcm); 1975 pDelete(&Lp.sig); 1976 Lp.lcm=NULL; 1977 pDelete (&m1); 1978 pDelete (&m2); 1979 return; 1980 } 1981 // in any case Lp is checked up to the next strat->P which is added 1982 // to S right after this critical pair creation. 1983 // NOTE: this even holds if the 2nd generator gives the bigger signature 1984 // moreover, this improves rewCriterion, 1985 // i.e. strat->checked > strat->from if and only if the 2nd generator 1986 // gives the bigger signature. 1987 Lp.checked = strat->sl+1; 1988 // at this point it is clear that the pair will be added to L, since it has 1989 // passed all tests up to now 1990 1991 // adds buchberger's first criterion 1992 if (pLmCmp(m2,pHead(p)) == 0) { 1993 Lp.prod_crit = TRUE; // Product Criterion 1994 #if 0 1995 int pos = posInSyz(strat, Lp.sig); 1996 enterSyz(Lp, strat, pos); 1997 Lp.lcm=NULL; 1998 pDelete (&m1); 1999 pDelete (&m2); 2000 return; 2001 #endif 2002 } 2003 pDelete (&m1); 2004 pDelete (&m2); 2005 #if DEBUGF5 2006 printf("SIGNATURE OF PAIR: "); 2007 pWrite(Lp.sig); 2008 #endif 2007 2009 /*- the pair (S[i],p) enters B -*/ 2008 2010 Lp.p1 = strat->S[i]; … … 4631 4633 } 4632 4634 4635 // for sba, sorting syzygies 4636 int posInSyz (const kStrategy strat, poly sig) 4637 { 4638 if (strat->syzl==0) return 0; 4639 if (pLmCmp(strat->syz[strat->syzl-1],sig) != currRing->OrdSgn) 4640 return strat->syzl; 4641 int i; 4642 int an = 0; 4643 int en= strat->syzl-1; 4644 loop 4645 { 4646 if (an >= en-1) 4647 { 4648 if (pLmCmp(strat->syz[an],sig) != currRing->OrdSgn) return en; 4649 return an; 4650 } 4651 i=(an+en) / 2; 4652 if (pLmCmp(strat->syz[i],sig) != currRing->OrdSgn) an=i; 4653 else en=i; 4654 /*aend. fuer lazy == in !=- machen */ 4655 } 4656 } 4657 4633 4658 /*2 4634 4659 * … … 5065 5090 for (int k=0; k<strat->syzl; k++) 5066 5091 { 5092 //printf("-%d",k); 5067 5093 //#if 1 5068 5094 #ifdef DEBUGF5 5069 Print("checking with: %d -- ",k);5095 Print("checking with: %d / %d -- \n",k,strat->syzl); 5070 5096 pWrite(pHead(strat->syz[k])); 5071 5097 #endif … … 5076 5102 printf("DELETE!\n"); 5077 5103 #endif 5104 //printf("- T -\n\n"); 5078 5105 return TRUE; 5079 5106 } 5080 5107 } 5108 //printf("- F -\n\n"); 5081 5109 return FALSE; 5082 5110 } … … 5089 5117 //#if 1 5090 5118 #ifdef DEBUGF5 5091 Print(" syzygy criterion checks: ");5119 Print("--- syzygy criterion checks: "); 5092 5120 pWrite(sig); 5093 5121 #endif … … 5112 5140 for (int k=min; k<max; k++) 5113 5141 { 5114 #ifdef DEBUGF55142 #ifdef F5DEBUG 5115 5143 printf("COMP %d/%d - MIN %d - MAX %d - SYZL %ld\n",comp,strat->currIdx,min,max,strat->syzl); 5116 5144 Print("checking with: %d -- ",k); … … 5127 5155 * REWRITTEN CRITERION for signature-based standard basis algorithms 5128 5156 */ 5129 BOOLEAN faugereRewCriterion(poly sig, unsigned long not_sevSig, kStrategy strat, int start=0)5157 BOOLEAN faugereRewCriterion(poly sig, unsigned long not_sevSig, poly /*lm*/, kStrategy strat, int start=0) 5130 5158 { 5131 5159 //printf("Faugere Rewritten Criterion\n"); … … 5135 5163 pWrite(sig); 5136 5164 #endif 5137 //for(int k = start; k<strat->sl+1; k++)5138 5165 for(int k = strat->sl; k>start; k--) 5139 5166 { … … 5145 5172 #endif 5146 5173 if (p_LmShortDivisibleBy(strat->sig[k], strat->sevSig[k], sig, not_sevSig, currRing)) 5147 //if (p_LmEqual(strat->sig[k], sig, currRing))5148 5174 { 5149 5175 //#if 1 … … 5153 5179 return TRUE; 5154 5180 } 5181 //k--; 5155 5182 } 5156 5183 #ifdef DEBUGF5 … … 5184 5211 // critical pair. In this situation we can discard the critical pair 5185 5212 // completely. 5186 BOOLEAN arriRewCriterion(poly /*sig*/, unsigned long /*not_sevSig*/, kStrategy strat, int /*start=0*/) 5187 { 5188 //printf("Arri Rewritten Criterion\n"); 5189 while (strat->Ll > 0 && pLmEqual(strat->L[strat->Ll].sig,strat->P.sig)) 5190 { 5191 // deletes the short spoly 5192 #ifdef HAVE_RINGS 5193 if (rField_is_Ring(currRing)) 5194 pLmDelete(strat->L[strat->Ll].p); 5195 else 5196 #endif 5197 pLmFree(strat->L[strat->Ll].p); 5198 5199 // TODO: needs some masking 5200 // TODO: masking needs to vanish once the signature 5201 // sutff is completely implemented 5202 strat->L[strat->Ll].p = NULL; 5203 poly m1 = NULL, m2 = NULL; 5204 5205 // check that spoly creation is ok 5206 while (strat->tailRing != currRing && 5207 !kCheckSpolyCreation(&(strat->L[strat->Ll]), strat, m1, m2)) 5208 { 5209 assume(m1 == NULL && m2 == NULL); 5210 // if not, change to a ring where exponents are at least 5211 // large enough 5212 if (!kStratChangeTailRing(strat)) 5213 { 5214 WerrorS("OVERFLOW..."); 5215 break; 5216 } 5217 } 5218 // create the real one 5219 ksCreateSpoly(&(strat->L[strat->Ll]), NULL, strat->use_buckets, 5220 strat->tailRing, m1, m2, strat->R); 5221 if (strat->P.GetLmCurrRing() == NULL) 5222 { 5223 deleteInL(strat->L,&strat->Ll,strat->Ll,strat); 5224 } 5225 if (strat->L[strat->Ll].GetLmCurrRing() == NULL) 5226 { 5227 strat->P.Delete(); 5228 strat->P = strat->L[strat->Ll]; 5229 strat->Ll--; 5230 } 5231 5232 if ((strat->P.GetLmCurrRing() != NULL) 5233 && (strat->L[strat->Ll].GetLmCurrRing() != NULL)) 5234 { 5235 if (pLmCmp(strat->P.GetLmCurrRing(),strat->L[strat->Ll].GetLmCurrRing()) == -1) 5236 { 5237 deleteInL(strat->L,&strat->Ll,strat->Ll,strat); 5238 } 5239 else 5240 { 5241 strat->P.Delete(); 5242 strat->P = strat->L[strat->Ll]; 5243 strat->Ll--; 5244 } 5245 } 5246 } 5213 BOOLEAN arriRewCriterion(poly /*sig*/, unsigned long /*not_sevSig*/, poly /*lm*/, kStrategy strat, int start=0) 5214 { 5215 poly p1 = pOne(); 5216 poly p2 = pOne(); 5217 for (int ii=strat->sl; ii>start; ii--) 5218 { 5219 if (p_LmShortDivisibleBy(strat->sig[ii], strat->sevSig[ii], strat->P.sig, ~strat->P.sevSig, currRing)) 5220 { 5221 p_ExpVectorSum(p1,strat->P.sig,strat->S[ii],currRing); 5222 p_ExpVectorSum(p2,strat->sig[ii],strat->P.p,currRing); 5223 if (!(pLmCmp(p1,p2) == 1)) 5224 { 5225 pDelete(&p1); 5226 pDelete(&p2); 5227 return TRUE; 5228 } 5229 } 5230 } 5231 pDelete(&p1); 5232 pDelete(&p2); 5233 return FALSE; 5234 } 5235 5236 BOOLEAN arriRewCriterionPre(poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int /*start=0*/) 5237 { 5238 int found = -1; 5239 for (int i=strat->Bl; i>-1; i--) { 5240 if (pLmEqual(strat->B[i].sig,sig)) { 5241 found = i; 5242 break; 5243 } 5244 } 5245 if (found != -1) { 5246 if (pLmCmp(lm,strat->B[found].GetLmCurrRing()) == -1) { 5247 deleteInL(strat->B,&strat->Bl,found,strat); 5248 } else { 5249 return TRUE; 5250 } 5251 } 5252 poly p1 = pOne(); 5253 poly p2 = pOne(); 5247 5254 for (int ii=strat->sl; ii>-1; ii--) 5248 5255 { 5249 if (p_LmShortDivisibleBy(strat->sig[ii], strat->sevSig[ii], strat->P.sig, ~strat->P.sevSig, currRing)) 5250 { 5251 if (!(pLmCmp(ppMult_mm(strat->P.sig,pHead(strat->S[ii])),ppMult_mm(strat->sig[ii],strat->P.GetLmCurrRing())) == 1)) 5252 { 5253 strat->P.Delete(); 5256 if (p_LmShortDivisibleBy(strat->sig[ii], strat->sevSig[ii], sig, not_sevSig, currRing)) 5257 { 5258 p_ExpVectorSum(p1,sig,strat->S[ii],currRing); 5259 p_ExpVectorSum(p2,strat->sig[ii],lm,currRing); 5260 if (!(pLmCmp(p1,p2) == 1)) 5261 { 5262 pDelete(&p1); 5263 pDelete(&p2); 5254 5264 return TRUE; 5255 5265 } 5256 5266 } 5257 5267 } 5268 pDelete(&p1); 5269 pDelete(&p2); 5258 5270 return FALSE; 5259 5271 } … … 5929 5941 else i=setmaxT; 5930 5942 strat->ecartS = initec(i); 5931 strat->fromS = initec(i);5932 5943 strat->sevS = initsevS(i); 5933 5944 strat->sevSig = initsevS(i); … … 5937 5948 strat->S = strat->Shdl->m; 5938 5949 strat->sig = (poly *)omAlloc0(i*sizeof(poly)); 5939 if ( !strat->incremental)5950 if (strat->sbaOrder != 1) 5940 5951 { 5941 5952 strat->syz = (poly *)omAlloc0(i*sizeof(poly)); … … 5999 6010 // => we can keep the underlying monomial order and get a Schreyer 6000 6011 // order without any bigger overhead 6001 if ( !strat->incremental)6012 if (strat->sbaOrder == 0 || strat->sbaOrder == 3) 6002 6013 { 6003 6014 p_ExpVectorAdd (h.sig,F->m[i],currRing); … … 6036 6047 } 6037 6048 /* 6038 if ( !strat->incremental)6049 if (strat->sbaOrder != 1) 6039 6050 { 6040 6051 for(j=0;j<i;j++) … … 6103 6114 strat->sevSyz = initsevS(ps); 6104 6115 strat->syz = (poly *)omAlloc(ps*sizeof(poly)); 6105 strat->syzl = strat->syzmax = ps; 6116 strat->syzmax = ps; 6117 strat->syzl = 0; 6106 6118 strat->syzidxmax = comp; 6107 6119 #if defined(DEBUGF5) || defined(DEBUGF51) … … 6141 6153 strat->syzIdx[j] = ctr; 6142 6154 j++; 6155 LObject Q; 6156 int pos; 6143 6157 for (k = 0; k<i; k++) 6144 6158 { 6145 poly p= pOne();6146 p Lcm(strat->S[k],strat->S[i],p);6147 strat->syz[ctr] = p;6148 p _SetCompP (strat->syz[ctr], comp,currRing);6149 p oly q = p_Copy(p,currRing);6159 Q.sig = pOne(); 6160 p_ExpVectorCopy(Q.sig,strat->S[k],currRing); 6161 p_SetCompP (Q.sig, comp, currRing); 6162 poly q = p_One(currRing); 6163 p_ExpVectorCopy(q,strat->S[i],currRing); 6150 6164 q = p_Neg (q, currRing); 6151 6165 p_SetCompP (q, p_GetComp(strat->sig[k], currRing), currRing); 6152 strat->syz[ctr] = p_Add_q (strat->syz[ctr], q, currRing); 6153 #if defined(DEBUGF5) || defined(DEBUGF51) 6154 pWrite(strat->syz[ctr]); 6155 #endif 6156 strat->sevSyz[ctr] = p_GetShortExpVector(strat->syz[ctr],currRing); 6166 Q.sig = p_Add_q (Q.sig, q, currRing); 6167 Q.sevSig = p_GetShortExpVector(Q.sig,currRing); 6168 pos = posInSyz(strat, Q.sig); 6169 enterSyz(Q, strat, pos); 6157 6170 ctr++; 6158 6171 } … … 6180 6193 } 6181 6194 strat->syzIdx[j] = ctr; 6195 LObject Q; 6196 int pos; 6182 6197 for (k = 0; k<strat->sl+1; k++) 6183 6198 { 6184 strat->syz[ctr] = p_Copy (pHead(strat->S[k]), currRing); 6185 p_SetCompP (strat->syz[ctr], comp, currRing); 6186 poly q = p_Copy (pHead(strat->L[strat->Ll].p), currRing); 6199 Q.sig = pOne(); 6200 p_ExpVectorCopy(Q.sig,strat->S[k],currRing); 6201 p_SetCompP (Q.sig, comp, currRing); 6202 poly q = p_One(currRing); 6203 p_ExpVectorCopy(q,strat->L[strat->Ll].p,currRing); 6187 6204 q = p_Neg (q, currRing); 6188 6205 p_SetCompP (q, p_GetComp(strat->sig[k], currRing), currRing); 6189 strat->syz[ctr] = p_Add_q (strat->syz[ctr], q, currRing); 6190 //#if 1 6191 #if DEBUGF5 || DEBUGF51 6192 printf(".."); 6193 pWrite(strat->syz[ctr]); 6194 #endif 6195 strat->sevSyz[ctr] = p_GetShortExpVector(strat->syz[ctr],currRing); 6206 Q.sig = p_Add_q (Q.sig, q, currRing); 6207 Q.sevSig = p_GetShortExpVector(Q.sig,currRing); 6208 pos = posInSyz(strat, Q.sig); 6209 enterSyz(Q, strat, pos); 6196 6210 ctr++; 6197 6211 } … … 6199 6213 #ifdef DEBUGF5 6200 6214 Print("Principal syzygies:\n"); 6215 printf("syzl %d\n",strat->syzl); 6216 printf("syzmax %d\n",strat->syzmax); 6217 printf("ps %d\n",ps); 6201 6218 Print("--------------------------------\n"); 6202 for(i=0;i<=ps-1;i++) 6203 { 6219 for(i=0;i<=strat->syzl-1;i++) 6220 { 6221 printf("%d - ",i); 6204 6222 pWrite(strat->syz[i]); 6223 } 6224 for(i=0;i<strat->currIdx;i++) 6225 { 6226 printf("%d - %d\n",i,strat->syzIdx[i]); 6205 6227 } 6206 6228 Print("--------------------------------\n"); … … 6367 6389 else i=setmaxT; 6368 6390 i=((i+IDELEMS(F)+IDELEMS(P)+15)/16)*16; 6369 strat->fromS=initec(i);6370 6391 strat->sevS=initsevS(i); 6371 6392 strat->sevSig=initsevS(i); … … 7024 7045 (IDELEMS(strat->Shdl)+setmaxTinc) 7025 7046 *sizeof(int)); 7026 strat->fromS = (intset)omReallocSize(strat->fromS,7027 IDELEMS(strat->Shdl)*sizeof(int),7028 (IDELEMS(strat->Shdl)+setmaxTinc)7029 *sizeof(int));7030 7047 strat->S_2_R = (int*) omRealloc0Size(strat->S_2_R, 7031 7048 IDELEMS(strat->Shdl)*sizeof(int), … … 7064 7081 memmove(&(strat->ecartS[atS+1]), &(strat->ecartS[atS]), 7065 7082 (strat->sl - atS + 1)*sizeof(int)); 7066 memmove(&(strat->fromS[atS+1]), &(strat->fromS[atS]),7067 (strat->sl - atS + 1)*sizeof(int));7068 7083 memmove(&(strat->sevS[atS+1]), &(strat->sevS[atS]), 7069 7084 (strat->sl - atS + 1)*sizeof(unsigned long)); … … 7081 7096 strat->S[i] = strat->S[i-1]; 7082 7097 strat->ecartS[i] = strat->ecartS[i-1]; 7083 strat->fromS[i] = strat->fromS[i-1];7084 7098 strat->sevS[i] = strat->sevS[i-1]; 7085 7099 strat->S_2_R[i] = strat->S_2_R[i-1]; … … 7128 7142 } 7129 7143 strat->ecartS[atS] = p.ecart; 7130 strat->fromS[atS] = p.from;7131 7144 strat->S_2_R[atS] = atR; 7132 7145 strat->sl++; … … 7214 7227 } 7215 7228 7229 7216 7230 /*2 7217 7231 * puts signature p.sig to the set syz 7218 7232 */ 7219 void enterSyz(LObject p, kStrategy strat) 7220 { 7221 int i = strat->syzl; 7222 7233 void enterSyz(LObject p, kStrategy strat, int atT) 7234 { 7235 int i; 7223 7236 strat->newt = TRUE; 7224 if (strat->syzl == strat->syzmax )7237 if (strat->syzl == strat->syzmax-1) 7225 7238 { 7226 7239 pEnlargeSet(&strat->syz,strat->syzmax,setmaxTinc); … … 7231 7244 strat->syzmax += setmaxTinc; 7232 7245 } 7233 strat->syz[i] = p.sig; 7234 strat->sevSyz[i] = p.sevSig; 7246 if (atT < strat->syzl) 7247 { 7248 #ifdef ENTER_USE_MEMMOVE 7249 memmove(&(strat->syz[atT+1]), &(strat->syz[atT]), 7250 (strat->syzl-atT+1)*sizeof(poly)); 7251 memmove(&(strat->sevSyz[atT+1]), &(strat->sevSyz[atT]), 7252 (strat->syzl-atT+1)*sizeof(unsigned long)); 7253 #endif 7254 for (i=strat->syzl; i>=atT+1; i--) 7255 { 7256 #ifndef ENTER_USE_MEMMOVE 7257 strat->syz[i] = strat->syz[i-1]; 7258 strat->sevSyz[i] = strat->sevSyz[i-1]; 7259 #endif 7260 } 7261 } 7262 //i = strat->syzl; 7263 i = atT; 7264 strat->syz[atT] = p.sig; 7265 strat->sevSyz[atT] = p.sevSig; 7235 7266 strat->syzl++; 7236 #if def DEBUGF57237 Print(" last element in strat->syz: %d--%d ",i+1,strat->syzmax);7238 pWrite(strat->syz[ i]);7267 #if F5DEBUG 7268 Print("element in strat->syz: %d--%d ",atT+1,strat->syzmax); 7269 pWrite(strat->syz[atT]); 7239 7270 #endif 7240 7271 // recheck pairs in strat->L with new rule and delete correspondingly … … 7242 7273 while (cc>-1) 7243 7274 { 7244 if (p_LmShortDivisibleBy( strat->syz[ strat->syzl-1], strat->sevSyz[strat->syzl-1],7275 if (p_LmShortDivisibleBy( strat->syz[atT], strat->sevSyz[atT], 7245 7276 strat->L[cc].sig, ~strat->L[cc].sevSig, currRing)) 7246 7277 { … … 7249 7280 cc--; 7250 7281 } 7251 7282 //#if 1 7283 #ifdef DEBUGF5 7284 Print("--- Syzygies ---\n"); 7285 printf("syzl %d\n",strat->syzl); 7286 printf("syzmax %d\n",strat->syzmax); 7287 Print("--------------------------------\n"); 7288 for(i=0;i<=strat->syzl-1;i++) 7289 { 7290 printf("%d - ",i); 7291 pWrite(strat->syz[i]); 7292 } 7293 Print("--------------------------------\n"); 7294 #endif 7252 7295 } 7253 7296 … … 7341 7384 *****************************************/ 7342 7385 //strat->rewCrit1 = faugereRewCriterion; 7343 if (strat-> incremental)7386 if (strat->sbaOrder == 1) 7344 7387 { 7345 7388 strat->syzCrit = syzCriterionInc; … … 7760 7803 omFreeSize(strat->sevT, (strat->tmax)*sizeof(unsigned long)); 7761 7804 omFreeSize(strat->ecartS,IDELEMS(strat->Shdl)*sizeof(int)); 7762 omFreeSize(strat->fromS,IDELEMS(strat->Shdl)*sizeof(int));7763 7805 omFreeSize((ADDRESS)strat->sevS,IDELEMS(strat->Shdl)*sizeof(unsigned long)); 7764 7806 omFreeSize((ADDRESS)strat->sevSig,IDELEMS(strat->Shdl)*sizeof(unsigned long)); 7765 7807 omFreeSize((ADDRESS)strat->syz,(strat->syzmax)*sizeof(poly)); 7766 7808 omFreeSize((ADDRESS)strat->sevSyz,(strat->syzmax)*sizeof(unsigned long)); 7767 if (strat-> incremental)7809 if (strat->sbaOrder == 1) 7768 7810 { 7769 7811 omFreeSize(strat->syzIdx,(strat->syzidxmax)*sizeof(int)); … … 8254 8296 { 8255 8297 int n = rBlocks(r); // Including trailing zero! 8256 // if incremental=> use (C,monomial order from r)8257 if (strat-> incremental)8298 // if sbaOrder == 1 => use (C,monomial order from r) 8299 if (strat->sbaOrder == 1) 8258 8300 { 8259 8301 if (r->order[0] == ringorder_C || r->order[0] == ringorder_c) … … 8261 8303 return r; 8262 8304 } 8263 ring res = rCopy0(r, FALSE, TRUE); 8264 for (int i=1; i<n-1; i++) 8265 { 8266 res->order[i] = res->order[i-1]; 8267 res->block0[i] = res->block0[i-1]; 8268 res->block1[i] = res->block1[i-1]; 8269 res->wvhdl[i] = res->wvhdl[i-1]; 8305 ring res = rCopy0(r, TRUE, FALSE); 8306 res->order = (int *)omAlloc0((n+1)*sizeof(int)); 8307 res->block0 = (int *)omAlloc0((n+1)*sizeof(int)); 8308 res->block1 = (int *)omAlloc0((n+1)*sizeof(int)); 8309 int **wvhdl = (int **)omAlloc0((n+1)*sizeof(int*)); 8310 res->wvhdl = wvhdl; 8311 for (int i=1; i<n; i++) 8312 { 8313 res->order[i] = r->order[i-1]; 8314 res->block0[i] = r->block0[i-1]; 8315 res->block1[i] = r->block1[i-1]; 8316 res->wvhdl[i] = r->wvhdl[i-1]; 8270 8317 } 8271 8318 8272 8319 // new 1st block 8273 8320 res->order[0] = ringorder_C; // Prefix 8274 res->block0[0] = 1; 8275 res->block1[0] = res->N; 8276 //res->wvhdl[j] = NULL; 8277 // res->order [j] = 0; // The End! 8321 // removes useless secondary component order if defined in old ring 8322 for (int i=rBlocks(res); i>0; --i) { 8323 if (res->order[i] == ringorder_C || res->order[i] == ringorder_c) { 8324 res->order[i] = 0; 8325 break; 8326 } 8327 } 8278 8328 rComplete(res, 1); 8279 8329 #ifdef HAVE_PLURAL … … 8297 8347 return (res); 8298 8348 } 8299 8300 // not incremental => use Schreyer order 8349 // if sbaOrder == 3 => degree - position - ring order 8350 if (strat->sbaOrder == 3) 8351 { 8352 ring res = rCopy0(r, TRUE, FALSE); 8353 res->order = (int *)omAlloc0((n+2)*sizeof(int)); 8354 res->block0 = (int *)omAlloc0((n+2)*sizeof(int)); 8355 res->block1 = (int *)omAlloc0((n+2)*sizeof(int)); 8356 int **wvhdl = (int **)omAlloc0((n+2)*sizeof(int*)); 8357 res->wvhdl = wvhdl; 8358 for (int i=2; i<n+2; i++) 8359 { 8360 res->order[i] = r->order[i-2]; 8361 res->block0[i] = r->block0[i-2]; 8362 res->block1[i] = r->block1[i-2]; 8363 res->wvhdl[i] = r->wvhdl[i-2]; 8364 } 8365 8366 // new 1st block 8367 res->order[0] = ringorder_a; // Prefix 8368 res->block0[0] = 1; 8369 res->wvhdl[0] = (int *)omAlloc(res->N*sizeof(int)); 8370 for (int i=0; i<res->N; ++i) 8371 res->wvhdl[0][i] = 1; 8372 res->block1[0] = si_min(res->N, rVar(res)); 8373 // new 2nd block 8374 res->order[1] = ringorder_C; // Prefix 8375 res->wvhdl[1] = NULL; 8376 // removes useless secondary component order if defined in old ring 8377 for (int i=rBlocks(res); i>0; --i) { 8378 if (res->order[i] == ringorder_C || res->order[i] == ringorder_c) { 8379 res->order[i] = 0; 8380 break; 8381 } 8382 } 8383 rComplete(res, 1); 8384 #ifdef HAVE_PLURAL 8385 if (rIsPluralRing(r)) 8386 { 8387 if ( nc_rComplete(r, res, false) ) // no qideal! 8388 { 8389 #ifndef NDEBUG 8390 WarnS("error in nc_rComplete"); 8391 #endif 8392 // cleanup? 8393 8394 // rDelete(res); 8395 // return r; 8396 8397 // just go on.. 8398 } 8399 } 8400 #endif 8401 strat->tailRing = res; 8402 return (res); 8403 } 8404 8405 // not sbaOrder == 1 => use Schreyer order 8301 8406 // this is done by a trick when initializing the signatures 8302 8407 // in initSLSba(): -
kernel/kutil.h
r74c4462 r6e3023a 180 180 public: 181 181 unsigned long sev; 182 unsigned long from; // from which polynomial it comes from 183 // this is important for signature-based 184 // algorithms 182 unsigned long from; // index in sig up to which the correspongin LObject was already checked 185 183 unsigned long checked; // this is the index of S up to which 186 184 // the corresponding LObject was already checked in … … 188 186 // reduction process it is enough to start a second 189 187 // rewritten criterion check from checked+1 onwards 190 // NOTE: If checked = 3 then the corresponding pair is 188 BOOLEAN prod_crit; 189 // NOTE: If prod_crit = TRUE then the corresponding pair is 191 190 // detected by Buchberger's Product Criterion and can be 192 191 // deleted … … 289 288 void (*chainCrit) (poly p,int ecart,kStrategy strat); 290 289 BOOLEAN (*syzCrit) (poly sig, unsigned long not_sevSig, kStrategy strat); 291 BOOLEAN (*rewCrit1) (poly sig, unsigned long not_sevSig, kStrategy strat, int start /*= 0*/); 292 BOOLEAN (*rewCrit2) (poly sig, unsigned long not_sevSig, kStrategy strat, int start /*= 0*/); 290 BOOLEAN (*rewCrit1) (poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int start /*= 0*/); 291 BOOLEAN (*rewCrit2) (poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int start /*= 0*/); 292 BOOLEAN (*rewCrit3) (poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int start /*= 0*/); 293 293 pFDegProc pOrigFDeg; 294 294 pLDegProc pOrigLDeg; … … 310 310 // syzygy of component i comes up 311 311 // important for signature-based algorithms 312 BOOLEAN incremental;312 unsigned sbaOrder; 313 313 unsigned long currIdx; 314 314 int max_lower_index; … … 446 446 int posInLSig (const LSet set, const int length, 447 447 LObject* L,const kStrategy strat); 448 int posInSyz (const kStrategy strat, const poly sig); 448 449 int posInL0 (const LSet set, const int length, 449 450 LObject* L,const kStrategy strat); … … 467 468 poly redtailBba (LObject *L, int pos,kStrategy strat, 468 469 BOOLEAN withT = FALSE,BOOLEAN normalize=FALSE); 470 poly redtailSba (LObject *L, int pos,kStrategy strat, 471 BOOLEAN withT = FALSE,BOOLEAN normalize=FALSE); 469 472 poly redtailBba (TObject *T, int pos,kStrategy strat); 470 473 poly redtail (poly p,int pos,kStrategy strat); … … 518 521 void initSyzRules (kStrategy strat); 519 522 void updateS(BOOLEAN toT,kStrategy strat); 520 void enterSyz (LObject p,kStrategy strat );523 void enterSyz (LObject p,kStrategy strat, int atT); 521 524 void enterT (LObject p,kStrategy strat, int atT = -1); 522 525 void cancelunit (LObject* p,BOOLEAN inNF=FALSE); … … 542 545 BOOLEAN syzCriterion(poly sig, unsigned long not_sevSig, kStrategy strat); 543 546 BOOLEAN syzCriterionInc(poly sig, unsigned long not_sevSig, kStrategy strat); 544 KINLINE BOOLEAN arriRewDummy(poly sig, unsigned long not_sevSig, kStrategy strat, int start); 545 BOOLEAN arriRewCriterion(poly sig, unsigned long not_sevSig, kStrategy strat, int start); 546 BOOLEAN faugereRewCriterion(poly sig, unsigned long not_sevSig, kStrategy strat, int start); 547 KINLINE BOOLEAN arriRewDummy(poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int start); 548 BOOLEAN arriRewCriterion(poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int start); 549 BOOLEAN arriRewCriterionPre(poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int start); 550 BOOLEAN faugereRewCriterion(poly sig, unsigned long not_sevSig, poly lm, kStrategy strat, int start); 547 551 BOOLEAN findMinLMPair(poly sig, unsigned long not_sevSig, kStrategy strat, int start); 548 552 // returns index of p in TSet, or -1 if not found
Note: See TracChangeset
for help on using the changeset viewer.