[6fe9a5] | 1 | ////////////////////////////////////////////////////////////////////////////// |
---|
| 2 | version="$Id$"; |
---|
| 3 | category="Noncommutative"; |
---|
| 4 | info=" |
---|
| 5 | LIBRARY: bfun.lib Algorithms for b-functions and Bernstein-Sato polynomial |
---|
| 6 | AUTHORS: Daniel Andres, daniel.andres@math.rwth-aachen.de |
---|
| 7 | @* Viktor Levandovskyy, levandov@math.rwth-aachen.de |
---|
| 8 | |
---|
| 9 | OVERVIEW: |
---|
| 10 | Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R, |
---|
| 11 | one is interested in the global b-function (also known as Bernstein-Sato |
---|
| 12 | polynomial) b(s) in K[s], defined to be the non-zero monic polynomial of minimal |
---|
| 13 | degree, satisfying a functional identity L * F^{s+1} = b(s) F^s, |
---|
| 14 | for some operator L in D[s] (* stands for the action of differential operator)@* |
---|
| 15 | By D one denotes the n-th Weyl algebra |
---|
| 16 | K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>. |
---|
| 17 | One is interested in the following data:@* |
---|
| 18 | - Bernstein-Sato polynomial b(s) in K[s],@* |
---|
| 19 | - the list of its roots, which are known to be rational@* |
---|
| 20 | - the multiplicities of the roots.@* |
---|
| 21 | |
---|
| 22 | There is a constructive definition of a b-function of a holonomic ideal I in D |
---|
| 23 | (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module) |
---|
| 24 | with respect to the given weight vector w: For a polynomial p in D, its initial |
---|
| 25 | form w.r.t. (-w,w) is defined as the sum of all terms of p which have |
---|
| 26 | maximal weighted total degree where the weight of x_i is -w_i and the weight |
---|
| 27 | of d_i is w_i. Let J be the initial ideal of I w.r.t. (-w,w), i.e. the |
---|
| 28 | K-vector space generated by all initial forms w.r.t (-w,w) of elements of I. |
---|
| 29 | Put s = w_1 x_1 d_1 + ... + w_n x_n d_n. Then the monic generator b_w(s) of |
---|
| 30 | the intersection of J with the PID K[s] is called the b-function of I w.r.t. w. |
---|
| 31 | Unlike Bernstein-Sato polynomial, general b-function with respect to |
---|
| 32 | arbitrary weights need not have rational roots at all. However, b-function |
---|
| 33 | of a holonomic ideal is known to be non-zero as well. |
---|
| 34 | |
---|
| 35 | REFERENCES: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of |
---|
| 36 | Hypergeometric Differential Equations (2000),@* |
---|
| 37 | Noro: An Efficient Modular Algorithm for Computing the Global b-function, |
---|
| 38 | (2002). |
---|
| 39 | |
---|
| 40 | |
---|
| 41 | PROCEDURES: |
---|
| 42 | bfct(f[,s,t,v]); compute the BS polynomial of f with linear reductions |
---|
| 43 | bfctSyz(f[,r,s,t,u,v]); compute the BS polynomial of f with syzygy-solver |
---|
| 44 | bfctAnn(f[,s]); compute the BS polynomial of f via Ann f^s + f |
---|
| 45 | bfctOneGB(f[,s,t]); compute the BS polynomial of f by just one GB computation |
---|
| 46 | bfctIdeal(I,w[,s,t]); compute the b-function of ideal w.r.t. weight |
---|
| 47 | pIntersect(f,I[,s]); intersection of ideal with subalgebra K[f] for a polynomial f |
---|
| 48 | pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] with syz-solver |
---|
| 49 | linReduce(f,I[,s]); reduce a polynomial by linear reductions w.r.t. ideal |
---|
| 50 | linReduceIdeal(I,[s,t]); interreduce generators of ideal by linear reductions |
---|
| 51 | linSyzSolve(I[,s]); compute a linear dependency of elements of ideal I |
---|
| 52 | |
---|
| 53 | allPositive(v); checks whether all entries of an intvec are positive |
---|
| 54 | scalarProd(v,w); computes the standard scalar product of two intvecs |
---|
| 55 | vec2poly(v[,i]); constructs an univariate polynomial with given coefficients |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | SEE ALSO: dmod_lib, dmodapp_lib, dmodvar_lib, gmssing_lib |
---|
| 59 | |
---|
| 60 | KEYWORDS: D-module; global Bernstein-Sato polynomial; Bernstein-Sato polynomial; b-function; |
---|
| 61 | graded Weyl algebra; initial ideal; initial form; principal intersection; linear interreduction; |
---|
| 62 | initial ideal approach |
---|
| 63 | "; |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | |
---|
| 67 | /* |
---|
| 68 | CHANGELOG |
---|
| 69 | |
---|
| 70 | 03.03.11: |
---|
| 71 | - simplified scalarProd |
---|
| 72 | - fixed bug in bfct when user used vars t,Dt |
---|
| 73 | - now bFactor is used by bfct, bfctAnn, i.e. the static procs |
---|
| 74 | addRoot, listofroots are superfluous |
---|
| 75 | - fixed printlevel/debug message issue in bfct, bfctAnn |
---|
| 76 | - fixed small issue for zero ideal input in linReduceIdeal |
---|
| 77 | |
---|
| 78 | 16.03.11 |
---|
| 79 | - fixed bug in linReduceIdeal when ideal contained unlucky constellation |
---|
| 80 | of zeros |
---|
| 81 | - fixed printlevel/debug message issue in linReduceIdeal |
---|
| 82 | */ |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | LIB "qhmoduli.lib"; // for Max |
---|
| 87 | LIB "dmod.lib"; // for SannfsBFCT etc |
---|
| 88 | LIB "dmodapp.lib"; // for initialIdealW etc |
---|
| 89 | LIB "nctools.lib"; // for isWeyl etc |
---|
| 90 | LIB "presolve.lib"; // for valvars |
---|
| 91 | |
---|
| 92 | //--------------- auxiliary procedures ---------------------------------------- |
---|
| 93 | |
---|
| 94 | /* |
---|
| 95 | static proc gradedWeyl (intvec u,intvec v) |
---|
| 96 | "USAGE: gradedWeyl(u,v); u,v intvecs |
---|
| 97 | RETURN: a ring, the associated graded ring of the basering w.r.t. u and v |
---|
| 98 | PURPOSE: compute the associated graded ring of the basering w.r.t. u and v |
---|
| 99 | ASSUME: basering is a Weyl algebra |
---|
| 100 | EXAMPLE: example gradedWeyl; shows examples |
---|
| 101 | NOTE: u[i] is the weight of x(i), v[i] the weight of D(i). |
---|
| 102 | @* u+v has to be a non-negative intvec. |
---|
| 103 | " |
---|
| 104 | { |
---|
| 105 | // todo check assumption |
---|
| 106 | int i; |
---|
| 107 | def save = basering; |
---|
| 108 | int n = nvars(save)/2; |
---|
| 109 | if (nrows(u)<>n || nrows(v)<>n) |
---|
| 110 | { |
---|
| 111 | ERROR("weight vectors have wrong dimension"); |
---|
| 112 | } |
---|
| 113 | intvec uv,gr; |
---|
| 114 | uv = u+v; |
---|
| 115 | for (i=1; i<=n; i++) |
---|
| 116 | { |
---|
| 117 | if (uv[i]>=0) |
---|
| 118 | { |
---|
| 119 | if (uv[i]==0) |
---|
| 120 | { |
---|
| 121 | gr[i] = 0; |
---|
| 122 | } |
---|
| 123 | else |
---|
| 124 | { |
---|
| 125 | gr[i] = 1; |
---|
| 126 | } |
---|
| 127 | } |
---|
| 128 | else |
---|
| 129 | { |
---|
| 130 | ERROR("the sum of the weight vectors has to be a non-negative intvec"); |
---|
| 131 | } |
---|
| 132 | } |
---|
| 133 | list l = ringlist(save); |
---|
| 134 | list l2 = l[2]; |
---|
| 135 | matrix l6 = l[6]; |
---|
| 136 | for (i=1; i<=n; i++) |
---|
| 137 | { |
---|
| 138 | if (gr[i] == 1) |
---|
| 139 | { |
---|
| 140 | l2[n+i] = "xi("+string(i)+")"; |
---|
| 141 | l6[i,n+i] = 0; |
---|
| 142 | } |
---|
| 143 | } |
---|
| 144 | l[2] = l2; |
---|
| 145 | l[6] = l6; |
---|
| 146 | def G = ring(l); |
---|
| 147 | return(G); |
---|
| 148 | } |
---|
| 149 | example |
---|
| 150 | { |
---|
| 151 | "EXAMPLE:"; echo = 2; |
---|
| 152 | ring @D = 0,(x,y,z,Dx,Dy,Dz),dp; |
---|
| 153 | def D = Weyl(); |
---|
| 154 | setring D; |
---|
| 155 | intvec u = -1,-1,1; intvec v = 2,1,1; |
---|
| 156 | def G = gradedWeyl(u,v); |
---|
| 157 | setring G; G; |
---|
| 158 | } |
---|
| 159 | */ |
---|
| 160 | |
---|
| 161 | static proc safeVarName (string s) |
---|
| 162 | { |
---|
| 163 | string S = "," + charstr(basering) + "," + varstr(basering) + ","; |
---|
| 164 | s = "," + s + ","; |
---|
| 165 | while (find(S,s) <> 0) |
---|
| 166 | { |
---|
| 167 | s[1] = "@"; |
---|
| 168 | s = "," + s; |
---|
| 169 | } |
---|
| 170 | s = s[2..size(s)-1]; |
---|
| 171 | return(s) |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | proc allPositive (intvec v) |
---|
| 175 | "USAGE: allPositive(v); v an intvec |
---|
| 176 | RETURN: int, 1 if all components of v are positive, or 0 otherwise |
---|
| 177 | PURPOSE: check whether all components of an intvec are positive |
---|
| 178 | EXAMPLE: example allPositive; shows an example |
---|
| 179 | " |
---|
| 180 | { |
---|
| 181 | int i; |
---|
| 182 | for (i=1; i<=size(v); i++) |
---|
| 183 | { |
---|
| 184 | if (v[i]<=0) |
---|
| 185 | { |
---|
| 186 | return(0); |
---|
| 187 | break; |
---|
| 188 | } |
---|
| 189 | } |
---|
| 190 | return(1); |
---|
| 191 | } |
---|
| 192 | example |
---|
| 193 | { |
---|
| 194 | "EXAMPLE:"; echo = 2; |
---|
| 195 | intvec v = 1,2,3; |
---|
| 196 | allPositive(v); |
---|
| 197 | intvec w = 1,-2,3; |
---|
| 198 | allPositive(w); |
---|
| 199 | } |
---|
| 200 | |
---|
| 201 | static proc findFirst (list l, i) |
---|
| 202 | "USAGE: findFirst(l,i); l a list, i an argument of any type |
---|
| 203 | RETURN: int, the position of the first appearance of i in l, |
---|
| 204 | @* or 0 if i is not a member of l |
---|
| 205 | PURPOSE: check whether the second argument is a member of a list |
---|
| 206 | EXAMPLE: example findFirst; shows an example |
---|
| 207 | " |
---|
| 208 | { |
---|
| 209 | int j; |
---|
| 210 | for (j=1; j<=size(l); j++) |
---|
| 211 | { |
---|
| 212 | if (l[j]==i) |
---|
| 213 | { |
---|
| 214 | return(j); |
---|
| 215 | break; |
---|
| 216 | } |
---|
| 217 | } |
---|
| 218 | return(0); |
---|
| 219 | } |
---|
| 220 | // findFirst(list(1, 2, list(1)),2); // seems to be a bit buggy, |
---|
| 221 | // findFirst(list(1, 2, list(1)),3); // but works for the purposes |
---|
| 222 | // findFirst(list(1, 2, list(1)),list(1)); // of this library |
---|
| 223 | // findFirst(l,list(2)); |
---|
| 224 | example |
---|
| 225 | { |
---|
| 226 | "EXAMPLE:"; echo = 2; |
---|
| 227 | ring r = 0,x,dp; |
---|
| 228 | list l = 1,2,3; |
---|
| 229 | findFirst(l,4); |
---|
| 230 | findFirst(l,2); |
---|
| 231 | } |
---|
| 232 | |
---|
| 233 | proc scalarProd (intvec v, intvec w) |
---|
| 234 | "USAGE: scalarProd(v,w); v,w intvecs |
---|
| 235 | RETURN: int, the standard scalar product of v and w |
---|
| 236 | PURPOSE: computes the scalar product of two intvecs |
---|
| 237 | ASSUME: the arguments are of the same size |
---|
| 238 | EXAMPLE: example scalarProd; shows examples |
---|
| 239 | " |
---|
| 240 | { |
---|
| 241 | if (size(v)!=size(w)) |
---|
| 242 | { |
---|
| 243 | ERROR("non-matching dimensions"); |
---|
| 244 | } |
---|
| 245 | else |
---|
| 246 | { |
---|
| 247 | intvec u = transpose(v)*w; |
---|
| 248 | } |
---|
| 249 | return(u[1]); |
---|
| 250 | } |
---|
| 251 | example |
---|
| 252 | { |
---|
| 253 | "EXAMPLE:"; echo = 2; |
---|
| 254 | intvec v = 1,2,3; |
---|
| 255 | intvec w = 4,5,6; |
---|
| 256 | scalarProd(v,w); |
---|
| 257 | } |
---|
| 258 | |
---|
| 259 | //-------------- main procedures ------------------------------------------------------- |
---|
| 260 | |
---|
| 261 | proc linReduceIdeal(ideal I, list #) |
---|
| 262 | "USAGE: linReduceIdeal(I [,s,t,u]); I an ideal, s,t,u optional ints |
---|
| 263 | RETURN: ideal or list, linear reductum (over field) of I by its elements |
---|
| 264 | PURPOSE: reduces a list of polys only by linear reductions (no monomial |
---|
| 265 | @* multiplications) |
---|
| 266 | NOTE: If s<>0, a list consisting of the reduced ideal and the coefficient |
---|
| 267 | @* vectors of the used reductions given as module is returned. |
---|
| 268 | @* Otherwise (and by default), only the reduced ideal is returned. |
---|
| 269 | @* If t<>0 (and by default) all monomials are reduced (if possible), |
---|
| 270 | @* otherwise, only leading monomials are reduced. |
---|
| 271 | @* If u<>0 (and by default), the ideal is first sorted in increasing order. |
---|
| 272 | @* If u is set to 0 and the given ideal is not sorted in the way described, |
---|
| 273 | @* the result might not be as expected. |
---|
| 274 | DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, |
---|
| 275 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 276 | EXAMPLE: example linReduceIdeal; shows examples |
---|
| 277 | " |
---|
| 278 | { |
---|
| 279 | // #[1] = remembercoeffs |
---|
| 280 | // #[2] = redtail |
---|
| 281 | // #[3] = sortideal |
---|
| 282 | int ppl = printlevel - voice + 2; |
---|
| 283 | int remembercoeffs = 0; // default |
---|
| 284 | int redtail = 1; // default |
---|
| 285 | int sortideal = 1; // default |
---|
| 286 | if (size(#)>0) |
---|
| 287 | { |
---|
| 288 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 289 | { |
---|
| 290 | remembercoeffs = #[1]; |
---|
| 291 | } |
---|
| 292 | if (size(#)>1) |
---|
| 293 | { |
---|
| 294 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 295 | { |
---|
| 296 | redtail = #[2]; |
---|
| 297 | } |
---|
| 298 | if (size(#)>2) |
---|
| 299 | { |
---|
| 300 | if (typeof(#[3])=="int" || typeof(#[3])=="number") |
---|
| 301 | { |
---|
| 302 | sortideal = #[3]; |
---|
| 303 | } |
---|
| 304 | } |
---|
| 305 | } |
---|
| 306 | } |
---|
| 307 | int sI = ncols(I); |
---|
| 308 | int sZeros = sI - size(I); |
---|
| 309 | int i,j,k; |
---|
| 310 | ideal J,lmJ,ordJ; |
---|
| 311 | list lJ = sort(I); |
---|
| 312 | intvec iv,iv2; //todo |
---|
| 313 | module M; // for the coefficients |
---|
| 314 | // step 1: prepare, e.g. sort I |
---|
| 315 | if (sortideal <> 0) |
---|
| 316 | { |
---|
| 317 | if (sZeros > 0) // I contains zeros |
---|
| 318 | { |
---|
| 319 | if (remembercoeffs <> 0) |
---|
| 320 | { |
---|
| 321 | j = 0; |
---|
| 322 | k = 0; |
---|
| 323 | intvec posNonZero; |
---|
| 324 | for (i=1; i<=sI; i++) |
---|
| 325 | { |
---|
| 326 | if (I[i] == 0) |
---|
| 327 | { |
---|
| 328 | j++; |
---|
| 329 | J[j] = 0; |
---|
| 330 | ordJ[j] = -1; |
---|
| 331 | M[j] = gen(i); |
---|
| 332 | } |
---|
| 333 | else |
---|
| 334 | { |
---|
| 335 | k++; |
---|
| 336 | M[k+sZeros] = gen(lJ[2][k]); |
---|
| 337 | posNonZero = posNonZero,i; |
---|
| 338 | } |
---|
| 339 | } |
---|
| 340 | posNonZero = posNonZero[2..nrows(posNonZero)]; |
---|
| 341 | posNonZero = posNonZero[lJ[2]]; |
---|
| 342 | for (i=1; i<=size(lJ[1]); i++) |
---|
| 343 | { |
---|
| 344 | M[i+sZeros] = gen(posNonZero[i]); |
---|
| 345 | } |
---|
| 346 | kill posNonZero; |
---|
| 347 | } |
---|
| 348 | else |
---|
| 349 | { |
---|
| 350 | for (i=1; i<=sZeros; i++) |
---|
| 351 | { |
---|
| 352 | J[i] = 0; |
---|
| 353 | ordJ[i] = -1; |
---|
| 354 | } |
---|
| 355 | } |
---|
| 356 | I = J,lJ[1]; |
---|
| 357 | } |
---|
| 358 | else // I contains no zeros |
---|
| 359 | { |
---|
| 360 | I = lJ[1]; |
---|
| 361 | if (remembercoeffs <> 0) |
---|
| 362 | { |
---|
| 363 | for (i=1; i<=size(lJ[1]); i++) { M[i] = gen(lJ[2][i]); } |
---|
| 364 | } |
---|
| 365 | } |
---|
| 366 | } |
---|
| 367 | else // assume I is already sorted |
---|
| 368 | { |
---|
| 369 | if (remembercoeffs <> 0) |
---|
| 370 | { |
---|
| 371 | for (i=1; i<=ncols(I); i++) { M[i] = gen(i); } |
---|
| 372 | } |
---|
| 373 | } |
---|
| 374 | dbprint(ppl,"// initially sorted ideal:", I); |
---|
| 375 | if (remembercoeffs <> 0) { dbprint(ppl,"// used permutations:", M); } |
---|
| 376 | // step 2: reduce leading monomials by linear reductions |
---|
| 377 | poly lm,c,redpoly,maxlmJ; |
---|
| 378 | J[sZeros+1] = I[sZeros+1]; |
---|
| 379 | lm = I[sZeros+1]; |
---|
| 380 | maxlmJ = leadmonom(lm); |
---|
| 381 | lmJ[sZeros+1] = maxlmJ; |
---|
| 382 | int ordlm,reduction; |
---|
| 383 | ordJ[sZeros+1] = ord(lm); |
---|
| 384 | vector v; |
---|
| 385 | int noRedPast; |
---|
| 386 | for (i=sZeros+2; i<=sI; i++) |
---|
| 387 | { |
---|
| 388 | redpoly = I[i]; |
---|
| 389 | lm = leadmonom(redpoly); |
---|
| 390 | ordlm = ord(lm); |
---|
| 391 | if (remembercoeffs <> 0) { v = M[i]; } |
---|
| 392 | reduction = 1; |
---|
| 393 | while (reduction == 1) // while there was a reduction |
---|
| 394 | { |
---|
| 395 | noRedPast = i; |
---|
| 396 | reduction = 0; |
---|
| 397 | for (j=sZeros+1; j<noRedPast; j++) |
---|
| 398 | { |
---|
| 399 | if (lm == 0) { break; } // nothing more to reduce |
---|
| 400 | if (lm > maxlmJ) { break; } //lm is bigger than maximal monomial to reduce with |
---|
| 401 | if (ordlm == ordJ[j]) |
---|
| 402 | { |
---|
| 403 | if (lm == lmJ[j]) |
---|
| 404 | { |
---|
| 405 | dbprint(ppl-1,"// reducing " + string(redpoly)); |
---|
| 406 | dbprint(ppl-1,"// with " + string(J[j])); |
---|
| 407 | c = leadcoef(redpoly)/leadcoef(J[j]); |
---|
| 408 | redpoly = redpoly - c*J[j]; |
---|
| 409 | dbprint(ppl-1,"// to " + string(redpoly)); |
---|
| 410 | lm = leadmonom(redpoly); |
---|
| 411 | ordlm = ord(lm); |
---|
| 412 | if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; } |
---|
| 413 | noRedPast = j; |
---|
| 414 | reduction = 1; |
---|
| 415 | } |
---|
| 416 | } |
---|
| 417 | } |
---|
| 418 | } |
---|
| 419 | for (j=sZeros+1; j<i; j++) |
---|
| 420 | { |
---|
| 421 | if (redpoly < J[j]) { break; } |
---|
| 422 | } |
---|
| 423 | J = insertGenerator(J,redpoly,j); |
---|
| 424 | lmJ = insertGenerator(lmJ,lm,j); |
---|
| 425 | ordJ = insertGenerator(ordJ,poly(ordlm),j); |
---|
| 426 | if (remembercoeffs <> 0) |
---|
| 427 | { |
---|
| 428 | v = M[i]; |
---|
| 429 | M = deleteGenerator(M,i); |
---|
| 430 | M = insertGenerator(M,v,j); |
---|
| 431 | } |
---|
| 432 | maxlmJ = lmJ[i]; |
---|
| 433 | } |
---|
| 434 | // step 3: reduce tails by linear reductions as well |
---|
| 435 | if (redtail <> 0) |
---|
| 436 | { |
---|
| 437 | dbprint(ppl,"// finished reducing leading monomials"); |
---|
| 438 | dbprint(ppl-1,string(J)); |
---|
| 439 | if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(M)); } |
---|
| 440 | for (i=sZeros+1; i<=sI; i++) |
---|
| 441 | { |
---|
| 442 | lm = lmJ[i]; |
---|
| 443 | for (j=i+1; j<=sI; j++) |
---|
| 444 | { |
---|
| 445 | for (k=2; k<=size(J[j]); k++) // run over all terms in J[j] |
---|
| 446 | { |
---|
| 447 | if (ordJ[i] == ord(J[j][k])) |
---|
| 448 | { |
---|
| 449 | if (lm == normalize(J[j][k])) |
---|
| 450 | { |
---|
| 451 | c = leadcoef(J[j][k])/leadcoef(J[i]); |
---|
| 452 | dbprint(ppl-1,"// reducing " + string(J[j])); |
---|
| 453 | dbprint(ppl-1,"// with " + string(J[i])); |
---|
| 454 | J[j] = J[j] - c*J[i]; |
---|
| 455 | dbprint(ppl-1,"// to " + string(J[j])); |
---|
| 456 | if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; } |
---|
| 457 | } |
---|
| 458 | } |
---|
| 459 | } |
---|
| 460 | } |
---|
| 461 | } |
---|
| 462 | } |
---|
| 463 | if (sI == sZeros) // if I=0,0,...,0, we now have one too much by construction due to sort |
---|
| 464 | { |
---|
| 465 | J = J[1..sZeros]; |
---|
| 466 | } |
---|
| 467 | if (remembercoeffs <> 0) { return(list(J,M)); } |
---|
| 468 | else { return(J); } |
---|
| 469 | } |
---|
| 470 | example |
---|
| 471 | { |
---|
| 472 | "EXAMPLE:"; echo = 2; |
---|
| 473 | ring r = 0,(x,y),dp; |
---|
| 474 | ideal I = 3,x+9,y4+5x,2y4+7x+2; |
---|
| 475 | linReduceIdeal(I); // reduces tails |
---|
| 476 | linReduceIdeal(I,0,0); // no reductions of tails |
---|
| 477 | list l = linReduceIdeal(I,1); // reduces tails and shows reductions used |
---|
| 478 | l; |
---|
| 479 | module M = I; |
---|
| 480 | matrix(l[1]) - M*l[2]; |
---|
| 481 | } |
---|
| 482 | |
---|
| 483 | proc linReduce(poly f, ideal I, list #) |
---|
| 484 | "USAGE: linReduce(f, I [,s,t,u]); f a poly, I an ideal, s,t,u optional ints |
---|
| 485 | RETURN: poly or list, linear reductum (over field) of f by elements from I |
---|
| 486 | PURPOSE: reduce a polynomial only by linear reductions (no monomial multiplications) |
---|
| 487 | NOTE: If s<>0, a list consisting of the reduced polynomial and the coefficient |
---|
| 488 | @* vector of the used reductions is returned, otherwise (and by default) |
---|
| 489 | @* only reduced polynomial is returned. |
---|
| 490 | @* If t<>0 (and by default) all monomials are reduced (if possible), |
---|
| 491 | @* otherwise, only leading monomials are reduced. |
---|
| 492 | @* If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. |
---|
| 493 | @* instead of the given ideal, the output of @code{linReduceIdeal} is used. |
---|
| 494 | @* If u is set to 0 and the given ideal does not equal the output of |
---|
| 495 | @* @code{linReduceIdeal}, the result might not be as expected. |
---|
| 496 | DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, |
---|
| 497 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 498 | EXAMPLE: example linReduce; shows examples |
---|
| 499 | " |
---|
| 500 | { |
---|
| 501 | int ppl = printlevel - voice + 2; |
---|
| 502 | int remembercoeffs = 0; // default |
---|
| 503 | int redtail = 1; // default |
---|
| 504 | int prepareideal = 1; // default |
---|
| 505 | if (size(#)>0) |
---|
| 506 | { |
---|
| 507 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 508 | { |
---|
| 509 | remembercoeffs = #[1]; |
---|
| 510 | } |
---|
| 511 | if (size(#)>1) |
---|
| 512 | { |
---|
| 513 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 514 | { |
---|
| 515 | redtail = #[2]; |
---|
| 516 | } |
---|
| 517 | if (size(#)>2) |
---|
| 518 | { |
---|
| 519 | if (typeof(#[3])=="int" || typeof(#[3])=="number") |
---|
| 520 | { |
---|
| 521 | prepareideal = #[3]; |
---|
| 522 | } |
---|
| 523 | } |
---|
| 524 | } |
---|
| 525 | } |
---|
| 526 | int i,j,k; |
---|
| 527 | int sI = ncols(I); |
---|
| 528 | // pre-reduce I: |
---|
| 529 | module M; |
---|
| 530 | if (prepareideal <> 0) |
---|
| 531 | { |
---|
| 532 | if (remembercoeffs <> 0) |
---|
| 533 | { |
---|
| 534 | list sortedI = linReduceIdeal(I,1,redtail); |
---|
| 535 | I = sortedI[1]; |
---|
| 536 | M = sortedI[2]; |
---|
| 537 | dbprint(ppl-1,"// prepared ideal:" +string(I)); |
---|
| 538 | dbprint(ppl-1,"// with operations:" +string(M)); |
---|
| 539 | } |
---|
| 540 | else { I = linReduceIdeal(I,0,redtail); } |
---|
| 541 | } |
---|
| 542 | else |
---|
| 543 | { |
---|
| 544 | if (remembercoeffs <> 0) |
---|
| 545 | { |
---|
| 546 | for (i=1; i<=sI; i++) { M[i] = gen(i); } |
---|
| 547 | } |
---|
| 548 | } |
---|
| 549 | ideal lmI,lcI,ordI; |
---|
| 550 | for (i=1; i<=sI; i++) |
---|
| 551 | { |
---|
| 552 | lmI[i] = leadmonom(I[i]); |
---|
| 553 | lcI[i] = leadcoef(I[i]); |
---|
| 554 | ordI[i] = ord(lmI[i]); |
---|
| 555 | } |
---|
| 556 | vector v; |
---|
| 557 | poly c; |
---|
| 558 | // === reduce leading monomials === |
---|
| 559 | poly lm = leadmonom(f); |
---|
| 560 | int ordf = ord(lm); |
---|
| 561 | for (i=sI; i>=1; i--) // I is sorted: smallest lm's on top |
---|
| 562 | { |
---|
| 563 | if (lm == 0) { break; } |
---|
| 564 | else |
---|
| 565 | { |
---|
| 566 | if (ordf == ordI[i]) |
---|
| 567 | { |
---|
| 568 | if (lm == lmI[i]) |
---|
| 569 | { |
---|
| 570 | c = leadcoef(f)/lcI[i]; |
---|
| 571 | f = f - c*I[i]; |
---|
| 572 | lm = leadmonom(f); |
---|
| 573 | ordf = ord(lm); |
---|
| 574 | if (remembercoeffs <> 0) { v = v - c * M[i]; } |
---|
| 575 | } |
---|
| 576 | } |
---|
| 577 | } |
---|
| 578 | } |
---|
| 579 | // === reduce tails as well === |
---|
| 580 | if (redtail <> 0) |
---|
| 581 | { |
---|
| 582 | dbprint(ppl,"// finished reducing leading monomials"); |
---|
| 583 | dbprint(ppl-1,string(f)); |
---|
| 584 | if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(v)); } |
---|
| 585 | for (i=1; i<=sI; i++) |
---|
| 586 | { |
---|
| 587 | dbprint(ppl,"// testing ideal entry "+string(i)); |
---|
| 588 | for (j=1; j<=size(f); j++) |
---|
| 589 | { |
---|
| 590 | if (ord(f[j]) == ordI[i]) |
---|
| 591 | { |
---|
| 592 | if (normalize(f[j]) == lmI[i]) |
---|
| 593 | { |
---|
| 594 | c = leadcoef(f[j])/lcI[i]; |
---|
| 595 | f = f - c*I[i]; |
---|
| 596 | dbprint(ppl-1,"// reducing with " + string(I[i])); |
---|
| 597 | dbprint(ppl-1,"// to " + string(f)); |
---|
| 598 | if (remembercoeffs <> 0) { v = v - c * M[i]; } |
---|
| 599 | } |
---|
| 600 | } |
---|
| 601 | } |
---|
| 602 | } |
---|
| 603 | } |
---|
| 604 | if (remembercoeffs <> 0) |
---|
| 605 | { |
---|
| 606 | list l = f,v; |
---|
| 607 | return(l); |
---|
| 608 | } |
---|
| 609 | else { return(f); } |
---|
| 610 | } |
---|
| 611 | example |
---|
| 612 | { |
---|
| 613 | "EXAMPLE:"; echo = 2; |
---|
| 614 | ring r = 0,(x,y),dp; |
---|
| 615 | ideal I = 1,y,xy; |
---|
| 616 | poly f = 5xy+7y+3; |
---|
| 617 | poly g = 7x+5y+3; |
---|
| 618 | linReduce(g,I); // reduces tails |
---|
| 619 | linReduce(g,I,0,0); // no reductions of tails |
---|
| 620 | linReduce(f,I,1); // reduces tails and shows reductions used |
---|
| 621 | f = x3+y2+x2+y+x; |
---|
| 622 | I = x3-y3, y3-x2,x3-y2,x2-y,y2-x; |
---|
| 623 | list l = linReduce(f,I,1); |
---|
| 624 | l; |
---|
| 625 | module M = I; |
---|
| 626 | f - (l[1]-(M*l[2])[1,1]); |
---|
| 627 | } |
---|
| 628 | |
---|
| 629 | proc linSyzSolve (ideal I, list #) |
---|
| 630 | "USAGE: linSyzSolve(I[,s]); I an ideal, s an optional int |
---|
| 631 | RETURN: vector, coefficient vector of linear combination of 0 in elements of I |
---|
| 632 | PURPOSE: compute a linear dependency between the elements of an ideal |
---|
| 633 | @* if such one exists |
---|
| 634 | NOTE: If s<>0, @code{std} is used for Groebner basis computations, |
---|
| 635 | @* otherwise, @code{slimgb} is used. |
---|
| 636 | @* By default, @code{slimgb} is used in char 0 and @code{std} in char >0. |
---|
| 637 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 638 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 639 | EXAMPLE: example linSyzSolve; shows examples |
---|
| 640 | " |
---|
| 641 | { |
---|
| 642 | int whichengine = 0; // default |
---|
| 643 | int enginespecified = 0; // default |
---|
| 644 | if (size(#)>0) |
---|
| 645 | { |
---|
| 646 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 647 | { |
---|
| 648 | whichengine = int( #[1]); |
---|
| 649 | enginespecified = 1; |
---|
| 650 | } |
---|
| 651 | } |
---|
| 652 | int ppl = printlevel - voice +2; |
---|
| 653 | int sI = ncols(I); |
---|
| 654 | // check if we are done |
---|
| 655 | if (I[sI]==0) |
---|
| 656 | { |
---|
| 657 | vector v = gen(sI); |
---|
| 658 | } |
---|
| 659 | else |
---|
| 660 | { |
---|
| 661 | // ------- 1. introduce undefined coeffs ------------------ |
---|
| 662 | def save = basering; |
---|
| 663 | list RL = ringlist(save); |
---|
| 664 | int nv = nvars(save); |
---|
| 665 | int np = npars(save); |
---|
| 666 | int p = char(save); |
---|
| 667 | string cs = "(" + charstr(save) + ")"; |
---|
| 668 | if (enginespecified == 0) |
---|
| 669 | { |
---|
| 670 | if (p <> 0) |
---|
| 671 | { |
---|
| 672 | whichengine = 1; |
---|
| 673 | } |
---|
| 674 | } |
---|
| 675 | int i; |
---|
| 676 | list Lvar; |
---|
| 677 | for (i=1; i<=sI; i++) |
---|
| 678 | { |
---|
| 679 | Lvar[i] = safeVarName("@a(" + string(i) + ")"); |
---|
| 680 | } |
---|
| 681 | list L@A = RL[1..4]; |
---|
| 682 | L@A[2] = Lvar; |
---|
| 683 | L@A[3] = list(list("lp",intvec(1:sI)),list("C",intvec(0))); |
---|
| 684 | def @A = ring(L@A); |
---|
| 685 | L@A[2] = list(safeVarName("@z")); |
---|
| 686 | L@A[3][1] = list("dp",intvec(1)); |
---|
| 687 | if (size(L@A[1])>1) |
---|
| 688 | { |
---|
| 689 | L@A[1][2] = L@A[1][2] + Lvar; |
---|
| 690 | L@A[1][3][size(L@A[1][3])+1] = list("lp",intvec(1:sI)); |
---|
| 691 | } |
---|
| 692 | else |
---|
| 693 | { |
---|
| 694 | L@A[1] = list(p,Lvar,list(list("lp",intvec(1:sI))),ideal(0)); |
---|
| 695 | } |
---|
| 696 | def @aA = ring(L@A); |
---|
| 697 | def @B = save + @aA; |
---|
| 698 | setring @B; |
---|
| 699 | ideal I = imap(save,I); |
---|
| 700 | // ------- 2. form the linear system for the undef coeffs --- |
---|
| 701 | poly W; ideal QQ; |
---|
| 702 | for (i=1; i<=sI; i++) |
---|
| 703 | { |
---|
| 704 | W = W + par(np+i)*I[i]; |
---|
| 705 | } |
---|
| 706 | while (W!=0) |
---|
| 707 | { |
---|
| 708 | QQ = QQ,leadcoef(W); |
---|
| 709 | W = W - lead(W); |
---|
| 710 | } |
---|
| 711 | // QQ consists of polynomial expressions in @a(i) of type number |
---|
| 712 | setring @A; |
---|
| 713 | ideal QQ = imap(@B,QQ); |
---|
| 714 | // ------- 3. this QQ is a polynomial ideal, so "solve" the system ----- |
---|
| 715 | dbprint(ppl, "// linSyzSolve: starting Groebner basis computation with engine:", whichengine); |
---|
| 716 | QQ = engine(QQ,whichengine); |
---|
| 717 | dbprint(ppl, "// QQ after engine:", QQ); |
---|
| 718 | if (dim(QQ) == -1) |
---|
| 719 | { |
---|
| 720 | dbprint(ppl+1, "// no solutions by linSyzSolve"); |
---|
| 721 | // output zeroes |
---|
| 722 | setring save; |
---|
| 723 | kill @A,@aA,@B; |
---|
| 724 | return(v); |
---|
| 725 | } |
---|
| 726 | // ------- 4. in order to get the numeric values ------- |
---|
| 727 | matrix AA = matrix(maxideal(1)); |
---|
| 728 | module MQQ = std(module(QQ)); |
---|
| 729 | AA = NF(AA,MQQ); // todo: we still receive NF warnings |
---|
| 730 | dbprint(ppl, "// AA after NF:",AA); |
---|
| 731 | // "AA after NF:"; print(AA); |
---|
| 732 | for(i=1; i<=sI; i++) |
---|
| 733 | { |
---|
| 734 | AA = subst(AA,var(i),1); |
---|
| 735 | } |
---|
| 736 | dbprint(ppl, "// AA after subst:",AA); |
---|
| 737 | // "AA after subst: "; print(AA); |
---|
| 738 | vector v = (module(transpose(AA)))[1]; |
---|
| 739 | setring save; |
---|
| 740 | vector v = imap(@A,v); |
---|
| 741 | kill @A,@aA,@B; |
---|
| 742 | } |
---|
| 743 | return(v); |
---|
| 744 | } |
---|
| 745 | example |
---|
| 746 | { |
---|
| 747 | "EXAMPLE:"; echo = 2; |
---|
| 748 | ring r = 0,x,dp; |
---|
| 749 | ideal I = x,2x; |
---|
| 750 | linSyzSolve(I); |
---|
| 751 | ideal J = x,x2; |
---|
| 752 | linSyzSolve(J); |
---|
| 753 | } |
---|
| 754 | |
---|
| 755 | proc pIntersect (poly s, ideal I, list #) |
---|
| 756 | "USAGE: pIntersect(f, I [,s]); f a poly, I an ideal, s an optional int |
---|
| 757 | RETURN: vector, coefficient vector of the monic polynomial |
---|
| 758 | PURPOSE: compute the intersection of ideal I with the subalgebra K[f] |
---|
| 759 | ASSUME: I is given as Groebner basis, basering is not a qring. |
---|
| 760 | NOTE: If the intersection is zero, this proc might not terminate. |
---|
| 761 | @* If s>0 is given, it is searched for the generator of the intersection |
---|
| 762 | @* only up to degree s. Otherwise (and by default), no bound is assumed. |
---|
| 763 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 764 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 765 | EXAMPLE: example pIntersect; shows examples |
---|
| 766 | " |
---|
| 767 | { |
---|
| 768 | def save = basering; |
---|
| 769 | int n = nvars(save); |
---|
| 770 | list RL = ringlist(save); |
---|
| 771 | int i,j,k; |
---|
| 772 | if (RL[4]<>0) |
---|
| 773 | { |
---|
| 774 | ERROR ("basering must not be a qring"); |
---|
| 775 | } |
---|
| 776 | // assume I is given in Groebner basis |
---|
| 777 | if (attrib(I,"isSB") <> 1) |
---|
| 778 | { |
---|
| 779 | print("// WARNING: The input has no SB attribute!"); |
---|
| 780 | print("// Treating it as if it were a Groebner basis and proceeding..."); |
---|
| 781 | attrib(I,"isSB",1); // set attribute for suppressing NF messages |
---|
| 782 | } |
---|
| 783 | int bound = 0; // default |
---|
| 784 | if (size(#)>0) |
---|
| 785 | { |
---|
| 786 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 787 | { |
---|
| 788 | bound = #[1]; |
---|
| 789 | } |
---|
| 790 | } |
---|
| 791 | int ppl = printlevel-voice+2; |
---|
| 792 | // ---case 1: I = basering--- |
---|
| 793 | if (size(I) == 1) |
---|
| 794 | { |
---|
| 795 | if (simplify(I,2)[1] == 1) |
---|
| 796 | { |
---|
| 797 | return(gen(1)); // = 1 |
---|
| 798 | } |
---|
| 799 | } |
---|
| 800 | // ---case 2: intersection is zero--- |
---|
| 801 | intvec degs = leadexp(s); |
---|
| 802 | intvec possdegbounds; |
---|
| 803 | list degI; |
---|
| 804 | i = 1; |
---|
| 805 | while (i <= ncols(I)) |
---|
| 806 | { |
---|
| 807 | if (i == ncols(I)+1) { break; } |
---|
| 808 | degI[i] = leadexp(I[i]); |
---|
| 809 | for (j=1; j<=n; j++) |
---|
| 810 | { |
---|
| 811 | if (degs[j] == 0) |
---|
| 812 | { |
---|
| 813 | if (degI[i][j] <> 0) { break; } |
---|
| 814 | } |
---|
| 815 | if (j == n) |
---|
| 816 | { |
---|
| 817 | k++; |
---|
| 818 | possdegbounds[k] = Max(degI[i]); |
---|
| 819 | } |
---|
| 820 | } |
---|
| 821 | i++; |
---|
| 822 | } |
---|
| 823 | int degbound = Min(possdegbounds); |
---|
| 824 | if (bound>0 && bound<degbound) // given bound is too small |
---|
| 825 | { |
---|
| 826 | print("// Try a bound of at least " + string(degbound)); |
---|
| 827 | return(vector(0)); |
---|
| 828 | } |
---|
| 829 | dbprint(ppl,"// lower bound for the degree of the insection is " +string(degbound)); |
---|
| 830 | if (degbound == 0) // lm(s) does not appear in lm(I) |
---|
| 831 | { |
---|
| 832 | print("// Intersection is zero"); |
---|
| 833 | return(vector(0)); |
---|
| 834 | } |
---|
| 835 | // ---case 3: intersection is non-trivial--- |
---|
| 836 | ideal redNI = 1; |
---|
| 837 | vector v; |
---|
| 838 | list l,ll; |
---|
| 839 | l[1] = vector(0); |
---|
| 840 | poly toNF,tobracket,newNF,rednewNF,oldNF,secNF; |
---|
| 841 | i = 1; |
---|
| 842 | while (1) |
---|
| 843 | { |
---|
| 844 | if (bound>0 && i>bound) { return(vector(0)); } |
---|
| 845 | dbprint(ppl,"// Testing degree "+string(i)); |
---|
| 846 | if (i>1) |
---|
| 847 | { |
---|
| 848 | oldNF = newNF; |
---|
| 849 | tobracket = s^(i-1) - oldNF; |
---|
| 850 | if (tobracket==0) // todo bug in bracket? |
---|
| 851 | { |
---|
| 852 | toNF = 0; |
---|
| 853 | } |
---|
| 854 | else |
---|
| 855 | { |
---|
| 856 | toNF = bracket(tobracket,secNF); |
---|
| 857 | } |
---|
| 858 | newNF = NF(toNF+oldNF*secNF,I); // = NF(s^i,I) |
---|
| 859 | } |
---|
| 860 | else |
---|
| 861 | { |
---|
| 862 | newNF = NF(s,I); |
---|
| 863 | secNF = newNF; |
---|
| 864 | } |
---|
| 865 | ll = linReduce(newNF,redNI,1); |
---|
| 866 | rednewNF = ll[1]; |
---|
| 867 | l[i+1] = ll[2]; |
---|
| 868 | dbprint(ppl-1,"// newNF is: "+string(newNF)); |
---|
| 869 | dbprint(ppl-1,"// rednewNF is: "+string(rednewNF)); |
---|
| 870 | if (rednewNF != 0) // no linear dependency |
---|
| 871 | { |
---|
| 872 | redNI[i+1] = rednewNF; |
---|
| 873 | i++; |
---|
| 874 | } |
---|
| 875 | else // there is a linear dependency, hence we are done |
---|
| 876 | { |
---|
| 877 | dbprint(ppl,"// degree of the generator of the intersection is: "+string(i)); |
---|
| 878 | break; |
---|
| 879 | } |
---|
| 880 | } |
---|
| 881 | dbprint(ppl-1,"// used linear reductions:", l); |
---|
| 882 | // we obtain the coefficients of the generator by the used reductions: |
---|
| 883 | list Lvar; |
---|
| 884 | for (j=1; j<=i+1; j++) |
---|
| 885 | { |
---|
| 886 | Lvar[j] = safeVarName("a("+string(j)+")"); |
---|
| 887 | } |
---|
| 888 | list Lord = list("dp",intvec(1:(i+1))),list("C",intvec(0)); |
---|
| 889 | list L@R = RL[1..4]; |
---|
| 890 | L@R[2] = Lvar; L@R[3] = Lord; |
---|
| 891 | def @R = ring(L@R); setring @R; |
---|
| 892 | list l = imap(save,l); |
---|
| 893 | ideal C; |
---|
| 894 | for (j=1;j<=i+1;j++) |
---|
| 895 | { |
---|
| 896 | C[j] = 0; |
---|
| 897 | for (k=1;k<=j;k++) |
---|
| 898 | { |
---|
| 899 | C[j] = C[j]+l[j][k]*var(k); |
---|
| 900 | } |
---|
| 901 | } |
---|
| 902 | for (j=i;j>=1;j--) |
---|
| 903 | { |
---|
| 904 | C[i+1]= subst(C[i+1],var(j),var(j)+C[j]); |
---|
| 905 | } |
---|
| 906 | matrix m = coeffs(C[i+1],maxideal(1)); |
---|
| 907 | vector v = gen(i+1); |
---|
| 908 | for (j=1;j<=i+1;j++) |
---|
| 909 | { |
---|
| 910 | v = v + m[j,1]*gen(j); |
---|
| 911 | } |
---|
| 912 | setring save; |
---|
| 913 | v = imap(@R,v); |
---|
| 914 | kill @R; |
---|
| 915 | return(v); |
---|
| 916 | } |
---|
| 917 | example |
---|
| 918 | { |
---|
| 919 | "EXAMPLE:"; echo = 2; |
---|
| 920 | ring r = 0,(x,y),dp; |
---|
| 921 | poly f = x^2+y^3+x*y^2; |
---|
| 922 | def D = initialMalgrange(f); |
---|
| 923 | setring D; |
---|
| 924 | inF; |
---|
| 925 | pIntersect(t*Dt,inF); |
---|
| 926 | pIntersect(t*Dt,inF,1); |
---|
| 927 | } |
---|
| 928 | |
---|
| 929 | proc pIntersectSyz (poly s, ideal I, list #) |
---|
| 930 | "USAGE: pIntersectSyz(f, I [,p,s,t]); f poly, I ideal, p,t optial ints, p prime |
---|
| 931 | RETURN: vector, coefficient vector of the monic polynomial |
---|
| 932 | PURPOSE: compute the intersection of an ideal I with the subalgebra K[f] |
---|
| 933 | ASSUME: I is given as Groebner basis. |
---|
| 934 | NOTE: If the intersection is zero, this procedure might not terminate. |
---|
| 935 | @* If p>0 is given, this proc computes the generator of the intersection in |
---|
| 936 | @* char p first and then only searches for a generator of the obtained |
---|
| 937 | @* degree in the basering. Otherwise, it searches for all degrees by |
---|
| 938 | @* computing syzygies. |
---|
| 939 | @* If s<>0, @code{std} is used for Groebner basis computations in char 0, |
---|
| 940 | @* otherwise, and by default, @code{slimgb} is used. |
---|
| 941 | @* If t<>0 and by default, @code{std} is used for Groebner basis |
---|
| 942 | @* computations in char >0, otherwise, @code{slimgb} is used. |
---|
| 943 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 944 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 945 | EXAMPLE: example pIntersectSyz; shows examples |
---|
| 946 | " |
---|
| 947 | { |
---|
| 948 | // assume I is given as Groebner basis |
---|
| 949 | if (attrib(I,"isSB") <> 1) |
---|
| 950 | { |
---|
| 951 | print("// WARNING: The input has no SB attribute!"); |
---|
| 952 | print("// Treating it as if it were a Groebner basis and proceeding..."); |
---|
| 953 | attrib(I,"isSB",1); // set attribute for suppressing NF messages |
---|
| 954 | } |
---|
| 955 | int ppl = printlevel-voice+2; |
---|
| 956 | int whichengine = 0; // default |
---|
| 957 | int modengine = 1; // default |
---|
| 958 | int solveincharp = 0; // default |
---|
| 959 | def save = basering; |
---|
| 960 | int n = nvars(save); |
---|
| 961 | if (size(#)>0) |
---|
| 962 | { |
---|
| 963 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 964 | { |
---|
| 965 | solveincharp = int(#[1]); |
---|
| 966 | } |
---|
| 967 | if (size(#)>1) |
---|
| 968 | { |
---|
| 969 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 970 | { |
---|
| 971 | whichengine = int(#[2]); |
---|
| 972 | } |
---|
| 973 | if (size(#)>2) |
---|
| 974 | { |
---|
| 975 | if (typeof(#[3])=="int" || typeof(#[3])=="number") |
---|
| 976 | { |
---|
| 977 | modengine = int(#[3]); |
---|
| 978 | } |
---|
| 979 | } |
---|
| 980 | } |
---|
| 981 | } |
---|
| 982 | int i,j; |
---|
| 983 | vector v; |
---|
| 984 | poly tobracket,toNF,newNF,p; |
---|
| 985 | ideal NI = 1; |
---|
| 986 | newNF = NF(s,I); |
---|
| 987 | NI[2] = newNF; |
---|
| 988 | list RL = ringlist(save); |
---|
| 989 | if (solveincharp) |
---|
| 990 | { |
---|
| 991 | int psolveincharp = prime(solveincharp); |
---|
| 992 | if (solveincharp <> psolveincharp) |
---|
| 993 | { |
---|
| 994 | print("// " + string(solveincharp) + " is invalid characteristic of ground field."); |
---|
| 995 | print("// " + string(psolveincharp) + " is used."); |
---|
| 996 | solveincharp = psolveincharp; |
---|
| 997 | kill psolveincharp; |
---|
| 998 | } |
---|
| 999 | list RLp = RL[1..4]; |
---|
| 1000 | if (typeof(RL[1]) == "int") { RLp[1] = solveincharp; } |
---|
| 1001 | else { RLp[1][1] = solveincharp; } // parameters |
---|
| 1002 | def @Rp = ring(RLp); |
---|
| 1003 | setring @Rp; |
---|
| 1004 | number c; |
---|
| 1005 | setring save; |
---|
| 1006 | int shortSave = short; // workaround for maps Q(a_i) -> Z/p(a_i) |
---|
| 1007 | short = 0; |
---|
| 1008 | string str; |
---|
| 1009 | int badprime; |
---|
| 1010 | i=1; |
---|
| 1011 | while (badprime == 0 && i<=size(s)) // detect bad primes |
---|
| 1012 | { |
---|
| 1013 | str = string(denominator(leadcoef(s[i]))); |
---|
| 1014 | str = "c = " + str + ";"; |
---|
| 1015 | setring @Rp; |
---|
| 1016 | execute(str); |
---|
| 1017 | if (c == 0) { badprime = 1; } |
---|
| 1018 | setring save; |
---|
| 1019 | i++; |
---|
| 1020 | } |
---|
| 1021 | str = "poly s = " + string(s) + ";"; |
---|
| 1022 | if (size(RL) > 4) // basering is NC-algebra |
---|
| 1023 | { |
---|
| 1024 | string RL5 = "@C = " + string(RL[5]) + ";"; |
---|
| 1025 | string RL6 = "@D = " + string(RL[6]) + ";"; |
---|
| 1026 | setring @Rp; |
---|
| 1027 | matrix @C[n][n]; matrix @D[n][n]; |
---|
| 1028 | execute(RL5); execute(RL6); |
---|
| 1029 | def Rp = nc_algebra(@C,@D); |
---|
| 1030 | } |
---|
| 1031 | else { def Rp = @Rp; } |
---|
| 1032 | setring Rp; |
---|
| 1033 | kill @Rp; |
---|
| 1034 | dbprint(ppl-1,"// solving in ring ", Rp); |
---|
| 1035 | execute(str); |
---|
| 1036 | vector v; |
---|
| 1037 | number c; |
---|
| 1038 | ideal NI = 1; |
---|
| 1039 | setring save; |
---|
| 1040 | i=1; |
---|
| 1041 | while (badprime == 0 && i<=size(I)) // detect bad primes |
---|
| 1042 | { |
---|
| 1043 | str = string(leadcoef(cleardenom(I[i]))); |
---|
| 1044 | str = "c = " + str + ";"; |
---|
| 1045 | setring Rp; |
---|
| 1046 | execute(str); |
---|
| 1047 | if (c == 0) { badprime = 1; } |
---|
| 1048 | setring save; |
---|
| 1049 | i++; |
---|
| 1050 | } |
---|
| 1051 | if (badprime == 1) |
---|
| 1052 | { |
---|
| 1053 | print("// WARNING: bad prime"); |
---|
| 1054 | short = shortSave; |
---|
| 1055 | return(vector(0)); |
---|
| 1056 | } |
---|
| 1057 | } |
---|
| 1058 | i = 1; |
---|
| 1059 | dbprint(ppl,"// pIntersectSyz starts..."); |
---|
| 1060 | dbprint(ppl-1,"// with ideal I=", I); |
---|
| 1061 | while (1) |
---|
| 1062 | { |
---|
| 1063 | dbprint(ppl,"// testing degree: "+string(i)); |
---|
| 1064 | if (i>1) |
---|
| 1065 | { |
---|
| 1066 | tobracket = s^(i-1)-NI[i]; |
---|
| 1067 | if (tobracket!=0) |
---|
| 1068 | { |
---|
| 1069 | toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2]; |
---|
| 1070 | } |
---|
| 1071 | else |
---|
| 1072 | { |
---|
| 1073 | toNF = NI[i]*NI[2]; |
---|
| 1074 | } |
---|
| 1075 | newNF = NF(toNF,I); |
---|
| 1076 | NI[i+1] = newNF; |
---|
| 1077 | } |
---|
| 1078 | // look for a solution |
---|
| 1079 | dbprint(ppl-1,"// linSyzSolve starts with: "+string(matrix(NI))); |
---|
| 1080 | if (solveincharp) // modular method |
---|
| 1081 | { |
---|
| 1082 | for (j=1; j<=size(newNF); j++) |
---|
| 1083 | { |
---|
| 1084 | str = string(denominator(leadcoef(s[i]))); |
---|
| 1085 | str = "c = " + str + ";"; |
---|
| 1086 | setring Rp; |
---|
| 1087 | execute(str); |
---|
| 1088 | if (c == 0) |
---|
| 1089 | { |
---|
| 1090 | print("// WARNING: bad prime"); |
---|
| 1091 | setring save; |
---|
| 1092 | short = shortSave; |
---|
| 1093 | return(vector(0)); |
---|
| 1094 | } |
---|
| 1095 | setring save; |
---|
| 1096 | } |
---|
| 1097 | str = "NI[" + string(i) + "+1] = " + string(newNF) + ";"; |
---|
| 1098 | setring Rp; |
---|
| 1099 | execute(str); // NI[i+1] = [newNF]_{solveincharp} |
---|
| 1100 | v = linSyzSolve(NI,modengine); |
---|
| 1101 | if (v!=0) // there is a modular solution |
---|
| 1102 | { |
---|
| 1103 | dbprint(ppl,"// got solution in char "+string(solveincharp)+" of degree "+string(i)); |
---|
| 1104 | setring save; |
---|
| 1105 | v = linSyzSolve(NI,whichengine); |
---|
| 1106 | if (v==0) |
---|
| 1107 | { |
---|
| 1108 | break; |
---|
| 1109 | } |
---|
| 1110 | } |
---|
| 1111 | else // no modular solution |
---|
| 1112 | { |
---|
| 1113 | setring save; |
---|
| 1114 | v = 0; |
---|
| 1115 | } |
---|
| 1116 | } |
---|
| 1117 | else // non-modular method |
---|
| 1118 | { |
---|
| 1119 | v = linSyzSolve(NI,whichengine); |
---|
| 1120 | } |
---|
| 1121 | matrix MM[1][nrows(v)] = matrix(v); |
---|
| 1122 | dbprint(ppl-1,"// linSyzSolve ready with: "+string(MM)); |
---|
| 1123 | kill MM; |
---|
| 1124 | // "linSyzSolve ready with"; print(v); |
---|
| 1125 | if (v!=0) |
---|
| 1126 | { |
---|
| 1127 | // a solution: |
---|
| 1128 | //check for the reality of the solution |
---|
| 1129 | p = 0; |
---|
| 1130 | for (j=1; j<=i+1; j++) |
---|
| 1131 | { |
---|
| 1132 | p = p + v[j]*NI[j]; |
---|
| 1133 | } |
---|
| 1134 | if (p!=0) |
---|
| 1135 | { |
---|
| 1136 | dbprint(ppl,"// linSyzSolve: bad solution!"); |
---|
| 1137 | } |
---|
| 1138 | else |
---|
| 1139 | { |
---|
| 1140 | dbprint(ppl,"// linSyzSolve: got solution!"); |
---|
| 1141 | // "got solution!"; |
---|
| 1142 | break; |
---|
| 1143 | } |
---|
| 1144 | } |
---|
| 1145 | // no solution: |
---|
| 1146 | i++; |
---|
| 1147 | } |
---|
| 1148 | dbprint(ppl,"// pIntersectSyz finished"); |
---|
| 1149 | if (solveincharp) { short = shortSave; } |
---|
| 1150 | return(v); |
---|
| 1151 | } |
---|
| 1152 | example |
---|
| 1153 | { |
---|
| 1154 | "EXAMPLE:"; echo = 2; |
---|
| 1155 | ring r = 0,(x,y),dp; |
---|
| 1156 | poly f = x^2+y^3+x*y^2; |
---|
| 1157 | def D = initialMalgrange(f); |
---|
| 1158 | setring D; |
---|
| 1159 | inF; |
---|
| 1160 | poly s = t*Dt; |
---|
| 1161 | pIntersectSyz(s,inF); |
---|
| 1162 | int p = prime(20000); |
---|
| 1163 | pIntersectSyz(s,inF,p,0,0); |
---|
| 1164 | } |
---|
| 1165 | |
---|
| 1166 | proc vec2poly (list #) |
---|
| 1167 | "USAGE: vec2poly(v [,i]); v a vector or an intvec, i an optional int |
---|
| 1168 | RETURN: poly, an univariate polynomial in i-th variable with coefficients given by v |
---|
| 1169 | PURPOSE: constructs an univariate polynomial in K[var(i)] with given coefficients, |
---|
| 1170 | @* such that the coefficient at var(i)^{j-1} is v[j]. |
---|
| 1171 | NOTE: The optional argument i must be positive, by default i is 1. |
---|
| 1172 | EXAMPLE: example vec2poly; shows examples |
---|
| 1173 | " |
---|
| 1174 | { |
---|
| 1175 | def save = basering; |
---|
| 1176 | int i,ringvar; |
---|
| 1177 | ringvar = 1; // default |
---|
| 1178 | if (size(#) > 0) |
---|
| 1179 | { |
---|
| 1180 | if (typeof(#[1])=="vector" || typeof(#[1])=="intvec") |
---|
| 1181 | { |
---|
| 1182 | def v = #[1]; |
---|
| 1183 | } |
---|
| 1184 | else |
---|
| 1185 | { |
---|
| 1186 | ERROR("wrong input: expected vector/intvec expression"); |
---|
| 1187 | } |
---|
| 1188 | if (size(#) > 1) |
---|
| 1189 | { |
---|
| 1190 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 1191 | { |
---|
| 1192 | ringvar = int(#[2]); |
---|
| 1193 | } |
---|
| 1194 | } |
---|
| 1195 | } |
---|
| 1196 | if (ringvar > nvars(save)) |
---|
| 1197 | { |
---|
| 1198 | ERROR("var out of range"); |
---|
| 1199 | } |
---|
| 1200 | poly p; |
---|
| 1201 | for (i=1; i<=nrows(v); i++) |
---|
| 1202 | { |
---|
| 1203 | p = p + v[i]*(var(ringvar))^(i-1); |
---|
| 1204 | } |
---|
| 1205 | return(p); |
---|
| 1206 | } |
---|
| 1207 | example |
---|
| 1208 | { |
---|
| 1209 | "EXAMPLE:"; echo = 2; |
---|
| 1210 | ring r = 0,(x,y),dp; |
---|
| 1211 | vector v = gen(1) + 3*gen(3) + 22/9*gen(4); |
---|
| 1212 | intvec iv = 3,2,1; |
---|
| 1213 | vec2poly(v,2); |
---|
| 1214 | vec2poly(iv); |
---|
| 1215 | } |
---|
| 1216 | |
---|
| 1217 | /* |
---|
| 1218 | // // listofroots and addRoots aren't needed anymore due to some modifications |
---|
| 1219 | // |
---|
| 1220 | // static proc listofroots (list #) |
---|
| 1221 | // { |
---|
| 1222 | // def save = basering; |
---|
| 1223 | // int n = nvars(save); |
---|
| 1224 | // int i; |
---|
| 1225 | // poly p; |
---|
| 1226 | // if (typeof(#[1])=="vector") |
---|
| 1227 | // { |
---|
| 1228 | // vector b = #[1]; |
---|
| 1229 | // for (i=1; i<=nrows(b); i++) |
---|
| 1230 | // { |
---|
| 1231 | // p = p + b[i]*(var(1))^(i-1); |
---|
| 1232 | // } |
---|
| 1233 | // } |
---|
| 1234 | // else |
---|
| 1235 | // { |
---|
| 1236 | // p = #[1]; |
---|
| 1237 | // } |
---|
| 1238 | // int substitution = int(#[2]); |
---|
| 1239 | // string s = safeVarName("s"); |
---|
| 1240 | // list RL = ringlist(save); RL = RL[1..4]; |
---|
| 1241 | // RL[2] = list(s); RL[3] = list(list("dp",intvec(1)),list("C",0)); |
---|
| 1242 | // def S = ring(RL); setring S; |
---|
| 1243 | // ideal J; |
---|
| 1244 | // for (i=1; i<=n; i++) |
---|
| 1245 | // { |
---|
| 1246 | // J[i] = var(1); |
---|
| 1247 | // } |
---|
| 1248 | // map @m = save,J; |
---|
| 1249 | // poly p = @m(p); |
---|
| 1250 | // if (substitution == 1) |
---|
| 1251 | // { |
---|
| 1252 | // p = subst(p,var(1),-var(1)-1); |
---|
| 1253 | // } |
---|
| 1254 | // // the rest of this proc is nicked from bernsteinBM from dmod.lib |
---|
| 1255 | // list P = factorize(p);//with constants and multiplicities |
---|
| 1256 | // ideal bs; intvec m; //the BS polynomial is monic, so we are not interested in constants |
---|
| 1257 | // for (i=2; i<= size(P[1]); i++) //we delete P[1][1] and P[2][1] |
---|
| 1258 | // { |
---|
| 1259 | // bs[i-1] = P[1][i]; |
---|
| 1260 | // m[i-1] = P[2][i]; |
---|
| 1261 | // } |
---|
| 1262 | // bs = normalize(bs); |
---|
| 1263 | // bs = -subst(bs,var(1),0); |
---|
| 1264 | // setring save; |
---|
| 1265 | // ideal bs = imap(S,bs); |
---|
| 1266 | // kill S; |
---|
| 1267 | // list BS = bs,m; |
---|
| 1268 | // return(BS); |
---|
| 1269 | // } |
---|
| 1270 | // |
---|
| 1271 | // |
---|
| 1272 | // static proc addRoot(number q, list L) |
---|
| 1273 | // { |
---|
| 1274 | // // add root to list in bFactor format |
---|
| 1275 | // int i,qInL; |
---|
| 1276 | // ideal I = L[1]; |
---|
| 1277 | // intvec v = L[2]; |
---|
| 1278 | // list LL; |
---|
| 1279 | // if (v == 0) |
---|
| 1280 | // { |
---|
| 1281 | // I = poly(q); |
---|
| 1282 | // v = 1; |
---|
| 1283 | // LL = I,v; |
---|
| 1284 | // } |
---|
| 1285 | // else |
---|
| 1286 | // { |
---|
| 1287 | // for (i=1; i<=ncols(I); i++) |
---|
| 1288 | // { |
---|
| 1289 | // if (I[i] == q) |
---|
| 1290 | // { |
---|
| 1291 | // qInL = i; |
---|
| 1292 | // break; |
---|
| 1293 | // } |
---|
| 1294 | // } |
---|
| 1295 | // if (qInL) |
---|
| 1296 | // { |
---|
| 1297 | // v[qInL] = v[qInL] + 1; |
---|
| 1298 | // } |
---|
| 1299 | // else |
---|
| 1300 | // { |
---|
| 1301 | // I = q,I; |
---|
| 1302 | // v = 1,v; |
---|
| 1303 | // } |
---|
| 1304 | // } |
---|
| 1305 | // LL = I,v; |
---|
| 1306 | // if (size(L) == 3) // irreducible factor |
---|
| 1307 | // { |
---|
| 1308 | // if (L[3] <> "0" && L[3] <> "1") |
---|
| 1309 | // { |
---|
| 1310 | // LL = LL + list(L[3]); |
---|
| 1311 | // } |
---|
| 1312 | // } |
---|
| 1313 | // return(LL); |
---|
| 1314 | // } |
---|
| 1315 | */ |
---|
| 1316 | |
---|
| 1317 | static proc bfctengine (poly f, int inorann, int whichengine, int addPD, int stdsum, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0) |
---|
| 1318 | { |
---|
| 1319 | int printlevelsave = printlevel; |
---|
| 1320 | printlevel = printlevel + 1; |
---|
| 1321 | int ppl = printlevel - voice + 2; |
---|
| 1322 | int i; |
---|
| 1323 | def save = basering; |
---|
| 1324 | int n = nvars(save); |
---|
| 1325 | if (isCommutative() == 0) { ERROR("basering must be commutative"); } |
---|
| 1326 | if (char(save) <> 0) { ERROR("characteristic of basering has to be 0"); } |
---|
| 1327 | list L = ringlist(save); |
---|
| 1328 | int qr; |
---|
| 1329 | if (L[4] <> 0) // qring |
---|
| 1330 | { |
---|
| 1331 | print("// basering is qring:"); |
---|
| 1332 | print("// discarding the quotient and proceeding..."); |
---|
| 1333 | L[4] = ideal(0); |
---|
| 1334 | qr = 1; |
---|
| 1335 | def save2 = ring(L); setring save2; |
---|
| 1336 | poly f = imap(save,f); |
---|
| 1337 | } |
---|
| 1338 | if (inorann == 0) // bfct using initial ideal |
---|
| 1339 | { |
---|
| 1340 | // list L = ringlist(basering); |
---|
| 1341 | intvec iv = valvars(f)[1]; // heuristacally better ordering for initialMalgrange |
---|
| 1342 | list varL = L[2]; |
---|
| 1343 | varL = varL[iv]; |
---|
| 1344 | L[2] = varL; |
---|
| 1345 | if (u0 <> 0) |
---|
| 1346 | { |
---|
| 1347 | u0 = u0[iv]; |
---|
| 1348 | } |
---|
| 1349 | def newr = ring(L); |
---|
| 1350 | kill varL,iv,L; |
---|
| 1351 | setring newr; |
---|
| 1352 | poly f = imap(save,f); |
---|
| 1353 | dbprint(ppl,"// starting computation of the initial ideal of the Malgrange ideal..."); |
---|
| 1354 | def D = initialMalgrange(f,whichengine,methodord,u0); |
---|
| 1355 | dbprint(ppl,"// ...done"); |
---|
| 1356 | setring D; |
---|
| 1357 | ideal J = inF; |
---|
| 1358 | kill inF; |
---|
| 1359 | poly s = var(1)*var(n+2); |
---|
| 1360 | } |
---|
| 1361 | else // bfct using Ann(f^s) |
---|
| 1362 | { |
---|
| 1363 | dbprint(ppl,"// starting computation of the s-parametric annihilator..."); |
---|
| 1364 | def D = SannfsBFCT(f,addPD,whichengine,stdsum); |
---|
| 1365 | dbprint(ppl,"// ...done"); |
---|
| 1366 | setring D; |
---|
| 1367 | ideal J = LD; |
---|
| 1368 | kill LD; |
---|
| 1369 | poly s = var(1); |
---|
| 1370 | } |
---|
| 1371 | vector b; |
---|
| 1372 | dbprint(ppl,"// starting to intersect with subalgebra..."); |
---|
| 1373 | // try it modular |
---|
| 1374 | if (methodpIntersect <> 0) // pIntersectSyz |
---|
| 1375 | { |
---|
| 1376 | if (pIntersectchar == 0) // pIntersectSyz::modular |
---|
| 1377 | { |
---|
| 1378 | list L = ringlist(D); |
---|
| 1379 | int lb = 10000; int ub = 536870909; // bounds for random primes |
---|
| 1380 | list usedprimes; |
---|
| 1381 | int sJ = size(J); |
---|
| 1382 | int sLJq; |
---|
| 1383 | ideal LJ; |
---|
| 1384 | for (i=1; i<=sJ; i++) |
---|
| 1385 | { |
---|
| 1386 | LJ[i] = leadcoef(cleardenom(J[i])); |
---|
| 1387 | } |
---|
| 1388 | int short_save = short; // workaround for map Q(a_i) -> Z/q(a_i) |
---|
| 1389 | short = 0; |
---|
| 1390 | string strLJq = "ideal LJq = " + string(LJ) + ";"; |
---|
| 1391 | int nD = nvars(D); |
---|
| 1392 | string L5 = "matrix @C[nD][nD]; @C = " + string(L[5]) + ";"; |
---|
| 1393 | string L6 = "matrix @D[nD][nD]; @D = " + string(L[6]) + ";"; |
---|
| 1394 | L = L[1..4]; |
---|
| 1395 | i = 1; |
---|
| 1396 | while (b == 0) |
---|
| 1397 | { |
---|
| 1398 | dbprint(ppl,"// number of run in the loop: "+string(i)); |
---|
| 1399 | int q = prime(random(lb,ub)); |
---|
| 1400 | if (findFirst(usedprimes,q)==0) // if q has not been used already |
---|
| 1401 | { |
---|
| 1402 | usedprimes = usedprimes,q; |
---|
| 1403 | dbprint(ppl,"// using prime: "+string(q)); |
---|
| 1404 | if (typeof(L[1]) == "int") { L[1] = q; } |
---|
| 1405 | else { L[1][1] = q; } // parameters |
---|
| 1406 | def @Rq = ring(L); setring @Rq; |
---|
| 1407 | execute(L5); execute(L6); |
---|
| 1408 | def Rq = nc_algebra(@C,@D); // def Rq = nc_algebra(1,@D); |
---|
| 1409 | setring Rq; kill @Rq; |
---|
| 1410 | execute(strLJq); |
---|
| 1411 | sLJq = size(LJq); |
---|
| 1412 | setring D; kill Rq; |
---|
| 1413 | if (sLJq <> sJ ) // detect unlucky prime |
---|
| 1414 | { |
---|
| 1415 | dbprint(ppl,"// " +string(q) + " is unlucky"); |
---|
| 1416 | b = 0; |
---|
| 1417 | } |
---|
| 1418 | else |
---|
| 1419 | { |
---|
| 1420 | b = pIntersectSyz(s,J,q,whichengine,modengine); |
---|
| 1421 | } |
---|
| 1422 | } |
---|
| 1423 | i++; |
---|
| 1424 | } |
---|
| 1425 | short = short_save; |
---|
| 1426 | } |
---|
| 1427 | else // pIntersectSyz::non-modular |
---|
| 1428 | { |
---|
| 1429 | b = pIntersectSyz(s,J,0,whichengine); |
---|
| 1430 | } |
---|
| 1431 | } |
---|
| 1432 | else // pIntersect: linReduce |
---|
| 1433 | { |
---|
| 1434 | b = pIntersect(s,J); |
---|
| 1435 | } |
---|
| 1436 | dbprint(ppl,"// ...done"); // with the intersection |
---|
| 1437 | poly pb = vec2poly(b); |
---|
| 1438 | if (inorann == 0) |
---|
| 1439 | { |
---|
| 1440 | pb = subst(pb,var(1),-var(1)-1); |
---|
| 1441 | } |
---|
| 1442 | else // bfctAnn |
---|
| 1443 | { |
---|
| 1444 | if (addPD) |
---|
| 1445 | { |
---|
| 1446 | pb = pb*(var(1)+1); |
---|
| 1447 | } |
---|
| 1448 | } |
---|
| 1449 | list l = bFactor(pb); |
---|
| 1450 | setring save; |
---|
| 1451 | list l = imap(D,l); |
---|
| 1452 | printlevel = printlevelsave; |
---|
| 1453 | return(l); |
---|
| 1454 | } |
---|
| 1455 | |
---|
| 1456 | proc bfct (poly f, list #) |
---|
| 1457 | "USAGE: bfct(f [,s,t,v]); f a poly, s,t optional ints, v an optional intvec |
---|
| 1458 | RETURN: list of ideal and intvec |
---|
| 1459 | PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) |
---|
| 1460 | @* for the hypersurface defined by f. |
---|
| 1461 | ASSUME: The basering is commutative and of characteristic 0. |
---|
| 1462 | BACKGROUND: In this proc, the initial Malgrange ideal is computed according to |
---|
| 1463 | @* the algorithm by Masayuki Noro and then a system of linear equations is |
---|
| 1464 | @* solved by linear reductions. |
---|
| 1465 | NOTE: In the output list, the ideal contains all the roots |
---|
| 1466 | @* and the intvec their multiplicities. |
---|
| 1467 | @* If s<>0, @code{std} is used for GB computations, |
---|
| 1468 | @* otherwise, and by default, @code{slimgb} is used. |
---|
| 1469 | @* If t<>0, a matrix ordering is used for Groebner basis computations, |
---|
| 1470 | @* otherwise, and by default, a block ordering is used. |
---|
| 1471 | @* If v is a positive weight vector, v is used for homogenization |
---|
| 1472 | @* computations, otherwise and by default, no weights are used. |
---|
| 1473 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 1474 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 1475 | EXAMPLE: example bfct; shows examples |
---|
| 1476 | " |
---|
| 1477 | { |
---|
| 1478 | int ppl = printlevel - voice +2; |
---|
| 1479 | int i; |
---|
| 1480 | int n = nvars(basering); |
---|
| 1481 | // in # we have two switches: |
---|
| 1482 | // one for the engine used for Groebner basis computations, |
---|
| 1483 | // one for M() ordering or its realization |
---|
| 1484 | // in # can also be the optional weight vector |
---|
| 1485 | int whichengine = 0; // default |
---|
| 1486 | int methodord = 0; // default |
---|
| 1487 | intvec u0 = 0; // default |
---|
| 1488 | if (size(#)>0) |
---|
| 1489 | { |
---|
| 1490 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 1491 | { |
---|
| 1492 | whichengine = int(#[1]); |
---|
| 1493 | } |
---|
| 1494 | if (size(#)>1) |
---|
| 1495 | { |
---|
| 1496 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 1497 | { |
---|
| 1498 | methodord = int(#[2]); |
---|
| 1499 | } |
---|
| 1500 | if (size(#)>2) |
---|
| 1501 | { |
---|
| 1502 | if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1) |
---|
| 1503 | { |
---|
| 1504 | u0 = #[3]; |
---|
| 1505 | } |
---|
| 1506 | } |
---|
| 1507 | } |
---|
| 1508 | } |
---|
| 1509 | list b = bfctengine(f,0,whichengine,0,0,methodord,0,0,0,u0); |
---|
| 1510 | return(b); |
---|
| 1511 | } |
---|
| 1512 | example |
---|
| 1513 | { |
---|
| 1514 | "EXAMPLE:"; echo = 2; |
---|
| 1515 | ring r = 0,(x,y),dp; |
---|
| 1516 | poly f = x^2+y^3+x*y^2; |
---|
| 1517 | bfct(f); |
---|
| 1518 | intvec v = 3,2; |
---|
| 1519 | bfct(f,1,0,v); |
---|
| 1520 | } |
---|
| 1521 | |
---|
| 1522 | proc bfctSyz (poly f, list #) |
---|
| 1523 | "USAGE: bfctSyz(f [,r,s,t,u,v]); f poly, r,s,t,u optional ints, v opt. intvec |
---|
| 1524 | RETURN: list of ideal and intvec |
---|
| 1525 | PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) |
---|
| 1526 | @* for the hypersurface defined by f |
---|
| 1527 | ASSUME: The basering is commutative and of characteristic 0. |
---|
| 1528 | BACKGROUND: In this proc, the initial Malgrange ideal is computed according to |
---|
| 1529 | @* the algorithm by Masayuki Noro and then a system of linear equations is |
---|
| 1530 | @* solved by computing syzygies. |
---|
| 1531 | NOTE: In the output list, the ideal contains all the roots and the intvec |
---|
| 1532 | @* their multiplicities. |
---|
| 1533 | @* If r<>0, @code{std} is used for GB computations in characteristic 0, |
---|
| 1534 | @* otherwise, and by default, @code{slimgb} is used. |
---|
| 1535 | @* If s<>0, a matrix ordering is used for GB computations, otherwise, |
---|
| 1536 | @* and by default, a block ordering is used. |
---|
| 1537 | @* If t<>0, the computation of the intersection is solely performed over |
---|
| 1538 | @* charasteristic 0, otherwise and by default, a modular method is used. |
---|
| 1539 | @* If u<>0 and by default, @code{std} is used for GB computations in |
---|
| 1540 | @* characteristic >0, otherwise, @code{slimgb} is used. |
---|
| 1541 | @* If v is a positive weight vector, v is used for homogenization |
---|
| 1542 | @* computations, otherwise and by default, no weights are used. |
---|
| 1543 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 1544 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 1545 | EXAMPLE: example bfctSyz; shows examples |
---|
| 1546 | " |
---|
| 1547 | { |
---|
| 1548 | int ppl = printlevel - voice +2; |
---|
| 1549 | int i; |
---|
| 1550 | // in # we have four switches: |
---|
| 1551 | // one for the engine used for Groebner basis computations in char 0, |
---|
| 1552 | // one for M() ordering or its realization |
---|
| 1553 | // one for a modular method when computing the intersection |
---|
| 1554 | // and one for the engine used for Groebner basis computations in char >0 |
---|
| 1555 | // in # can also be the optional weight vector |
---|
| 1556 | int n = nvars(basering); |
---|
| 1557 | int whichengine = 0; // default |
---|
| 1558 | int methodord = 0; // default |
---|
| 1559 | int pIntersectchar = 0; // default |
---|
| 1560 | int modengine = 1; // default |
---|
| 1561 | intvec u0 = 0; // default |
---|
| 1562 | if (size(#)>0) |
---|
| 1563 | { |
---|
| 1564 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 1565 | { |
---|
| 1566 | whichengine = int(#[1]); |
---|
| 1567 | } |
---|
| 1568 | if (size(#)>1) |
---|
| 1569 | { |
---|
| 1570 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 1571 | { |
---|
| 1572 | methodord = int(#[2]); |
---|
| 1573 | } |
---|
| 1574 | if (size(#)>2) |
---|
| 1575 | { |
---|
| 1576 | if (typeof(#[3])=="int" || typeof(#[3])=="number") |
---|
| 1577 | { |
---|
| 1578 | pIntersectchar = int(#[3]); |
---|
| 1579 | } |
---|
| 1580 | if (size(#)>3) |
---|
| 1581 | { |
---|
| 1582 | if (typeof(#[4])=="int" || typeof(#[4])=="number") |
---|
| 1583 | { |
---|
| 1584 | modengine = int(#[4]); |
---|
| 1585 | } |
---|
| 1586 | if (size(#)>4) |
---|
| 1587 | { |
---|
| 1588 | if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1) |
---|
| 1589 | { |
---|
| 1590 | u0 = #[5]; |
---|
| 1591 | } |
---|
| 1592 | } |
---|
| 1593 | } |
---|
| 1594 | } |
---|
| 1595 | } |
---|
| 1596 | } |
---|
| 1597 | list b = bfctengine(f,0,whichengine,0,0,methodord,1,pIntersectchar,modengine,u0); |
---|
| 1598 | return(b); |
---|
| 1599 | } |
---|
| 1600 | example |
---|
| 1601 | { |
---|
| 1602 | "EXAMPLE:"; echo = 2; |
---|
| 1603 | ring r = 0,(x,y),dp; |
---|
| 1604 | poly f = x^2+y^3+x*y^2; |
---|
| 1605 | bfctSyz(f); |
---|
| 1606 | intvec v = 3,2; |
---|
| 1607 | bfctSyz(f,0,1,1,0,v); |
---|
| 1608 | } |
---|
| 1609 | |
---|
| 1610 | proc bfctIdeal (ideal I, intvec w, list #) |
---|
| 1611 | "USAGE: bfctIdeal(I,w[,s,t]); I an ideal, w an intvec, s,t optional ints |
---|
| 1612 | RETURN: list of ideal and intvec |
---|
| 1613 | PURPOSE: computes the roots of the global b-function of I w.r.t. the weight (-w,w). |
---|
| 1614 | ASSUME: The basering is the n-th Weyl algebra in characteristic 0 and for all |
---|
| 1615 | @* 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the |
---|
| 1616 | @* sequence of variables is given by x(1),...,x(n),D(1),...,D(n), |
---|
| 1617 | @* where D(i) is the differential operator belonging to x(i). |
---|
| 1618 | @* Further we assume that I is holonomic. |
---|
| 1619 | BACKGROUND: In this proc, the initial ideal of I is computed according to the |
---|
| 1620 | @* algorithm by Masayuki Noro and then a system of linear equations is |
---|
| 1621 | @* solved by linear reductions. |
---|
| 1622 | NOTE: In the output list, say L, |
---|
| 1623 | @* - L[1] of type ideal contains all the rational roots of a b-function, |
---|
| 1624 | @* - L[2] of type intvec contains the multiplicities of above roots, |
---|
| 1625 | @* - optional L[3] of type string is the part of b-function without rational roots. |
---|
| 1626 | @* Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and |
---|
| 1627 | @* L[3] is 1 (for nonzero constant) or 0 (for zero b-function). |
---|
| 1628 | @* If s<>0, @code{std} is used for GB computations in characteristic 0, |
---|
| 1629 | @* otherwise, and by default, @code{slimgb} is used. |
---|
| 1630 | @* If t<>0, a matrix ordering is used for GB computations, otherwise, |
---|
| 1631 | @* and by default, a block ordering is used. |
---|
| 1632 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 1633 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 1634 | EXAMPLE: example bfctIdeal; shows examples |
---|
| 1635 | " |
---|
| 1636 | { |
---|
| 1637 | int ppl = printlevel - voice +2; |
---|
| 1638 | int i; |
---|
| 1639 | def save = basering; |
---|
| 1640 | int n = nvars(save)/2; |
---|
| 1641 | int whichengine = 0; // default |
---|
| 1642 | int methodord = 0; // default |
---|
| 1643 | if (size(#)>0) |
---|
| 1644 | { |
---|
| 1645 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 1646 | { |
---|
| 1647 | whichengine = int(#[1]); |
---|
| 1648 | } |
---|
| 1649 | if (size(#)>1) |
---|
| 1650 | { |
---|
| 1651 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 1652 | { |
---|
| 1653 | methodord = int(#[2]); |
---|
| 1654 | } |
---|
| 1655 | } |
---|
| 1656 | } |
---|
| 1657 | if (isWeyl()==0) { ERROR("basering is not a Weyl algebra"); } |
---|
| 1658 | for (i=1; i<=n; i++) |
---|
| 1659 | { |
---|
| 1660 | if (bracket(var(i+n),var(i))<>1) |
---|
| 1661 | { |
---|
| 1662 | ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i))); |
---|
| 1663 | } |
---|
| 1664 | } |
---|
| 1665 | int isH = isHolonomic(I); |
---|
| 1666 | if (isH<>1) |
---|
| 1667 | { |
---|
| 1668 | print("WARNING: given ideal is not holonomic"); |
---|
| 1669 | print("... setting bound for degree of b-function to 10 and proceeding"); |
---|
| 1670 | isH = 10; |
---|
| 1671 | } |
---|
| 1672 | else { isH = 0; } // no degree bound for pIntersect |
---|
| 1673 | ideal J = initialIdealW(I,-w,w,whichengine,methodord); |
---|
| 1674 | poly s; |
---|
| 1675 | for (i=1; i<=n; i++) |
---|
| 1676 | { |
---|
| 1677 | s = s + w[i]*var(i)*var(n+i); |
---|
| 1678 | } |
---|
| 1679 | vector b = pIntersect(s,J,isH); |
---|
| 1680 | list RL = ringlist(save); RL = RL[1..4]; |
---|
| 1681 | RL[2] = list(safeVarName("s")); |
---|
| 1682 | RL[3] = list(list("dp",intvec(1)),list("C",intvec(0))); |
---|
| 1683 | def @S = ring(RL); setring @S; |
---|
| 1684 | vector b = imap(save,b); |
---|
| 1685 | poly bs = vec2poly(b); |
---|
| 1686 | list l = bFactor(bs); |
---|
| 1687 | setring save; |
---|
| 1688 | list l = imap(@S,l); |
---|
| 1689 | return(l); |
---|
| 1690 | } |
---|
| 1691 | example |
---|
| 1692 | { |
---|
| 1693 | "EXAMPLE:"; echo = 2; |
---|
| 1694 | ring @D = 0,(x,y,Dx,Dy),dp; |
---|
| 1695 | def D = Weyl(); |
---|
| 1696 | setring D; |
---|
| 1697 | ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6; I = std(I); |
---|
| 1698 | intvec w1 = 0,1; |
---|
| 1699 | intvec w2 = 2,3; |
---|
| 1700 | bfctIdeal(I,w1); |
---|
| 1701 | bfctIdeal(I,w2,0,1); |
---|
| 1702 | ideal J = I[size(I)]; // J is not holonomic by construction |
---|
| 1703 | bfctIdeal(J,w1); // b-function of D/J w.r.t. w1 is non-zero |
---|
| 1704 | bfctIdeal(J,w2); // b-function of D/J w.r.t. w2 is zero |
---|
| 1705 | } |
---|
| 1706 | |
---|
| 1707 | proc bfctOneGB (poly f,list #) |
---|
| 1708 | "USAGE: bfctOneGB(f [,s,t]); f a poly, s,t optional ints |
---|
| 1709 | RETURN: list of ideal and intvec |
---|
| 1710 | PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the |
---|
| 1711 | @* hypersurface defined by f, using only one GB computation |
---|
| 1712 | ASSUME: The basering is commutative and of characteristic 0. |
---|
| 1713 | BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the |
---|
| 1714 | @* algorithm by Masayuki Noro and combined with an elimination ordering. |
---|
| 1715 | NOTE: In the output list, the ideal contains all the roots and the intvec |
---|
| 1716 | @* their multiplicities. |
---|
| 1717 | @* If s<>0, @code{std} is used for the GB computation, otherwise, |
---|
| 1718 | @* and by default, @code{slimgb} is used. |
---|
| 1719 | @* If t<>0, a matrix ordering is used for GB computations, |
---|
| 1720 | @* otherwise, and by default, a block ordering is used. |
---|
| 1721 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 1722 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 1723 | EXAMPLE: example bfctOneGB; shows examples |
---|
| 1724 | " |
---|
| 1725 | { |
---|
| 1726 | int ppl = printlevel - voice +2; |
---|
| 1727 | if (!isCommutative()) { ERROR("Basering must be commutative"); } |
---|
| 1728 | def save = basering; |
---|
| 1729 | int n = nvars(save); |
---|
| 1730 | if (char(save) <> 0) |
---|
| 1731 | { |
---|
| 1732 | ERROR("characteristic of basering has to be 0"); |
---|
| 1733 | } |
---|
| 1734 | list L = ringlist(save); |
---|
| 1735 | int qr; |
---|
| 1736 | if (L[4] <> 0) // qring? |
---|
| 1737 | { |
---|
| 1738 | print("// basering is qring:"); |
---|
| 1739 | print("// discarding the quotient and proceeding..."); |
---|
| 1740 | L[4] = ideal(0); |
---|
| 1741 | qr = 1; |
---|
| 1742 | def save2 = ring(L); |
---|
| 1743 | setring save2; |
---|
| 1744 | poly f = imap(save,f); |
---|
| 1745 | } |
---|
| 1746 | int N = 2*n+4; |
---|
| 1747 | int i; |
---|
| 1748 | int whichengine = 0; // default |
---|
| 1749 | int methodord = 0; // default |
---|
| 1750 | if (size(#)>0) |
---|
| 1751 | { |
---|
| 1752 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 1753 | { |
---|
| 1754 | whichengine = int(#[1]); |
---|
| 1755 | } |
---|
| 1756 | if (size(#)>1) |
---|
| 1757 | { |
---|
| 1758 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 1759 | { |
---|
| 1760 | methodord = int(#[2]); |
---|
| 1761 | } |
---|
| 1762 | } |
---|
| 1763 | } |
---|
| 1764 | // creating the homogenized extended Weyl algebra |
---|
| 1765 | // create names for vars |
---|
| 1766 | list Lvar; |
---|
| 1767 | Lvar[1] = safeVarName("t"); |
---|
| 1768 | Lvar[2] = safeVarName("s"); |
---|
| 1769 | Lvar[n+3] = safeVarName("D"+Lvar[1]); |
---|
| 1770 | Lvar[N] = safeVarName("h"); |
---|
| 1771 | for (i=1; i<=n; i++) |
---|
| 1772 | { |
---|
| 1773 | Lvar[i+2] = string(var(i)); |
---|
| 1774 | Lvar[i+n+3] = safeVarName("D" + string(var(i))); |
---|
| 1775 | } |
---|
| 1776 | // create ordering |
---|
| 1777 | intvec uv = -1; uv[n+3] = 1; uv[N] = 0; |
---|
| 1778 | intvec @a = 1:N; @a[2] = 2; |
---|
| 1779 | intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0; |
---|
| 1780 | list Lord; |
---|
| 1781 | Lord[1] = list("a",@a); Lord[2] = list("a",@a2); |
---|
| 1782 | if (methodord == 0) // default: block ordering |
---|
| 1783 | { |
---|
| 1784 | //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(N-1),lp(1)); |
---|
| 1785 | Lord[3] = list("a",uv); |
---|
| 1786 | Lord[4] = list("dp",intvec(1:(N-1))); |
---|
| 1787 | Lord[5] = list("lp",intvec(1)); |
---|
| 1788 | Lord[6] = list("C",intvec(0)); |
---|
| 1789 | } |
---|
| 1790 | else // M() ordering |
---|
| 1791 | { |
---|
| 1792 | intmat @Ord[N][N]; |
---|
| 1793 | @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1); |
---|
| 1794 | for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; } |
---|
| 1795 | dbprint(ppl,"// weights for ordering: "+string(transpose(@a))); |
---|
| 1796 | dbprint(ppl,"// the ordering matrix:",@Ord); |
---|
| 1797 | //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord)); |
---|
| 1798 | Lord[3] = list("M",intvec(@Ord)); |
---|
| 1799 | Lord[4] = list("C",intvec(0)); |
---|
| 1800 | } |
---|
| 1801 | // create commutative ring |
---|
| 1802 | list L@Dh = ringlist(basering); |
---|
| 1803 | L@Dh = L@Dh[1..4]; // if basering is commutative nc_algebra |
---|
| 1804 | L@Dh[2] = Lvar; L@Dh[3] = Lord; |
---|
| 1805 | def @Dh = ring(L@Dh); setring @Dh; |
---|
| 1806 | dbprint(ppl,"// the ring @Dh:",@Dh); |
---|
| 1807 | // create non-commutative relations |
---|
| 1808 | matrix @relD[N][N]; |
---|
| 1809 | @relD[1,2] = var(1)*var(N)^2; // s*t = t*s + t*h^2 |
---|
| 1810 | @relD[2,n+3] = var(n+3)*var(N)^2; // Dt*s = s*Dt+Dt*h^2 |
---|
| 1811 | @relD[1,n+3] = var(N)^2; |
---|
| 1812 | for (i=1; i<=n; i++) |
---|
| 1813 | { |
---|
| 1814 | @relD[i+2,n+3+i] = var(N)^2; |
---|
| 1815 | } |
---|
| 1816 | dbprint(ppl,"// nc relations:",@relD); |
---|
| 1817 | def Dh = nc_algebra(1,@relD); |
---|
| 1818 | setring Dh; kill @Dh; |
---|
| 1819 | dbprint(ppl,"// computing in ring",Dh); |
---|
| 1820 | poly f = imap(save,f); |
---|
| 1821 | f = homog(f,h); |
---|
| 1822 | // create the Malgrange ideal |
---|
| 1823 | ideal I = var(1) - f, var(1)*var(n+3) - var(2); |
---|
| 1824 | for (i=1; i<=n; i++) |
---|
| 1825 | { |
---|
| 1826 | I[3+i] = var(i+n+3)+diff(f,var(i+2))*var(n+3); |
---|
| 1827 | } |
---|
| 1828 | dbprint(ppl-1, "// the Malgrange ideal: " +string(I)); |
---|
| 1829 | // the hard part: Groebner basis computation |
---|
| 1830 | dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine)); |
---|
| 1831 | I = engine(I, whichengine); |
---|
| 1832 | dbprint(ppl, "// finished Groebner basis computation"); |
---|
| 1833 | I = subst(I,h,1); // dehomogenization |
---|
| 1834 | dbprint(ppl-1,string(I)); |
---|
| 1835 | // 3.3 the initial form |
---|
| 1836 | I = inForm(I,uv); |
---|
| 1837 | dbprint(ppl, "// the initial ideal:", string(matrix(I))); |
---|
| 1838 | // read off the solution |
---|
| 1839 | intvec tonselect = 1; |
---|
| 1840 | for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; } |
---|
| 1841 | I = nselect(I,tonselect); |
---|
| 1842 | dbprint(ppl, "// generators containing only s:", string(matrix(I))); |
---|
| 1843 | I = engine(I, whichengine); // is now a principal ideal; |
---|
| 1844 | if (qr == 1) { setring save2; } |
---|
| 1845 | else { setring save; } |
---|
| 1846 | ideal J; J[2] = var(1); |
---|
| 1847 | map @m = Dh,J; |
---|
| 1848 | ideal I = @m(I); |
---|
| 1849 | poly p = I[1]; |
---|
| 1850 | p = subst(p,var(1),-var(1)-1); |
---|
| 1851 | list l = bFactor(p); |
---|
| 1852 | if (qr == 1) |
---|
| 1853 | { |
---|
| 1854 | setring save; |
---|
| 1855 | list l = imap(save2,l); |
---|
| 1856 | } |
---|
| 1857 | return(l); |
---|
| 1858 | } |
---|
| 1859 | example |
---|
| 1860 | { |
---|
| 1861 | "EXAMPLE:"; echo = 2; |
---|
| 1862 | ring r = 0,(x,y),dp; |
---|
| 1863 | poly f = x^2+y^3+x*y^2; |
---|
| 1864 | bfctOneGB(f); |
---|
| 1865 | bfctOneGB(f,1,1); |
---|
| 1866 | } |
---|
| 1867 | |
---|
| 1868 | |
---|
| 1869 | |
---|
| 1870 | proc bfctAnn (poly f, list #) |
---|
| 1871 | "USAGE: bfctAnn(f [,a,b,c]); f a poly, a, b, c optional ints |
---|
| 1872 | RETURN: list of ideal and intvec |
---|
| 1873 | PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the |
---|
| 1874 | @* hypersurface defined by f. |
---|
| 1875 | ASSUME: The basering is commutative and of characteristic 0. |
---|
| 1876 | BACKGROUND: In this proc, Ann(f^s) is computed and then a system of linear |
---|
| 1877 | @* equations is solved by linear reductions. |
---|
| 1878 | NOTE: In the output list, the ideal contains all the roots and the intvec |
---|
| 1879 | @* their multiplicities. |
---|
| 1880 | @* If a<>0, only f is appended to Ann(f^s), otherwise, and by default, |
---|
| 1881 | @* f and all its partial derivatives are appended. |
---|
| 1882 | @* If b<>0, @code{std} is used for GB computations, otherwise, and by |
---|
| 1883 | @* default, @code{slimgb} is used. |
---|
| 1884 | @* If c<>0, @code{std} is used for Groebner basis computations of ideals |
---|
| 1885 | @* <I+J> when I is already a Groebner basis of <I>. |
---|
| 1886 | @* Otherwise, and by default the engine determined by the switch b is used. |
---|
| 1887 | @* Note that in the case c<>0, the choice for b will be overwritten only |
---|
| 1888 | @* for the types of ideals mentioned above. |
---|
| 1889 | @* This means that if b<>0, specifying c has no effect. |
---|
| 1890 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
| 1891 | @* if printlevel>=2, all the debug messages will be printed. |
---|
| 1892 | EXAMPLE: example bfctAnn; shows examples |
---|
| 1893 | " |
---|
| 1894 | { |
---|
| 1895 | def save = basering; |
---|
| 1896 | int ppl = printlevel - voice + 2; |
---|
| 1897 | int addPD = 1; // default |
---|
| 1898 | int whichengine = 0; // default |
---|
| 1899 | int stdsum = 0; // default |
---|
| 1900 | if (size(#)>0) |
---|
| 1901 | { |
---|
| 1902 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 1903 | { |
---|
| 1904 | addPD = int(#[1]); |
---|
| 1905 | } |
---|
| 1906 | if (size(#)>1) |
---|
| 1907 | { |
---|
| 1908 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 1909 | { |
---|
| 1910 | whichengine = int(#[2]); |
---|
| 1911 | } |
---|
| 1912 | if (size(#)>2) |
---|
| 1913 | { |
---|
| 1914 | if (typeof(#[3])=="int" || typeof(#[3])=="number") |
---|
| 1915 | { |
---|
| 1916 | stdsum = int(#[3]); |
---|
| 1917 | } |
---|
| 1918 | } |
---|
| 1919 | } |
---|
| 1920 | } |
---|
| 1921 | list b = bfctengine(f,1,whichengine,addPD,stdsum,0,0,0,0,0); |
---|
| 1922 | return(b); |
---|
| 1923 | } |
---|
| 1924 | example |
---|
| 1925 | { |
---|
| 1926 | "EXAMPLE:"; echo = 2; |
---|
| 1927 | ring r = 0,(x,y),dp; |
---|
| 1928 | poly f = x^2+y^3+x*y^2; |
---|
| 1929 | bfctAnn(f); |
---|
| 1930 | def R = reiffen(4,5); setring R; |
---|
| 1931 | RC; // the Reiffen curve in 4,5 |
---|
| 1932 | bfctAnn(RC,0,1); |
---|
| 1933 | } |
---|
| 1934 | |
---|
| 1935 | /* |
---|
| 1936 | static proc hardexamples () |
---|
| 1937 | { |
---|
| 1938 | // some hard examples |
---|
| 1939 | ring r1 = 0,(x,y,z,w),dp; |
---|
| 1940 | // ab34 |
---|
| 1941 | poly ab34 = (z3+w4)*(3z2x+4w3y); |
---|
| 1942 | bfct(ab34); |
---|
| 1943 | // ha3 |
---|
| 1944 | poly ha3 = xyzw*(x+y)*(x+z)*(x+w)*(y+z+w); |
---|
| 1945 | bfct(ha3); |
---|
| 1946 | // ha4 |
---|
| 1947 | poly ha4 = xyzw*(x+y)*(x+z)*(x+w)*(y+z)*(y+w); |
---|
| 1948 | bfct(ha4); |
---|
| 1949 | // chal4: reiffen(4,5)*reiffen(5,4) |
---|
| 1950 | ring r2 = 0,(x,y),dp; |
---|
| 1951 | poly chal4 = (x4+xy4+y5)*(x5+x4y+y4); |
---|
| 1952 | bfct(chal4); |
---|
| 1953 | // (xy+z)*reiffen(4,5) |
---|
| 1954 | ring r3 = 0,(x,y,z),dp; |
---|
| 1955 | poly xyzreiffen45 = (xy+z)*(y4+yz4+z5); |
---|
| 1956 | bfct(xyzreiffen45); |
---|
| 1957 | |
---|
| 1958 | // sparse ideal as suggested by Alex; gives 1 as std |
---|
| 1959 | ideal I1 = 28191*y^2+14628*x*Dy, 24865*x^2+24072*x*Dx+17756*Dy^2; |
---|
| 1960 | std(I1); |
---|
| 1961 | } |
---|
| 1962 | */ |
---|