[380a17b] | 1 | ///////////////////////////////////////////////////////////////////// |
---|
[3686937] | 2 | version="version dmodloc.lib 4.0.0.0 Jun_2013 "; // $Id$ |
---|
[1e1ec4] | 3 | category="Noncommutative"; |
---|
| 4 | info=" |
---|
| 5 | LIBRARY: dmodloc.lib Localization of algebraic D-modules and applications |
---|
| 6 | AUTHOR: Daniel Andres, daniel.andres@math.rwth-aachen.de |
---|
| 7 | |
---|
| 8 | Support: DFG Graduiertenkolleg 1632 `Experimentelle und konstruktive Algebra' |
---|
| 9 | |
---|
| 10 | |
---|
| 11 | OVERVIEW: |
---|
| 12 | Let I be a left ideal in the n-th polynomial Weyl algebra D=K[x]<d> and |
---|
| 13 | let f be a polynomial in K[x]. |
---|
| 14 | |
---|
| 15 | If D/I is a holonomic module over D, it is known that the localization of D/I |
---|
| 16 | at f is also holonomic. The procedure @code{Dlocalization} computes an ideal |
---|
| 17 | J in D such that this localization is isomorphic to D/J as D-modules. |
---|
| 18 | |
---|
| 19 | If one regards I as an ideal in the rational Weyl algebra as above, K(x)<d>*I, |
---|
| 20 | and intersects with K[x]<d>, the result is called the Weyl closure of I. |
---|
| 21 | The procedures @code{WeylClosure} (if I has finite holonomic rank) and |
---|
| 22 | @code{WeylClosure1} (if I is in the first Weyl algebra) can be used for |
---|
| 23 | computations. |
---|
| 24 | |
---|
| 25 | As an application of the Weyl closure, the procedure @code{annRatSyz} computes |
---|
| 26 | a holonomic part of the annihilator of a rational function by computing certain |
---|
| 27 | syzygies. The full annihilator can be obtained by taking the Weyl closure of |
---|
| 28 | the result. |
---|
| 29 | |
---|
| 30 | If one regards the left ideal I as system of linear PDEs, one can find its |
---|
| 31 | polynomial solutions with @code{polSol} (if I is holonomic) or |
---|
| 32 | @code{polSolFiniteRank} (if I is of finite holonomic rank). Rational solutions |
---|
| 33 | can be obtained with @code{ratSol}. |
---|
| 34 | |
---|
| 35 | The procedure @code{bfctBound} computes a possible multiple of the b-function |
---|
| 36 | for f^s*u at a generic root of f. Here, u stands for [1] in D/I. |
---|
| 37 | |
---|
| 38 | This library also offers the procedures @code{holonomicRank} and |
---|
| 39 | @code{DsingularLocus} to compute the holonomic rank and the singular locus |
---|
| 40 | of the D-module D/I. |
---|
| 41 | |
---|
| 42 | |
---|
| 43 | REFERENCES: |
---|
| 44 | (OT) T. Oaku, N. Takayama: `Algorithms for D-modules', |
---|
| 45 | Journal of Pure and Applied Algebra, 1998. |
---|
| 46 | @* (OTT) T. Oaku, N. Takayama, H. Tsai: `Polynomial and rational solutions |
---|
| 47 | of holonomic systems', Journal of Pure and Applied Algebra, 2001. |
---|
| 48 | @* (OTW) T. Oaku, N. Takayama, U. Walther: `A Localization Algorithm for |
---|
| 49 | D-modules', Journal of Symbolic Computation, 2000. |
---|
| 50 | @* (Tsa) H. Tsai: `Algorithms for algebraic analysis', PhD thesis, 2000. |
---|
| 51 | |
---|
| 52 | |
---|
| 53 | PROCEDURES: |
---|
| 54 | Dlocalization(I,f[,k,e]); computes the localization of a D-module |
---|
| 55 | WeylClosure(I); computes the Weyl closure of an ideal in the Weyl algebra |
---|
| 56 | WeylClosure1(L); computes the Weyl closure of operator in first Weyl algebra |
---|
| 57 | holonomicRank(I); computes the holonomic rank of I |
---|
| 58 | DsingularLocus(I); computes the singular locus of a D-module |
---|
| 59 | polSol(I[,w,m]); computes basis of polynomial solutions to the given system |
---|
| 60 | polSolFiniteRank(I[,w]); computes basis of polynomial solutions to given system |
---|
| 61 | ratSol(I); computes basis of rational solutions to the given system |
---|
| 62 | bfctBound(I,f[,primdec]); computes multiple of b-function for f^s*u |
---|
| 63 | annRatSyz(f,g[,db,eng]); computes part of annihilator of rational function g/f |
---|
| 64 | |
---|
| 65 | dmodGeneralAssumptionCheck(); check general assumptions |
---|
| 66 | safeVarName(s); finds a free name to use for a new variable |
---|
| 67 | extendWeyl(S); extends basering (Weyl algebra) by given vars |
---|
| 68 | polyVars(f,v); checks whether f contains only variables indexed by v |
---|
| 69 | monomialInIdeal(I); computes all monomials appearing in generators of ideal |
---|
| 70 | vars2pars(v); converts variables specified by v into parameters |
---|
| 71 | minIntRoot2(L); finds minimal integer root in a list of roots |
---|
| 72 | maxIntRoot(L); finds maximal integer root in a list of roots |
---|
| 73 | dmodAction(id,f[,v]); computes the natural action of a D-module on K[x] |
---|
| 74 | dmodActionRat(id,w); computes the natural action of a D-module on K(x) |
---|
| 75 | simplifyRat(v); simplifies rational function |
---|
| 76 | addRat(v,w); adds rational functions |
---|
| 77 | multRat(v,w); multiplies rational functions |
---|
| 78 | diffRat(v,j); derives rational function |
---|
| 79 | commRing(); deletes non-commutative relations from ring |
---|
| 80 | rightNFWeyl(id,k); computes right NF wrt right ideal (var(k)) in Weyl algebra |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | KEYWORDS: D-module; holonomic rank; singular locus of D-module; |
---|
| 84 | D-localization; localization of D-module; characteristic variety; |
---|
| 85 | Weyl closure; polynomial solutions; rational solutions; |
---|
| 86 | annihilator of rational function |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | SEE ALSO: bfun_lib, dmod_lib, dmodapp_lib, dmodvar_lib, gmssing_lib |
---|
| 90 | "; |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | /* |
---|
| 94 | CHANGELOG: |
---|
| 95 | 12.11.12: bugfixes, updated docu |
---|
| 96 | 17.12.12: updated docu, removed redundant procedure killTerms |
---|
| 97 | */ |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | LIB "bfun.lib"; // for pIntersect etc |
---|
| 101 | LIB "dmodapp.lib"; // for GBWeight, charVariety etc |
---|
| 102 | LIB "nctools.lib"; // for Weyl, isWeyl etc |
---|
| 103 | // TODO uncomment this once chern.lib is ready |
---|
| 104 | // LIB "chern.lib"; // for orderedPartition |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | // testing for consistency of the library ///////////////////////////////////// |
---|
| 108 | |
---|
| 109 | static proc testdmodloc() |
---|
| 110 | { |
---|
| 111 | example dmodGeneralAssumptionCheck; |
---|
| 112 | example safeVarName; |
---|
| 113 | example extendWeyl; |
---|
| 114 | example polyVars; |
---|
| 115 | example monomialInIdeal; |
---|
| 116 | example vars2pars; |
---|
| 117 | example minIntRoot2; |
---|
| 118 | example maxIntRoot; |
---|
| 119 | example dmodAction; |
---|
| 120 | example dmodActionRat; |
---|
| 121 | example simplifyRat; |
---|
| 122 | example addRat; |
---|
| 123 | example multRat; |
---|
| 124 | example diffRat; |
---|
| 125 | example commRing; |
---|
| 126 | example holonomicRank; |
---|
| 127 | example DsingularLocus; |
---|
| 128 | example rightNFWeyl; |
---|
| 129 | example Dlocalization; |
---|
| 130 | example WeylClosure1; |
---|
| 131 | example WeylClosure; |
---|
| 132 | example polSol; |
---|
| 133 | example polSolFiniteRank; |
---|
| 134 | example ratSol; |
---|
| 135 | example bfctBound; |
---|
| 136 | example annRatSyz; |
---|
| 137 | } |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | // tools ////////////////////////////////////////////////////////////////////// |
---|
| 141 | |
---|
| 142 | proc dmodGeneralAssumptionCheck () |
---|
| 143 | " |
---|
| 144 | USAGE: dmodGeneralAssumptionCheck(); |
---|
| 145 | RETURN: nothing, but checks general assumptions on the basering |
---|
| 146 | NOTE: This procedure checks the following conditions on the basering R |
---|
| 147 | and prints an error message if any of them is violated: |
---|
| 148 | @* - R is the n-th Weyl algebra over a field of characteristic 0, |
---|
| 149 | @* - R is not a qring, |
---|
| 150 | @* - for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 |
---|
| 151 | holds, i.e. the sequence of variables is given by |
---|
| 152 | x(1),...,x(n),D(1),...,D(n), where D(i) is the differential |
---|
| 153 | operator belonging to x(i). |
---|
| 154 | EXAMPLE: example dmodGeneralAssumptionCheck; shows examples |
---|
| 155 | " |
---|
| 156 | { |
---|
| 157 | // char K <> 0, qring |
---|
| 158 | if ( (size(ideal(basering)) >0) || (char(basering) >0) ) |
---|
| 159 | { |
---|
| 160 | ERROR("Basering is inappropriate: characteristic>0 or qring present"); |
---|
| 161 | } |
---|
| 162 | // no Weyl algebra |
---|
| 163 | if (isWeyl() == 0) |
---|
| 164 | { |
---|
| 165 | ERROR("Basering is not a Weyl algebra"); |
---|
| 166 | } |
---|
| 167 | // wrong sequence of vars |
---|
| 168 | int i,n; |
---|
| 169 | n = nvars(basering) div 2; |
---|
| 170 | for (i=1; i<=n; i++) |
---|
| 171 | { |
---|
| 172 | if (bracket(var(i+n),var(i))<>1) |
---|
| 173 | { |
---|
| 174 | ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i))); |
---|
| 175 | } |
---|
| 176 | } |
---|
| 177 | return(); |
---|
| 178 | } |
---|
| 179 | example |
---|
| 180 | { |
---|
| 181 | "EXAMPLE"; echo=2; |
---|
| 182 | ring r = 0,(x,D),dp; |
---|
| 183 | dmodGeneralAssumptionCheck(); // prints error message |
---|
| 184 | def W = Weyl(); |
---|
| 185 | setring W; |
---|
| 186 | dmodGeneralAssumptionCheck(); // returns nothing |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | |
---|
| 190 | proc safeVarName (string s) |
---|
| 191 | " |
---|
| 192 | USAGE: safeVarName(s); s string |
---|
| 193 | RETURN: string, returns s if s is not the name of a par/var of basering |
---|
| 194 | and `@' + s otherwise |
---|
| 195 | EXAMPLE: example safeVarName; shows examples |
---|
| 196 | " |
---|
| 197 | { |
---|
| 198 | string S = "," + charstr(basering) + "," + varstr(basering) + ","; |
---|
| 199 | s = "," + s + ","; |
---|
| 200 | while (find(S,s) <> 0) |
---|
| 201 | { |
---|
| 202 | s[1] = "@"; |
---|
| 203 | s = "," + s; |
---|
| 204 | } |
---|
| 205 | s = s[2..size(s)-1]; |
---|
| 206 | return(s); |
---|
| 207 | } |
---|
| 208 | example |
---|
| 209 | { |
---|
| 210 | "EXAMPLE:"; echo = 2; |
---|
| 211 | ring r = (0,a),(w,@w,x,y),dp; |
---|
| 212 | safeVarName("a"); |
---|
| 213 | safeVarName("x"); |
---|
| 214 | safeVarName("z"); |
---|
| 215 | safeVarName("w"); |
---|
| 216 | } |
---|
| 217 | |
---|
| 218 | |
---|
| 219 | proc extendWeyl (def newVars) |
---|
| 220 | " |
---|
| 221 | USAGE: extendWeyl(S); S string or list of strings |
---|
| 222 | ASSUME: The basering is the n-th Weyl algebra over a field of |
---|
| 223 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 224 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 225 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) |
---|
| 226 | is the differential operator belonging to x(i). |
---|
| 227 | RETURN: ring, Weyl algebra extended by vars given by S |
---|
| 228 | EXAMPLE: example extendWeyl; shows examples |
---|
| 229 | " |
---|
| 230 | { |
---|
| 231 | dmodGeneralAssumptionCheck(); |
---|
| 232 | int i,s; |
---|
| 233 | string inpt = typeof(newVars); |
---|
| 234 | list L; |
---|
| 235 | if (inpt=="string") |
---|
| 236 | { |
---|
| 237 | s = 1; |
---|
| 238 | L = newVars; |
---|
| 239 | } |
---|
| 240 | else |
---|
| 241 | { |
---|
| 242 | if (inpt=="list") |
---|
| 243 | { |
---|
| 244 | s = size(newVars); |
---|
| 245 | if (s<1) |
---|
| 246 | { |
---|
| 247 | ERROR("No new variables specified."); |
---|
| 248 | } |
---|
| 249 | for (i=1; i<=s; i++) |
---|
| 250 | { |
---|
| 251 | if (typeof(newVars[i]) <> "string") |
---|
| 252 | { |
---|
| 253 | ERROR("Entries of input list must be of type string."); |
---|
| 254 | } |
---|
| 255 | } |
---|
| 256 | L = newVars; |
---|
| 257 | } |
---|
| 258 | else |
---|
| 259 | { |
---|
| 260 | ERROR("Expected string or list of strings as input."); |
---|
| 261 | } |
---|
| 262 | } |
---|
| 263 | def save = basering; |
---|
| 264 | int n = nvars(save) div 2; |
---|
| 265 | list RL = ringlist(save); |
---|
| 266 | RL = RL[1..4]; |
---|
| 267 | list Ltemp = L; |
---|
| 268 | for (i=s; i>0; i--) |
---|
| 269 | { |
---|
| 270 | Ltemp[n+s+i] = "D" + newVars[i]; |
---|
| 271 | } |
---|
| 272 | for (i=n; i>0; i--) |
---|
| 273 | { |
---|
| 274 | Ltemp[s+i] = RL[2][i]; |
---|
| 275 | Ltemp[n+2*s+i] = RL[2][n+i]; |
---|
| 276 | } |
---|
| 277 | RL[2] = Ltemp; |
---|
| 278 | Ltemp = list(); |
---|
| 279 | Ltemp[1] = list("dp",intvec(1:(2*n+2*s))); |
---|
| 280 | Ltemp[2] = list("C",intvec(0)); |
---|
| 281 | RL[3] = Ltemp; |
---|
| 282 | kill Ltemp; |
---|
| 283 | def @Dv = ring(RL); |
---|
| 284 | setring @Dv; |
---|
| 285 | def Dv = Weyl(); |
---|
| 286 | setring save; |
---|
| 287 | return(Dv); |
---|
| 288 | } |
---|
| 289 | example |
---|
| 290 | { |
---|
| 291 | "EXAMPLE:"; echo = 2; |
---|
| 292 | ring @D2 = 0,(x,y,Dx,Dy),dp; |
---|
| 293 | def D2 = Weyl(); |
---|
| 294 | setring D2; |
---|
| 295 | def D3 = extendWeyl("t"); |
---|
| 296 | setring D3; D3; |
---|
| 297 | list L = "u","v"; |
---|
| 298 | def D5 = extendWeyl(L); |
---|
| 299 | setring D5; |
---|
| 300 | D5; |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | |
---|
| 304 | proc polyVars (poly f, intvec v) |
---|
| 305 | " |
---|
| 306 | USAGE: polyVars(f,v); f poly, v intvec |
---|
| 307 | RETURN: int, 1 if f contains only variables indexed by v, 0 otherwise |
---|
| 308 | EXAMPLE: example polyVars; shows examples |
---|
| 309 | " |
---|
| 310 | { |
---|
| 311 | ideal varsf = variables(f); // vars contained in f |
---|
| 312 | ideal V; |
---|
| 313 | int i; |
---|
| 314 | int n = nvars(basering); |
---|
| 315 | for (i=1; i<=nrows(v); i++) |
---|
| 316 | { |
---|
| 317 | if ( (v[i]<0) || (v[i]>n) ) |
---|
| 318 | { |
---|
| 319 | ERROR("var(" + string(v[i]) + ") out of range"); |
---|
| 320 | } |
---|
| 321 | V[i] = var(v[i]); |
---|
| 322 | } |
---|
| 323 | attrib(V,"isSB",1); |
---|
| 324 | ideal N = NF(varsf,V); |
---|
| 325 | N = simplify(N,2); |
---|
| 326 | if (N[1]==0) |
---|
| 327 | { |
---|
| 328 | return(1); |
---|
| 329 | } |
---|
| 330 | else |
---|
| 331 | { |
---|
| 332 | return(0); |
---|
| 333 | } |
---|
| 334 | } |
---|
| 335 | example |
---|
| 336 | { |
---|
| 337 | "EXAMPLE:"; echo = 2; |
---|
| 338 | ring r = 0,(x,y,z),dp; |
---|
| 339 | poly f = y^2+zy; |
---|
| 340 | intvec v = 1,2; |
---|
| 341 | polyVars(f,v); // does f depend only on x,y? |
---|
| 342 | v = 2,3; |
---|
| 343 | polyVars(f,v); // does f depend only on y,z? |
---|
| 344 | } |
---|
| 345 | |
---|
| 346 | |
---|
| 347 | proc monomialInIdeal (ideal I) |
---|
| 348 | " |
---|
| 349 | USAGE: monomialInIdeal(I); I ideal |
---|
| 350 | RETURN: ideal consisting of all monomials appearing in generators of ideal |
---|
| 351 | EXAMLPE: example monomialInIdeal; shows examples |
---|
| 352 | " |
---|
| 353 | { |
---|
| 354 | // returns ideal consisting of all monomials appearing in generators of ideal |
---|
| 355 | I = simplify(I,2+8); |
---|
| 356 | int i; |
---|
| 357 | poly p; |
---|
| 358 | ideal M; |
---|
| 359 | for (i=1; i<=size(I); i++) |
---|
| 360 | { |
---|
| 361 | p = I[i]; |
---|
| 362 | while (p<>0) |
---|
| 363 | { |
---|
| 364 | M[size(M)+1] = leadmonom(p); |
---|
| 365 | p = p - lead(p); |
---|
| 366 | } |
---|
| 367 | } |
---|
| 368 | M = simplify(M,4+2); |
---|
| 369 | return(M); |
---|
| 370 | } |
---|
| 371 | example |
---|
| 372 | { |
---|
| 373 | "EXAMPLE"; echo=2; |
---|
| 374 | ring r = 0,(x,y),dp; |
---|
| 375 | ideal I = x2+5x3y7, x-x2-6xy; |
---|
| 376 | monomialInIdeal(I); |
---|
| 377 | } |
---|
| 378 | |
---|
| 379 | |
---|
| 380 | proc vars2pars (intvec v) |
---|
| 381 | " |
---|
| 382 | USAGE: vars2pars(v); v intvec |
---|
| 383 | ASSUME: The basering is commutative. |
---|
| 384 | RETURN: ring with variables specified by v converted into parameters |
---|
| 385 | EXAMPLE: example vars2pars; shows examples |
---|
| 386 | " |
---|
| 387 | { |
---|
| 388 | if (isCommutative() == 0) |
---|
| 389 | { |
---|
| 390 | ERROR("The basering must be commutative."); |
---|
| 391 | } |
---|
| 392 | v = sortIntvec(v)[1]; |
---|
| 393 | int sv = size(v); |
---|
| 394 | if ( (v[1]<1) || (v[sv]<1) ) |
---|
| 395 | { |
---|
| 396 | ERROR("Expected entries of intvec in the range 1.."+string(n)); |
---|
| 397 | } |
---|
| 398 | def save = basering; |
---|
| 399 | int i,j,n; |
---|
| 400 | n = nvars(save); |
---|
| 401 | list RL = ringlist(save); |
---|
| 402 | list Lp,Lv,L1; |
---|
| 403 | if (typeof(RL[1]) == "list") |
---|
| 404 | { |
---|
| 405 | L1 = RL[1]; |
---|
| 406 | Lp = L1[2]; |
---|
| 407 | } |
---|
| 408 | else |
---|
| 409 | { |
---|
| 410 | L1[1] = RL[1]; |
---|
| 411 | L1[4] = ideal(0); |
---|
| 412 | } |
---|
| 413 | j = sv; |
---|
| 414 | for (i=1; i<=n; i++) |
---|
| 415 | { |
---|
| 416 | if (j>0) |
---|
| 417 | { |
---|
| 418 | if (v[j]==i) |
---|
| 419 | { |
---|
| 420 | Lp[size(Lp)+1] = string(var(i)); |
---|
| 421 | j--; |
---|
| 422 | } |
---|
| 423 | else |
---|
| 424 | { |
---|
| 425 | Lv[size(Lv)+1] = string(var(i)); |
---|
| 426 | } |
---|
| 427 | } |
---|
| 428 | else |
---|
| 429 | { |
---|
| 430 | Lv[size(Lv)+1] = string(var(i)); |
---|
| 431 | } |
---|
| 432 | } |
---|
| 433 | RL[2] = Lv; |
---|
| 434 | L1[2] = Lp; |
---|
| 435 | L1[3] = list(list("lp",intvec(1:size(Lp)))); |
---|
| 436 | RL[1] = L1; |
---|
| 437 | L1 = list(); |
---|
| 438 | L1[1] = list("dp",intvec(1:sv)); |
---|
| 439 | L1[2] = list("C",intvec(0)); |
---|
| 440 | RL[3] = L1; |
---|
| 441 | // RL; |
---|
| 442 | def R = ring(RL); |
---|
| 443 | return(R); |
---|
| 444 | } |
---|
| 445 | example |
---|
| 446 | { |
---|
| 447 | "EXAMPLE:"; echo = 2; |
---|
| 448 | ring r = 0,(x,y,z,a,b,c),dp; |
---|
| 449 | intvec v = 4,5,6; |
---|
| 450 | def R = vars2pars(v); |
---|
| 451 | setring R; |
---|
| 452 | R; |
---|
| 453 | v = 1,2; |
---|
| 454 | def RR = vars2pars(v); |
---|
| 455 | setring RR; |
---|
| 456 | RR; |
---|
| 457 | } |
---|
| 458 | |
---|
| 459 | |
---|
| 460 | static proc minMaxIntRoot (list L, string minmax) |
---|
| 461 | { |
---|
| 462 | int win; |
---|
| 463 | if (size(L)>1) |
---|
| 464 | { |
---|
| 465 | if ( (typeof(L[1])<>"ideal") || (typeof(L[2])<>"intvec") ) |
---|
| 466 | { |
---|
| 467 | win = 1; |
---|
| 468 | } |
---|
| 469 | } |
---|
| 470 | else |
---|
| 471 | { |
---|
| 472 | win = 1; |
---|
| 473 | } |
---|
| 474 | if (win) |
---|
| 475 | { |
---|
| 476 | ERROR("Expected list in the format of bFactor"); |
---|
| 477 | } |
---|
| 478 | if (size(L)>2) |
---|
| 479 | { |
---|
| 480 | if ( (L[3]=="1") || (L[3]=="0") ) |
---|
| 481 | { |
---|
| 482 | print("// Warning: Constant poly. Returning 0."); |
---|
| 483 | return(int(0)); |
---|
| 484 | } |
---|
| 485 | } |
---|
| 486 | ideal I = L[1]; |
---|
| 487 | int i,k,b; |
---|
| 488 | if (minmax=="min") |
---|
| 489 | { |
---|
| 490 | i = ncols(I); |
---|
| 491 | k = -1; |
---|
| 492 | b = 0; |
---|
| 493 | } |
---|
| 494 | else // minmax=="max" |
---|
| 495 | { |
---|
| 496 | i = 1; |
---|
| 497 | k = 1; |
---|
| 498 | b = ncols(I); |
---|
| 499 | } |
---|
| 500 | for (; k*i<k*b; i=i+k) |
---|
| 501 | { |
---|
| 502 | if (isInt(leadcoef(I[i]))) |
---|
| 503 | { |
---|
| 504 | return(int(leadcoef(I[i]))); |
---|
| 505 | } |
---|
| 506 | } |
---|
| 507 | print("// Warning: No integer root found. Returning 0."); |
---|
| 508 | return(int(0)); |
---|
| 509 | } |
---|
| 510 | |
---|
| 511 | |
---|
| 512 | //TODO rename? minIntRoot is name of proc in dmod.lib |
---|
| 513 | proc minIntRoot2 (list L) |
---|
| 514 | " |
---|
| 515 | USAGE: minIntRoot2(L); L list |
---|
| 516 | ASSUME: L is the output of bFactor. |
---|
| 517 | RETURN: int, the minimal integer root in a list of roots |
---|
| 518 | SEE ALSO: minIntRoot, maxIntRoot, bFactor |
---|
| 519 | EXAMPLE: example minIntRoot2; shows examples |
---|
| 520 | " |
---|
| 521 | { |
---|
| 522 | return(minMaxIntRoot(L,"min")); |
---|
| 523 | } |
---|
| 524 | example |
---|
| 525 | { |
---|
| 526 | "EXAMPLE"; echo=2; |
---|
| 527 | ring r = 0,x,dp; |
---|
| 528 | poly f = x*(x+1)*(x-2)*(x-5/2)*(x+5/2); |
---|
| 529 | list L = bFactor(f); |
---|
| 530 | minIntRoot2(L); |
---|
| 531 | } |
---|
| 532 | |
---|
| 533 | |
---|
| 534 | proc maxIntRoot (list L) |
---|
| 535 | " |
---|
| 536 | USAGE: maxIntRoot(L); L list |
---|
| 537 | ASSUME: L is the output of bFactor. |
---|
| 538 | RETURN: int, the maximal integer root in a list of roots |
---|
| 539 | SEE ALSO: minIntRoot2, bFactor |
---|
| 540 | EXAMPLE: example maxIntRoot; shows examples |
---|
| 541 | " |
---|
| 542 | { |
---|
| 543 | return(minMaxIntRoot(L,"max")); |
---|
| 544 | } |
---|
| 545 | example |
---|
| 546 | { |
---|
| 547 | "EXAMPLE"; echo=2; |
---|
| 548 | ring r = 0,x,dp; |
---|
| 549 | poly f = x*(x+1)*(x-2)*(x-5/2)*(x+5/2); |
---|
| 550 | list L = bFactor(f); |
---|
| 551 | maxIntRoot(L); |
---|
| 552 | } |
---|
| 553 | |
---|
| 554 | |
---|
| 555 | proc dmodAction (def id, poly f, list #) |
---|
| 556 | " |
---|
| 557 | USAGE: dmodAction(id,f[,v]); id ideal or poly, f poly, v optional intvec |
---|
| 558 | ASSUME: If v is not given, the basering is the n-th Weyl algebra W over a |
---|
| 559 | field of characteristic 0 and for all 1<=i<=n the identity |
---|
| 560 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 561 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the |
---|
| 562 | differential operator belonging to x(i). |
---|
| 563 | Otherwise, v is assumed to specify positions of variables, which form |
---|
| 564 | a Weyl algebra as a subalgebra of the basering: |
---|
| 565 | If size(v) equals 2*n, then bracket(var(v[i]),var(v[j])) must equal |
---|
| 566 | 1 if and only if j equals i+n, and 0 otherwise, for all 1<=i,j<=n. |
---|
| 567 | @* Further, assume that f does not contain any D(i). |
---|
| 568 | RETURN: same type as id, the result of the natural D-module action of id on f |
---|
| 569 | NOTE: The assumptions made are not checked. |
---|
| 570 | EXAMPLE: example dmodAction; shows examples |
---|
| 571 | " |
---|
| 572 | { |
---|
| 573 | string inp1 = typeof(id); |
---|
| 574 | if ((inp1<>"poly") && (inp1<>"ideal")) |
---|
| 575 | { |
---|
| 576 | ERROR("Expected first argument to be poly or ideal but received "+inp1); |
---|
| 577 | } |
---|
| 578 | intvec posXD = 1..nvars(basering); |
---|
| 579 | if (size(#)>0) |
---|
| 580 | { |
---|
| 581 | if (typeof(#[1])=="intvec") |
---|
| 582 | { |
---|
| 583 | posXD = #[1]; |
---|
| 584 | } |
---|
| 585 | } |
---|
| 586 | if ((size(posXD) mod 2)<>0) |
---|
| 587 | { |
---|
| 588 | ERROR("Even number of variables expected.") |
---|
| 589 | } |
---|
| 590 | int n = (size(posXD)) div 2; |
---|
| 591 | int i,j,k,l; |
---|
| 592 | ideal resI = id; |
---|
| 593 | int sid = ncols(resI); |
---|
| 594 | intvec v; |
---|
| 595 | poly P,h; |
---|
| 596 | for (l=1; l<=sid; l++) |
---|
| 597 | { |
---|
| 598 | P = resI[l]; |
---|
| 599 | resI[l] = 0; |
---|
| 600 | for (i=1; i<=size(P); i++) |
---|
| 601 | { |
---|
| 602 | v = leadexp(P[i]); |
---|
| 603 | h = f; |
---|
| 604 | for (j=1; j<=n; j++) |
---|
| 605 | { |
---|
| 606 | for (k=1; k<=v[posXD[j+n]]; k++) // action of Dx |
---|
| 607 | { |
---|
| 608 | h = diff(h,var(posXD[j])); |
---|
| 609 | } |
---|
| 610 | h = h*var(posXD[j])^v[posXD[j]]; // action of x |
---|
| 611 | } |
---|
| 612 | h = leadcoef(P[i])*h; |
---|
| 613 | resI[l] = resI[l] + h; |
---|
| 614 | } |
---|
| 615 | } |
---|
| 616 | if (inp1 == "ideal") |
---|
| 617 | { |
---|
| 618 | return(resI); |
---|
| 619 | } |
---|
| 620 | else |
---|
| 621 | { |
---|
| 622 | return(resI[1]); |
---|
| 623 | } |
---|
| 624 | } |
---|
| 625 | example |
---|
| 626 | { |
---|
| 627 | ring r = 0,(x,y,z),dp; |
---|
| 628 | poly f = x^2*z - y^3; |
---|
| 629 | def A = annPoly(f); |
---|
| 630 | setring A; |
---|
| 631 | poly f = imap(r,f); |
---|
| 632 | dmodAction(LD,f); |
---|
| 633 | poly P = y*Dy+3*z*Dz-3; |
---|
| 634 | dmodAction(P,f); |
---|
| 635 | dmodAction(P[1],f); |
---|
| 636 | } |
---|
| 637 | |
---|
| 638 | |
---|
| 639 | static proc checkRatInput (vector I) |
---|
| 640 | { |
---|
| 641 | // check for valid input |
---|
| 642 | int wrginpt; |
---|
| 643 | if (nrows(I)<>2) |
---|
| 644 | { |
---|
| 645 | wrginpt = 1; |
---|
| 646 | } |
---|
| 647 | else |
---|
| 648 | { |
---|
| 649 | if (I[2] == 0) |
---|
| 650 | { |
---|
| 651 | wrginpt = 1; |
---|
| 652 | } |
---|
| 653 | } |
---|
| 654 | if (wrginpt) |
---|
| 655 | { |
---|
| 656 | ERROR("Vector must consist of exactly two components, second one not 0"); |
---|
| 657 | } |
---|
| 658 | return(); |
---|
| 659 | } |
---|
| 660 | |
---|
| 661 | |
---|
| 662 | proc dmodActionRat(def id, vector w) |
---|
| 663 | " |
---|
| 664 | USAGE: dmodActionRat(id,w); id ideal or poly, f vector |
---|
| 665 | ASSUME: The basering is the n-th Weyl algebra W over a field of |
---|
| 666 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 667 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 668 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the |
---|
| 669 | differential operator belonging to x(i). |
---|
| 670 | @* Further, assume that w has exactly two components, second one not 0, |
---|
| 671 | and that w does not contain any D(i). |
---|
| 672 | RETURN: same type as id, the result of the natural D-module action of id on |
---|
| 673 | the rational function w[1]/w[2] |
---|
| 674 | EXAMPLE: example dmodActionRat; shows examples |
---|
| 675 | " |
---|
| 676 | { |
---|
| 677 | string inp1 = typeof(id); |
---|
| 678 | if ( (inp1<>"poly") && (inp1<>"ideal") ) |
---|
| 679 | { |
---|
| 680 | ERROR("Expected first argument to be poly or ideal but received " + inp1); |
---|
| 681 | } |
---|
| 682 | checkRatInput(w); |
---|
| 683 | poly f = w[1]; |
---|
| 684 | finKx(f); |
---|
| 685 | f = w[2]; |
---|
| 686 | finKx(f); |
---|
| 687 | def save = basering; |
---|
| 688 | def r = commRing(); |
---|
| 689 | setring r; |
---|
| 690 | ideal I = imap(save,id); |
---|
| 691 | vector w = imap(save,w); |
---|
| 692 | int i,j,k,l; |
---|
| 693 | int n = nvars(basering) div 2; |
---|
| 694 | int sid = ncols(I); |
---|
| 695 | intvec v; |
---|
| 696 | poly P; |
---|
| 697 | vector h,resT; |
---|
| 698 | module resL; |
---|
| 699 | for (l=1; l<=sid; l++) |
---|
| 700 | { |
---|
| 701 | P = I[l]; |
---|
| 702 | resT = [0,1]; |
---|
| 703 | for (i=1; i<=size(P); i++) |
---|
| 704 | { |
---|
| 705 | v = leadexp(P[i]); |
---|
| 706 | h = w; |
---|
| 707 | for (j=1; j<=n; j++) |
---|
| 708 | { |
---|
| 709 | for (k=1; k<=v[j+n]; k++) // action of Dx |
---|
| 710 | { |
---|
| 711 | h = diffRat(h,j); |
---|
| 712 | } |
---|
| 713 | h = h + h[1]*(var(j)^v[j]-1)*gen(1); // action of x |
---|
| 714 | } |
---|
| 715 | h = h + (leadcoef(P[i])-1)*h[1]*gen(1); |
---|
| 716 | resT = addRat(resT,h); |
---|
| 717 | } |
---|
| 718 | resL[l] = resT; |
---|
| 719 | } |
---|
| 720 | setring save; |
---|
| 721 | module resL = imap(r,resL); |
---|
| 722 | return(resL); |
---|
| 723 | } |
---|
| 724 | example |
---|
| 725 | { |
---|
| 726 | "EXAMPLE:"; echo = 2; |
---|
| 727 | ring r = 0,(x,y),dp; |
---|
| 728 | poly f = 2*x; poly g = y; |
---|
| 729 | def A = annRat(f,g); setring A; |
---|
| 730 | poly f = imap(r,f); poly g = imap(r,g); |
---|
| 731 | vector v = [f,g]; // represents f/g |
---|
| 732 | // x and y act by multiplication |
---|
| 733 | dmodActionRat(x,v); |
---|
| 734 | dmodActionRat(y,v); |
---|
| 735 | // Dx and Dy act by partial derivation |
---|
| 736 | dmodActionRat(Dx,v); |
---|
| 737 | dmodActionRat(Dy,v); |
---|
| 738 | dmodActionRat(x*Dx+y*Dy,v); |
---|
| 739 | setring r; |
---|
| 740 | f = 2*x*y; g = x^2 - y^3; |
---|
| 741 | def B = annRat(f,g); setring B; |
---|
| 742 | poly f = imap(r,f); poly g = imap(r,g); |
---|
| 743 | vector v = [f,g]; |
---|
| 744 | dmodActionRat(LD,v); // hence LD is indeed the annihilator of f/g |
---|
| 745 | } |
---|
| 746 | |
---|
| 747 | |
---|
| 748 | static proc arithmeticRat (vector I, vector J, string op, list #) |
---|
| 749 | { |
---|
| 750 | // op = "+": return I+J |
---|
| 751 | // op = "*": return I*J |
---|
| 752 | // op = "s": return simplified I |
---|
| 753 | // op = "d": return diff(I,var(#[1])) |
---|
| 754 | int isComm = isCommutative(); |
---|
| 755 | if (!isComm) |
---|
| 756 | { |
---|
| 757 | def save = basering; |
---|
| 758 | def r = commRing(); |
---|
| 759 | setring r; |
---|
| 760 | ideal m = maxideal(1); |
---|
| 761 | map f = save,m; |
---|
| 762 | vector I = f(I); |
---|
| 763 | vector J = f(J); |
---|
| 764 | } |
---|
| 765 | vector K; |
---|
| 766 | poly p; |
---|
| 767 | if (op == "s") |
---|
| 768 | { |
---|
| 769 | p = gcd(I[1],I[2]); |
---|
| 770 | K = (I[1]/p)*gen(1) + (I[2]/p)*gen(2); |
---|
| 771 | } |
---|
| 772 | else |
---|
| 773 | { |
---|
| 774 | if (op == "+") |
---|
| 775 | { |
---|
| 776 | I = arithmeticRat(I,vector(0),"s"); |
---|
| 777 | J = arithmeticRat(J,vector(0),"s"); |
---|
| 778 | p = lcm(I[2],J[2]); |
---|
| 779 | K = (I[1]*p/I[2] + J[1]*p/J[2])*gen(1) + p*gen(2); |
---|
| 780 | } |
---|
| 781 | else |
---|
| 782 | { |
---|
| 783 | if (op == "*") |
---|
| 784 | { |
---|
| 785 | K = (I[1]*J[1])*gen(1) + (I[2]*J[2])*gen(2); |
---|
| 786 | } |
---|
| 787 | else |
---|
| 788 | { |
---|
| 789 | if (op == "d") |
---|
| 790 | { |
---|
| 791 | int j = #[1]; |
---|
| 792 | K = (diff(I[1],var(j))*I[2] - I[1]*diff(I[2],var(j)))*gen(1)+ (I[2]^2)*gen(2); |
---|
| 793 | } |
---|
| 794 | } |
---|
| 795 | } |
---|
| 796 | K = arithmeticRat(K,vector(0),"s"); |
---|
| 797 | } |
---|
| 798 | if (!isComm) |
---|
| 799 | { |
---|
| 800 | setring save; |
---|
| 801 | vector K = imap(r,K); |
---|
| 802 | } |
---|
| 803 | return(K); |
---|
| 804 | } |
---|
| 805 | |
---|
| 806 | |
---|
| 807 | proc simplifyRat (vector J) |
---|
| 808 | " |
---|
| 809 | USAGE: simplifyRat(v); v vector |
---|
| 810 | ASSUME: Assume that v has exactly two components, second one not 0. |
---|
| 811 | RETURN: vector, representing simplified rational function v[1]/v[2] |
---|
| 812 | NOTE: Possibly present non-commutative relations of the basering are |
---|
| 813 | ignored. |
---|
| 814 | EXAMPLE: example simplifyRat; shows examples |
---|
| 815 | " |
---|
| 816 | { |
---|
| 817 | checkRatInput(J); |
---|
| 818 | return(arithmeticRat(J,vector(0),"s")); |
---|
| 819 | } |
---|
| 820 | example |
---|
| 821 | { |
---|
| 822 | "EXAMPLE:"; echo = 2; |
---|
| 823 | ring r = 0,(x,y),dp; |
---|
| 824 | vector v = [x2-1,x+1]; |
---|
| 825 | simplifyRat(v); |
---|
| 826 | simplifyRat(v) - [x-1,1]; |
---|
| 827 | } |
---|
| 828 | |
---|
| 829 | |
---|
| 830 | proc addRat (vector I, vector J) |
---|
| 831 | " |
---|
| 832 | USAGE: addRat(v,w); v,w vectors |
---|
| 833 | ASSUME: Assume that v,w have exactly two components, second ones not 0. |
---|
| 834 | RETURN: vector, representing rational function (v[1]/v[2])+(w[1]/w[2]) |
---|
| 835 | NOTE: Possibly present non-commutative relations of the basering are |
---|
| 836 | ignored. |
---|
| 837 | EXAMPLE: example addRat; shows examples |
---|
| 838 | " |
---|
| 839 | { |
---|
| 840 | checkRatInput(I); |
---|
| 841 | checkRatInput(J); |
---|
| 842 | return(arithmeticRat(I,J,"+")); |
---|
| 843 | } |
---|
| 844 | example |
---|
| 845 | { |
---|
| 846 | "EXAMPLE:"; echo = 2; |
---|
| 847 | ring r = 0,(x,y),dp; |
---|
| 848 | vector v = [x,y]; |
---|
| 849 | vector w = [y,x]; |
---|
| 850 | addRat(v,w); |
---|
| 851 | addRat(v,w) - [x2+y2,xy]; |
---|
| 852 | } |
---|
| 853 | |
---|
| 854 | |
---|
| 855 | proc multRat (vector I, vector J) |
---|
| 856 | " |
---|
| 857 | USAGE: multRat(v,w); v,w vectors |
---|
| 858 | ASSUME: Assume that v,w have exactly two components, second ones not 0. |
---|
| 859 | RETURN: vector, representing rational function (v[1]/v[2])*(w[1]/w[2]) |
---|
| 860 | NOTE: Possibly present non-commutative relations of the basering are |
---|
| 861 | ignored. |
---|
| 862 | EXAMPLE: example multRat; shows examples |
---|
| 863 | " |
---|
| 864 | { |
---|
| 865 | checkRatInput(I); |
---|
| 866 | checkRatInput(J); |
---|
| 867 | return(arithmeticRat(I,J,"*")); |
---|
| 868 | } |
---|
| 869 | example |
---|
| 870 | { |
---|
| 871 | "EXAMPLE:"; echo = 2; |
---|
| 872 | ring r = 0,(x,y),dp; |
---|
| 873 | vector v = [x,y]; |
---|
| 874 | vector w = [y,x]; |
---|
| 875 | multRat(v,w); |
---|
| 876 | multRat(v,w) - [1,1]; |
---|
| 877 | } |
---|
| 878 | |
---|
| 879 | |
---|
| 880 | proc diffRat (vector I, int j) |
---|
| 881 | " |
---|
| 882 | USAGE: diffRat(v,j); v vector, j int |
---|
| 883 | ASSUME: Assume that v has exactly two components, second one not 0. |
---|
| 884 | RETURN: vector, representing rational function derivative of rational |
---|
| 885 | function (v[1]/v[2]) w.r.t. var(j) |
---|
| 886 | NOTE: Possibly present non-commutative relations of the basering are |
---|
| 887 | ignored. |
---|
| 888 | EXAMPLE: example diffRat; shows examples |
---|
| 889 | " |
---|
| 890 | { |
---|
| 891 | checkRatInput(I); |
---|
| 892 | if ( (j<1) || (j>nvars(basering)) ) |
---|
| 893 | { |
---|
| 894 | ERROR("Second argument must be in the range 1.."+string(nvars(basering))); |
---|
| 895 | } |
---|
| 896 | return(arithmeticRat(I,vector(0),"d",j)); |
---|
| 897 | } |
---|
| 898 | example |
---|
| 899 | { |
---|
| 900 | "EXAMPLE:"; echo = 2; |
---|
| 901 | ring r = 0,(x,y),dp; |
---|
| 902 | vector v = [x,y]; |
---|
| 903 | diffRat(v,1); |
---|
| 904 | diffRat(v,1) - [1,y]; |
---|
| 905 | diffRat(v,2); |
---|
| 906 | diffRat(v,2) - [-x,y2]; |
---|
| 907 | } |
---|
| 908 | |
---|
| 909 | |
---|
| 910 | proc commRing () |
---|
| 911 | " |
---|
| 912 | USAGE: commRing(); |
---|
| 913 | RETURN: ring, basering without non-commutative relations |
---|
| 914 | EXAMPLE: example commRing; shows examples |
---|
| 915 | " |
---|
| 916 | { |
---|
| 917 | list RL = ringlist(basering); |
---|
| 918 | if (size(RL)<=4) |
---|
| 919 | { |
---|
| 920 | return(basering); |
---|
| 921 | } |
---|
| 922 | RL = RL[1..4]; |
---|
| 923 | def r = ring(RL); |
---|
| 924 | return(r); |
---|
| 925 | } |
---|
| 926 | example |
---|
| 927 | { |
---|
| 928 | "EXAMPLE:"; echo = 2; |
---|
| 929 | def W = makeWeyl(3); |
---|
| 930 | setring W; W; |
---|
| 931 | def W2 = commRing(); |
---|
| 932 | setring W2; W2; |
---|
| 933 | ring r = 0,(x,y),dp; |
---|
| 934 | def r2 = commRing(); // same as r |
---|
| 935 | setring r2; r2; |
---|
| 936 | } |
---|
| 937 | |
---|
| 938 | |
---|
| 939 | // TODO remove this proc once chern.lib is ready |
---|
| 940 | static proc orderedPartition(int n, list #) |
---|
| 941 | " |
---|
| 942 | USUAGE: orderedPartition(n,a); n,a positive ints |
---|
| 943 | orderedPartition(n,w); n positive int, w positive intvec |
---|
| 944 | RETURN: list of intvecs |
---|
| 945 | PURPOSE: Computes all partitions of n of length a, if the second |
---|
| 946 | argument is an int, or computes all weighted partitions |
---|
| 947 | w.r.t. w of n of length size(w) if the second argument |
---|
| 948 | is an intvec. |
---|
| 949 | In both cases, zero parts are included. |
---|
| 950 | EXAMPLE: example orderedPartition; shows an example |
---|
| 951 | " |
---|
| 952 | { |
---|
| 953 | int a,wrongInpt,intInpt; |
---|
| 954 | intvec w = 1; |
---|
| 955 | if (size(#)>0) |
---|
| 956 | { |
---|
| 957 | if (typeof(#[1]) == "int") |
---|
| 958 | { |
---|
| 959 | a = #[1]; |
---|
| 960 | intInpt = 1; |
---|
| 961 | } |
---|
| 962 | else |
---|
| 963 | { |
---|
| 964 | if (typeof(#[1]) == "intvec") |
---|
| 965 | { |
---|
| 966 | w = #[1]; |
---|
| 967 | a = size(w); |
---|
| 968 | } |
---|
| 969 | else |
---|
| 970 | { |
---|
| 971 | wrongInpt = 1; |
---|
| 972 | } |
---|
| 973 | } |
---|
| 974 | } |
---|
| 975 | else |
---|
| 976 | { |
---|
| 977 | wrongInpt = 1; |
---|
| 978 | } |
---|
| 979 | if (wrongInpt) |
---|
| 980 | { |
---|
| 981 | ERROR("Expected second argument of type int or intvec."); |
---|
| 982 | } |
---|
| 983 | kill wrongInpt; |
---|
| 984 | if (n==0 && a>0) |
---|
| 985 | { |
---|
| 986 | return(list(0:a)); |
---|
| 987 | } |
---|
| 988 | if (n<=0 || a<=0 || allPositive(w)==0) |
---|
| 989 | { |
---|
| 990 | ERROR("Positive arguments expected."); |
---|
| 991 | } |
---|
| 992 | int baseringdef; |
---|
| 993 | if (defined(basering)) // if a basering is defined, it should be saved for later use |
---|
| 994 | { |
---|
| 995 | def save = basering; |
---|
| 996 | baseringdef = 1; |
---|
| 997 | } |
---|
| 998 | ring r = 0,(x(1..a)),dp; // all variables for partition of length a |
---|
| 999 | ideal M; |
---|
| 1000 | if (intInpt) |
---|
| 1001 | { |
---|
| 1002 | M = maxideal(n); // all monomials of total degree n |
---|
| 1003 | } |
---|
| 1004 | else |
---|
| 1005 | { |
---|
| 1006 | M = weightKB(ideal(0),n,w); // all monomials of total weighted degree n |
---|
| 1007 | } |
---|
| 1008 | list L; |
---|
| 1009 | int i; |
---|
| 1010 | for (i = 1; i <= ncols(M); i++) {L = insert(L,leadexp(M[i]));} |
---|
| 1011 | // the leadexp corresponds to a partition |
---|
| 1012 | if (baseringdef) // sets the old ring as basering again |
---|
| 1013 | { |
---|
| 1014 | setring save; |
---|
| 1015 | } |
---|
| 1016 | return(L); //returns the list of partitions |
---|
| 1017 | } |
---|
| 1018 | example |
---|
| 1019 | { |
---|
| 1020 | "EXAMPLE"; echo = 2; |
---|
| 1021 | orderedPartition(4,2); |
---|
| 1022 | orderedPartition(5,3); |
---|
| 1023 | orderedPartition(2,4); |
---|
| 1024 | orderedPartition(8,intvec(2,3)); |
---|
| 1025 | orderedPartition(7,intvec(2,2)); // no such partition |
---|
| 1026 | } |
---|
| 1027 | |
---|
| 1028 | |
---|
| 1029 | // applications of characteristic variety ///////////////////////////////////// |
---|
| 1030 | |
---|
| 1031 | proc holonomicRank (ideal I, list #) |
---|
| 1032 | " |
---|
| 1033 | USAGE: holonomicRank(I[,e]); I ideal, e optional int |
---|
| 1034 | ASSUME: The basering is the n-th Weyl algebra over a field of |
---|
| 1035 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 1036 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 1037 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) |
---|
| 1038 | is the differential operator belonging to x(i). |
---|
| 1039 | RETURN: int, the holonomic rank of I |
---|
| 1040 | REMARKS: The holonomic rank of I is defined to be the K(x(1..n))-dimension of |
---|
| 1041 | the module W/WI, where W is the rational Weyl algebra |
---|
| 1042 | K(x(1..n))<D(1..n)>. |
---|
| 1043 | If this dimension is infinite, -1 is returned. |
---|
| 1044 | NOTE: If e<>0, @code{std} is used for Groebner basis computations, |
---|
| 1045 | otherwise (and by default) @code{slimgb} is used. |
---|
| 1046 | @* If printlevel=1, progress debug messages will be printed, |
---|
| 1047 | if printlevel>=2, all the debug messages will be printed. |
---|
| 1048 | EXAMPLE: example holonomicRank; shows examples |
---|
| 1049 | " |
---|
| 1050 | { |
---|
| 1051 | // assumption check is done by charVariety |
---|
| 1052 | int ppl = printlevel - voice + 2; |
---|
| 1053 | int eng; |
---|
| 1054 | if (size(#)>0) |
---|
| 1055 | { |
---|
| 1056 | if(typeof(#[1])=="int") |
---|
| 1057 | { |
---|
| 1058 | eng = #[1]; |
---|
| 1059 | } |
---|
| 1060 | } |
---|
| 1061 | def save = basering; |
---|
| 1062 | dbprint(ppl ,"// Computing characteristic variety..."); |
---|
| 1063 | def grD = charVariety(I); |
---|
| 1064 | setring grD; // commutative ring |
---|
| 1065 | dbprint(ppl ,"// ...done."); |
---|
| 1066 | dbprint(ppl-1,"// " + string(charVar)); |
---|
| 1067 | int n = nvars(save) div 2; |
---|
| 1068 | intvec v = 1..n; |
---|
| 1069 | def R = vars2pars(v); |
---|
| 1070 | setring R; |
---|
| 1071 | ideal J = imap(grD,charVar); |
---|
| 1072 | dbprint(ppl ,"// Starting GB computation..."); |
---|
| 1073 | J = engine(J,0); // use slimgb |
---|
| 1074 | dbprint(ppl ,"// ...done."); |
---|
| 1075 | dbprint(ppl-1,"// " + string(J)); |
---|
| 1076 | int d = vdim(J); |
---|
| 1077 | setring save; |
---|
| 1078 | return(d); |
---|
| 1079 | } |
---|
| 1080 | example |
---|
| 1081 | { |
---|
| 1082 | "EXAMPLE:"; echo = 2; |
---|
| 1083 | // (OTW), Example 8 |
---|
| 1084 | ring r3 = 0,(x,y,z,Dx,Dy,Dz),dp; |
---|
| 1085 | def D3 = Weyl(); |
---|
| 1086 | setring D3; |
---|
| 1087 | poly f = x^3-y^2*z^2; |
---|
| 1088 | ideal I = f^2*Dx+3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z; |
---|
| 1089 | // I annihilates exp(1/f) |
---|
| 1090 | holonomicRank(I); |
---|
| 1091 | } |
---|
| 1092 | |
---|
| 1093 | |
---|
| 1094 | proc DsingularLocus (ideal I) |
---|
| 1095 | " |
---|
| 1096 | USAGE: DsingularLocus(I); I ideal |
---|
| 1097 | ASSUME: The basering is the n-th Weyl algebra over a field of |
---|
| 1098 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 1099 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 1100 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) |
---|
| 1101 | is the differential operator belonging to x(i). |
---|
| 1102 | RETURN: ideal, describing the singular locus of the D-module D/I |
---|
| 1103 | NOTE: If printlevel>=1, progress debug messages will be printed, |
---|
| 1104 | if printlevel>=2, all the debug messages will be printed |
---|
| 1105 | EXAMPLE: example DsingularLocus; shows examples |
---|
| 1106 | " |
---|
| 1107 | { |
---|
| 1108 | // assumption check is done by charVariety |
---|
| 1109 | int ppl = printlevel - voice + 2; |
---|
| 1110 | def save = basering; |
---|
| 1111 | dbprint(ppl ,"// Computing characteristic variety..."); |
---|
| 1112 | def grD = charVariety(I); |
---|
| 1113 | setring grD; |
---|
| 1114 | dbprint(ppl ,"// ...done"); |
---|
| 1115 | dbprint(ppl-1,"// " + string(charVar)); |
---|
| 1116 | poly pDD = 1; |
---|
| 1117 | ideal IDD; |
---|
| 1118 | int i; |
---|
| 1119 | int n = nvars(basering) div 2; |
---|
| 1120 | for (i=1; i<=n; i++) |
---|
| 1121 | { |
---|
| 1122 | pDD = pDD*var(i+n); |
---|
| 1123 | IDD[i] = var(i+n); |
---|
| 1124 | } |
---|
| 1125 | dbprint(ppl ,"// Computing saturation..."); |
---|
| 1126 | ideal S = sat(charVar,IDD)[1]; |
---|
| 1127 | dbprint(ppl ,"// ...done"); |
---|
| 1128 | dbprint(ppl-1,"// " + string(S)); |
---|
| 1129 | dbprint(ppl ,"// Computing elimination..."); |
---|
| 1130 | S = eliminate(S,pDD); |
---|
| 1131 | dbprint(ppl ,"// ...done"); |
---|
| 1132 | dbprint(ppl-1,"// " + string(S)); |
---|
| 1133 | dbprint(ppl ,"// Computing radical..."); |
---|
| 1134 | S = radical(S); |
---|
| 1135 | dbprint(ppl ,"// ...done"); |
---|
| 1136 | dbprint(ppl-1,"// " + string(S)); |
---|
| 1137 | setring save; |
---|
| 1138 | ideal S = imap(grD,S); |
---|
| 1139 | return(S); |
---|
| 1140 | } |
---|
| 1141 | example |
---|
| 1142 | { |
---|
| 1143 | "EXAMPLE:"; echo = 2; |
---|
| 1144 | // (OTW), Example 8 |
---|
| 1145 | ring @D3 = 0,(x,y,z,Dx,Dy,Dz),dp; |
---|
| 1146 | def D3 = Weyl(); |
---|
| 1147 | setring D3; |
---|
| 1148 | poly f = x^3-y^2*z^2; |
---|
| 1149 | ideal I = f^2*Dx + 3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z; |
---|
| 1150 | // I annihilates exp(1/f) |
---|
| 1151 | DsingularLocus(I); |
---|
| 1152 | } |
---|
| 1153 | |
---|
| 1154 | |
---|
| 1155 | // localization /////////////////////////////////////////////////////////////// |
---|
| 1156 | |
---|
| 1157 | static proc finKx(poly f) |
---|
| 1158 | { |
---|
| 1159 | int n = nvars(basering) div 2; |
---|
| 1160 | intvec iv = 1..n; |
---|
| 1161 | if (polyVars(f,iv) == 0) |
---|
| 1162 | { |
---|
| 1163 | ERROR("Given poly may not contain differential operators."); |
---|
| 1164 | } |
---|
| 1165 | return(); |
---|
| 1166 | } |
---|
| 1167 | |
---|
| 1168 | |
---|
| 1169 | proc rightNFWeyl (def id, int k) |
---|
| 1170 | " |
---|
| 1171 | USAGE: rightNFWeyl(id,k); id ideal or poly, k int |
---|
| 1172 | ASSUME: The basering is the n-th Weyl algebra over a field of |
---|
| 1173 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 1174 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 1175 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) |
---|
| 1176 | is the differential operator belonging to x(i). |
---|
| 1177 | RETURN: same type as id, the right normal form of id with respect to the |
---|
| 1178 | principal right ideal generated by the k-th variable |
---|
| 1179 | NOTE: No Groebner basis computation is used. |
---|
| 1180 | EXAMPLE: example rightNFWeyl; shows examples. |
---|
| 1181 | " |
---|
| 1182 | { |
---|
| 1183 | dmodGeneralAssumptionCheck(); |
---|
| 1184 | string inpt = typeof(id); |
---|
| 1185 | if (inpt=="ideal" || inpt=="poly") |
---|
| 1186 | { |
---|
| 1187 | ideal I = id; |
---|
| 1188 | } |
---|
| 1189 | else |
---|
| 1190 | { |
---|
| 1191 | ERROR("Expected first input to be of type ideal or poly."); |
---|
| 1192 | } |
---|
| 1193 | def save = basering; |
---|
| 1194 | int n = nvars(save) div 2; |
---|
| 1195 | if (0>k || k>2*n) |
---|
| 1196 | { |
---|
| 1197 | ERROR("Expected second input to be in the range 1.."+string(2*n)+"."); |
---|
| 1198 | } |
---|
| 1199 | int i,j; |
---|
| 1200 | if (k>n) // var(k) = Dx(k-n) |
---|
| 1201 | { |
---|
| 1202 | // switch var(k),var(k-n) |
---|
| 1203 | list RL = ringlist(save); |
---|
| 1204 | matrix rel = RL[6]; |
---|
| 1205 | rel[k-n,k] = -1; |
---|
| 1206 | RL = RL[1..4]; |
---|
| 1207 | list L = RL[2]; |
---|
| 1208 | string str = L[k-n]; |
---|
| 1209 | L[k-n] = L[k]; |
---|
| 1210 | L[k] = str; |
---|
| 1211 | RL[2] = L; |
---|
| 1212 | def @W = ring(RL); |
---|
| 1213 | kill L,RL,str; |
---|
| 1214 | ideal @mm = maxideal(1); |
---|
| 1215 | setring @W; |
---|
| 1216 | matrix rel = imap(save,rel); |
---|
| 1217 | def W = nc_algebra(1,rel); |
---|
| 1218 | setring W; |
---|
| 1219 | ideal @mm = imap(save,@mm); |
---|
| 1220 | map mm = save,@mm; |
---|
| 1221 | ideal I = mm(I); |
---|
| 1222 | i = k-n; |
---|
| 1223 | } |
---|
| 1224 | else // var(k) = x(k) |
---|
| 1225 | { |
---|
| 1226 | def W = save; |
---|
| 1227 | i = k; |
---|
| 1228 | } |
---|
| 1229 | for (j=1; j<=ncols(I); j++) |
---|
| 1230 | { |
---|
| 1231 | I[j] = subst(I[j],var(i),0); |
---|
| 1232 | } |
---|
| 1233 | setring save; |
---|
| 1234 | I = imap(W,I); |
---|
| 1235 | if (inpt=="poly") |
---|
| 1236 | { |
---|
| 1237 | return(I[1]); |
---|
| 1238 | } |
---|
| 1239 | else |
---|
| 1240 | { |
---|
| 1241 | return(I); |
---|
| 1242 | } |
---|
| 1243 | } |
---|
| 1244 | example |
---|
| 1245 | { |
---|
| 1246 | "EXAMPLE:"; echo = 2; |
---|
| 1247 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
| 1248 | def W = Weyl(); |
---|
| 1249 | setring W; |
---|
| 1250 | ideal I = x^3*Dx^3, y^2*Dy^2, x*Dy, y*Dx; |
---|
| 1251 | rightNFWeyl(I,1); // right NF wrt principal right ideal x*W |
---|
| 1252 | rightNFWeyl(I,3); // right NF wrt principal right ideal Dx*W |
---|
| 1253 | rightNFWeyl(I,2); // right NF wrt principal right ideal y*W |
---|
| 1254 | rightNFWeyl(I,4); // right NF wrt principal right ideal Dy*W |
---|
| 1255 | poly p = x*Dx+1; |
---|
| 1256 | rightNFWeyl(p,1); // right NF wrt principal right ideal x*W |
---|
| 1257 | } |
---|
| 1258 | |
---|
| 1259 | |
---|
| 1260 | // TODO check OTW for assumptions on holonomicity |
---|
| 1261 | proc Dlocalization (ideal J, poly f, list #) |
---|
| 1262 | " |
---|
| 1263 | USAGE: Dlocalization(I,f[,k,e]); I ideal, f poly, k,e optional ints |
---|
| 1264 | ASSUME: The basering is the n-th Weyl algebra over a field of |
---|
| 1265 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 1266 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 1267 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) |
---|
| 1268 | is the differential operator belonging to x(i). |
---|
| 1269 | @* Further, assume that f does not contain any D(i) and that I is |
---|
| 1270 | holonomic on K^n\V(f). |
---|
| 1271 | RETURN: ideal or list, computes an ideal J such that D/J is isomorphic |
---|
| 1272 | to D/I localized at f as D-modules. |
---|
| 1273 | If k<>0, a list consisting of J and an integer m is returned, |
---|
| 1274 | such that f^m represents the natural map from D/I to D/J. |
---|
| 1275 | Otherwise (and by default), only the ideal J is returned. |
---|
| 1276 | REMARKS: It is known that a localization at f of a holonomic D-module is |
---|
| 1277 | again a holonomic D-module. |
---|
| 1278 | @* Reference: (OTW) |
---|
| 1279 | NOTE: If e<>0, @code{std} is used for Groebner basis computations, |
---|
| 1280 | otherwise (and by default) @code{slimgb} is used. |
---|
| 1281 | @* If printlevel=1, progress debug messages will be printed, |
---|
| 1282 | if printlevel>=2, all the debug messages will be printed. |
---|
| 1283 | SEE ALSO: DLoc, SDLoc, DLoc0 |
---|
| 1284 | EXAMPLE: example Dlocalization; shows examples |
---|
| 1285 | " |
---|
| 1286 | { |
---|
| 1287 | dmodGeneralAssumptionCheck(); |
---|
| 1288 | finKx(f); |
---|
| 1289 | int ppl = printlevel - voice + 2; |
---|
| 1290 | int outList,eng; |
---|
| 1291 | if (size(#)>0) |
---|
| 1292 | { |
---|
| 1293 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
| 1294 | { |
---|
| 1295 | outList = int(#[1]); |
---|
| 1296 | } |
---|
| 1297 | if (size(#)>1) |
---|
| 1298 | { |
---|
| 1299 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
| 1300 | { |
---|
| 1301 | eng = int(#[2]); |
---|
| 1302 | } |
---|
| 1303 | } |
---|
| 1304 | } |
---|
| 1305 | int i,j; |
---|
| 1306 | def save = basering; |
---|
| 1307 | int n = nvars(save) div 2; |
---|
| 1308 | def Dv = extendWeyl(safeVarName("v")); |
---|
| 1309 | setring Dv; |
---|
| 1310 | poly f = imap(save,f); |
---|
| 1311 | ideal phiI; |
---|
| 1312 | for (i=n; i>0; i--) |
---|
| 1313 | { |
---|
| 1314 | phiI[i+n] = var(i+n+2)-var(1)^2*bracket(var(i+n+2),f)*var(n+2); |
---|
| 1315 | phiI[i] = var(i+1); |
---|
| 1316 | } |
---|
| 1317 | map phi = save,phiI; |
---|
| 1318 | ideal J = phi(J); |
---|
| 1319 | J = J, 1-f*var(1); |
---|
| 1320 | // TODO original J has to be holonomic only on K^n\V(f), not on all of K^n |
---|
| 1321 | // does is suffice to show that new J is holonomic on Dv?? |
---|
| 1322 | if (isHolonomic(J) == 0) |
---|
| 1323 | { |
---|
| 1324 | ERROR("Module is not holonomic."); |
---|
| 1325 | } |
---|
| 1326 | intvec w = 1; w[n+1]=0; |
---|
| 1327 | ideal G = GBWeight(J,w,-w,eng); |
---|
| 1328 | dbprint(ppl ,"// found GB wrt weight " +string(-w)); |
---|
| 1329 | dbprint(ppl-1,"// " + string(G)); |
---|
| 1330 | intvec ww = w,-w; |
---|
| 1331 | ideal inG = inForm(G,ww); |
---|
| 1332 | inG = engine(inG,eng); |
---|
| 1333 | poly s = var(1)*var(n+2); // s=v*Dv |
---|
| 1334 | vector intersecvec = pIntersect(s,inG); |
---|
| 1335 | s = vec2poly(intersecvec); |
---|
| 1336 | s = subst(s,var(1),-var(1)-1); |
---|
| 1337 | list L = bFactor(s); |
---|
| 1338 | dbprint(ppl ,"// found b-function"); |
---|
| 1339 | dbprint(ppl-1,"// roots: "+string(L[1])); |
---|
| 1340 | dbprint(ppl-1,"// multiplicities: "+string(L[2])); |
---|
| 1341 | kill inG,intersecvec,s; |
---|
| 1342 | // TODO: use maxIntRoot |
---|
| 1343 | L = intRoots(L); // integral roots of b-function |
---|
| 1344 | if (L[2]==0:size(L[2])) // no integral roots |
---|
| 1345 | { |
---|
| 1346 | setring save; |
---|
| 1347 | return(ideal(1)); |
---|
| 1348 | } |
---|
| 1349 | intvec iv; |
---|
| 1350 | for (i=1; i<=ncols(L[1]); i++) |
---|
| 1351 | { |
---|
| 1352 | iv[i] = int(L[1][i]); |
---|
| 1353 | } |
---|
| 1354 | int l0 = Max(iv); |
---|
| 1355 | dbprint(ppl,"// maximal integral root is " +string(l0)); |
---|
| 1356 | kill L,iv; |
---|
| 1357 | intvec degG; |
---|
| 1358 | ideal Gk; |
---|
| 1359 | for (j=1; j<=ncols(G); j++) |
---|
| 1360 | { |
---|
| 1361 | degG[j] = deg(G[j],ww); |
---|
| 1362 | for (i=0; i<=l0-degG[j]; i++) |
---|
| 1363 | { |
---|
| 1364 | Gk[ncols(Gk)+1] = var(1)^i*G[j]; |
---|
| 1365 | } |
---|
| 1366 | } |
---|
| 1367 | Gk = rightNFWeyl(Gk,n+2); |
---|
| 1368 | dbprint(ppl,"// found right normalforms"); |
---|
| 1369 | module M = coeffs(Gk,var(1)); |
---|
| 1370 | setring save; |
---|
| 1371 | def mer = makeModElimRing(save); |
---|
| 1372 | setring mer; |
---|
| 1373 | module M = imap(Dv,M); |
---|
| 1374 | kill Dv; |
---|
| 1375 | M = engine(M,eng); |
---|
| 1376 | dbprint(ppl ,"// found GB of free module of rank " + string(l0+1)); |
---|
| 1377 | dbprint(ppl-1,"// " + string(M)); |
---|
| 1378 | M = prune(M); |
---|
| 1379 | setring save; |
---|
| 1380 | matrix M = imap(mer,M); |
---|
| 1381 | kill mer; |
---|
| 1382 | int ro = nrows(M); |
---|
| 1383 | int co = ncols(M); |
---|
| 1384 | ideal I; |
---|
| 1385 | if (ro == 1) // nothing to do |
---|
| 1386 | { |
---|
| 1387 | I = M; |
---|
| 1388 | } |
---|
| 1389 | else |
---|
| 1390 | { |
---|
| 1391 | matrix zm[ro-1][1]; // zero matrix |
---|
| 1392 | matrix v[ro-1][1]; |
---|
| 1393 | for (i=1; i<=co; i++) |
---|
| 1394 | { |
---|
| 1395 | v = M[1..ro-1,i]; |
---|
| 1396 | if (v == zm) |
---|
| 1397 | { |
---|
| 1398 | I[size(I)+1] = M[ro,i]; |
---|
| 1399 | } |
---|
| 1400 | } |
---|
| 1401 | } |
---|
| 1402 | if (outList<>0) |
---|
| 1403 | { |
---|
| 1404 | return(list(I,l0+2)); |
---|
| 1405 | } |
---|
| 1406 | else |
---|
| 1407 | { |
---|
| 1408 | return(I); |
---|
| 1409 | } |
---|
| 1410 | } |
---|
| 1411 | example |
---|
| 1412 | { |
---|
| 1413 | "EXAMPLE:"; echo = 2; |
---|
| 1414 | // (OTW), Example 8 |
---|
| 1415 | ring r = 0,(x,y,z,Dx,Dy,Dz),dp; |
---|
| 1416 | def W = Weyl(); |
---|
| 1417 | setring W; |
---|
| 1418 | poly f = x^3-y^2*z^2; |
---|
| 1419 | ideal I = f^2*Dx+3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z; |
---|
| 1420 | // I annihilates exp(1/f) |
---|
| 1421 | ideal J = Dlocalization(I,f); |
---|
| 1422 | J; |
---|
| 1423 | Dlocalization(I,f,1); // The natural map D/I -> D/J is given by 1/f^2 |
---|
| 1424 | } |
---|
| 1425 | |
---|
| 1426 | |
---|
| 1427 | |
---|
| 1428 | // Weyl closure /////////////////////////////////////////////////////////////// |
---|
| 1429 | |
---|
| 1430 | static proc orderFiltrationD1 (poly f) |
---|
| 1431 | { |
---|
| 1432 | // returns list of ideal and intvec |
---|
| 1433 | // ideal contains x-parts, intvec corresponding degree in Dx |
---|
| 1434 | poly g,h; |
---|
| 1435 | g = f; |
---|
| 1436 | ideal I; |
---|
| 1437 | intvec v,w,u; |
---|
| 1438 | w = 0,1; |
---|
| 1439 | int i,j; |
---|
| 1440 | i = 1; |
---|
| 1441 | while (g<>0) |
---|
| 1442 | { |
---|
| 1443 | h = inForm(g,w); |
---|
| 1444 | I[i] = 0; |
---|
| 1445 | for (j=1; j<=size(h); j++) |
---|
| 1446 | { |
---|
| 1447 | v = leadexp(h[j]); |
---|
| 1448 | u[i] = v[2]; |
---|
| 1449 | v[2] = 0; |
---|
| 1450 | I[i] = I[i] + leadcoef(h[j])*monomial(v); |
---|
| 1451 | } |
---|
| 1452 | g = g-h; |
---|
| 1453 | i++; |
---|
| 1454 | } |
---|
| 1455 | return(list(I,u)); |
---|
| 1456 | } |
---|
| 1457 | |
---|
| 1458 | |
---|
| 1459 | static proc kerLinMapD1 (ideal W, poly L, poly p) |
---|
| 1460 | { |
---|
| 1461 | // computes kernel of right multiplication with L viewed |
---|
| 1462 | // as homomorphism of K-vector spaces span(W) -> D1/p*D1 |
---|
| 1463 | // assume p in K[x], basering is K<x,Dx> |
---|
| 1464 | ideal G,K; |
---|
| 1465 | G = std(p); |
---|
| 1466 | list l; |
---|
| 1467 | int i,j; |
---|
| 1468 | // first, compute the image of span(W) |
---|
| 1469 | if (simplify(W,2)[1] == 0) |
---|
| 1470 | { |
---|
| 1471 | return(K); // = 0 |
---|
| 1472 | } |
---|
| 1473 | for (i=1; i<=size(W); i++) |
---|
| 1474 | { |
---|
| 1475 | l = orderFiltrationD1(W[i]*L); |
---|
| 1476 | K[i] = 0; |
---|
| 1477 | for (j=1; j<=size(l[1]); j++) |
---|
| 1478 | { |
---|
| 1479 | K[i] = K[i] + NF(l[1][j],G)*var(2)^(l[2][j]); |
---|
| 1480 | } |
---|
| 1481 | } |
---|
| 1482 | // now, we get the kernel by linear algebra |
---|
| 1483 | l = linReduceIdeal(K,1); |
---|
| 1484 | i = ncols(l[1]) - size(l[1]); |
---|
| 1485 | if (i<>0) |
---|
| 1486 | { |
---|
| 1487 | K = module(W)*l[2]; |
---|
| 1488 | K = K[1..i]; |
---|
| 1489 | } |
---|
| 1490 | else |
---|
| 1491 | { |
---|
| 1492 | K = 0; |
---|
| 1493 | } |
---|
| 1494 | return(K); |
---|
| 1495 | } |
---|
| 1496 | |
---|
| 1497 | |
---|
| 1498 | static proc leftDivisionKxD1 (poly p, poly L) |
---|
| 1499 | { |
---|
| 1500 | // basering is D1 = K<x,Dx> |
---|
| 1501 | // p in K[x] |
---|
| 1502 | // compute p^(-1)*L if p is a left divisor of L |
---|
| 1503 | // if (rightNF(L,ideal(p))<>0) |
---|
| 1504 | // { |
---|
| 1505 | // ERROR("First poly is not a right factor of second poly"); |
---|
| 1506 | // } |
---|
| 1507 | def save = basering; |
---|
| 1508 | list l = orderFiltrationD1(L); |
---|
| 1509 | ideal l1 = l[1]; |
---|
| 1510 | ring r = 0,x,dp; |
---|
| 1511 | ideal l1 = fetch(save,l1); |
---|
| 1512 | poly p = fetch(save,p); |
---|
| 1513 | int i; |
---|
| 1514 | for (i=1; i<=ncols(l1); i++) |
---|
| 1515 | { |
---|
| 1516 | l1[i] = division(l1[i],p)[1][1,1]; |
---|
| 1517 | } |
---|
| 1518 | setring save; |
---|
| 1519 | ideal I = fetch(r,l1); |
---|
| 1520 | poly f; |
---|
| 1521 | for (i=1; i<=ncols(I); i++) |
---|
| 1522 | { |
---|
| 1523 | f = f + I[i]*var(2)^(l[2][i]); |
---|
| 1524 | } |
---|
| 1525 | return(f); |
---|
| 1526 | } |
---|
| 1527 | |
---|
| 1528 | |
---|
| 1529 | proc WeylClosure1 (poly L) |
---|
| 1530 | " |
---|
| 1531 | USAGE: WeylClosure1(L); L a poly |
---|
| 1532 | ASSUME: The basering is the first Weyl algebra D=K<x,d|dx=xd+1> over a field |
---|
| 1533 | K of characteristic 0. |
---|
| 1534 | RETURN: ideal, the Weyl closure of the principal left ideal generated by L |
---|
| 1535 | REMARKS: The Weyl closure of a left ideal I in the Weyl algebra D is defined |
---|
| 1536 | to be the intersection of I regarded as left ideal in the rational |
---|
| 1537 | Weyl algebra K(x)<d> with the polynomial Weyl algebra D. |
---|
| 1538 | @* Reference: (Tsa), Algorithm 1.2.4 |
---|
| 1539 | NOTE: If printlevel=1, progress debug messages will be printed, |
---|
| 1540 | if printlevel>=2, all the debug messages will be printed. |
---|
| 1541 | SEE ALSO: WeylClosure |
---|
| 1542 | EXAMPLE: example WeylClosure1; shows examples |
---|
| 1543 | " |
---|
| 1544 | { |
---|
| 1545 | dmodGeneralAssumptionCheck(); // assumption check |
---|
| 1546 | int ppl = printlevel - voice + 2; |
---|
| 1547 | def save = basering; |
---|
| 1548 | intvec w = 0,1; // for order filtration |
---|
| 1549 | poly p = inForm(L,w); |
---|
| 1550 | ring @R = 0,var(1),dp; |
---|
| 1551 | ideal mm = var(1),1; |
---|
| 1552 | map m = save,mm; |
---|
| 1553 | ideal @p = m(p); |
---|
| 1554 | poly p = @p[1]; |
---|
| 1555 | poly g = gcd(p,diff(p,var(1))); |
---|
| 1556 | if (g == 1) |
---|
| 1557 | { |
---|
| 1558 | g = p; |
---|
| 1559 | } |
---|
| 1560 | ideal facp = factorize(g,1); // g is squarefree, constants aren't interesting |
---|
| 1561 | dbprint(ppl-1, |
---|
| 1562 | "// squarefree part of highest coefficient w.r.t. order filtration:"); |
---|
| 1563 | dbprint(ppl-1, "// " + string(facp)); |
---|
| 1564 | setring save; |
---|
| 1565 | p = imap(@R,p); |
---|
| 1566 | // 1-1 extend basering by parameter and introduce new var t=x*d |
---|
| 1567 | list RL = ringlist(save); |
---|
| 1568 | RL = RL[1..4]; |
---|
| 1569 | list l; |
---|
| 1570 | l[1] = int(0); |
---|
| 1571 | l[2] = list(safeVarName("a")); |
---|
| 1572 | l[3] = list(list("lp",intvec(1))); |
---|
| 1573 | l[4] = ideal(0); |
---|
| 1574 | RL[1] = l; |
---|
| 1575 | l = RL[2] + list(safeVarName("t")); |
---|
| 1576 | RL[2] = l; |
---|
| 1577 | l = list(); |
---|
| 1578 | l[1] = list("dp",intvec(1,1)); |
---|
| 1579 | l[2] = list("dp",intvec(1)); |
---|
| 1580 | l[3] = list("C",intvec(0)); |
---|
| 1581 | RL[3] = l; |
---|
| 1582 | def @Wat = ring(RL); |
---|
| 1583 | kill RL,l; |
---|
| 1584 | setring @Wat; |
---|
| 1585 | matrix relD[3][3]; |
---|
| 1586 | relD[1,2] = 1; |
---|
| 1587 | relD[1,3] = var(1); |
---|
| 1588 | relD[2,3] = -var(2); |
---|
| 1589 | def Wat = nc_algebra(1,relD); |
---|
| 1590 | setring Wat; |
---|
| 1591 | kill @Wat; |
---|
| 1592 | // 1-2 rewrite L using Euler operators |
---|
| 1593 | ideal mm = var(1)+par(1),var(2); |
---|
| 1594 | map m = save,mm; |
---|
| 1595 | poly L = m(L); |
---|
| 1596 | w = -1,1,0; // for Bernstein filtration |
---|
| 1597 | int i = 1; |
---|
| 1598 | ideal Q; |
---|
| 1599 | poly p = L; |
---|
| 1600 | intvec d; |
---|
| 1601 | while (p<>0) |
---|
| 1602 | { |
---|
| 1603 | Q[i] = inForm(p,w); |
---|
| 1604 | p = p - Q[i]; |
---|
| 1605 | d[i] = -deg(Q[i],w); |
---|
| 1606 | i++; |
---|
| 1607 | } |
---|
| 1608 | ideal S = std(var(1)*var(2)-var(3)); |
---|
| 1609 | Q = NF(Q,S); |
---|
| 1610 | dbprint(ppl, "// found Euler representation of operator"); |
---|
| 1611 | dbprint(ppl-1,"// " + string(Q)); |
---|
| 1612 | Q = subst(Q,var(1),1); |
---|
| 1613 | Q = subst(Q,var(2),1); |
---|
| 1614 | // 1-3 prepare for algebraic extensions with minpoly = facp[i] |
---|
| 1615 | list RL = ringlist(Wat); |
---|
| 1616 | RL = RL[1..4]; |
---|
| 1617 | list l; |
---|
| 1618 | l = string(var(3)); |
---|
| 1619 | RL[2] = l; |
---|
| 1620 | l = list(); |
---|
| 1621 | l[1] = list("dp",intvec(1)); |
---|
| 1622 | l[2] = list("C",intvec(0)); |
---|
| 1623 | RL[3] = l; |
---|
| 1624 | mm = par(1); |
---|
| 1625 | m = @R,par(1); |
---|
| 1626 | ideal facp = m(facp); |
---|
| 1627 | kill @R,m,mm,l,S; |
---|
| 1628 | intvec maxroots,testroots; |
---|
| 1629 | int sq = size(Q); |
---|
| 1630 | string strQ = "ideal Q = " + string(Q) + ";"; |
---|
| 1631 | // TODO do it without string workaround when issue with maps from |
---|
| 1632 | // transcendental to algebraic extension fields is fixed |
---|
| 1633 | int j,maxr; |
---|
| 1634 | // 2-1 get max int root of lowest nonzero entry of Q in algebraic extension |
---|
| 1635 | for (i=1; i<=size(facp); i++) |
---|
| 1636 | { |
---|
| 1637 | testroots = 0; |
---|
| 1638 | def Ra = ring(RL); |
---|
| 1639 | setring Ra; |
---|
| 1640 | ideal mm = 1,1,var(1); |
---|
| 1641 | map m = Wat,mm; |
---|
| 1642 | ideal facp = m(facp); |
---|
| 1643 | minpoly = leadcoef(facp[i]); |
---|
| 1644 | execute(strQ); |
---|
| 1645 | if (simplify(Q,2)[1] == poly(0)) |
---|
| 1646 | { |
---|
| 1647 | break; |
---|
| 1648 | } |
---|
| 1649 | j = 1; |
---|
| 1650 | while (j<sq) |
---|
| 1651 | { |
---|
| 1652 | if (Q[j]==0) |
---|
| 1653 | { |
---|
| 1654 | j++; |
---|
| 1655 | } |
---|
| 1656 | else |
---|
| 1657 | { |
---|
| 1658 | break; |
---|
| 1659 | } |
---|
| 1660 | } |
---|
| 1661 | maxroots[i] = d[j]; // d[j] = r_k |
---|
| 1662 | list LR = bFactor(Q[j]); |
---|
| 1663 | LR = intRoots(LR); |
---|
| 1664 | if (LR[2]<>0:size(LR[2])) // there are integral roots |
---|
| 1665 | { |
---|
| 1666 | for (j=1; j<=ncols(LR[1]); j++) |
---|
| 1667 | { |
---|
| 1668 | testroots[j] = int(LR[1][j]); |
---|
| 1669 | } |
---|
| 1670 | maxr = Max(testroots); |
---|
| 1671 | if(maxr<0) |
---|
| 1672 | { |
---|
| 1673 | maxr = 0; |
---|
| 1674 | } |
---|
| 1675 | maxroots[i] = maxroots[i] + maxr; |
---|
| 1676 | } |
---|
| 1677 | kill LR; |
---|
| 1678 | setring Wat; |
---|
| 1679 | kill Ra; |
---|
| 1680 | } |
---|
| 1681 | maxr = Max(maxroots); |
---|
| 1682 | // 3-1 build basis of vectorspace |
---|
| 1683 | setring save; |
---|
| 1684 | ideal KB; |
---|
| 1685 | for (i=0; i<deg(p); i++) // it's really <, not <= |
---|
| 1686 | { |
---|
| 1687 | for (j=0; j<=maxr; j++) // it's really <=, not < |
---|
| 1688 | { |
---|
| 1689 | KB[size(KB)+1] = monomial(intvec(i,j)); |
---|
| 1690 | } |
---|
| 1691 | } |
---|
| 1692 | dbprint(ppl,"// got vector space basis"); |
---|
| 1693 | dbprint(ppl-1, "// " + string(KB)); |
---|
| 1694 | // 3-2 get kernel of *L: span(KB)->D/pD |
---|
| 1695 | KB = kerLinMapD1(KB,L,p); |
---|
| 1696 | dbprint(ppl,"// got kernel"); |
---|
| 1697 | dbprint(ppl-1, "// " + string(KB)); |
---|
| 1698 | // 4-1 get (1/p)*f*L where f in KB |
---|
| 1699 | for (i=1; i<=ncols(KB); i++) |
---|
| 1700 | { |
---|
| 1701 | KB[i] = leftDivisionKxD1(p,KB[i]*L); |
---|
| 1702 | } |
---|
| 1703 | KB = L,KB; |
---|
| 1704 | // 4-2 done |
---|
| 1705 | return(KB); |
---|
| 1706 | } |
---|
| 1707 | example |
---|
| 1708 | { |
---|
| 1709 | "EXAMPLE:"; echo = 2; |
---|
| 1710 | ring r = 0,(x,Dx),dp; |
---|
| 1711 | def W = Weyl(); |
---|
| 1712 | setring W; |
---|
| 1713 | poly L = (x^3+2)*Dx-3*x^2; |
---|
| 1714 | WeylClosure1(L); |
---|
| 1715 | L = (x^4-4*x^3+3*x^2)*Dx^2+(-6*x^3+20*x^2-12*x)*Dx+(12*x^2-32*x+12); |
---|
| 1716 | WeylClosure1(L); |
---|
| 1717 | } |
---|
| 1718 | |
---|
| 1719 | |
---|
| 1720 | proc WeylClosure (ideal I) |
---|
| 1721 | " |
---|
| 1722 | USAGE: WeylClosure(I); I an ideal |
---|
| 1723 | ASSUME: The basering is the n-th Weyl algebra W over a field of |
---|
| 1724 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 1725 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 1726 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the |
---|
| 1727 | differential operator belonging to x(i). |
---|
| 1728 | @* Moreover, assume that the holonomic rank of W/I is finite. |
---|
| 1729 | RETURN: ideal, the Weyl closure of I |
---|
| 1730 | REMARKS: The Weyl closure of a left ideal I in the Weyl algebra W is defined to |
---|
| 1731 | be the intersection of I regarded as left ideal in the rational Weyl |
---|
| 1732 | algebra K(x(1..n))<D(1..n)> with the polynomial Weyl algebra W. |
---|
| 1733 | @* Reference: (Tsa), Algorithm 2.2.4 |
---|
| 1734 | NOTE: If printlevel=1, progress debug messages will be printed, |
---|
| 1735 | if printlevel>=2, all the debug messages will be printed. |
---|
| 1736 | SEE ALSO: WeylClosure1 |
---|
| 1737 | EXAMPLE: example WeylClosure; shows examples |
---|
| 1738 | " |
---|
| 1739 | { |
---|
| 1740 | // assumption check |
---|
| 1741 | dmodGeneralAssumptionCheck(); |
---|
| 1742 | if (holonomicRank(I)==-1) |
---|
| 1743 | { |
---|
| 1744 | ERROR("Input is not of finite holonomic rank."); |
---|
| 1745 | } |
---|
| 1746 | int ppl = printlevel - voice + 2; |
---|
| 1747 | int eng = 0; // engine |
---|
| 1748 | def save = basering; |
---|
| 1749 | dbprint(ppl ,"// Starting to compute singular locus..."); |
---|
| 1750 | ideal sl = DsingularLocus(I); |
---|
| 1751 | sl = simplify(sl,2); |
---|
| 1752 | dbprint(ppl ,"// ...done."); |
---|
| 1753 | dbprint(ppl-1,"// " + string(sl)); |
---|
| 1754 | if (sl[1] == 0) // can never get here |
---|
| 1755 | { |
---|
| 1756 | ERROR("Can't find polynomial in K[x] vanishing on singular locus."); |
---|
| 1757 | } |
---|
| 1758 | poly f = sl[1]; |
---|
| 1759 | dbprint(ppl ,"// Found poly vanishing on singular locus: " + string(f)); |
---|
| 1760 | dbprint(ppl ,"// Starting to compute localization..."); |
---|
| 1761 | list L = Dlocalization(I,f,1); |
---|
| 1762 | ideal G = L[1]; |
---|
| 1763 | dbprint(ppl ,"// ...done."); |
---|
| 1764 | dbprint(ppl-1,"// " + string(G)); |
---|
| 1765 | dbprint(ppl ,"// Starting to compute kernel of localization map..."); |
---|
| 1766 | if (eng == 0) |
---|
| 1767 | { |
---|
| 1768 | G = moduloSlim(f^L[2],G); |
---|
| 1769 | } |
---|
| 1770 | else |
---|
| 1771 | { |
---|
| 1772 | G = modulo(f^L[2],G); |
---|
| 1773 | } |
---|
| 1774 | dbprint(ppl ,"// ...done."); |
---|
| 1775 | return(G); |
---|
| 1776 | } |
---|
| 1777 | example |
---|
| 1778 | { |
---|
| 1779 | "EXAMPLE:"; echo = 2; |
---|
| 1780 | // (OTW), Example 8 |
---|
| 1781 | ring r = 0,(x,y,z,Dx,Dy,Dz),dp; |
---|
| 1782 | def D3 = Weyl(); |
---|
| 1783 | setring D3; |
---|
| 1784 | poly f = x^3-y^2*z^2; |
---|
| 1785 | ideal I = f^2*Dx + 3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z; |
---|
| 1786 | // I annihilates exp(1/f) |
---|
| 1787 | WeylClosure(I); |
---|
| 1788 | } |
---|
| 1789 | |
---|
| 1790 | |
---|
| 1791 | |
---|
| 1792 | // solutions to systems of PDEs /////////////////////////////////////////////// |
---|
| 1793 | |
---|
| 1794 | proc polSol (ideal I, list #) |
---|
| 1795 | " |
---|
| 1796 | USAGE: polSol(I[,w,m]); I ideal, w optional intvec, m optional int |
---|
| 1797 | ASSUME: The basering is the n-th Weyl algebra W over a field of |
---|
| 1798 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 1799 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 1800 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the |
---|
| 1801 | differential operator belonging to x(i). |
---|
| 1802 | @* Moreover, assume that I is holonomic. |
---|
| 1803 | RETURN: ideal, a basis of the polynomial solutions to the given system of |
---|
| 1804 | linear PDEs with polynomial coefficients, encoded via I |
---|
| 1805 | REMARKS: If w is given, w should consist of n strictly negative entries. |
---|
| 1806 | Otherwise and by default, w is set to -1:n. |
---|
| 1807 | In this case, w is used as weight vector for the computation of a |
---|
| 1808 | b-function. |
---|
| 1809 | @* If m is given, m is assumed to be the minimal integer root of the |
---|
| 1810 | b-function of I w.r.t. w. Note that this assumption is not checked. |
---|
| 1811 | @* Reference: (OTT), Algorithm 2.4 |
---|
| 1812 | NOTE: If printlevel=1, progress debug messages will be printed, |
---|
| 1813 | if printlevel>=2, all the debug messages will be printed. |
---|
| 1814 | SEE ALSO: polSolFiniteRank, ratSol |
---|
| 1815 | EXAMPLE: example polSol; shows examples |
---|
| 1816 | " |
---|
| 1817 | { |
---|
| 1818 | dmodGeneralAssumptionCheck(); |
---|
| 1819 | int ppl = printlevel - voice + 2; |
---|
| 1820 | int mr,mrgiven; |
---|
| 1821 | def save = basering; |
---|
| 1822 | int n = nvars(save); |
---|
| 1823 | intvec w = -1:(n div 2); |
---|
| 1824 | if (size(#)>0) |
---|
| 1825 | { |
---|
| 1826 | if (typeof(#[1])=="intvec") |
---|
| 1827 | { |
---|
| 1828 | if (allPositive(-#[1])) |
---|
| 1829 | { |
---|
| 1830 | w = #[1]; |
---|
| 1831 | } |
---|
| 1832 | } |
---|
| 1833 | if (size(#)>1) |
---|
| 1834 | { |
---|
| 1835 | if (typeof(#[2])=="int") |
---|
| 1836 | { |
---|
| 1837 | mr = #[2]; |
---|
| 1838 | mrgiven = 1; |
---|
| 1839 | } |
---|
| 1840 | } |
---|
| 1841 | } |
---|
| 1842 | // Step 1: the b-function |
---|
| 1843 | list L; |
---|
| 1844 | if (!mrgiven) |
---|
| 1845 | { |
---|
| 1846 | if (!isHolonomic(I)) |
---|
| 1847 | { |
---|
| 1848 | ERROR("Ideal is not holonomic. Try polSolFiniteRank."); |
---|
| 1849 | } |
---|
| 1850 | dbprint(ppl,"// Computing b-function..."); |
---|
| 1851 | L = bfctIdeal(I,w); |
---|
| 1852 | dbprint(ppl,"// ...done."); |
---|
| 1853 | dbprint(ppl-1,"// Roots: " + string(L[1])); |
---|
| 1854 | dbprint(ppl-1,"// Multiplicities: " + string(L[2])); |
---|
| 1855 | mr = minIntRoot2(L); |
---|
| 1856 | dbprint(ppl,"// Minimal integer root is " + string(mr) + "."); |
---|
| 1857 | } |
---|
| 1858 | if (mr>0) |
---|
| 1859 | { |
---|
| 1860 | return(ideal(0)); |
---|
| 1861 | } |
---|
| 1862 | // Step 2: get the form of a solution f |
---|
| 1863 | int i; |
---|
| 1864 | L = list(); |
---|
| 1865 | for (i=0; i<=-mr; i++) |
---|
| 1866 | { |
---|
| 1867 | L = L + orderedPartition(i,-w); |
---|
| 1868 | } |
---|
| 1869 | ideal mons; |
---|
| 1870 | for (i=1; i<=size(L); i++) |
---|
| 1871 | { |
---|
| 1872 | mons[i] = monomial(L[i]); |
---|
| 1873 | } |
---|
| 1874 | kill L; |
---|
| 1875 | mons = simplify(mons,2+4); // L might contain lots of 0s by construction |
---|
| 1876 | ring @C = (0,@c(1..size(mons))),dummyvar,dp; |
---|
| 1877 | def WC = save + @C; |
---|
| 1878 | setring WC; |
---|
| 1879 | ideal mons = imap(save,mons); |
---|
| 1880 | poly f; |
---|
| 1881 | for (i=1; i<=size(mons); i++) |
---|
| 1882 | { |
---|
| 1883 | f = f + par(i)*mons[i]; |
---|
| 1884 | } |
---|
| 1885 | // Step 3: determine values of @c(i) by equating coefficients |
---|
| 1886 | ideal I = imap(save,I); |
---|
| 1887 | I = dmodAction(I,f,1..n); |
---|
| 1888 | ideal M = monomialInIdeal(I); |
---|
| 1889 | matrix CC = coeffs(I,M); |
---|
| 1890 | int j; |
---|
| 1891 | ideal C; |
---|
| 1892 | for (i=1; i<=nrows(CC); i++) |
---|
| 1893 | { |
---|
| 1894 | f = 0; |
---|
| 1895 | for (j=1; j<=ncols(CC); j++) |
---|
| 1896 | { |
---|
| 1897 | f = f + CC[i,j]; |
---|
| 1898 | } |
---|
| 1899 | C[size(C)+1] = f; |
---|
| 1900 | } |
---|
| 1901 | // Step 3.1: solve a linear system |
---|
| 1902 | ring rC = 0,(@c(1..size(mons))),dp; |
---|
| 1903 | ideal C = imap(WC,C); |
---|
| 1904 | matrix M = coeffs(C,maxideal(1)); |
---|
| 1905 | module MM = leftKernel(M); |
---|
| 1906 | setring WC; |
---|
| 1907 | module MM = imap(rC,MM); |
---|
| 1908 | // Step 3.2: return the solution |
---|
| 1909 | ideal F = ideal(MM*transpose(mons)); |
---|
| 1910 | setring save; |
---|
| 1911 | ideal F = imap(WC,F); |
---|
| 1912 | return(F); |
---|
| 1913 | } |
---|
| 1914 | example |
---|
| 1915 | { |
---|
| 1916 | "EXAMPLE:"; echo=2; |
---|
| 1917 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
| 1918 | def W = Weyl(); |
---|
| 1919 | setring W; |
---|
| 1920 | poly tx,ty = x*Dx, y*Dy; |
---|
| 1921 | ideal I = // Appel F1 with parameters (2,-3,-2,5) |
---|
| 1922 | tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3), |
---|
| 1923 | ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2), |
---|
| 1924 | (x-y)*Dx*Dy+2*Dx-3*Dy; |
---|
| 1925 | intvec w = -1,-1; |
---|
| 1926 | polSol(I,w); |
---|
| 1927 | } |
---|
| 1928 | |
---|
| 1929 | |
---|
| 1930 | static proc ex_polSol() |
---|
| 1931 | { ring r = 0,(x,y,Dx,Dy),dp; |
---|
| 1932 | def W = Weyl(); |
---|
| 1933 | setring W; |
---|
| 1934 | poly tx,ty = x*Dx, y*Dy; |
---|
| 1935 | ideal I = // Appel F1 with parameters (2,-3,-2,5) |
---|
| 1936 | tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3), |
---|
| 1937 | ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2), |
---|
| 1938 | (x-y)*Dx*Dy+2*Dx-3*Dy; |
---|
| 1939 | intvec w = -5,-7; |
---|
| 1940 | // the following gives a bug |
---|
| 1941 | polSol(I,w); |
---|
| 1942 | // this is due to a bug in weightKB, see ticket #339 |
---|
| 1943 | // http://www.singular.uni-kl.de:8002/trac/ticket/339 |
---|
| 1944 | } |
---|
| 1945 | |
---|
| 1946 | |
---|
| 1947 | proc polSolFiniteRank (ideal I, list #) |
---|
| 1948 | " |
---|
| 1949 | USAGE: polSolFiniteRank(I[,w]); I ideal, w optional intvec |
---|
| 1950 | ASSUME: The basering is the n-th Weyl algebra W over a field of |
---|
| 1951 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 1952 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 1953 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the |
---|
| 1954 | differential operator belonging to x(i). |
---|
| 1955 | @* Moreover, assume that I is of finite holonomic rank. |
---|
| 1956 | RETURN: ideal, a basis of the polynomial solutions to the given system of |
---|
| 1957 | linear PDEs with polynomial coefficients, encoded via I |
---|
| 1958 | REMARKS: If w is given, w should consist of n strictly negative entries. |
---|
| 1959 | Otherwise and by default, w is set to -1:n. |
---|
| 1960 | In this case, w is used as weight vector for the computation of a |
---|
| 1961 | b-function. |
---|
| 1962 | @* Reference: (OTT), Algorithm 2.6 |
---|
| 1963 | NOTE: If printlevel=1, progress debug messages will be printed, |
---|
| 1964 | if printlevel>=2, all the debug messages will be printed. |
---|
| 1965 | SEE ALSO: polSol, ratSol |
---|
| 1966 | EXAMPLE: example polSolFiniteRank; shows examples |
---|
| 1967 | " |
---|
| 1968 | { |
---|
| 1969 | dmodGeneralAssumptionCheck(); |
---|
| 1970 | if (holonomicRank(I)==-1) |
---|
| 1971 | { |
---|
| 1972 | ERROR("Ideal is not of finite holonomic rank."); |
---|
| 1973 | } |
---|
| 1974 | int ppl = printlevel - voice + 2; |
---|
| 1975 | int n = nvars(basering) div 2; |
---|
| 1976 | int eng; |
---|
| 1977 | intvec w = -1:(n div 2); |
---|
| 1978 | if (size(#)>0) |
---|
| 1979 | { |
---|
| 1980 | if (typeof(#[1])=="intvec") |
---|
| 1981 | { |
---|
| 1982 | if (allPositive(-#[1])) |
---|
| 1983 | { |
---|
| 1984 | w = #[1]; |
---|
| 1985 | } |
---|
| 1986 | } |
---|
| 1987 | } |
---|
| 1988 | dbprint(ppl,"// Computing initial ideal..."); |
---|
| 1989 | ideal J = initialIdealW(I,-w,w); |
---|
| 1990 | dbprint(ppl,"// ...done."); |
---|
| 1991 | dbprint(ppl,"// Computing Weyl closure..."); |
---|
| 1992 | J = WeylClosure(J); |
---|
| 1993 | J = engine(J,eng); |
---|
| 1994 | dbprint(ppl,"// ...done."); |
---|
| 1995 | poly s; |
---|
| 1996 | int i; |
---|
| 1997 | for (i=1; i<=n; i++) |
---|
| 1998 | { |
---|
| 1999 | s = s + w[i]*var(i)*var(i+n); |
---|
| 2000 | } |
---|
| 2001 | dbprint(ppl,"// Computing intersection..."); |
---|
| 2002 | vector v = pIntersect(s,J); |
---|
| 2003 | list L = bFactor(vec2poly(v)); |
---|
| 2004 | dbprint(ppl-1,"// roots: " + string(L[1])); |
---|
| 2005 | dbprint(ppl-1,"// multiplicities: " + string(L[2])); |
---|
| 2006 | dbprint(ppl,"// ...done."); |
---|
| 2007 | int mr = minIntRoot2(L); |
---|
| 2008 | int pl = printlevel; |
---|
| 2009 | printlevel = printlevel + 1; |
---|
| 2010 | ideal K = polSol(I,w,mr); |
---|
| 2011 | printlevel = printlevel - 1; |
---|
| 2012 | return(K); |
---|
| 2013 | } |
---|
| 2014 | example |
---|
| 2015 | { |
---|
| 2016 | "EXAMPLE:"; echo=2; |
---|
| 2017 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
| 2018 | def W = Weyl(); |
---|
| 2019 | setring W; |
---|
| 2020 | poly tx,ty = x*Dx, y*Dy; |
---|
| 2021 | ideal I = // Appel F1 with parameters (2,-3,-2,5) |
---|
| 2022 | tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3), |
---|
| 2023 | ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2), |
---|
| 2024 | (x-y)*Dx*Dy+2*Dx-3*Dy; |
---|
| 2025 | intvec w = -1,-1; |
---|
| 2026 | polSolFiniteRank(I,w); |
---|
| 2027 | } |
---|
| 2028 | |
---|
| 2029 | |
---|
| 2030 | static proc twistedIdeal(ideal I, poly f, intvec k, ideal F) |
---|
| 2031 | { |
---|
| 2032 | // I subset D_n, f in K[x], F = factorize(f,1), size(k) = size(F), k[i]>0 |
---|
| 2033 | def save = basering; |
---|
| 2034 | int n = nvars(save) div 2; |
---|
| 2035 | int i,j; |
---|
| 2036 | intvec a,v,w; |
---|
| 2037 | w = (0:n),(1:n); |
---|
| 2038 | for (i=1; i<=size(I); i++) |
---|
| 2039 | { |
---|
| 2040 | a[i] = deg(I[i],w); |
---|
| 2041 | } |
---|
| 2042 | ring FD = 0,(fd(1..n)),dp; |
---|
| 2043 | def @@WFD = save + FD; |
---|
| 2044 | setring @@WFD; |
---|
| 2045 | poly f = imap(save,f); |
---|
| 2046 | list RL = ringlist(basering); |
---|
| 2047 | RL = RL[1..4]; |
---|
| 2048 | list L = RL[3]; |
---|
| 2049 | v = (1:(2*n)),((deg(f)+1):n); |
---|
| 2050 | L = insert(L,list("a",v)); |
---|
| 2051 | RL[3] = L; |
---|
| 2052 | def @WFD = ring(RL); |
---|
| 2053 | setring @WFD; |
---|
| 2054 | poly f = imap(save,f); |
---|
| 2055 | matrix Drel[3*n][3*n]; |
---|
| 2056 | for (i=1; i<=n; i++) |
---|
| 2057 | { |
---|
| 2058 | Drel[i,i+n] = 1; // [D,x] |
---|
| 2059 | Drel[i,i+2*n] = f; // [fD,x] |
---|
| 2060 | for (j=1; j<=n; j++) |
---|
| 2061 | { |
---|
| 2062 | Drel[i+n,j+2*n] = -diff(f,var(i))*var(j+n); // [fD,D] |
---|
| 2063 | Drel[j+2*n,i+2*n] = diff(f,var(i))*var(j+2*n) - diff(f,var(j))*var(i+2*n); |
---|
| 2064 | // [fD,fD] |
---|
| 2065 | } |
---|
| 2066 | } |
---|
| 2067 | def WFD = nc_algebra(1,Drel); |
---|
| 2068 | setring WFD; |
---|
| 2069 | kill @WFD,@@WFD; |
---|
| 2070 | ideal I = imap(save,I); |
---|
| 2071 | poly f = imap(save,f); |
---|
| 2072 | for (i=1; i<=size(I); i++) |
---|
| 2073 | { |
---|
| 2074 | I[i] = f^(a[i])*I[i]; |
---|
| 2075 | } |
---|
| 2076 | ideal S; |
---|
| 2077 | for (i=1; i<=n; i++) |
---|
| 2078 | { |
---|
| 2079 | S[size(S)+1] = var(i+2*n) - f*var(i+n); |
---|
| 2080 | } |
---|
| 2081 | S = slimgb(S); |
---|
| 2082 | I = NF(I,S); |
---|
| 2083 | if (select1(I,intvec((n+1)..2*n))[1] <> 0) |
---|
| 2084 | { |
---|
| 2085 | // should never get here |
---|
| 2086 | ERROR("Something's wrong. Please inform the author."); |
---|
| 2087 | } |
---|
| 2088 | setring save; |
---|
| 2089 | ideal mm = maxideal(1); |
---|
| 2090 | poly s; |
---|
| 2091 | for (i=1; i<=n; i++) |
---|
| 2092 | { |
---|
| 2093 | s = f*var(i+n); |
---|
| 2094 | for (j=1; j<=size(F); j++) |
---|
| 2095 | { |
---|
| 2096 | s = s + k[j]*(f/F[j])*bracket(var(i+n),F[j]); |
---|
| 2097 | } |
---|
| 2098 | mm[i+2*n] = s; |
---|
| 2099 | } |
---|
| 2100 | map m = WFD,mm; |
---|
| 2101 | ideal J = m(I); |
---|
| 2102 | return(J); |
---|
| 2103 | } |
---|
| 2104 | example |
---|
| 2105 | { |
---|
| 2106 | "EXAMPLE"; echo=2; |
---|
| 2107 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
| 2108 | def W = Weyl(); |
---|
| 2109 | setring W; |
---|
| 2110 | poly tx,ty = x*Dx, y*Dy; |
---|
| 2111 | ideal I = // Appel F1 with parameters (3,-1,1,1) is a solution |
---|
| 2112 | tx*(tx+ty)-x*(tx+ty+3)*(tx-1), |
---|
| 2113 | ty*(tx+ty)-y*(tx+ty+3)*(ty+1); |
---|
| 2114 | kill tx,ty; |
---|
| 2115 | poly f = x^3*y^2-x^2*y^3-x^3*y+x*y^3+x^2*y-x*y^2; |
---|
| 2116 | ideal F = x-1,x,-x+y,y-1,y; |
---|
| 2117 | intvec k = -1,-1,-1,-3,-1; |
---|
| 2118 | ideal T = twistedIdeal(I,f,k,F); |
---|
| 2119 | // TODO change the ordering of WFD |
---|
| 2120 | // introduce new var f?? |
---|
| 2121 | //paper: |
---|
| 2122 | poly fx = diff(f,x); |
---|
| 2123 | poly fy = diff(f,y); |
---|
| 2124 | poly fDx = f*Dx; |
---|
| 2125 | poly fDy = f*Dy; |
---|
| 2126 | poly fd(1) = fDx; |
---|
| 2127 | poly fd(2) = fDy; |
---|
| 2128 | ideal K= |
---|
| 2129 | (x^2-x^3)*(fDx)^2+x*((1-3*x)*f-(1-x)*y*fy-(1-x)*x*fx)*(fDx) |
---|
| 2130 | +x*(1-x)*y*(fDy)*(fDx)+x*y*f*(fDy)+3*x*f^2, |
---|
| 2131 | (y^2-y^3)*(fDy)^2+y*((1-5*y)*f-(1-y)*x*fx-(1-y)*y*fy)*(fDy) |
---|
| 2132 | +y*(1-y)*x*(fDx)*(fDy)-y*x*f*(fDx)-3*y*f^2; |
---|
| 2133 | } |
---|
| 2134 | |
---|
| 2135 | |
---|
| 2136 | proc ratSol (ideal I) |
---|
| 2137 | " |
---|
| 2138 | USAGE: ratSol(I); I ideal |
---|
| 2139 | ASSUME: The basering is the n-th Weyl algebra W over a field of |
---|
| 2140 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 2141 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 2142 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the |
---|
| 2143 | differential operator belonging to x(i). |
---|
| 2144 | @* Moreover, assume that I is holonomic. |
---|
| 2145 | RETURN: module, a basis of the rational solutions to the given system of |
---|
| 2146 | linear PDEs with polynomial coefficients, encoded via I |
---|
| 2147 | Note that each entry has two components, the first one standing for |
---|
| 2148 | the enumerator, the second one for the denominator. |
---|
| 2149 | REMARKS: Reference: (OTT), Algorithm 3.10 |
---|
| 2150 | NOTE: If printlevel=1, progress debug messages will be printed, |
---|
| 2151 | if printlevel>=2, all the debug messages will be printed. |
---|
| 2152 | SEE ALSO: polSol, polSolFiniteRank |
---|
| 2153 | EXAMPLE: example ratSol; shows examples |
---|
| 2154 | " |
---|
| 2155 | { |
---|
| 2156 | dmodGeneralAssumptionCheck(); |
---|
| 2157 | if (!isHolonomic(I)) |
---|
| 2158 | { |
---|
| 2159 | ERROR("Ideal is not holonomic."); |
---|
| 2160 | } |
---|
| 2161 | int ppl = printlevel - voice + 2; |
---|
| 2162 | def save = basering; |
---|
| 2163 | dbprint(ppl,"// computing singular locus..."); |
---|
| 2164 | ideal S = DsingularLocus(I); |
---|
| 2165 | dbprint(ppl,"// ...done."); |
---|
| 2166 | poly f = S[1]; |
---|
| 2167 | dbprint(ppl,"// considering poly " + string(f)); |
---|
| 2168 | int n = nvars(save) div 2; |
---|
| 2169 | list RL = ringlist(save); |
---|
| 2170 | RL = RL[1..4]; |
---|
| 2171 | list L = RL[2]; |
---|
| 2172 | L = list(L[1..n]); |
---|
| 2173 | RL[2] = L; |
---|
| 2174 | L = list(); |
---|
| 2175 | L[1] = list("dp",intvec(1:n)); |
---|
| 2176 | L[2] = list("C",intvec(0)); |
---|
| 2177 | RL[3] = L; |
---|
| 2178 | def rr = ring(RL); |
---|
| 2179 | setring rr; |
---|
| 2180 | poly f = imap(save,f); |
---|
| 2181 | ideal F = factorize(f,1); // not interested in multiplicities |
---|
| 2182 | dbprint(ppl,"// with irreducible factors " + string(F)); |
---|
| 2183 | setring save; |
---|
| 2184 | ideal F = imap(rr,F); |
---|
| 2185 | kill rr,RL; |
---|
| 2186 | int i; |
---|
| 2187 | intvec k; |
---|
| 2188 | ideal FF = 1,1; |
---|
| 2189 | dbprint(ppl,"// computing b-functions of irreducible factors..."); |
---|
| 2190 | for (i=1; i<=size(F); i++) |
---|
| 2191 | { |
---|
| 2192 | dbprint(ppl,"// considering " + string(F[i]) + "..."); |
---|
| 2193 | L = bfctBound(I,F[i]); |
---|
| 2194 | if (size(L) == 3) // bfct is constant |
---|
| 2195 | { |
---|
| 2196 | dbprint(ppl,"// ...got " + string(L[3])); |
---|
| 2197 | if (L[3] == "1") |
---|
| 2198 | { |
---|
| 2199 | return(0); // TODO type // no rational solutions |
---|
| 2200 | } |
---|
| 2201 | else // should never get here |
---|
| 2202 | { |
---|
| 2203 | ERROR("Oops, something went wrong. Please inform the author."); |
---|
| 2204 | } |
---|
| 2205 | } |
---|
| 2206 | else |
---|
| 2207 | { |
---|
| 2208 | dbprint(ppl,"// ...got roots " + string(L[1])); |
---|
| 2209 | dbprint(ppl,"// with multiplicities " + string(L[2])); |
---|
| 2210 | k[i] = -maxIntRoot(L)-1; |
---|
| 2211 | if (k[i] < 0) |
---|
| 2212 | { |
---|
| 2213 | FF[2] = FF[2]*F[i]^(-k[i]); |
---|
| 2214 | } |
---|
| 2215 | else |
---|
| 2216 | { |
---|
| 2217 | FF[1] = FF[1]*F[i]^(k[i]); |
---|
| 2218 | } |
---|
| 2219 | } |
---|
| 2220 | } |
---|
| 2221 | vector v = FF[1]*gen(1) + FF[2]*gen(2); |
---|
| 2222 | kill FF; |
---|
| 2223 | dbprint(ppl,"// ...done"); |
---|
| 2224 | ideal twI = twistedIdeal(I,f,k,F); |
---|
| 2225 | intvec w = -1:n; |
---|
| 2226 | dbprint(ppl,"// computing polynomial solutions of twisted system..."); |
---|
| 2227 | if (isHolonomic(twI)) |
---|
| 2228 | { |
---|
| 2229 | ideal P = polSol(twI,w); |
---|
| 2230 | } |
---|
| 2231 | else |
---|
| 2232 | { |
---|
| 2233 | ideal P = polSolFiniteRank(twI,w); |
---|
| 2234 | } |
---|
| 2235 | module M; |
---|
| 2236 | vector vv; |
---|
| 2237 | for (i=1; i<=ncols(P); i++) |
---|
| 2238 | { |
---|
| 2239 | vv = P[i]*gen(1) + 1*gen(2); |
---|
| 2240 | M[i] = multRat(v,vv); |
---|
| 2241 | } |
---|
| 2242 | dbprint(ppl,"// ...done"); |
---|
| 2243 | return (M); |
---|
| 2244 | } |
---|
| 2245 | example |
---|
| 2246 | { |
---|
| 2247 | "EXAMPLE"; echo=2; |
---|
| 2248 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
| 2249 | def W = Weyl(); |
---|
| 2250 | setring W; |
---|
| 2251 | poly tx,ty = x*Dx, y*Dy; |
---|
| 2252 | ideal I = // Appel F1 with parameters (3,-1,1,1) is a solution |
---|
| 2253 | tx*(tx+ty)-x*(tx+ty+3)*(tx-1), |
---|
| 2254 | ty*(tx+ty)-y*(tx+ty+3)*(ty+1); |
---|
| 2255 | module M = ratSol(I); |
---|
| 2256 | // We obtain a basis of the rational solutions to I represented by a |
---|
| 2257 | // module / matrix with two rows. |
---|
| 2258 | // Each column of the matrix represents a rational function, where |
---|
| 2259 | // the first row correspond to the enumerator and the second row to |
---|
| 2260 | // the denominator. |
---|
| 2261 | print(M); |
---|
| 2262 | } |
---|
| 2263 | |
---|
| 2264 | |
---|
| 2265 | proc bfctBound (ideal I, poly f, list #) |
---|
| 2266 | " |
---|
| 2267 | USAGE: bfctBound (I,f[,primdec]); I ideal, f poly, primdec optional string |
---|
| 2268 | ASSUME: The basering is the n-th Weyl algebra W over a field of |
---|
| 2269 | characteristic 0 and for all 1<=i<=n the identity |
---|
| 2270 | var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of |
---|
| 2271 | variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the |
---|
| 2272 | differential operator belonging to x(i). |
---|
| 2273 | @* Moreover, assume that I is holonomic. |
---|
| 2274 | RETURN: list of roots (of type ideal) and multiplicities (of type intvec) of |
---|
| 2275 | a multiple of the b-function for f^s*u at a generic root of f. |
---|
| 2276 | Here, u stands for [1] in D/I. |
---|
| 2277 | REMARKS: Reference: (OTT), Algorithm 3.4 |
---|
| 2278 | NOTE: This procedure requires to compute a primary decomposition in a |
---|
| 2279 | commutative ring. The optional string primdec can be used to specify |
---|
| 2280 | the algorithm to do so. It may either be `GTZ' (Gianni, Trager, |
---|
| 2281 | Zacharias) or `SY' (Shimoyama, Yokoyama). By default, `GTZ' is used. |
---|
| 2282 | @* If printlevel=1, progress debug messages will be printed, |
---|
| 2283 | if printlevel>=2, all the debug messages will be printed. |
---|
| 2284 | SEE ALSO: bernstein, bfct, bfctAnn |
---|
| 2285 | EXAMPLE: example bfctBound; shows examples |
---|
| 2286 | " |
---|
| 2287 | { |
---|
| 2288 | dmodGeneralAssumptionCheck(); |
---|
| 2289 | finKx(f); |
---|
| 2290 | if (!isHolonomic(I)) |
---|
| 2291 | { |
---|
| 2292 | ERROR("Ideal is not holonomic."); |
---|
| 2293 | } |
---|
| 2294 | int ppl = printlevel - voice + 2; |
---|
| 2295 | string primdec = "GTZ"; |
---|
| 2296 | if (size(#)>1) |
---|
| 2297 | { |
---|
| 2298 | if (typeof(#[1])=="string") |
---|
| 2299 | { |
---|
| 2300 | if ( (#[1]=="SY") || (#[1]=="sy") || (#[1]=="Sy") ) |
---|
| 2301 | { |
---|
| 2302 | primdec = "SY"; |
---|
| 2303 | } |
---|
| 2304 | else |
---|
| 2305 | { |
---|
| 2306 | if ( (#[1]<>"GTZ") && (#[1]<>"gtz") && (#[1]<>"Gtz") ) |
---|
| 2307 | { |
---|
| 2308 | print("// Warning: optional string may either be `GTZ' or `SY',"); |
---|
| 2309 | print("// proceeding with `GTZ'."); |
---|
| 2310 | primdec = "GTZ"; |
---|
| 2311 | } |
---|
| 2312 | } |
---|
| 2313 | } |
---|
| 2314 | } |
---|
| 2315 | def save = basering; |
---|
| 2316 | int n = nvars(save) div 2; |
---|
| 2317 | // step 1 |
---|
| 2318 | ideal mm = maxideal(1); |
---|
| 2319 | def Wt = extendWeyl(safeVarName("t")); |
---|
| 2320 | setring Wt; |
---|
| 2321 | poly f = imap(save,f); |
---|
| 2322 | ideal mm = imap(save,mm); |
---|
| 2323 | int i; |
---|
| 2324 | for (i=1; i<=n; i++) |
---|
| 2325 | { |
---|
| 2326 | mm[i+n] = var(i+n+2) + bracket(var(i+n+2),f)*var(n+2); |
---|
| 2327 | } |
---|
| 2328 | map m = save,mm; |
---|
| 2329 | ideal I = m(I); |
---|
| 2330 | I = I, var(1)-f; |
---|
| 2331 | // step 2 |
---|
| 2332 | intvec w = 1,(0:n); |
---|
| 2333 | dbprint(ppl ,"// Computing initial ideal..."); |
---|
| 2334 | I = initialIdealW(I,-w,w); |
---|
| 2335 | dbprint(ppl ,"// ...done."); |
---|
| 2336 | dbprint(ppl-1,"// " + string(I)); |
---|
| 2337 | // step 3: rewrite I using Euler operator t*Dt |
---|
| 2338 | list RL = ringlist(Wt); |
---|
| 2339 | RL = RL[1..4]; |
---|
| 2340 | list L = RL[2] + list(safeVarName("s")); // s=t*Dt |
---|
| 2341 | RL[2] = L; |
---|
| 2342 | L = list(); |
---|
| 2343 | L[1] = list("dp",intvec(1:(2*n+2))); |
---|
| 2344 | L[2] = list("dp",intvec(1)); |
---|
| 2345 | L[3] = list("C",intvec(0)); |
---|
| 2346 | RL[3] = L; |
---|
| 2347 | def @Wts = ring(RL); |
---|
| 2348 | kill L,RL; |
---|
| 2349 | setring @Wts; |
---|
| 2350 | matrix relD[2*n+3][2*n+3]; |
---|
| 2351 | relD[1,2*n+3] = var(1); |
---|
| 2352 | relD[n+2,2*n+3] = -var(n+2); |
---|
| 2353 | for (i=1; i<=n+1; i++) |
---|
| 2354 | { |
---|
| 2355 | relD[i,n+i+1] = 1; |
---|
| 2356 | } |
---|
| 2357 | def Wts = nc_algebra(1,relD); |
---|
| 2358 | setring Wts; |
---|
| 2359 | ideal I = imap(Wt,I); |
---|
| 2360 | kill Wt,@Wts; |
---|
| 2361 | ideal S = var(1)*var(n+2)-var(2*n+3); |
---|
| 2362 | attrib(S,"isSB",1); |
---|
| 2363 | dbprint(ppl ,"// Computing Euler representation..."); |
---|
| 2364 | // I = NF(I,S); |
---|
| 2365 | int d; |
---|
| 2366 | intvec ww = 0:(2*2+2); ww[1] = -1; ww[n+2] = 1; |
---|
| 2367 | for (i=1; i<=size(I); i++) |
---|
| 2368 | { |
---|
| 2369 | d = deg(I[i],ww); |
---|
| 2370 | if (d>0) |
---|
| 2371 | { |
---|
| 2372 | I[i] = var(1)^d*I[i]; |
---|
| 2373 | } |
---|
| 2374 | if (d<0) |
---|
| 2375 | { |
---|
| 2376 | d = -d; |
---|
| 2377 | I[i] = var(n+2)^d*I[i]; |
---|
| 2378 | } |
---|
| 2379 | } |
---|
| 2380 | I = NF(I,S); // now there are no t,Dt in I, only s |
---|
| 2381 | dbprint(ppl ,"// ...done."); |
---|
| 2382 | I = subst(I,var(2*n+3),-var(2*n+3)-1); |
---|
| 2383 | ring Ks = 0,s,dp; |
---|
| 2384 | def Ws = save + Ks; |
---|
| 2385 | setring Ws; |
---|
| 2386 | ideal I = imap(Wts,I); |
---|
| 2387 | kill Wts; |
---|
| 2388 | poly DD = 1; |
---|
| 2389 | for (i=1; i<=n; i++) |
---|
| 2390 | { |
---|
| 2391 | DD = DD * var(n+i); |
---|
| 2392 | } |
---|
| 2393 | dbprint(ppl ,"// Eliminating differential operators..."); |
---|
| 2394 | ideal J = eliminate(I,DD); // J subset K[x,s] |
---|
| 2395 | dbprint(ppl ,"// ...done."); |
---|
| 2396 | dbprint(ppl-1,"// " + string(J)); |
---|
| 2397 | list RL = ringlist(Ws); |
---|
| 2398 | RL = RL[1..4]; |
---|
| 2399 | list L = RL[2]; |
---|
| 2400 | L = list(L[1..n]) + list(L[2*n+1]); |
---|
| 2401 | RL[2] = L; |
---|
| 2402 | L = list(); |
---|
| 2403 | L[1] = list("dp",intvec(1:(n+1))); |
---|
| 2404 | L[2] = list("C",intvec(0)); |
---|
| 2405 | RL[3] = L; |
---|
| 2406 | def Kxs = ring(RL); |
---|
| 2407 | setring Kxs; |
---|
| 2408 | ideal J = imap(Ws,J); |
---|
| 2409 | dbprint(ppl ,"// Computing primary decomposition with engine " |
---|
| 2410 | + primdec + "..."); |
---|
| 2411 | if (primdec == "GTZ") |
---|
| 2412 | { |
---|
| 2413 | list P = primdecGTZ(J); |
---|
| 2414 | } |
---|
| 2415 | else |
---|
| 2416 | { |
---|
| 2417 | list P = primdecSY(J); |
---|
| 2418 | } |
---|
| 2419 | dbprint(ppl ,"// ...done."); |
---|
| 2420 | dbprint(ppl-1,"// " + string(P)); |
---|
| 2421 | ideal GP,Qix,rad,B; |
---|
| 2422 | poly f = imap(save,f); |
---|
| 2423 | vector v; |
---|
| 2424 | for (i=1; i<=size(P); i++) |
---|
| 2425 | { |
---|
| 2426 | dbprint(ppl ,"// Considering primary component " + string(i) |
---|
| 2427 | + " of " + string(size(P)) + "..."); |
---|
| 2428 | dbprint(ppl ,"// Intersecting with K[x] and computing radical..."); |
---|
| 2429 | GP = std(P[i][1]); |
---|
| 2430 | Qix = eliminate(GP,var(n+1)); // subset K[x] |
---|
| 2431 | rad = radical(Qix); |
---|
| 2432 | rad = std(rad); |
---|
| 2433 | dbprint(ppl ,"// ...done."); |
---|
| 2434 | dbprint(ppl-1,"// " + string(rad)); |
---|
| 2435 | if (rad[1]==0 || NF(f,rad)==0) |
---|
| 2436 | { |
---|
| 2437 | dbprint(ppl ,"// Intersecting with K[s]..."); |
---|
| 2438 | v = pIntersect(var(n+1),GP); |
---|
| 2439 | B[size(B)+1] = vec2poly(v,n+1); |
---|
| 2440 | dbprint(ppl ,"// ...done."); |
---|
| 2441 | dbprint(ppl-1,"// " + string(B[size(B)])); |
---|
| 2442 | } |
---|
| 2443 | dbprint(ppl ,"// ...done."); |
---|
| 2444 | } |
---|
| 2445 | f = lcm(B); // =lcm(B[1],...,B[size(B)]) |
---|
| 2446 | list bb = bFactor(f); |
---|
| 2447 | setring save; |
---|
| 2448 | list bb = imap(Kxs,bb); |
---|
| 2449 | return(bb); |
---|
| 2450 | } |
---|
| 2451 | example |
---|
| 2452 | { |
---|
| 2453 | "EXAMPLE"; echo=2; |
---|
| 2454 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
| 2455 | def W = Weyl(); |
---|
| 2456 | setring W; |
---|
| 2457 | poly tx,ty = x*Dx, y*Dy; |
---|
| 2458 | ideal I = // Appel F1 with parameters (2,-3,-2,5) |
---|
| 2459 | tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3), |
---|
| 2460 | ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2), |
---|
| 2461 | (x-y)*Dx*Dy+2*Dx-3*Dy; |
---|
| 2462 | kill tx,ty; |
---|
| 2463 | poly f = x-1; |
---|
| 2464 | bfctBound(I,f); |
---|
| 2465 | } |
---|
| 2466 | |
---|
| 2467 | |
---|
| 2468 | //TODO check f/g or g/f, check Weyl closure of result |
---|
| 2469 | proc annRatSyz (poly f, poly g, list #) |
---|
| 2470 | " |
---|
| 2471 | USAGE: annRatSyz(f,g[,db,eng]); f, g polynomials, db,eng optional integers |
---|
| 2472 | ASSUME: The basering is commutative and over a field of characteristic 0. |
---|
| 2473 | RETURN: ring (a Weyl algebra) containing an ideal `LD', which is (part of) |
---|
| 2474 | the annihilator of the rational function g/f in the corresponding |
---|
| 2475 | Weyl algebra |
---|
| 2476 | REMARKS: This procedure uses the computation of certain syzygies. |
---|
| 2477 | One can obtain the full annihilator by computing the Weyl closure of |
---|
| 2478 | the ideal LD. |
---|
| 2479 | NOTE: Activate the output ring with the @code{setring} command. |
---|
| 2480 | In the output ring, the ideal `LD' (in Groebner basis) is (part of) |
---|
| 2481 | the annihilator of g/f. |
---|
| 2482 | @* If db>0 is given, operators of order up to db are considered, |
---|
| 2483 | otherwise, and by default, a minimal holonomic solution is computed. |
---|
| 2484 | @* If eng<>0, @code{std} is used for Groebner basis computations, |
---|
| 2485 | otherwise, and by default, @code{slimgb} is used. |
---|
| 2486 | @* If printlevel =1, progress debug messages will be printed, |
---|
| 2487 | if printlevel>=2, all the debug messages will be printed. |
---|
| 2488 | SEE ALSO: annRat, annPoly |
---|
| 2489 | EXAMPLE: example annRatSyz; shows examples |
---|
| 2490 | " |
---|
| 2491 | { |
---|
| 2492 | // check assumptions |
---|
| 2493 | if (!isCommutative()) |
---|
| 2494 | { |
---|
| 2495 | ERROR("Basering must be commutative."); |
---|
| 2496 | } |
---|
| 2497 | if ( (size(ideal(basering)) >0) || (char(basering) >0) ) |
---|
| 2498 | { |
---|
| 2499 | ERROR("Basering is inappropriate: characteristic>0 or qring present."); |
---|
| 2500 | } |
---|
| 2501 | if (g == 0) |
---|
| 2502 | { |
---|
| 2503 | ERROR("Second polynomial must not be zero."); |
---|
| 2504 | } |
---|
| 2505 | int db,eng; |
---|
| 2506 | if (size(#)>0) |
---|
| 2507 | { |
---|
| 2508 | if (typeof(#[1]) == "int") |
---|
| 2509 | { |
---|
| 2510 | db = int(#[1]); |
---|
| 2511 | } |
---|
| 2512 | if (size(#)>1) |
---|
| 2513 | { |
---|
| 2514 | if (typeof(#[2]) == "int") |
---|
| 2515 | { |
---|
| 2516 | eng = int(#[1]); |
---|
| 2517 | } |
---|
| 2518 | } |
---|
| 2519 | } |
---|
| 2520 | int ppl = printlevel - voice + 2; |
---|
| 2521 | vector I = f*gen(1)+g*gen(2); |
---|
| 2522 | checkRatInput(I); |
---|
| 2523 | int i,j; |
---|
| 2524 | def R = basering; |
---|
| 2525 | int n = nvars(R); |
---|
| 2526 | list RL = ringlist(R); |
---|
| 2527 | RL = RL[1..4]; |
---|
| 2528 | list Ltmp = RL[2]; |
---|
| 2529 | for (i=1; i<=n; i++) |
---|
| 2530 | { |
---|
| 2531 | Ltmp[i+n] = safeVarName("D" + Ltmp[i]); |
---|
| 2532 | } |
---|
| 2533 | RL[2] = Ltmp; |
---|
| 2534 | Ltmp = list(); |
---|
| 2535 | Ltmp[1] = list("dp",intvec(1:2*n)); |
---|
| 2536 | Ltmp[2] = list("C",intvec(0)); |
---|
| 2537 | RL[3] = Ltmp; |
---|
| 2538 | kill Ltmp; |
---|
| 2539 | def @D = ring(RL); |
---|
| 2540 | setring @D; |
---|
| 2541 | def D = Weyl(); |
---|
| 2542 | setring D; |
---|
| 2543 | ideal DD = 1; |
---|
| 2544 | ideal Dcd,Dnd,LD,tmp; |
---|
| 2545 | Dnd = 1; |
---|
| 2546 | module DS; |
---|
| 2547 | poly DJ; |
---|
| 2548 | kill @D; |
---|
| 2549 | setring R; |
---|
| 2550 | module Rnd,Rcd; |
---|
| 2551 | Rnd[1] = I; |
---|
| 2552 | vector RJ; |
---|
| 2553 | ideal L = I[1]; |
---|
| 2554 | module RS; |
---|
| 2555 | poly p,pnew; |
---|
| 2556 | pnew = I[2]; |
---|
| 2557 | int k,c; |
---|
| 2558 | while(1) |
---|
| 2559 | { |
---|
| 2560 | k++; |
---|
| 2561 | setring R; |
---|
| 2562 | dbprint(ppl,"// Testing order: " + string(k)); |
---|
| 2563 | Rcd = Rnd; |
---|
| 2564 | Rnd = 0; |
---|
| 2565 | setring D; |
---|
| 2566 | Dcd = Dnd; |
---|
| 2567 | Dnd = 0; |
---|
| 2568 | dbprint(ppl-1,"// Current members of the annihilator: " + string(LD)); |
---|
| 2569 | setring R; |
---|
| 2570 | c = size(Rcd); |
---|
| 2571 | p = pnew; |
---|
| 2572 | for (i=1; i<=c; i++) |
---|
| 2573 | { |
---|
| 2574 | for (j=1; j<=n; j++) |
---|
| 2575 | { |
---|
| 2576 | RJ = diffRat(Rcd[i],j); |
---|
| 2577 | setring D; |
---|
| 2578 | DJ = Dcd[i]*var(n+j); |
---|
| 2579 | tmp = Dnd,DJ; |
---|
| 2580 | if (size(Dnd) <> size(simplify(tmp,4))) // new element |
---|
| 2581 | { |
---|
| 2582 | Dnd[size(Dnd)+1] = DJ; |
---|
| 2583 | setring R; |
---|
| 2584 | Rnd[size(Rnd)+1] = RJ; |
---|
| 2585 | pnew = lcm(pnew,RJ[2]); |
---|
| 2586 | } |
---|
| 2587 | else // already have DJ in Dnd |
---|
| 2588 | { |
---|
| 2589 | setring R; |
---|
| 2590 | } |
---|
| 2591 | } |
---|
| 2592 | } |
---|
| 2593 | p = pnew/p; |
---|
| 2594 | for (i=1; i<=size(L); i++) |
---|
| 2595 | { |
---|
| 2596 | L[i] = p*L[i]; |
---|
| 2597 | } |
---|
| 2598 | for (i=1; i<=size(Rnd); i++) |
---|
| 2599 | { |
---|
| 2600 | L[size(L)+1] = pnew/Rnd[i][2]*Rnd[i][1]; |
---|
| 2601 | } |
---|
| 2602 | RS = syz(L); |
---|
| 2603 | setring D; |
---|
| 2604 | DD = DD,Dnd; |
---|
| 2605 | setring R; |
---|
| 2606 | if (RS <> 0) |
---|
| 2607 | { |
---|
| 2608 | setring D; |
---|
| 2609 | DS = imap(R,RS); |
---|
| 2610 | LD = ideal(transpose(DS)*transpose(DD)); |
---|
| 2611 | } |
---|
| 2612 | else |
---|
| 2613 | { |
---|
| 2614 | setring D; |
---|
| 2615 | } |
---|
| 2616 | LD = engine(LD,eng); |
---|
| 2617 | // test if we're done |
---|
| 2618 | if (db<=0) |
---|
| 2619 | { |
---|
| 2620 | if (isHolonomic(LD)) { break; } |
---|
| 2621 | } |
---|
| 2622 | else |
---|
| 2623 | { |
---|
| 2624 | if (k==db) { break; } |
---|
| 2625 | } |
---|
| 2626 | } |
---|
| 2627 | export(LD); |
---|
| 2628 | setring R; |
---|
| 2629 | return(D); |
---|
| 2630 | } |
---|
| 2631 | example |
---|
| 2632 | { |
---|
| 2633 | "EXAMPLE:"; echo = 2; |
---|
| 2634 | // printlevel = 3; |
---|
| 2635 | ring r = 0,(x,y),dp; |
---|
| 2636 | poly f = 2*x*y; poly g = x^2 - y^3; |
---|
| 2637 | def A = annRatSyz(f,g); // compute a holonomic solution |
---|
| 2638 | setring A; A; |
---|
| 2639 | LD; |
---|
| 2640 | setring r; |
---|
| 2641 | def B = annRatSyz(f,g,5); // compute a solution up to degree 5 |
---|
| 2642 | setring B; |
---|
| 2643 | LD; // this is the full annihilator as we will check below |
---|
| 2644 | setring r; |
---|
| 2645 | def C = annRat(f,g); setring C; |
---|
| 2646 | LD; // the full annihilator |
---|
| 2647 | ideal BLD = imap(B,LD); |
---|
| 2648 | NF(LD,std(BLD)); |
---|
| 2649 | } |
---|
| 2650 | |
---|