1 | ////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id: dmodapp.lib,v 1.31 2009-07-07 20:46:51 levandov Exp $"; |
---|
3 | category="Noncommutative"; |
---|
4 | info=" |
---|
5 | LIBRARY: dmodapp.lib Applications of algebraic D-modules |
---|
6 | AUTHORS: Viktor Levandovskyy, levandov@math.rwth-aachen.de |
---|
7 | @* Daniel Andres, daniel.andres@math.rwth-aachen.de |
---|
8 | |
---|
9 | GUIDE: Let K be a field of characteristic 0, R = K[x1,..xN] and |
---|
10 | @* D be the Weyl algebra in variables x1,..xN,d1,..dN. |
---|
11 | @* In this library there are the following procedures for algebraic D-modules: |
---|
12 | @* - localization of a holonomic module D/I with respect to a mult. closed set |
---|
13 | @* of all powers of a given polynomial F from R. Our aim is to compute an |
---|
14 | @* ideal L in D, such that D/L is a presentation of a localized module. Such L |
---|
15 | @* always exists, since such localizations are known to be holonomic and thus |
---|
16 | @* cyclic modules. The procedures for the localization are DLoc,SDLoc and DLoc0. |
---|
17 | @* |
---|
18 | @* - annihilator in D of a given polynomial F from R as well as |
---|
19 | @* of a given rational function G/F from Quot(R). These can be computed via |
---|
20 | @* procedures annPoly resp. annRat. |
---|
21 | @* |
---|
22 | @* - initial form and initial ideals in Weyl algebras with respect to a given |
---|
23 | @* weight vector can be computed with inForm, initialMalgrange, initialIdealW. |
---|
24 | @* |
---|
25 | @* - appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebras, |
---|
26 | @* which annihilate corresponding Appel hypergeometric functions. |
---|
27 | |
---|
28 | REFERENCES: |
---|
29 | @* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric |
---|
30 | @* Differential Equations', Springer, 2000 |
---|
31 | @* (ONW) Oaku, Takayama, Walther 'A Localization Algorithm for D-modules', 2000 |
---|
32 | |
---|
33 | MAIN PROCEDURES: |
---|
34 | |
---|
35 | annPoly(f); annihilator of a polynomial f in the corr. Weyl algebra |
---|
36 | annRat(f,g); annihilator of a rational function f/g in the corr. Weyl algebra |
---|
37 | DLoc(I,F); presentation of the localization of D/I w.r.t. f^s |
---|
38 | SDLoc(I, F); a generic presentation of the localization of D/I w.r.t. f^s |
---|
39 | DLoc0(I, F); presentation of the localization of D/I w.r.t. f^s, based on SDLoc |
---|
40 | |
---|
41 | initialMalgrange(f[,s,t,v]); Groebner basis of the initial Malgrange ideal for f |
---|
42 | initialIdealW(I,u,v[,s,t]); initial ideal of a given ideal w.r.t. given weights |
---|
43 | inForm(f,w); initial form of a poly/ideal w.r.t. a given weight |
---|
44 | isFsat(I, F); check whether the ideal I is F-saturated |
---|
45 | |
---|
46 | AUXILIARY PROCEDURES: |
---|
47 | |
---|
48 | bFactor(F); computes the roots of irreducible factors of an univariate poly |
---|
49 | appelF1(); create an ideal annihilating Appel F1 function |
---|
50 | appelF2(); create an ideal annihilating Appel F2 function |
---|
51 | appelF4(); create an ideal annihilating Appel F4 function |
---|
52 | engine(I,i); computes a Groebner basis with the algorithm given by i |
---|
53 | poly2list(f); decompose a polynomial into a list of terms and exponents |
---|
54 | fl2poly(L,s); reconstruct a monic univariate polynomial from its factorization |
---|
55 | insertGenerator(id,p[,k]); insert an element into an ideal/module |
---|
56 | deleteGenerator(id,k); delete the k-th element from an ideal/module |
---|
57 | |
---|
58 | |
---|
59 | SEE ALSO: dmod_lib, gmssing_lib, bfun_lib |
---|
60 | "; |
---|
61 | |
---|
62 | LIB "poly.lib"; |
---|
63 | LIB "sing.lib"; |
---|
64 | LIB "primdec.lib"; |
---|
65 | LIB "dmod.lib"; // loads e.g. nctools.lib |
---|
66 | LIB "bfun.lib"; //formerly LIB "bfct.lib"; |
---|
67 | LIB "nctools.lib"; // for isWeyl etc |
---|
68 | LIB "gkdim.lib"; |
---|
69 | |
---|
70 | // todo: complete and include into above list |
---|
71 | // charVariety(I); compute the characteristic variety of the ideal I |
---|
72 | // charInfo(); ??? |
---|
73 | |
---|
74 | |
---|
75 | /////////////////////////////////////////////////////////////////////////////// |
---|
76 | // testing for consistency of the library: |
---|
77 | proc testdmodapp() |
---|
78 | { |
---|
79 | example initialIdealW; |
---|
80 | example initialMalgrange; |
---|
81 | example DLoc; |
---|
82 | example DLoc0; |
---|
83 | example SDLoc; |
---|
84 | example inForm; |
---|
85 | example isFsat; |
---|
86 | example annRat; |
---|
87 | example annPoly; |
---|
88 | example appelF1; |
---|
89 | example appelF2; |
---|
90 | example appelF4; |
---|
91 | example poly2list; |
---|
92 | example fl2poly; |
---|
93 | example insertGenerator; |
---|
94 | example deleteGenerator; |
---|
95 | example bFactor; |
---|
96 | } |
---|
97 | |
---|
98 | proc inForm (ideal I, intvec w) |
---|
99 | "USAGE: inForm(I,w); I ideal, w intvec |
---|
100 | RETURN: the initial form of I w.r.t. the weight vector w |
---|
101 | PURPOSE: computes the initial form of an ideal w.r.t. a given weight vector |
---|
102 | NOTE: the size of the weight vector must be equal to the number of variables |
---|
103 | @* of the basering. |
---|
104 | EXAMPLE: example inForm; shows examples |
---|
105 | " |
---|
106 | { |
---|
107 | if (size(w) != nvars(basering)) |
---|
108 | { |
---|
109 | ERROR("weight vector has wrong dimension"); |
---|
110 | } |
---|
111 | if (I == 0) |
---|
112 | { |
---|
113 | return(I); |
---|
114 | } |
---|
115 | int j,i,s,m; |
---|
116 | list l; |
---|
117 | poly g; |
---|
118 | ideal J; |
---|
119 | for (j=1; j<=ncols(I); j++) |
---|
120 | { |
---|
121 | l = poly2list(I[j]); |
---|
122 | m = scalarProd(w,l[1][1]); |
---|
123 | g = l[1][2]; |
---|
124 | for (i=2; i<=size(l); i++) |
---|
125 | { |
---|
126 | s = scalarProd(w,l[i][1]); |
---|
127 | if (s == m) |
---|
128 | { |
---|
129 | g = g + l[i][2]; |
---|
130 | } |
---|
131 | else |
---|
132 | { |
---|
133 | if (s > m) |
---|
134 | { |
---|
135 | m = s; |
---|
136 | g = l[i][2]; |
---|
137 | } |
---|
138 | } |
---|
139 | } |
---|
140 | J[j] = g; |
---|
141 | } |
---|
142 | return(J); |
---|
143 | } |
---|
144 | example |
---|
145 | { |
---|
146 | "EXAMPLE:"; echo = 2; |
---|
147 | ring @D = 0,(x,y,Dx,Dy),dp; |
---|
148 | def D = Weyl(); |
---|
149 | setring D; |
---|
150 | poly F = 3*x^2*Dy+2*y*Dx; |
---|
151 | poly G = 2*x*Dx+3*y*Dy+6; |
---|
152 | ideal I = F,G; |
---|
153 | intvec w1 = -1,-1,1,1; |
---|
154 | intvec w2 = -1,-2,1,2; |
---|
155 | intvec w3 = -2,-3,2,3; |
---|
156 | inForm(I,w1); |
---|
157 | inForm(I,w2); |
---|
158 | inForm(I,w3); |
---|
159 | } |
---|
160 | |
---|
161 | /* |
---|
162 | |
---|
163 | proc charVariety(ideal I) |
---|
164 | "USAGE: charVariety(I); I an ideal |
---|
165 | RETURN: ring |
---|
166 | PURPOSE: compute the characteristic variety of a D-module D/I |
---|
167 | STATUS: experimental, todo |
---|
168 | ASSUME: the ground ring is the Weyl algebra with x's before d's |
---|
169 | NOTE: activate the output ring with the @code{setring} command. |
---|
170 | @* In the output (in a commutative ring): |
---|
171 | @* - the ideal CV is the characteristic variety char(I) |
---|
172 | @* If @code{printlevel}=1, progress debug messages will be printed, |
---|
173 | @* if @code{printlevel}>=2, all the debug messages will be printed. |
---|
174 | EXAMPLE: example charVariety; shows examples |
---|
175 | " |
---|
176 | { |
---|
177 | // 1. introduce the weights 0, 1 |
---|
178 | def save = basering; |
---|
179 | list LL = ringlist(save); |
---|
180 | list L; |
---|
181 | int i; |
---|
182 | for(i=1;i<=4;i++) |
---|
183 | { |
---|
184 | L[i] = LL[i]; |
---|
185 | } |
---|
186 | list OLD = L[3]; |
---|
187 | list NEW; list tmp; |
---|
188 | tmp[1] = "a"; // string |
---|
189 | intvec iv; |
---|
190 | int N = nvars(basering); N = N div 2; |
---|
191 | for(i=N+1; i<=2*N; i++) |
---|
192 | { |
---|
193 | iv[i] = 1; |
---|
194 | } |
---|
195 | tmp[2] = iv; |
---|
196 | NEW[1] = tmp; |
---|
197 | for (i=2; i<=size(OLD);i++) |
---|
198 | { |
---|
199 | NEW[i] = OLD[i-1]; |
---|
200 | } |
---|
201 | L[3] = NEW; |
---|
202 | list ncr =ncRelations(save); |
---|
203 | matrix @C = ncr[1]; |
---|
204 | matrix @D = ncr[2]; |
---|
205 | def @U = ring(L); |
---|
206 | // 2. create the commutative ring |
---|
207 | setring save; |
---|
208 | list CL; |
---|
209 | for(i=1;i<=4;i++) |
---|
210 | { |
---|
211 | CL[i] = L[i]; |
---|
212 | } |
---|
213 | CL[3] = OLD; |
---|
214 | def @CU = ring(CL); |
---|
215 | // comm ring is ready |
---|
216 | setring @U; |
---|
217 | // make @U noncommutative |
---|
218 | matrix @C = imap(save,@C); |
---|
219 | matrix @D = imap(save,@D); |
---|
220 | def @@U = nc_algebra(@C,@D); |
---|
221 | setring @@U; kill @U; |
---|
222 | // 2. compute Groebner basis |
---|
223 | ideal I = imap(save,I); |
---|
224 | // I = groebner(I); |
---|
225 | I = slimgb(I); // a bug? |
---|
226 | setring @CU; |
---|
227 | ideal CV = imap(@@U,I); |
---|
228 | // CV = groebner(CV); // cosmetics |
---|
229 | CV = slimgb(CV); |
---|
230 | export CV; |
---|
231 | return(@CU); |
---|
232 | } |
---|
233 | example |
---|
234 | { |
---|
235 | "EXAMPLE:"; echo = 2; |
---|
236 | ring r = 0,(x,y),Dp; |
---|
237 | poly F = x3-y2; |
---|
238 | printlevel = 0; |
---|
239 | def A = annfs(F); |
---|
240 | setring A; // Weyl algebra |
---|
241 | LD; // the ideal |
---|
242 | def CA = charVariety(LD); |
---|
243 | setring CA; |
---|
244 | CV; |
---|
245 | dim(CV); |
---|
246 | } |
---|
247 | |
---|
248 | /* |
---|
249 | |
---|
250 | // TODO |
---|
251 | |
---|
252 | /* |
---|
253 | proc charInfo(ideal I) |
---|
254 | "USAGE: charInfo(I); I an ideal |
---|
255 | RETURN: ring |
---|
256 | STATUS: experimental, todo |
---|
257 | PURPOSE: compute the characteristic information for I |
---|
258 | ASSUME: the ground ring is the Weyl algebra with x's before d's |
---|
259 | NOTE: activate the output ring with the @code{setring} command. |
---|
260 | @* In the output (in a commutative ring): |
---|
261 | @* - the ideal CV is the characteristic variety char(I) |
---|
262 | @* - the ideal SL is the singular locus of char(I) |
---|
263 | @* - the list PD is the primary decomposition of char(I) |
---|
264 | @* If @code{printlevel}=1, progress debug messages will be printed, |
---|
265 | @* if @code{printlevel}>=2, all the debug messages will be printed. |
---|
266 | EXAMPLE: example annfs; shows examples |
---|
267 | " |
---|
268 | { |
---|
269 | def save = basering; |
---|
270 | def @A = charVariety(I); |
---|
271 | setring @A; |
---|
272 | // run slocus |
---|
273 | // run primdec |
---|
274 | } |
---|
275 | */ |
---|
276 | |
---|
277 | |
---|
278 | proc appelF1() |
---|
279 | "USAGE: appelF1(); |
---|
280 | RETURN: ring (and exports an ideal into it) |
---|
281 | PURPOSE: define the ideal in a parametric Weyl algebra, |
---|
282 | @* which annihilates Appel F1 hypergeometric function |
---|
283 | NOTE: the ideal called IAppel1 is exported to the output ring |
---|
284 | EXAMPLE: example appelF1; shows examples |
---|
285 | " |
---|
286 | { |
---|
287 | // Appel F1, d = b', SST p.48 |
---|
288 | ring @r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp); |
---|
289 | matrix @D[4][4]; |
---|
290 | @D[1,3]=1; @D[2,4]=1; |
---|
291 | def @S = nc_algebra(1,@D); |
---|
292 | setring @S; |
---|
293 | ideal IAppel1 = |
---|
294 | (x*Dx)*(x*Dx+y*Dy+c-1) - x*(x*Dx+y*Dy+a)*(x*Dx+b), |
---|
295 | (y*Dy)*(x*Dx+y*Dy+c-1) - y*(x*Dx+y*Dy+a)*(y*Dy+d), |
---|
296 | (x-y)*Dx*Dy - d*Dx + b*Dy; |
---|
297 | export IAppel1; |
---|
298 | kill @r; |
---|
299 | return(@S); |
---|
300 | } |
---|
301 | example |
---|
302 | { |
---|
303 | "EXAMPLE:"; echo = 2; |
---|
304 | def A = appelF1(); |
---|
305 | setring A; |
---|
306 | IAppel1; |
---|
307 | } |
---|
308 | |
---|
309 | proc appelF2() |
---|
310 | "USAGE: appelF2(); |
---|
311 | RETURN: ring (and exports an ideal into it) |
---|
312 | PURPOSE: define the ideal in a parametric Weyl algebra, |
---|
313 | @* which annihilates Appel F2 hypergeometric function |
---|
314 | NOTE: the ideal called IAppel2 is exported to the output ring |
---|
315 | EXAMPLE: example appelF2; shows examples |
---|
316 | " |
---|
317 | { |
---|
318 | // Appel F2, c = b', SST p.85 |
---|
319 | ring @r = (0,a,b,c),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp); |
---|
320 | matrix @D[4][4]; |
---|
321 | @D[1,3]=1; @D[2,4]=1; |
---|
322 | def @S = nc_algebra(1,@D); |
---|
323 | setring @S; |
---|
324 | ideal IAppel2 = |
---|
325 | (x*Dx)^2 - x*(x*Dx+y*Dy+a)*(x*Dx+b), |
---|
326 | (y*Dy)^2 - y*(x*Dx+y*Dy+a)*(y*Dy+c); |
---|
327 | export IAppel2; |
---|
328 | kill @r; |
---|
329 | return(@S); |
---|
330 | } |
---|
331 | example |
---|
332 | { |
---|
333 | "EXAMPLE:"; echo = 2; |
---|
334 | def A = appelF2(); |
---|
335 | setring A; |
---|
336 | IAppel2; |
---|
337 | } |
---|
338 | |
---|
339 | proc appelF4() |
---|
340 | "USAGE: appelF4(); |
---|
341 | RETURN: ring (and exports an ideal into it) |
---|
342 | PURPOSE: define the ideal in a parametric Weyl algebra, |
---|
343 | @* which annihilates Appel F4 hypergeometric function |
---|
344 | NOTE: the ideal called IAppel4 is exported to the output ring |
---|
345 | EXAMPLE: example appelF4; shows examples |
---|
346 | " |
---|
347 | { |
---|
348 | // Appel F4, d = c', SST, p. 39 |
---|
349 | ring @r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp); |
---|
350 | matrix @D[4][4]; |
---|
351 | @D[1,3]=1; @D[2,4]=1; |
---|
352 | def @S = nc_algebra(1,@D); |
---|
353 | setring @S; |
---|
354 | ideal IAppel4 = |
---|
355 | Dx*(x*Dx+c-1) - (x*Dx+y*Dy+a)*(x*Dx+y*Dy+b), |
---|
356 | Dy*(y*Dy+d-1) - (x*Dx+y*Dy+a)*(x*Dx+y*Dy+b); |
---|
357 | export IAppel4; |
---|
358 | kill @r; |
---|
359 | return(@S); |
---|
360 | } |
---|
361 | example |
---|
362 | { |
---|
363 | "EXAMPLE:"; echo = 2; |
---|
364 | def A = appelF4(); |
---|
365 | setring A; |
---|
366 | IAppel4; |
---|
367 | } |
---|
368 | |
---|
369 | proc poly2list (poly f) |
---|
370 | "USAGE: poly2list(f); f a poly |
---|
371 | RETURN: list of exponents and corresponding terms of f |
---|
372 | PURPOSE: convert a polynomial to a list of exponents and corresponding terms |
---|
373 | EXAMPLE: example poly2list; shows examples |
---|
374 | " |
---|
375 | { |
---|
376 | list l; |
---|
377 | int i = 1; |
---|
378 | if (f == 0) // just for the zero polynomial |
---|
379 | { |
---|
380 | l[1] = list(leadexp(f), lead(f)); |
---|
381 | } |
---|
382 | else { l[size(f)] = list(); } // memory pre-allocation |
---|
383 | while (f != 0) |
---|
384 | { |
---|
385 | l[i] = list(leadexp(f), lead(f)); |
---|
386 | f = f - lead(f); |
---|
387 | i++; |
---|
388 | } |
---|
389 | return(l); |
---|
390 | } |
---|
391 | example |
---|
392 | { |
---|
393 | "EXAMPLE:"; echo = 2; |
---|
394 | ring r = 0,x,dp; |
---|
395 | poly F = x; |
---|
396 | poly2list(F); |
---|
397 | ring r2 = 0,(x,y),dp; |
---|
398 | poly F = x2y+x*y2; |
---|
399 | poly2list(F); |
---|
400 | } |
---|
401 | |
---|
402 | proc isFsat(ideal I, poly F) |
---|
403 | "USAGE: isFsat(I, F); I an ideal, F a poly |
---|
404 | RETURN: int |
---|
405 | PURPOSE: check whether the ideal I is F-saturated |
---|
406 | NOTE: 1 is returned if I is F-saturated, otherwise 0 is returned. |
---|
407 | @* we check indeed that Ker(D --F--> D/I) is (0) |
---|
408 | EXAMPLE: example isFsat; shows examples |
---|
409 | " |
---|
410 | { |
---|
411 | /* checks whether I is F-saturated, that is Ke (D -F-> D/I) is 0 */ |
---|
412 | /* works in any algebra */ |
---|
413 | /* for simplicity : later check attrib */ |
---|
414 | /* returns -1 if true */ |
---|
415 | if (attrib(I,"isSB")!=1) |
---|
416 | { |
---|
417 | I = groebner(I); |
---|
418 | } |
---|
419 | matrix @M = matrix(I); |
---|
420 | matrix @F[1][1] = F; |
---|
421 | module S = modulo(@F,@M); |
---|
422 | S = NF(S,I); |
---|
423 | S = groebner(S); |
---|
424 | return( (gkdim(S) == -1) ); |
---|
425 | } |
---|
426 | example |
---|
427 | { |
---|
428 | "EXAMPLE:"; echo = 2; |
---|
429 | ring r = 0,(x,y),dp; |
---|
430 | poly G = x*(x-y)*y; |
---|
431 | def A = annfs(G); |
---|
432 | setring A; |
---|
433 | poly F = x3-y2; |
---|
434 | isFsat(LD,F); |
---|
435 | ideal J = LD*F; |
---|
436 | isFsat(J,F); |
---|
437 | } |
---|
438 | |
---|
439 | proc DLoc(ideal I, poly F) |
---|
440 | "USAGE: DLoc(I, F); I an ideal, F a poly |
---|
441 | RETURN: nothing (exports objects instead) |
---|
442 | ASSUME: the basering is a Weyl algebra |
---|
443 | PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s |
---|
444 | NOTE: In the basering, the following objects are exported: |
---|
445 | @* the ideal LD0 (in Groebner basis) is the presentation of the localization |
---|
446 | @* the list BS contains roots with multiplicities of Bernstein polynomial of (D/I)_f. |
---|
447 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
448 | @* if printlevel>=2, all the debug messages will be printed. |
---|
449 | EXAMPLE: example DLoc; shows examples |
---|
450 | " |
---|
451 | { |
---|
452 | /* runs SDLoc and DLoc0 */ |
---|
453 | /* assume: run from Weyl algebra */ |
---|
454 | if (dmodappassumeViolation()) |
---|
455 | { |
---|
456 | ERROR("Basering is inappropriate: characteristic>0 or qring present"); |
---|
457 | } |
---|
458 | if (!isWeyl()) |
---|
459 | { |
---|
460 | ERROR("Basering is not a Weyl algebra"); |
---|
461 | } |
---|
462 | if (defined(LD0) || defined(BS)) |
---|
463 | { |
---|
464 | ERROR("Reserved names LD0 and/or BS are used. Please rename the objects."); |
---|
465 | } |
---|
466 | int old_printlevel = printlevel; |
---|
467 | printlevel=printlevel+1; |
---|
468 | def @R = basering; |
---|
469 | def @R2 = SDLoc(I,F); |
---|
470 | setring @R2; |
---|
471 | poly F = imap(@R,F); |
---|
472 | def @R3 = DLoc0(LD,F); |
---|
473 | setring @R3; |
---|
474 | ideal bs = BS[1]; |
---|
475 | intvec m = BS[2]; |
---|
476 | setring @R; |
---|
477 | ideal LD0 = imap(@R3,LD0); |
---|
478 | export LD0; |
---|
479 | ideal bs = imap(@R3,bs); |
---|
480 | list BS; BS[1] = bs; BS[2] = m; |
---|
481 | export BS; |
---|
482 | kill @R3; |
---|
483 | printlevel = old_printlevel; |
---|
484 | } |
---|
485 | example; |
---|
486 | { |
---|
487 | "EXAMPLE:"; echo = 2; |
---|
488 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
489 | def R = Weyl(); setring R; // Weyl algebra in variables x,y,Dx,Dy |
---|
490 | poly F = x2-y3; |
---|
491 | ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F; |
---|
492 | // I is not holonomic, since its dimension is not 4/2=2 |
---|
493 | gkdim(I); |
---|
494 | DLoc(I, x2-y3); // exports LD0 and BS |
---|
495 | LD0; // localized module (R/I)_f is isomorphic to R/LD0 |
---|
496 | BS; // description of b-function for localization |
---|
497 | } |
---|
498 | |
---|
499 | proc DLoc0(ideal I, poly F) |
---|
500 | "USAGE: DLoc0(I, F); I an ideal, F a poly |
---|
501 | RETURN: ring |
---|
502 | PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s, |
---|
503 | @* where D is a Weyl Algebra, based on the output of procedure SDLoc |
---|
504 | ASSUME: the basering is similar to the output ring of SDLoc procedure |
---|
505 | NOTE: activate this ring with the @code{setring} command. In this ring, |
---|
506 | @* the ideal LD0 (in Groebner basis) is the presentation of the localization |
---|
507 | @* the list BS contains roots and multiplicities of Bernstein polynomial of (D/I)_f. |
---|
508 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
509 | @* if printlevel>=2, all the debug messages will be printed. |
---|
510 | EXAMPLE: example DLoc0; shows examples |
---|
511 | " |
---|
512 | { |
---|
513 | if (dmodappassumeViolation()) |
---|
514 | { |
---|
515 | ERROR("Basering is inappropriate: characteristic>0 or qring present"); |
---|
516 | } |
---|
517 | /* assume: to be run in the output ring of SDLoc */ |
---|
518 | /* doing: add F, eliminate vars*Dvars, factorize BS */ |
---|
519 | /* analogue to annfs0 */ |
---|
520 | def @R2 = basering; |
---|
521 | // we're in D_n[s], where the elim ord for s is set |
---|
522 | ideal J = NF(I,std(F)); |
---|
523 | // make leadcoeffs positive |
---|
524 | int i; |
---|
525 | for (i=1; i<= ncols(J); i++) |
---|
526 | { |
---|
527 | if (leadcoef(J[i]) <0 ) |
---|
528 | { |
---|
529 | J[i] = -J[i]; |
---|
530 | } |
---|
531 | } |
---|
532 | J = J,F; |
---|
533 | ideal M = groebner(J); |
---|
534 | int Nnew = nvars(@R2); |
---|
535 | ideal K2 = nselect(M,1..Nnew-1); |
---|
536 | int ppl = printlevel-voice+2; |
---|
537 | dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering"); |
---|
538 | dbprint(ppl-1, K2); |
---|
539 | // the ring @R3 and the search for minimal negative int s |
---|
540 | ring @R3 = 0,s,dp; |
---|
541 | dbprint(ppl,"// -2-1- the ring @R3 = K[s] is ready"); |
---|
542 | ideal K3 = imap(@R2,K2); |
---|
543 | poly p = K3[1]; |
---|
544 | dbprint(ppl,"// -2-2- attempt the factorization"); |
---|
545 | list PP = factorize(p); //with constants and multiplicities |
---|
546 | ideal bs; intvec m; //the Bernstein polynomial is monic, so we are not interested in constants |
---|
547 | for (i=2; i<= size(PP[1]); i++) //we delete P[1][1] and P[2][1] |
---|
548 | { |
---|
549 | bs[i-1] = PP[1][i]; |
---|
550 | m[i-1] = PP[2][i]; |
---|
551 | } |
---|
552 | ideal bbs; int srat=0; int HasRatRoots = 0; |
---|
553 | int sP; |
---|
554 | for (i=1; i<= size(bs); i++) |
---|
555 | { |
---|
556 | if (deg(bs[i]) == 1) |
---|
557 | { |
---|
558 | bbs = bbs,bs[i]; |
---|
559 | } |
---|
560 | } |
---|
561 | if (size(bbs)==0) |
---|
562 | { |
---|
563 | dbprint(ppl-1,"// -2-3- factorization: no rational roots"); |
---|
564 | // HasRatRoots = 0; |
---|
565 | HasRatRoots = 1; // s0 = -1 then |
---|
566 | sP = -1; |
---|
567 | // todo: return ideal with no subst and a b-function unfactorized |
---|
568 | } |
---|
569 | else |
---|
570 | { |
---|
571 | // exist rational roots |
---|
572 | dbprint(ppl-1,"// -2-3- factorization: rational roots found"); |
---|
573 | HasRatRoots = 1; |
---|
574 | // dbprint(ppl-1,bbs); |
---|
575 | bbs = bbs[2..ncols(bbs)]; |
---|
576 | ideal P = bbs; |
---|
577 | dbprint(ppl-1,P); |
---|
578 | srat = size(bs) - size(bbs); |
---|
579 | // define minIntRoot on linear factors or find out that it doesn't exist |
---|
580 | intvec vP; |
---|
581 | number nP; |
---|
582 | P = normalize(P); // now leadcoef = 1 |
---|
583 | P = lead(P)-P; |
---|
584 | sP = size(P); |
---|
585 | int cnt = 0; |
---|
586 | for (i=1; i<=sP; i++) |
---|
587 | { |
---|
588 | nP = leadcoef(P[i]); |
---|
589 | if ( (nP - int(nP)) == 0 ) |
---|
590 | { |
---|
591 | cnt++; |
---|
592 | vP[cnt] = int(nP); |
---|
593 | } |
---|
594 | } |
---|
595 | // if ( size(vP)>=2 ) |
---|
596 | // { |
---|
597 | // vP = vP[2..size(vP)]; |
---|
598 | // } |
---|
599 | if ( size(vP)==0 ) |
---|
600 | { |
---|
601 | // no roots! |
---|
602 | dbprint(ppl,"// -2-4- no integer root, setting s0 = -1"); |
---|
603 | sP = -1; |
---|
604 | // HasRatRoots = 0; // older stuff, here we do substitution |
---|
605 | HasRatRoots = 1; |
---|
606 | } |
---|
607 | else |
---|
608 | { |
---|
609 | HasRatRoots = 1; |
---|
610 | sP = -Max(-vP); |
---|
611 | dbprint(ppl,"// -2-4- minimal integer root found"); |
---|
612 | dbprint(ppl-1, sP); |
---|
613 | // int sP = minIntRoot(bbs,1); |
---|
614 | // P = normalize(P); |
---|
615 | // bs = -subst(bs,s,0); |
---|
616 | if (sP >=0) |
---|
617 | { |
---|
618 | dbprint(ppl,"// -2-5- nonnegative root, setting s0 = -1"); |
---|
619 | sP = -1; |
---|
620 | } |
---|
621 | else |
---|
622 | { |
---|
623 | dbprint(ppl,"// -2-5- the root is negative"); |
---|
624 | } |
---|
625 | } |
---|
626 | } |
---|
627 | |
---|
628 | if (HasRatRoots) |
---|
629 | { |
---|
630 | setring @R2; |
---|
631 | K2 = subst(I,s,sP); |
---|
632 | // IF min int root exists -> |
---|
633 | // create the ordinary Weyl algebra and put the result into it, |
---|
634 | // thus creating the ring @R5 |
---|
635 | // ELSE : return the same ring with new objects |
---|
636 | // keep: N, i,j,s, tmp, RL |
---|
637 | Nnew = Nnew - 1; // former 2*N; |
---|
638 | // list RL = ringlist(save); // is defined earlier |
---|
639 | // kill Lord, tmp, iv; |
---|
640 | list L = 0; |
---|
641 | list Lord, tmp; |
---|
642 | intvec iv; |
---|
643 | list RL = ringlist(basering); |
---|
644 | L[1] = RL[1]; |
---|
645 | L[4] = RL[4]; //char, minpoly |
---|
646 | // check whether vars have admissible names -> done earlier |
---|
647 | // list Name = RL[2]M |
---|
648 | // DName is defined earlier |
---|
649 | list NName; // = RL[2]; // skip the last var 's' |
---|
650 | for (i=1; i<=Nnew; i++) |
---|
651 | { |
---|
652 | NName[i] = RL[2][i]; |
---|
653 | } |
---|
654 | L[2] = NName; |
---|
655 | // dp ordering; |
---|
656 | string s = "iv="; |
---|
657 | for (i=1; i<=Nnew; i++) |
---|
658 | { |
---|
659 | s = s+"1,"; |
---|
660 | } |
---|
661 | s[size(s)] = ";"; |
---|
662 | execute(s); |
---|
663 | tmp = 0; |
---|
664 | tmp[1] = "dp"; // string |
---|
665 | tmp[2] = iv; // intvec |
---|
666 | Lord[1] = tmp; |
---|
667 | kill s; |
---|
668 | tmp[1] = "C"; |
---|
669 | iv = 0; |
---|
670 | tmp[2] = iv; |
---|
671 | Lord[2] = tmp; |
---|
672 | tmp = 0; |
---|
673 | L[3] = Lord; |
---|
674 | // we are done with the list |
---|
675 | // Add: Plural part |
---|
676 | def @R4@ = ring(L); |
---|
677 | setring @R4@; |
---|
678 | int N = Nnew/2; |
---|
679 | matrix @D[Nnew][Nnew]; |
---|
680 | for (i=1; i<=N; i++) |
---|
681 | { |
---|
682 | @D[i,N+i]=1; |
---|
683 | } |
---|
684 | def @R4 = nc_algebra(1,@D); |
---|
685 | setring @R4; |
---|
686 | kill @R4@; |
---|
687 | dbprint(ppl,"// -3-1- the ring @R4 is ready"); |
---|
688 | dbprint(ppl-1, @R4); |
---|
689 | ideal K4 = imap(@R2,K2); |
---|
690 | intvec vopt = option(get); |
---|
691 | option(redSB); |
---|
692 | dbprint(ppl,"// -3-2- the final cosmetic std"); |
---|
693 | K4 = groebner(K4); // std does the job too |
---|
694 | option(set,vopt); |
---|
695 | // total cleanup |
---|
696 | setring @R2; |
---|
697 | ideal bs = imap(@R3,bs); |
---|
698 | bs = -normalize(bs); // "-" for getting correct coeffs! |
---|
699 | bs = subst(bs,s,0); |
---|
700 | kill @R3; |
---|
701 | setring @R4; |
---|
702 | ideal bs = imap(@R2,bs); // only rationals are the entries |
---|
703 | list BS; BS[1] = bs; BS[2] = m; |
---|
704 | export BS; |
---|
705 | // list LBS = imap(@R3,LBS); |
---|
706 | // list BS; BS[1] = sbs; BS[2] = m; |
---|
707 | // BS; |
---|
708 | // export BS; |
---|
709 | ideal LD0 = K4; |
---|
710 | export LD0; |
---|
711 | return(@R4); |
---|
712 | } |
---|
713 | else |
---|
714 | { |
---|
715 | /* SHOULD NEVER GET THERE */ |
---|
716 | /* no rational/integer roots */ |
---|
717 | /* return objects in the copy of current ring */ |
---|
718 | setring @R2; |
---|
719 | ideal LD0 = I; |
---|
720 | poly BS = normalize(K2[1]); |
---|
721 | export LD0; |
---|
722 | export BS; |
---|
723 | return(@R2); |
---|
724 | } |
---|
725 | } |
---|
726 | example; |
---|
727 | { |
---|
728 | "EXAMPLE:"; echo = 2; |
---|
729 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
730 | def R = Weyl(); setring R; // Weyl algebra in variables x,y,Dx,Dy |
---|
731 | poly F = x2-y3; |
---|
732 | ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F; |
---|
733 | // moreover I is not holonomic, since its dimension is not 2 = 4/2 |
---|
734 | gkdim(I); // 3 |
---|
735 | def W = SDLoc(I,F); setring W; // creates ideal LD in W = R[s] |
---|
736 | def U = DLoc0(LD, x2-y3); setring U; // compute in R |
---|
737 | LD0; // Groebner basis of the presentation of localization |
---|
738 | BS; // description of b-function for localization |
---|
739 | } |
---|
740 | |
---|
741 | |
---|
742 | proc SDLoc(ideal I, poly F) |
---|
743 | "USAGE: SDLoc(I, F); I an ideal, F a poly |
---|
744 | RETURN: ring |
---|
745 | PURPOSE: compute a generic presentation of the localization of D/I w.r.t. f^s |
---|
746 | ASSUME: the basering D is a Weyl algebra |
---|
747 | NOTE: activate this ring with the @code{setring} command. In this ring, |
---|
748 | @* the ideal LD (in Groebner basis) is the presentation of the localization |
---|
749 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
750 | @* if printlevel>=2, all the debug messages will be printed. |
---|
751 | EXAMPLE: example SDLoc; shows examples |
---|
752 | " |
---|
753 | { |
---|
754 | /* analogue to Sannfs */ |
---|
755 | /* printlevel >=4 gives debug info */ |
---|
756 | /* assume: we're in the Weyl algebra D in x1,x2,...,d1,d2,... */ |
---|
757 | |
---|
758 | if (dmodappassumeViolation()) |
---|
759 | { |
---|
760 | ERROR("Basering is inappropriate: characteristic>0 or qring present"); |
---|
761 | } |
---|
762 | if (!isWeyl()) |
---|
763 | { |
---|
764 | ERROR("Basering is not a Weyl algebra"); |
---|
765 | } |
---|
766 | def save = basering; |
---|
767 | /* 1. create D <t, dt, s > as in LOT */ |
---|
768 | /* ordering: eliminate t,dt */ |
---|
769 | int ppl = printlevel-voice+2; |
---|
770 | int N = nvars(save); N = N div 2; |
---|
771 | int Nnew = 2*N + 3; // t,Dt,s |
---|
772 | int i,j; |
---|
773 | string s; |
---|
774 | list RL = ringlist(save); |
---|
775 | list L, Lord; |
---|
776 | list tmp; |
---|
777 | intvec iv; |
---|
778 | L[1] = RL[1]; // char |
---|
779 | L[4] = RL[4]; // char, minpoly |
---|
780 | // check whether vars have admissible names |
---|
781 | list Name = RL[2]; |
---|
782 | list RName; |
---|
783 | RName[1] = "@t"; |
---|
784 | RName[2] = "@Dt"; |
---|
785 | RName[3] = "@s"; |
---|
786 | for(i=1;i<=N;i++) |
---|
787 | { |
---|
788 | for(j=1; j<=size(RName);j++) |
---|
789 | { |
---|
790 | if (Name[i] == RName[j]) |
---|
791 | { |
---|
792 | ERROR("Variable names should not include @t,@Dt,@s"); |
---|
793 | } |
---|
794 | } |
---|
795 | } |
---|
796 | // now, create the names for new vars |
---|
797 | tmp = 0; |
---|
798 | tmp[1] = "@t"; |
---|
799 | tmp[2] = "@Dt"; |
---|
800 | list SName ; SName[1] = "@s"; |
---|
801 | list NName = tmp + Name + SName; |
---|
802 | L[2] = NName; |
---|
803 | tmp = 0; |
---|
804 | kill NName; |
---|
805 | // block ord (a(1,1),dp); |
---|
806 | tmp[1] = "a"; // string |
---|
807 | iv = 1,1; |
---|
808 | tmp[2] = iv; //intvec |
---|
809 | Lord[1] = tmp; |
---|
810 | // continue with dp 1,1,1,1... |
---|
811 | tmp[1] = "dp"; // string |
---|
812 | s = "iv="; |
---|
813 | for(i=1;i<=Nnew;i++) |
---|
814 | { |
---|
815 | s = s+"1,"; |
---|
816 | } |
---|
817 | s[size(s)]= ";"; |
---|
818 | execute(s); |
---|
819 | tmp[2] = iv; |
---|
820 | Lord[2] = tmp; |
---|
821 | tmp[1] = "C"; |
---|
822 | iv = 0; |
---|
823 | tmp[2] = iv; |
---|
824 | Lord[3] = tmp; |
---|
825 | tmp = 0; |
---|
826 | L[3] = Lord; |
---|
827 | // we are done with the list |
---|
828 | def @R@ = ring(L); |
---|
829 | setring @R@; |
---|
830 | matrix @D[Nnew][Nnew]; |
---|
831 | @D[1,2]=1; |
---|
832 | for(i=1; i<=N; i++) |
---|
833 | { |
---|
834 | @D[2+i,N+2+i]=1; |
---|
835 | } |
---|
836 | // ADD [s,t]=-t, [s,Dt]=Dt |
---|
837 | @D[1,Nnew] = -var(1); |
---|
838 | @D[2,Nnew] = var(2); |
---|
839 | def @R = nc_algebra(1,@D); |
---|
840 | setring @R; |
---|
841 | kill @R@; |
---|
842 | dbprint(ppl,"// -1-1- the ring @R(@t,@Dt,_x,_Dx,@s) is ready"); |
---|
843 | dbprint(ppl-1, @R); |
---|
844 | poly F = imap(save,F); |
---|
845 | ideal I = imap(save,I); |
---|
846 | dbprint(ppl-1, "the ideal after map:"); |
---|
847 | dbprint(ppl-1, I); |
---|
848 | poly p = 0; |
---|
849 | for(i=1; i<=N; i++) |
---|
850 | { |
---|
851 | p = diff(F,var(2+i))*@Dt + var(2+N+i); |
---|
852 | dbprint(ppl-1, p); |
---|
853 | I = subst(I,var(2+N+i),p); |
---|
854 | dbprint(ppl-1, var(2+N+i)); |
---|
855 | p = 0; |
---|
856 | } |
---|
857 | I = I, @t - F; |
---|
858 | // t*Dt + s +1 reduced with t-f gives f*Dt + s |
---|
859 | I = I, F*var(2) + var(Nnew); // @s |
---|
860 | // -------- the ideal I is ready ---------- |
---|
861 | dbprint(ppl,"// -1-2- starting the elimination of @t,@Dt in @R"); |
---|
862 | dbprint(ppl-1, I); |
---|
863 | // ideal J = engine(I,eng); |
---|
864 | ideal J = groebner(I); |
---|
865 | dbprint(ppl-1,"// -1-2-1- result of the elimination of @t,@Dt in @R"); |
---|
866 | dbprint(ppl-1, J);; |
---|
867 | ideal K = nselect(J,1..2); |
---|
868 | dbprint(ppl,"// -1-3- @t,@Dt are eliminated"); |
---|
869 | dbprint(ppl-1, K); // K is without t, Dt |
---|
870 | K = groebner(K); // std does the job too |
---|
871 | // now, we must change the ordering |
---|
872 | // and create a ring without t, Dt |
---|
873 | setring save; |
---|
874 | // ----------- the ring @R3 ------------ |
---|
875 | // _x, _Dx,s; elim.ord for _x,_Dx. |
---|
876 | // keep: N, i,j,s, tmp, RL |
---|
877 | Nnew = 2*N+1; |
---|
878 | kill Lord, tmp, iv, RName; |
---|
879 | list Lord, tmp; |
---|
880 | intvec iv; |
---|
881 | L[1] = RL[1]; |
---|
882 | L[4] = RL[4]; // char, minpoly |
---|
883 | // check whether vars hava admissible names -> done earlier |
---|
884 | // now, create the names for new var |
---|
885 | tmp[1] = "s"; |
---|
886 | list NName = Name + tmp; |
---|
887 | L[2] = NName; |
---|
888 | tmp = 0; |
---|
889 | // block ord (dp(N),dp); |
---|
890 | // string s is already defined |
---|
891 | s = "iv="; |
---|
892 | for (i=1; i<=Nnew-1; i++) |
---|
893 | { |
---|
894 | s = s+"1,"; |
---|
895 | } |
---|
896 | s[size(s)]=";"; |
---|
897 | execute(s); |
---|
898 | tmp[1] = "dp"; // string |
---|
899 | tmp[2] = iv; // intvec |
---|
900 | Lord[1] = tmp; |
---|
901 | // continue with dp 1,1,1,1... |
---|
902 | tmp[1] = "dp"; // string |
---|
903 | s[size(s)] = ","; |
---|
904 | s = s+"1;"; |
---|
905 | execute(s); |
---|
906 | kill s; |
---|
907 | kill NName; |
---|
908 | tmp[2] = iv; |
---|
909 | Lord[2] = tmp; |
---|
910 | tmp[1] = "C"; iv = 0; tmp[2]=iv; |
---|
911 | Lord[3] = tmp; tmp = 0; |
---|
912 | L[3] = Lord; |
---|
913 | // we are done with the list. Now add a Plural part |
---|
914 | def @R2@ = ring(L); |
---|
915 | setring @R2@; |
---|
916 | matrix @D[Nnew][Nnew]; |
---|
917 | for (i=1; i<=N; i++) |
---|
918 | { |
---|
919 | @D[i,N+i]=1; |
---|
920 | } |
---|
921 | def @R2 = nc_algebra(1,@D); |
---|
922 | setring @R2; |
---|
923 | kill @R2@; |
---|
924 | dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); |
---|
925 | dbprint(ppl-1, @R2); |
---|
926 | ideal MM = maxideal(1); |
---|
927 | MM = 0,s,MM; |
---|
928 | map R01 = @R, MM; |
---|
929 | ideal K = R01(K); |
---|
930 | // total cleanup |
---|
931 | ideal LD = K; |
---|
932 | // make leadcoeffs positive |
---|
933 | for (i=1; i<= ncols(LD); i++) |
---|
934 | { |
---|
935 | if (leadcoef(LD[i]) <0 ) |
---|
936 | { |
---|
937 | LD[i] = -LD[i]; |
---|
938 | } |
---|
939 | } |
---|
940 | export LD; |
---|
941 | kill @R; |
---|
942 | return(@R2); |
---|
943 | } |
---|
944 | example; |
---|
945 | { |
---|
946 | "EXAMPLE:"; echo = 2; |
---|
947 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
948 | def R = Weyl(); // Weyl algebra on the variables x,y,Dx,Dy |
---|
949 | setring R; |
---|
950 | poly F = x2-y3; |
---|
951 | ideal I = Dx*F, Dy*F; |
---|
952 | // note, that I is not holonomic, since it's dimension is not 2 |
---|
953 | gkdim(I); // 3, while dim R = 4 |
---|
954 | def W = SDLoc(I,F); |
---|
955 | setring W; // = R[s], where s is a new variable |
---|
956 | LD; // Groebner basis of s-parametric presentation |
---|
957 | } |
---|
958 | |
---|
959 | proc annRat(poly g, poly f) |
---|
960 | "USAGE: annRat(g,f); f, g polynomials |
---|
961 | RETURN: ring |
---|
962 | PURPOSE: compute the annihilator of the rational function g/f in the Weyl algebra D |
---|
963 | NOTE: activate the output ring with the @code{setring} command. |
---|
964 | @* In the output ring, the ideal LD (in Groebner basis) is the annihilator. |
---|
965 | @* The algorithm uses the computation of ann f^{-1} via D-modules. |
---|
966 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
967 | @* if printlevel>=2, all the debug messages will be printed. |
---|
968 | SEE ALSO: annPoly |
---|
969 | EXAMPLE: example annRat; shows examples |
---|
970 | " |
---|
971 | { |
---|
972 | |
---|
973 | if (dmodappassumeViolation()) |
---|
974 | { |
---|
975 | ERROR("Basering is inappropriate: characteristic>0 or qring present"); |
---|
976 | } |
---|
977 | |
---|
978 | // assumptions: f is not a constant |
---|
979 | if (f==0) { ERROR("Denominator cannot be zero"); } |
---|
980 | if (leadexp(f) == 0) |
---|
981 | { |
---|
982 | // f = const, so use annPoly |
---|
983 | g = g/f; |
---|
984 | def @R = annPoly(g); |
---|
985 | return(@R); |
---|
986 | } |
---|
987 | // computes the annihilator of g/f |
---|
988 | def save = basering; |
---|
989 | int ppl = printlevel-voice+2; |
---|
990 | dbprint(ppl,"// -1-[annRat] computing the ann f^s"); |
---|
991 | def @R1 = SannfsBM(f); |
---|
992 | setring @R1; |
---|
993 | poly f = imap(save,f); |
---|
994 | int i,mir; |
---|
995 | int isr = 0; // checkRoot1(LD,f,1); // roots are negative, have to enter positive int |
---|
996 | if (!isr) |
---|
997 | { |
---|
998 | // -1 is not the root |
---|
999 | // find the m.i.r iteratively |
---|
1000 | mir = 0; |
---|
1001 | for(i=nvars(save)+1; i>=1; i--) |
---|
1002 | { |
---|
1003 | isr = checkRoot1(LD,f,i); |
---|
1004 | if (isr) { mir =-i; break; } |
---|
1005 | } |
---|
1006 | if (mir ==0) |
---|
1007 | { |
---|
1008 | "No integer root found! Aborting computations, inform the authors!"; |
---|
1009 | return(0); |
---|
1010 | } |
---|
1011 | // now mir == i is m.i.r. |
---|
1012 | } |
---|
1013 | else |
---|
1014 | { |
---|
1015 | // -1 is the m.i.r |
---|
1016 | mir = -1; |
---|
1017 | } |
---|
1018 | dbprint(ppl,"// -2-[annRat] the minimal integer root is "); |
---|
1019 | dbprint(ppl-1, mir); |
---|
1020 | // use annfspecial |
---|
1021 | dbprint(ppl,"// -3-[annRat] running annfspecial "); |
---|
1022 | ideal AF = annfspecial(LD,f,mir,-1); // ann f^{-1} |
---|
1023 | // LD = subst(LD,s,j); |
---|
1024 | // LD = engine(LD,0); |
---|
1025 | // modify the ring: throw s away |
---|
1026 | // output ring comes from SannfsBM |
---|
1027 | list U = ringlist(@R1); |
---|
1028 | list tmp; // variables |
---|
1029 | for(i=1; i<=size(U[2])-1; i++) |
---|
1030 | { |
---|
1031 | tmp[i] = U[2][i]; |
---|
1032 | } |
---|
1033 | U[2] = tmp; |
---|
1034 | tmp = 0; |
---|
1035 | tmp[1] = U[3][1]; // x,Dx block |
---|
1036 | tmp[2] = U[3][3]; // module block |
---|
1037 | U[3] = tmp; |
---|
1038 | tmp = 0; |
---|
1039 | tmp = U[1],U[2],U[3],U[4]; |
---|
1040 | def @R2 = ring(tmp); |
---|
1041 | setring @R2; |
---|
1042 | // now supply with Weyl algebra relations |
---|
1043 | int N = nvars(@R2)/2; |
---|
1044 | matrix @D[2*N][2*N]; |
---|
1045 | for(i=1; i<=N; i++) |
---|
1046 | { |
---|
1047 | @D[i,N+i]=1; |
---|
1048 | } |
---|
1049 | def @R3 = nc_algebra(1,@D); |
---|
1050 | setring @R3; |
---|
1051 | dbprint(ppl,"// - -[annRat] ring without s is ready:"); |
---|
1052 | dbprint(ppl-1,@R3); |
---|
1053 | poly g = imap(save,g); |
---|
1054 | matrix G[1][1] = g; |
---|
1055 | matrix LL = matrix(imap(@R1,AF)); |
---|
1056 | kill @R1; kill @R2; |
---|
1057 | dbprint(ppl,"// -4-[annRat] running modulo"); |
---|
1058 | ideal LD = modulo(G,LL); |
---|
1059 | dbprint(ppl,"// -4-[annRat] running GB on the final result"); |
---|
1060 | LD = engine(LD,0); |
---|
1061 | export LD; |
---|
1062 | return(@R3); |
---|
1063 | } |
---|
1064 | example |
---|
1065 | { |
---|
1066 | "EXAMPLE:"; echo = 2; |
---|
1067 | ring r = 0,(x,y),dp; |
---|
1068 | poly g = 2*x*y; poly f = x^2 - y^3; |
---|
1069 | def B = annRat(g,f); |
---|
1070 | setring B; |
---|
1071 | LD; |
---|
1072 | // Now, compare with the output of Macaulay2: |
---|
1073 | ideal tst = 3*x*Dx + 2*y*Dy + 1, y^3*Dy^2 - x^2*Dy^2 + 6*y^2*Dy + 6*y, |
---|
1074 | 9*y^2*Dx^2*Dy-4*y*Dy^3+27*y*Dx^2+2*Dy^2, 9*y^3*Dx^2-4*y^2*Dy^2+10*y*Dy -10; |
---|
1075 | option(redSB); option(redTail); |
---|
1076 | LD = groebner(LD); |
---|
1077 | tst = groebner(tst); |
---|
1078 | print(matrix(NF(LD,tst))); print(matrix(NF(tst,LD))); |
---|
1079 | // So, these two answers are the same |
---|
1080 | } |
---|
1081 | |
---|
1082 | /* |
---|
1083 | //static proc ex_annRat() |
---|
1084 | { |
---|
1085 | // more complicated example for annRat |
---|
1086 | ring r = 0,(x,y,z),dp; |
---|
1087 | poly f = x3+y3+z3; // mir = -2 |
---|
1088 | poly g = x*y*z; |
---|
1089 | def A = annRat(g,f); |
---|
1090 | setring A; |
---|
1091 | } |
---|
1092 | */ |
---|
1093 | |
---|
1094 | proc annPoly(poly f) |
---|
1095 | "USAGE: annPoly(f); f a poly |
---|
1096 | RETURN: ring |
---|
1097 | PURPOSE: compute the complete annihilator ideal of f in the Weyl algebra D |
---|
1098 | NOTE: activate the output ring with the @code{setring} command. |
---|
1099 | @* In the output ring, the ideal LD (in Groebner basis) is the annihilator. |
---|
1100 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
1101 | @* if printlevel>=2, all the debug messages will be printed. |
---|
1102 | SEE ALSO: annRat |
---|
1103 | EXAMPLE: example annPoly; shows examples |
---|
1104 | " |
---|
1105 | { |
---|
1106 | // computes a system of linear PDEs with polynomial coeffs for f |
---|
1107 | def save = basering; |
---|
1108 | list L = ringlist(save); |
---|
1109 | list Name = L[2]; |
---|
1110 | int N = nvars(save); |
---|
1111 | int i; |
---|
1112 | for (i=1; i<=N; i++) |
---|
1113 | { |
---|
1114 | Name[N+i] = "D"+Name[i]; // concat |
---|
1115 | } |
---|
1116 | L[2] = Name; |
---|
1117 | def @R = ring(L); |
---|
1118 | setring @R; |
---|
1119 | def @@R = Weyl(); |
---|
1120 | setring @@R; |
---|
1121 | kill @R; |
---|
1122 | matrix M[1][N]; |
---|
1123 | for (i=1; i<=N; i++) |
---|
1124 | { |
---|
1125 | M[1,i] = var(N+i); |
---|
1126 | } |
---|
1127 | matrix F[1][1] = imap(save,f); |
---|
1128 | ideal I = modulo(F,M); |
---|
1129 | ideal LD = groebner(I); |
---|
1130 | export LD; |
---|
1131 | return(@@R); |
---|
1132 | } |
---|
1133 | example |
---|
1134 | { |
---|
1135 | "EXAMPLE:"; echo = 2; |
---|
1136 | ring r = 0,(x,y,z),dp; |
---|
1137 | poly f = x^2*z - y^3; |
---|
1138 | def A = annPoly(f); |
---|
1139 | setring A; // A is the 3rd Weyl algebra in 6 variables |
---|
1140 | LD; // the Groebner basis of annihilator |
---|
1141 | gkdim(LD); // must be 3 = 6/2, since A/LD is holonomic module |
---|
1142 | NF(Dy^4, LD); // must be 0 since Dy^4 clearly annihilates f |
---|
1143 | } |
---|
1144 | |
---|
1145 | /* DIFFERENT EXAMPLES |
---|
1146 | |
---|
1147 | //static proc exCusp() |
---|
1148 | { |
---|
1149 | "EXAMPLE:"; echo = 2; |
---|
1150 | ring r = 0,(x,y,Dx,Dy),dp; |
---|
1151 | def R = Weyl(); setring R; |
---|
1152 | poly F = x2-y3; |
---|
1153 | ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F; |
---|
1154 | def W = SDLoc(I,F); |
---|
1155 | setring W; |
---|
1156 | LD; |
---|
1157 | def U = DLoc0(LD,x2-y3); |
---|
1158 | setring U; |
---|
1159 | LD0; |
---|
1160 | BS; |
---|
1161 | // the same with DLoc: |
---|
1162 | setring R; |
---|
1163 | DLoc(I,F); |
---|
1164 | } |
---|
1165 | |
---|
1166 | //static proc exWalther1() |
---|
1167 | { |
---|
1168 | // p.18 Rem 3.10 |
---|
1169 | ring r = 0,(x,Dx),dp; |
---|
1170 | def R = nc_algebra(1,1); |
---|
1171 | setring R; |
---|
1172 | poly F = x; |
---|
1173 | ideal I = x*Dx+1; |
---|
1174 | def W = SDLoc(I,F); |
---|
1175 | setring W; |
---|
1176 | LD; |
---|
1177 | ideal J = LD, x; |
---|
1178 | eliminate(J,x*Dx); // must be [1]=s // agree! |
---|
1179 | // the same result with Dloc0: |
---|
1180 | def U = DLoc0(LD,x); |
---|
1181 | setring U; |
---|
1182 | LD0; |
---|
1183 | BS; |
---|
1184 | } |
---|
1185 | |
---|
1186 | //static proc exWalther2() |
---|
1187 | { |
---|
1188 | // p.19 Rem 3.10 cont'd |
---|
1189 | ring r = 0,(x,Dx),dp; |
---|
1190 | def R = nc_algebra(1,1); |
---|
1191 | setring R; |
---|
1192 | poly F = x; |
---|
1193 | ideal I = (x*Dx)^2+1; |
---|
1194 | def W = SDLoc(I,F); |
---|
1195 | setring W; |
---|
1196 | LD; |
---|
1197 | ideal J = LD, x; |
---|
1198 | eliminate(J,x*Dx); // must be [1]=s^2+2*s+2 // agree! |
---|
1199 | // the same result with Dloc0: |
---|
1200 | def U = DLoc0(LD,x); |
---|
1201 | setring U; |
---|
1202 | LD0; |
---|
1203 | BS; |
---|
1204 | // almost the same with DLoc |
---|
1205 | setring R; |
---|
1206 | DLoc(I,F); |
---|
1207 | LD0; BS; |
---|
1208 | } |
---|
1209 | |
---|
1210 | //static proc exWalther3() |
---|
1211 | { |
---|
1212 | // can check with annFs too :-) |
---|
1213 | // p.21 Ex 3.15 |
---|
1214 | LIB "nctools.lib"; |
---|
1215 | ring r = 0,(x,y,z,w,Dx,Dy,Dz,Dw),dp; |
---|
1216 | def R = Weyl(); |
---|
1217 | setring R; |
---|
1218 | poly F = x2+y2+z2+w2; |
---|
1219 | ideal I = Dx,Dy,Dz,Dw; |
---|
1220 | def W = SDLoc(I,F); |
---|
1221 | setring W; |
---|
1222 | LD; |
---|
1223 | ideal J = LD, x2+y2+z2+w2; |
---|
1224 | eliminate(J,x*y*z*w*Dx*Dy*Dz*Dw); // must be [1]=s^2+3*s+2 // agree |
---|
1225 | ring r2 = 0,(x,y,z,w),dp; |
---|
1226 | poly F = x2+y2+z2+w2; |
---|
1227 | def Z = annfs(F); |
---|
1228 | setring Z; |
---|
1229 | LD; |
---|
1230 | BS; |
---|
1231 | // the same result with Dloc0: |
---|
1232 | setring W; |
---|
1233 | def U = DLoc0(LD,x2+y2+z2+w2); |
---|
1234 | setring U; |
---|
1235 | LD0; BS; |
---|
1236 | // the same result with DLoc: |
---|
1237 | setring R; |
---|
1238 | DLoc(I,F); |
---|
1239 | LD0; BS; |
---|
1240 | } |
---|
1241 | |
---|
1242 | */ |
---|
1243 | |
---|
1244 | proc engine(def I, int i) |
---|
1245 | "USAGE: engine(I,i); I ideal/module/matrix, i an int |
---|
1246 | RETURN: the same type as I |
---|
1247 | PURPOSE: compute the Groebner basis of I with the algorithm, chosen via i |
---|
1248 | NOTE: By default and if i=0, slimgb is used; otherwise std does the job. |
---|
1249 | EXAMPLE: example engine; shows examples |
---|
1250 | " |
---|
1251 | { |
---|
1252 | /* std - slimgb mix */ |
---|
1253 | def J; |
---|
1254 | // ideal J; |
---|
1255 | if (i==0) |
---|
1256 | { |
---|
1257 | J = slimgb(I); |
---|
1258 | } |
---|
1259 | else |
---|
1260 | { |
---|
1261 | // without options -> strange! (ringlist?) |
---|
1262 | intvec v = option(get); |
---|
1263 | option(redSB); |
---|
1264 | option(redTail); |
---|
1265 | J = std(I); |
---|
1266 | option(set, v); |
---|
1267 | } |
---|
1268 | return(J); |
---|
1269 | } |
---|
1270 | example |
---|
1271 | { |
---|
1272 | "EXAMPLE:"; echo = 2; |
---|
1273 | ring r = 0,(x,y),Dp; |
---|
1274 | ideal I = y*(x3-y2),x*(x3-y2); |
---|
1275 | engine(I,0); // uses slimgb |
---|
1276 | engine(I,1); // uses std |
---|
1277 | } |
---|
1278 | |
---|
1279 | proc insertGenerator (list #) |
---|
1280 | "USAGE: insertGenerator(id,p[,k]); id an ideal/module, p a poly/vector, k an optional int |
---|
1281 | RETURN: same as id |
---|
1282 | PURPOSE: inserts p into the first argument at k-th index position and returns the enlarged object |
---|
1283 | NOTE: If k is given, p is inserted at position k, otherwise (and by default), |
---|
1284 | @* p is inserted at the beginning. |
---|
1285 | EXAMPLE: example insertGenerator; shows examples |
---|
1286 | " |
---|
1287 | { |
---|
1288 | if (size(#) < 2) |
---|
1289 | { |
---|
1290 | ERROR("insertGenerator has to be called with at least 2 arguments (ideal/module,poly/vector)"); |
---|
1291 | } |
---|
1292 | string inp1 = typeof(#[1]); |
---|
1293 | if (inp1 == "ideal" || inp1 == "module") |
---|
1294 | { |
---|
1295 | if (inp1 == "ideal") { ideal id = #[1]; } |
---|
1296 | else { module id = #[1]; } |
---|
1297 | } |
---|
1298 | else { ERROR("first argument has to be of type ideal or module"); } |
---|
1299 | string inp2 = typeof(#[2]); |
---|
1300 | if (inp2 == "poly" || inp2 == "vector") |
---|
1301 | { |
---|
1302 | if (inp2 =="poly") { poly f = #[2]; } |
---|
1303 | else |
---|
1304 | { |
---|
1305 | if (inp1 == "ideal") |
---|
1306 | { |
---|
1307 | ERROR("second argument has to be a polynomial if first argument is an ideal"); |
---|
1308 | } |
---|
1309 | else { vector f = #[2]; } |
---|
1310 | } |
---|
1311 | } |
---|
1312 | else { ERROR("second argument has to be of type poly/vector"); } |
---|
1313 | int n = ncols(id); |
---|
1314 | int k = 1; // default |
---|
1315 | if (size(#)>=3) |
---|
1316 | { |
---|
1317 | if (typeof(#[3]) == "int") |
---|
1318 | { |
---|
1319 | k = #[3]; |
---|
1320 | if (k<=0) |
---|
1321 | { |
---|
1322 | ERROR("third argument has to be positive"); |
---|
1323 | } |
---|
1324 | } |
---|
1325 | else { ERROR("third argument has to be of type int"); } |
---|
1326 | } |
---|
1327 | execute(inp1 +" J;"); |
---|
1328 | if (k == 1) { J = f,id; } |
---|
1329 | else |
---|
1330 | { |
---|
1331 | if (k>n) |
---|
1332 | { |
---|
1333 | J = id; |
---|
1334 | J[k] = f; |
---|
1335 | } |
---|
1336 | else // 1<k<=n |
---|
1337 | { |
---|
1338 | J[1..k-1] = id[1..k-1]; |
---|
1339 | J[k] = f; |
---|
1340 | J[k+1..n+1] = id[k..n]; |
---|
1341 | } |
---|
1342 | } |
---|
1343 | return(J); |
---|
1344 | } |
---|
1345 | example |
---|
1346 | { |
---|
1347 | "EXAMPLE:"; echo = 2; |
---|
1348 | ring r = 0,(x,y,z),dp; |
---|
1349 | ideal I = x^2,z^4; |
---|
1350 | insertGenerator(I,y^3); |
---|
1351 | insertGenerator(I,y^3,2); |
---|
1352 | module M = I; |
---|
1353 | insertGenerator(M,[x^3,y^2,z],2); |
---|
1354 | } |
---|
1355 | |
---|
1356 | proc deleteGenerator (list #) |
---|
1357 | "USAGE: deleteGenerator(id,k); id an ideal/module, k an int |
---|
1358 | RETURN: same as id |
---|
1359 | PURPOSE: deletes the k-th generator from the first argument and returns the altered object |
---|
1360 | EXAMPLE: example insertGenerator; shows examples |
---|
1361 | " |
---|
1362 | { |
---|
1363 | if (size(#) < 2) |
---|
1364 | { |
---|
1365 | ERROR("deleteGenerator has to be called with 2 arguments (ideal/module,int)"); |
---|
1366 | } |
---|
1367 | string inp1 = typeof(#[1]); |
---|
1368 | if (inp1 == "ideal" || inp1 == "module") |
---|
1369 | { |
---|
1370 | if (inp1 == "ideal") { ideal id = #[1]; } |
---|
1371 | else { module id = #[1]; } |
---|
1372 | } |
---|
1373 | else { ERROR("first argument has to be of type ideal or module"); } |
---|
1374 | string inp2 = typeof(#[2]); |
---|
1375 | if (inp2 == "int" || inp2 == "number") { int k = int(#[2]); } |
---|
1376 | else { ERROR("second argument has to be of type int"); } |
---|
1377 | int n = ncols(id); |
---|
1378 | if (n == 1) { ERROR(inp1+" must have more than one generator"); } |
---|
1379 | if (k<=0 || k>n) { ERROR("second argument has to be in the range 1,...,"+string(n)); } |
---|
1380 | execute(inp1 +" J;"); |
---|
1381 | if (k == 1) { J = id[2..n]; } |
---|
1382 | else |
---|
1383 | { |
---|
1384 | if (k == n) { J = id[1..n-1]; } |
---|
1385 | else |
---|
1386 | { |
---|
1387 | J[1..k-1] = id[1..k-1]; |
---|
1388 | J[k..n-1] = id[k+1..n]; |
---|
1389 | } |
---|
1390 | } |
---|
1391 | return(J); |
---|
1392 | } |
---|
1393 | example |
---|
1394 | { |
---|
1395 | "EXAMPLE:"; echo = 2; |
---|
1396 | ring r = 0,(x,y,z),dp; |
---|
1397 | ideal I = x^2,y^3,z^4; |
---|
1398 | deleteGenerator(I,2); |
---|
1399 | module M = [x,y,z],[x2,y2,z2],[x3,y3,z3];; |
---|
1400 | deleteGenerator(M,2); |
---|
1401 | } |
---|
1402 | |
---|
1403 | proc fl2poly(list L, string s) |
---|
1404 | "USAGE: fl2poly(L,s); L a list, s a string |
---|
1405 | RETURN: poly |
---|
1406 | PURPOSE: reconstruct a monic polynomial in one variable from its factorization |
---|
1407 | ASSUME: s is a string with the name of some variable and |
---|
1408 | @* L is supposed to consist of two entries: |
---|
1409 | @* L[1] of the type ideal with the roots of a polynomial |
---|
1410 | @* L[2] of the type intvec with the multiplicities of corr. roots |
---|
1411 | EXAMPLE: example fl2poly; shows examples |
---|
1412 | " |
---|
1413 | { |
---|
1414 | if (varNum(s)==0) |
---|
1415 | { |
---|
1416 | ERROR("no such variable found in the basering"); return(0); |
---|
1417 | } |
---|
1418 | poly x = var(varNum(s)); |
---|
1419 | poly P = 1; |
---|
1420 | int sl = size(L[1]); |
---|
1421 | ideal RR = L[1]; |
---|
1422 | intvec IV = L[2]; |
---|
1423 | for(int i=1; i<= sl; i++) |
---|
1424 | { |
---|
1425 | if (IV[i] > 0) |
---|
1426 | { |
---|
1427 | P = P*((x-RR[i])^IV[i]); |
---|
1428 | } |
---|
1429 | else |
---|
1430 | { |
---|
1431 | printf("Ignored the root with incorrect multiplicity %s",string(IV[i])); |
---|
1432 | } |
---|
1433 | } |
---|
1434 | return(P); |
---|
1435 | } |
---|
1436 | example |
---|
1437 | { |
---|
1438 | "EXAMPLE:"; echo = 2; |
---|
1439 | ring r = 0,(x,y,z,s),Dp; |
---|
1440 | ideal I = -1,-4/3,-5/3,-2; |
---|
1441 | intvec mI = 2,1,1,1; |
---|
1442 | list BS = I,mI; |
---|
1443 | poly p = fl2poly(BS,"s"); |
---|
1444 | p; |
---|
1445 | factorize(p,2); |
---|
1446 | } |
---|
1447 | |
---|
1448 | static proc safeVarName (string s, string cv) |
---|
1449 | { |
---|
1450 | string S; |
---|
1451 | if (cv == "v") { S = "," + "," + varstr(basering) + ","; } |
---|
1452 | if (cv == "c") { S = "," + "," + charstr(basering) + ","; } |
---|
1453 | if (cv == "cv") { S = "," + charstr(basering) + "," + varstr(basering) + ","; } |
---|
1454 | s = "," + s + ","; |
---|
1455 | while (find(S,s) <> 0) |
---|
1456 | { |
---|
1457 | s[1] = "@"; |
---|
1458 | s = "," + s; |
---|
1459 | } |
---|
1460 | s = s[2..size(s)-1]; |
---|
1461 | return(s) |
---|
1462 | } |
---|
1463 | |
---|
1464 | proc initialIdealW (ideal I, intvec u, intvec v, list #) |
---|
1465 | "USAGE: initialIdealW(I,u,v [,s,t,w]); I ideal, u,v intvecs, s,t optional ints, |
---|
1466 | @* w an optional intvec |
---|
1467 | RETURN: ideal, GB of initial ideal of the input ideal w.r.t. the weights u and v |
---|
1468 | ASSUME: The basering is the n-th Weyl algebra in characteristic 0 and for all |
---|
1469 | @* 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the |
---|
1470 | @* sequence of variables is given by x(1),...,x(n),D(1),...,D(n), |
---|
1471 | @* where D(i) is the differential operator belonging to x(i). |
---|
1472 | PURPOSE: computes the initial ideal with respect to given weights. |
---|
1473 | NOTE: u and v are understood as weight vectors for x(1..n) and D(1..n) |
---|
1474 | @* respectively. |
---|
1475 | @* If s<>0, @code{std} is used for Groebner basis computations, |
---|
1476 | @* otherwise, and by default, @code{slimgb} is used. |
---|
1477 | @* If t<>0, a matrix ordering is used for Groebner basis computations, |
---|
1478 | @* otherwise, and by default, a block ordering is used. |
---|
1479 | @* If w consist of 2n strictly positive entries, w is used for weighted |
---|
1480 | @* homogenization, otherwise, and by default, no weights are used. |
---|
1481 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
1482 | @* if printlevel>=2, all the debug messages will be printed. |
---|
1483 | EXAMPLE: example initialIdealW; shows examples |
---|
1484 | " |
---|
1485 | { |
---|
1486 | if (dmodappassumeViolation()) |
---|
1487 | { |
---|
1488 | ERROR("Basering is inappropriate: characteristic>0 or qring present"); |
---|
1489 | } |
---|
1490 | if (!isWeyl()) |
---|
1491 | { |
---|
1492 | ERROR("Basering is not a Weyl algebra."); |
---|
1493 | } |
---|
1494 | int ppl = printlevel - voice +2; |
---|
1495 | def save = basering; |
---|
1496 | int n = nvars(save)/2; |
---|
1497 | int N = 2*n+1; |
---|
1498 | list RL = ringlist(save); |
---|
1499 | int i; |
---|
1500 | int whichengine = 0; // default |
---|
1501 | int methodord = 0; // default |
---|
1502 | intvec homogweights = 1:(2*n); // default |
---|
1503 | if (size(#)>0) |
---|
1504 | { |
---|
1505 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
1506 | { |
---|
1507 | whichengine = int(#[1]); |
---|
1508 | } |
---|
1509 | if (size(#)>1) |
---|
1510 | { |
---|
1511 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
1512 | { |
---|
1513 | methodord = int(#[2]); |
---|
1514 | } |
---|
1515 | if (size(#)>2) |
---|
1516 | { |
---|
1517 | if (typeof(#[3])=="intvec") |
---|
1518 | { |
---|
1519 | if (size(#[3])==2*n && allPositive(#[3])==1) |
---|
1520 | { |
---|
1521 | homogweights = #[3]; |
---|
1522 | } |
---|
1523 | else |
---|
1524 | { |
---|
1525 | print("// Homogenization weight vector must consist of positive entries and be"); |
---|
1526 | print("// of size " + string(n) + ". Using no weights."); |
---|
1527 | } |
---|
1528 | } |
---|
1529 | } |
---|
1530 | } |
---|
1531 | } |
---|
1532 | for (i=1; i<=n; i++) |
---|
1533 | { |
---|
1534 | if (bracket(var(i+n),var(i))<>1) |
---|
1535 | { |
---|
1536 | ERROR(string(var(i+n)) + " is not a differential operator for " + string(var(i))); |
---|
1537 | } |
---|
1538 | } |
---|
1539 | // 1. create homogenized Weyl algebra |
---|
1540 | // 1.1 create ordering |
---|
1541 | intvec uv = u,v,0; |
---|
1542 | homogweights = homogweights,1; |
---|
1543 | list Lord = list(list("a",homogweights)); |
---|
1544 | list C0 = list("C",intvec(0)); |
---|
1545 | if (methodord == 0) // default: blockordering |
---|
1546 | { |
---|
1547 | Lord[2] = list("a",uv); |
---|
1548 | Lord[3] = list("dp",intvec(1:(N-1))); |
---|
1549 | Lord[4] = list("lp",intvec(1)); |
---|
1550 | Lord[5] = C0; |
---|
1551 | } |
---|
1552 | else // M() ordering |
---|
1553 | { |
---|
1554 | intmat @Ord[N][N]; |
---|
1555 | @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1); |
---|
1556 | for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; } |
---|
1557 | dbprint(ppl-1,"// the ordering matrix:",@Ord); |
---|
1558 | Lord[2] = list("M",intvec(@Ord)); |
---|
1559 | Lord[3] = C0; |
---|
1560 | } |
---|
1561 | // 1.2 the new var |
---|
1562 | list Lvar = RL[2]; Lvar[N] = safeVarName("h","cv"); |
---|
1563 | // 1.3 create commutative ring |
---|
1564 | list L@@Dh = RL; L@@Dh = L@@Dh[1..4]; |
---|
1565 | L@@Dh[2] = Lvar; L@@Dh[3] = Lord; |
---|
1566 | def @Dh = ring(L@@Dh); kill L@@Dh; |
---|
1567 | setring @Dh; |
---|
1568 | dbprint(ppl-1,"// the ring @Dh:",@Dh); |
---|
1569 | // 1.4 create non-commutative relations |
---|
1570 | matrix @relD[N][N]; |
---|
1571 | for (i=1; i<=n; i++) { @relD[i,n+i] = var(N)^(homogweights[i]+homogweights[n+i]); } |
---|
1572 | def Dh = nc_algebra(1,@relD); |
---|
1573 | setring Dh; kill @Dh; |
---|
1574 | dbprint(ppl-1,"// computing in ring",Dh); |
---|
1575 | // 2. Compute the initial ideal |
---|
1576 | ideal I = imap(save,I); |
---|
1577 | I = homog(I,h); |
---|
1578 | // 2.1 the hard part: Groebner basis computation |
---|
1579 | dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine)); |
---|
1580 | I = engine(I, whichengine); |
---|
1581 | dbprint(ppl, "// finished Groebner basis computation"); |
---|
1582 | dbprint(ppl-1, "// I before dehomogenization is " +string(I)); |
---|
1583 | I = subst(I,var(N),1); // dehomogenization |
---|
1584 | dbprint(ppl-1, "I after dehomogenization is " +string(I)); |
---|
1585 | // 2.2 the initial form |
---|
1586 | I = inForm(I,uv); |
---|
1587 | setring save; |
---|
1588 | I = imap(Dh,I); kill Dh; |
---|
1589 | // 2.3 the final GB |
---|
1590 | dbprint(ppl, "// starting cosmetic Groebner basis computation with engine: "+string(whichengine)); |
---|
1591 | I = engine(I, whichengine); |
---|
1592 | dbprint(ppl,"// finished cosmetic Groebner basis computation"); |
---|
1593 | return(I); |
---|
1594 | } |
---|
1595 | example |
---|
1596 | { |
---|
1597 | "EXAMPLE:"; echo = 2; |
---|
1598 | ring @D = 0,(x,Dx),dp; |
---|
1599 | def D = Weyl(); |
---|
1600 | setring D; |
---|
1601 | intvec u = -1; intvec v = 2; |
---|
1602 | ideal I = x^2*Dx^2,x*Dx^4; |
---|
1603 | ideal J = initialIdealW(I,u,v); J; |
---|
1604 | } |
---|
1605 | |
---|
1606 | proc initialMalgrange (poly f,list #) |
---|
1607 | "USAGE: initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec |
---|
1608 | RETURN: ring, Weyl algebra induced by basering, extended by two new vars t,Dt |
---|
1609 | PURPOSE: computes the initial Malgrange ideal of a given polynomial w.r.t. the weight |
---|
1610 | @* vector (-1,0...,0,1,0,...,0) such that the weight of t is -1 and the |
---|
1611 | @* weight of Dt is 1. |
---|
1612 | ASSUME: The basering is commutative and of characteristic 0. |
---|
1613 | NOTE: Activate the output ring with the @code{setring} command. |
---|
1614 | @* The returned ring contains the ideal \"inF\", being the initial ideal |
---|
1615 | @* of the Malgrange ideal of f. |
---|
1616 | @* Varnames of the basering should not include t and Dt. |
---|
1617 | @* If a<>0, @code{std} is used for Groebner basis computations, |
---|
1618 | @* otherwise, and by default, @code{slimgb} is used. |
---|
1619 | @* If b<>0, a matrix ordering is used for Groebner basis computations, |
---|
1620 | @* otherwise, and by default, a block ordering is used. |
---|
1621 | @* If a positive weight vector v is given, the weight |
---|
1622 | @* (d,v[1],...,v[n],1,d+1-v[1],...,d+1-v[n]) is used for homogenization |
---|
1623 | @* computations, where d denotes the weighted degree of f. |
---|
1624 | @* Otherwise and by default, v is set to (1,...,1). See Noro, 2002. |
---|
1625 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
1626 | @* if printlevel>=2, all the debug messages will be printed. |
---|
1627 | EXAMPLE: example initialMalgrange; shows examples |
---|
1628 | " |
---|
1629 | { |
---|
1630 | if (dmodappassumeViolation()) |
---|
1631 | { |
---|
1632 | ERROR("Basering is inappropriate: characteristic>0 or qring present"); |
---|
1633 | } |
---|
1634 | if (!isCommutative()) |
---|
1635 | { |
---|
1636 | ERROR("Basering must be commutative."); |
---|
1637 | } |
---|
1638 | int ppl = printlevel - voice +2; |
---|
1639 | def save = basering; |
---|
1640 | int n = nvars(save); |
---|
1641 | int i; |
---|
1642 | int whichengine = 0; // default |
---|
1643 | int methodord = 0; // default |
---|
1644 | intvec u0 = 1:n; // default |
---|
1645 | if (size(#)>0) |
---|
1646 | { |
---|
1647 | if (typeof(#[1])=="int" || typeof(#[1])=="number") |
---|
1648 | { |
---|
1649 | whichengine = int(#[1]); |
---|
1650 | } |
---|
1651 | if (size(#)>1) |
---|
1652 | { |
---|
1653 | if (typeof(#[2])=="int" || typeof(#[2])=="number") |
---|
1654 | { |
---|
1655 | methodord = int(#[2]); |
---|
1656 | } |
---|
1657 | if (size(#)>2) |
---|
1658 | { |
---|
1659 | if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1) |
---|
1660 | { |
---|
1661 | u0 = #[3]; |
---|
1662 | } |
---|
1663 | } |
---|
1664 | } |
---|
1665 | } |
---|
1666 | // creating the homogenized extended Weyl algebra |
---|
1667 | list RL = ringlist(save); |
---|
1668 | RL = RL[1..4]; // if basering is commutative nc_algebra |
---|
1669 | list C0 = list("C",intvec(0)); |
---|
1670 | // 1. create homogenization weights |
---|
1671 | // 1.1. get the weighted degree of f |
---|
1672 | list Lord = list(list("wp",u0),C0); |
---|
1673 | list L@r = RL; |
---|
1674 | L@r[3] = Lord; |
---|
1675 | def r = ring(L@r); kill L@r,Lord; |
---|
1676 | setring r; |
---|
1677 | poly f = imap(save,f); |
---|
1678 | int d = deg(f); |
---|
1679 | setring save; kill r; |
---|
1680 | // 1.2 the homogenization weights |
---|
1681 | intvec homogweights = d; |
---|
1682 | homogweights[n+2] = 1; |
---|
1683 | for (i=1; i<=n; i++) |
---|
1684 | { |
---|
1685 | homogweights[i+1] = u0[i]; |
---|
1686 | homogweights[n+2+i] = d+1-u0[i]; |
---|
1687 | } |
---|
1688 | // 2. create extended Weyl algebra |
---|
1689 | int N = 2*n+2; |
---|
1690 | // 2.1 create names for vars |
---|
1691 | string vart = safeVarName("t","cv"); |
---|
1692 | string varDt = safeVarName("D"+vart,"cv"); |
---|
1693 | while (varDt <> "D"+vart) |
---|
1694 | { |
---|
1695 | vart = safeVarName("@"+vart,"cv"); |
---|
1696 | varDt = safeVarName("D"+vart,"cv"); |
---|
1697 | } |
---|
1698 | list Lvar; |
---|
1699 | Lvar[1] = vart; Lvar[n+2] = varDt; |
---|
1700 | for (i=1; i<=n; i++) |
---|
1701 | { |
---|
1702 | Lvar[i+1] = string(var(i)); |
---|
1703 | Lvar[i+n+2] = safeVarName("D" + string(var(i)),"cv"); |
---|
1704 | } |
---|
1705 | // 2.2 create ordering |
---|
1706 | list Lord = list(list("dp",intvec(1:N)),C0); |
---|
1707 | // 2.3 create the (n+1)-th Weyl algebra |
---|
1708 | list L@D = RL; L@D[2] = Lvar; L@D[3] = Lord; |
---|
1709 | def @D = ring(L@D); kill L@D; |
---|
1710 | setring @D; |
---|
1711 | def D = Weyl(); |
---|
1712 | setring D; kill @D; |
---|
1713 | dbprint(ppl,"// the (n+1)-th Weyl algebra :" ,D); |
---|
1714 | // 3. compute the initial ideal |
---|
1715 | // 3.1 create the Malgrange ideal |
---|
1716 | poly f = imap(save,f); |
---|
1717 | ideal I = var(1)-f; |
---|
1718 | for (i=1; i<=n; i++) |
---|
1719 | { |
---|
1720 | I = I, var(n+2+i)+diff(f,var(i+1))*var(n+2); |
---|
1721 | } |
---|
1722 | // I = engine(I,whichengine); // todo is it efficient to compute a GB first wrt dp first? |
---|
1723 | // 3.2 computie the initial ideal |
---|
1724 | intvec w = 1,0:n; |
---|
1725 | I = initialIdealW(I,-w,w,whichengine,methodord,homogweights); |
---|
1726 | ideal inF = I; attrib(inF,"isSB",1); |
---|
1727 | export(inF); |
---|
1728 | setring save; |
---|
1729 | return(D); |
---|
1730 | } |
---|
1731 | example |
---|
1732 | { |
---|
1733 | "EXAMPLE:"; echo = 2; |
---|
1734 | ring r = 0,(x,y),dp; |
---|
1735 | poly f = x^2+y^3+x*y^2; |
---|
1736 | def D = initialMalgrange(f); |
---|
1737 | setring D; |
---|
1738 | inF; |
---|
1739 | setring r; |
---|
1740 | intvec v = 3,2; |
---|
1741 | def D2 = initialMalgrange(f,1,1,v); |
---|
1742 | setring D2; |
---|
1743 | inF; |
---|
1744 | } |
---|
1745 | |
---|
1746 | proc dmodappassumeViolation() |
---|
1747 | { |
---|
1748 | // returns Boolean : yes/no [for assume violation] |
---|
1749 | // char K = 0 |
---|
1750 | // no qring |
---|
1751 | if ( (size(ideal(basering)) >0) || (char(basering) >0) ) |
---|
1752 | { |
---|
1753 | // "ERROR: no qring is allowed"; |
---|
1754 | return(1); |
---|
1755 | } |
---|
1756 | return(0); |
---|
1757 | } |
---|
1758 | |
---|
1759 | proc bFactor (poly F) |
---|
1760 | "USAGE: bFactor(f); f poly |
---|
1761 | RETURN: list |
---|
1762 | PURPOSE: tries to compute the roots of a univariate poly f |
---|
1763 | NOTE: The output list consists of two or three entries: |
---|
1764 | @* roots of f as an ideal, their multiplicities as intvec, and, |
---|
1765 | @* if present, a third one being the product of all irreducible factors |
---|
1766 | @* of degree greater than one, given as string. |
---|
1767 | DISPLAY: If printlevel=1, progress debug messages will be printed, |
---|
1768 | @* if printlevel>=2, all the debug messages will be printed. |
---|
1769 | EXAMPLE: example bFactor; shows examples |
---|
1770 | " |
---|
1771 | { |
---|
1772 | int ppl = printlevel - voice +2; |
---|
1773 | def save = basering; |
---|
1774 | ideal LI = variables(F); |
---|
1775 | list L; |
---|
1776 | int i = size(LI); |
---|
1777 | if (i>1) { ERROR("poly has to be univariate")} |
---|
1778 | if (i == 0) |
---|
1779 | { |
---|
1780 | dbprint(ppl,"// poly is constant"); |
---|
1781 | L = list(ideal(0),intvec(0),string(F)); |
---|
1782 | return(L); |
---|
1783 | } |
---|
1784 | poly v = LI[1]; |
---|
1785 | L = ringlist(save); L = L[1..4]; |
---|
1786 | L[2] = list(string(v)); |
---|
1787 | L[3] = list(list("dp",intvec(1)),list("C",intvec(0))); |
---|
1788 | def @S = ring(L); |
---|
1789 | setring @S; |
---|
1790 | poly ir = 1; poly F = imap(save,F); |
---|
1791 | list L = factorize(F); |
---|
1792 | ideal I = L[1]; |
---|
1793 | intvec m = L[2]; |
---|
1794 | ideal II; intvec mm; |
---|
1795 | for (i=2; i<=ncols(I); i++) |
---|
1796 | { |
---|
1797 | if (deg(I[i]) > 1) |
---|
1798 | { |
---|
1799 | ir = ir * I[i]^m[i]; |
---|
1800 | } |
---|
1801 | else |
---|
1802 | { |
---|
1803 | II[size(II)+1] = I[i]; |
---|
1804 | mm[size(II)] = m[i]; |
---|
1805 | } |
---|
1806 | } |
---|
1807 | II = normalize(II); |
---|
1808 | II = subst(II,var(1),0); |
---|
1809 | II = -II; |
---|
1810 | if (size(II)>0) |
---|
1811 | { |
---|
1812 | dbprint(ppl,"// found roots"); |
---|
1813 | dbprint(ppl-1,string(II)); |
---|
1814 | } |
---|
1815 | else |
---|
1816 | { |
---|
1817 | dbprint(ppl,"// found no roots"); |
---|
1818 | } |
---|
1819 | L = list(II,mm); |
---|
1820 | if (ir <> 1) |
---|
1821 | { |
---|
1822 | dbprint(ppl,"// found irreducible factors"); |
---|
1823 | ir = cleardenom(ir); |
---|
1824 | dbprint(ppl-1,string(ir)); |
---|
1825 | L = L + list(string(ir)); |
---|
1826 | } |
---|
1827 | else |
---|
1828 | { |
---|
1829 | dbprint(ppl,"// no irreducible factors found"); |
---|
1830 | } |
---|
1831 | setring save; |
---|
1832 | L = imap(@S,L); |
---|
1833 | return(L); |
---|
1834 | } |
---|
1835 | example |
---|
1836 | { |
---|
1837 | "EXAMPLE:"; echo = 2; |
---|
1838 | ring r = 0,(x,y),dp; |
---|
1839 | bFactor((x^2-1)^2); |
---|
1840 | bFactor((x^2+1)^2); |
---|
1841 | bFactor((y^2+1/2)*(y+9)*(y-7)); |
---|
1842 | } |
---|