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