1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id: dmod.lib,v 1.3 2006-05-08 14:43:21 levandov Exp $"; |
---|
3 | category="Noncommutative"; |
---|
4 | info=" |
---|
5 | LIBRARY: dmod.lib Algorithms for algebraic D-modules |
---|
6 | AUTHORS: Viktor Levandovskyy, levandov@mathematik.uni-kl.de |
---|
7 | @* Jorge Martin Morales, jorge@unizar.es |
---|
8 | |
---|
9 | THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R, one is interested in the ring R[1/F^s] for a natural number s. |
---|
10 | @* In fact, the ring R[1/F^s] has a structure of a D(R)-module, where D(R) is a Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>. |
---|
11 | @* Constructively, one needs to find a left ideal I = I(F^s) in D(R), such |
---|
12 | that K[x_1,...,x_n,1/F^s] is isomorphic to D(R)/I as a D(R)-module. |
---|
13 | @* We provide two implementations: |
---|
14 | @* 1) the classical Ann F^s algorithm from Oaku and Takayama (J. Pure Applied Math., 1999) and |
---|
15 | @* 2) the newer Ann F^s algorithm by Briancon and Maisonobe. |
---|
16 | |
---|
17 | PROCEDURES: |
---|
18 | annfsOT(F[,eng]); compute Ann F^s for a poly F with the algorithm of Oaku-Takayama |
---|
19 | annfsBM(F[,eng]); compute Ann F^s for a poly F with the algorithm of Briancon-Maisonobe |
---|
20 | reiffen(p,q); create the polynomial, describing a Reiffen curve |
---|
21 | arrange(p); create the polynomial, describing a generic hyperplane arrangement |
---|
22 | isHolonomic(M); check whether a module is holonomic |
---|
23 | convloc(L); replace global orderings with local in the ringlist L |
---|
24 | "; |
---|
25 | // very experimental: |
---|
26 | // annfsgms(F[,eng]); compute Ann F^s for a poly F with the modified Oaku-Takayama |
---|
27 | // a new algorithm, based on the computation of a Bernstein polynomial |
---|
28 | // first (with the help of algorithms by Mathias Schulze, see gmssing.lib) and |
---|
29 | // consequent eliminations like in Oaku and Takayama. |
---|
30 | |
---|
31 | LIB "nctools.lib"; |
---|
32 | LIB "elim.lib"; |
---|
33 | LIB "qhmoduli.lib"; // for Max |
---|
34 | LIB "gkdim.lib"; |
---|
35 | |
---|
36 | proc engine(ideal I, int i) |
---|
37 | { |
---|
38 | /* std / slimgb mix */ |
---|
39 | ideal J; |
---|
40 | if (i==0) |
---|
41 | { |
---|
42 | J = slimgb(I); |
---|
43 | } |
---|
44 | else |
---|
45 | { |
---|
46 | // without options -> buggy!!! (ringlist?) |
---|
47 | option(redSB); |
---|
48 | option(redTail); |
---|
49 | J = std(I); |
---|
50 | } |
---|
51 | return(J); |
---|
52 | } |
---|
53 | |
---|
54 | proc annfsBM(poly F, list #) |
---|
55 | "USAGE: annfsBM(f [,eng]); f a poly, eng an optional int |
---|
56 | RETURN: ring |
---|
57 | PURPOSE: compute the D-module structure of basering[f^s], according |
---|
58 | to the algorithm by Briancon and Maisonobe |
---|
59 | NOTE: activate this ring with the @code{setring} command. In this ring, |
---|
60 | @* - the ideal LD is the needed D-mod structure, |
---|
61 | @* - the ideal BS is the list of roots of a Bernstein polynomial of f. |
---|
62 | @* If eng <>0, @code{std} is used for Groebner basis computations, |
---|
63 | @* otherwise (and by default) @code{slimgb} is used. |
---|
64 | @* If printlevel=1, progress debug messages will be printed, |
---|
65 | @* if printlevel>=2, all the debug messages will be printed. |
---|
66 | EXAMPLE: example annfsBM; shows examples |
---|
67 | " |
---|
68 | { |
---|
69 | int eng = 0; |
---|
70 | if ( size(#)>0 ) |
---|
71 | { |
---|
72 | if ( typeof(#[1]) == "int" ) |
---|
73 | { |
---|
74 | eng = int(#[1]); |
---|
75 | } |
---|
76 | } |
---|
77 | // returns a list with a ring and an ideal LD in it |
---|
78 | int ppl = printlevel-voice+2; |
---|
79 | // printf("plevel :%s, voice: %s",printlevel,voice); |
---|
80 | def save = basering; |
---|
81 | int N = nvars(basering); |
---|
82 | int Nnew = 2*N+2; |
---|
83 | int i,j; |
---|
84 | string s; |
---|
85 | list RL = ringlist(basering); |
---|
86 | list L, Lord; |
---|
87 | list tmp; |
---|
88 | intvec iv; |
---|
89 | L[1] = RL[1]; // char |
---|
90 | L[4] = RL[4]; // char, minpoly |
---|
91 | // check whether vars have admissible names |
---|
92 | list Name = RL[2]; |
---|
93 | list RName; |
---|
94 | RName[1] = "t"; |
---|
95 | RName[2] = "s"; |
---|
96 | for(i=1;i<=N;i++) |
---|
97 | { |
---|
98 | for(j=1; j<=size(RName);j++) |
---|
99 | { |
---|
100 | if (Name[i] == RName[j]) |
---|
101 | { |
---|
102 | "Variable names should not include t,s"; |
---|
103 | return(0); |
---|
104 | } |
---|
105 | } |
---|
106 | } |
---|
107 | // now, create the names for new vars |
---|
108 | list DName; |
---|
109 | for(i=1;i<=N;i++) |
---|
110 | { |
---|
111 | DName[i] = "D"+Name[i]; // concat |
---|
112 | } |
---|
113 | tmp[1] = "t"; |
---|
114 | tmp[2] = "s"; |
---|
115 | list NName = tmp + Name + DName; |
---|
116 | L[2] = NName; |
---|
117 | // Name, Dname will be used further |
---|
118 | kill NName; |
---|
119 | // block ord (lp(2),dp); |
---|
120 | tmp[1] = "lp"; // string |
---|
121 | iv = 1,1; |
---|
122 | tmp[2] = iv; //intvec |
---|
123 | Lord[1] = tmp; |
---|
124 | // continue with dp 1,1,1,1... |
---|
125 | tmp[1] = "dp"; // string |
---|
126 | s = "iv="; |
---|
127 | for(i=1;i<=Nnew;i++) |
---|
128 | { |
---|
129 | s = s+"1,"; |
---|
130 | } |
---|
131 | s[size(s)]= ";"; |
---|
132 | execute(s); |
---|
133 | kill s; |
---|
134 | tmp[2] = iv; |
---|
135 | Lord[2] = tmp; |
---|
136 | tmp[1] = "C"; |
---|
137 | iv = 0; |
---|
138 | tmp[2] = iv; |
---|
139 | Lord[3] = tmp; |
---|
140 | tmp = 0; |
---|
141 | L[3] = Lord; |
---|
142 | // we are done with the list |
---|
143 | def @R = ring(L); |
---|
144 | setring @R; |
---|
145 | matrix @D[Nnew][Nnew]; |
---|
146 | @D[1,2]=t; |
---|
147 | for(i=1; i<=N; i++) |
---|
148 | { |
---|
149 | @D[2+i,N+2+i]=1; |
---|
150 | } |
---|
151 | // L[5] = matrix(UpOneMatrix(Nnew)); |
---|
152 | // L[6] = @D; |
---|
153 | ncalgebra(1,@D); |
---|
154 | dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready"); |
---|
155 | dbprint(ppl-1, @R); |
---|
156 | // create the ideal I |
---|
157 | poly F = imap(save,F); |
---|
158 | ideal I = t*F+s; |
---|
159 | poly p; |
---|
160 | for(i=1; i<=N; i++) |
---|
161 | { |
---|
162 | p = t; // t |
---|
163 | p = diff(F,var(2+i))*p; |
---|
164 | I = I, var(N+2+i) + p; |
---|
165 | } |
---|
166 | // -------- the ideal I is ready ---------- |
---|
167 | dbprint(ppl,"// -1-2- starting the elimination of t in @R"); |
---|
168 | dbprint(ppl-1, I); |
---|
169 | ideal J = engine(I,eng); |
---|
170 | ideal K = nselect(J,1); |
---|
171 | dbprint(ppl,"// -1-3- t is eliminated"); |
---|
172 | dbprint(ppl-1, K); // K is without t |
---|
173 | setring save; |
---|
174 | // ----------- the ring @R3 ------------ |
---|
175 | // _x, _Dx,s; elim.ord for _x,_Dx. |
---|
176 | // keep: N, i,j,s, tmp, RL |
---|
177 | Nnew = 2*N+1; |
---|
178 | kill Lord, tmp, iv, RName; |
---|
179 | list Lord, tmp; |
---|
180 | intvec iv; |
---|
181 | L[1] = RL[1]; |
---|
182 | L[4] = RL[4]; // char, minpoly |
---|
183 | // check whether vars hava admissible names -> done earlier |
---|
184 | // now, create the names for new var |
---|
185 | tmp[1] = "s"; |
---|
186 | // DName is defined earlier |
---|
187 | list NName = Name + DName + tmp; |
---|
188 | L[2] = NName; |
---|
189 | tmp = 0; |
---|
190 | // block ord (dp(N),dp); |
---|
191 | string s = "iv="; |
---|
192 | for (i=1; i<=Nnew-1; i++) |
---|
193 | { |
---|
194 | s = s+"1,"; |
---|
195 | } |
---|
196 | s[size(s)]=";"; |
---|
197 | execute(s); |
---|
198 | tmp[1] = "dp"; // string |
---|
199 | tmp[2] = iv; // intvec |
---|
200 | Lord[1] = tmp; |
---|
201 | // continue with dp 1,1,1,1... |
---|
202 | tmp[1] = "dp"; // string |
---|
203 | s[size(s)] = ","; |
---|
204 | s = s+"1;"; |
---|
205 | execute(s); |
---|
206 | kill s; |
---|
207 | kill NName; |
---|
208 | tmp[2] = iv; |
---|
209 | Lord[2] = tmp; |
---|
210 | tmp[1] = "C"; iv = 0; tmp[2]=iv; |
---|
211 | Lord[3] = tmp; tmp = 0; |
---|
212 | L[3] = Lord; |
---|
213 | // we are done with the list. Now add a Plural part |
---|
214 | def @R2 = ring(L); |
---|
215 | setring @R2; |
---|
216 | matrix @D[Nnew][Nnew]; |
---|
217 | for (i=1; i<=N; i++) |
---|
218 | { |
---|
219 | @D[i,N+i]=1; |
---|
220 | } |
---|
221 | ncalgebra(1,@D); |
---|
222 | dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); |
---|
223 | dbprint(ppl-1, @R2); |
---|
224 | ideal MM = maxideal(1); |
---|
225 | MM = 0,s,MM; |
---|
226 | map R01 = @R, MM; |
---|
227 | ideal K = R01(K); |
---|
228 | poly F = imap(save,F); |
---|
229 | K = K,F; |
---|
230 | dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2"); |
---|
231 | dbprint(ppl-1, K); |
---|
232 | ideal M = engine(K,eng); |
---|
233 | ideal K2 = nselect(M,1,Nnew-1); |
---|
234 | dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2"); |
---|
235 | dbprint(ppl-1, K2); |
---|
236 | // the ring @R3 and the search for minimal negative int s |
---|
237 | ring @R3 = 0,s,dp; |
---|
238 | dbprint(ppl,"// -3-1- the ring @R3 is ready"); |
---|
239 | ideal K3 = imap(@R2,K2); |
---|
240 | poly p = K3[1]; |
---|
241 | dbprint(ppl,"// -3-2- factorization"); |
---|
242 | ideal P = factorize(p,1); //without constants and multiplicities |
---|
243 | // "--------- b-function factorizes into ---------"; P; |
---|
244 | int sP = minIntRoot(P,1); |
---|
245 | dbprint(ppl,"// -3-3- minimal interger root found"); |
---|
246 | dbprint(ppl-1, sP); |
---|
247 | // convert factors to a list of their roots |
---|
248 | // assume all factors are linear |
---|
249 | ideal BS = normalize(P); |
---|
250 | BS = subst(BS,s,0); |
---|
251 | BS = -BS; |
---|
252 | //TODO: sort BS! |
---|
253 | // --------- substitute s found in the ideal --------- |
---|
254 | // --------- going back to @R and substitute --------- |
---|
255 | setring @R; |
---|
256 | ideal K2 = subst(K,s,sP); |
---|
257 | // create the ordinary Weyl algebra and put the result into it, |
---|
258 | // thus creating the ring @R5 |
---|
259 | // keep: N, i,j,s, tmp, RL |
---|
260 | setring save; |
---|
261 | Nnew = 2*N; |
---|
262 | // list RL = ringlist(save); // is defined earlier |
---|
263 | kill Lord, tmp, iv; |
---|
264 | L = 0; |
---|
265 | list Lord, tmp; |
---|
266 | intvec iv; |
---|
267 | L[1] = RL[1]; |
---|
268 | L[4] = RL[4]; //char, minpoly |
---|
269 | // check whether vars have admissible names -> done earlier |
---|
270 | // list Name = RL[2]M |
---|
271 | // DName is defined earlier |
---|
272 | list NName = Name + DName; |
---|
273 | L[2] = NName; |
---|
274 | // dp ordering; |
---|
275 | string s = "iv="; |
---|
276 | for (i=1; i<=Nnew; i++) |
---|
277 | { |
---|
278 | s = s+"1,"; |
---|
279 | } |
---|
280 | s[size(s)] = ";"; |
---|
281 | execute(s); |
---|
282 | tmp = 0; |
---|
283 | tmp[1] = "dp"; // string |
---|
284 | tmp[2] = iv; // intvec |
---|
285 | Lord[1] = tmp; |
---|
286 | kill s; |
---|
287 | tmp[1] = "C"; |
---|
288 | iv = 0; |
---|
289 | tmp[2] = iv; |
---|
290 | Lord[2] = tmp; |
---|
291 | tmp = 0; |
---|
292 | L[3] = Lord; |
---|
293 | // we are done with the list |
---|
294 | // Add: Plural part |
---|
295 | def @R4 = ring(L); |
---|
296 | setring @R4; |
---|
297 | matrix @D[Nnew][Nnew]; |
---|
298 | for (i=1; i<=N; i++) |
---|
299 | { |
---|
300 | @D[i,N+i]=1; |
---|
301 | } |
---|
302 | ncalgebra(1,@D); |
---|
303 | dbprint(ppl,"// -4-1- the ring @R4 is ready"); |
---|
304 | dbprint(ppl-1, @R4); |
---|
305 | ideal K4 = imap(@R,K2); |
---|
306 | option(redSB); |
---|
307 | dbprint(ppl,"// -4-2- the final cosmetic std"); |
---|
308 | K4 = engine(K4,eng); // std does the job too |
---|
309 | // total cleanup |
---|
310 | kill @R; |
---|
311 | kill @R2; |
---|
312 | ideal BS = imap(@R3,BS); |
---|
313 | export BS; |
---|
314 | kill @R3; |
---|
315 | ideal LD = K4; |
---|
316 | export LD; |
---|
317 | return(@R4); |
---|
318 | } |
---|
319 | example |
---|
320 | { |
---|
321 | "EXAMPLE:"; echo = 2; |
---|
322 | ring r = 0,(x,y,z),Dp; |
---|
323 | poly F = x^2+y^3+z^5; |
---|
324 | printlevel = 0; |
---|
325 | def A = annfsBM(F); |
---|
326 | setring A; |
---|
327 | LD; |
---|
328 | print(matrix(BS)); |
---|
329 | } |
---|
330 | |
---|
331 | proc annfsOT(poly F, list #) |
---|
332 | "USAGE: annfsOT(f [,eng]); f a poly, eng an optional int |
---|
333 | RETURN: ring |
---|
334 | PURPOSE: compute the D-module structure of basering[f^s], according |
---|
335 | to the algorithm by Oaku and Takayama |
---|
336 | NOTE: activate this ring with the @code{setring} command. In this ring, |
---|
337 | @* - the ideal LD is the needed D-mod structure, |
---|
338 | @* - the ideal BS is the list of roots of a Bernstein polynomial of f. |
---|
339 | @* If eng <>0, @code{std} is used for Groebner basis computations, |
---|
340 | @* otherwise (and by default) @code{slimgb} is used. |
---|
341 | @* If printlevel=1, progress debug messages will be printed, |
---|
342 | @* if printlevel>=2, all the debug messages will be printed. |
---|
343 | EXAMPLE: example annfsOT; shows examples |
---|
344 | " |
---|
345 | { |
---|
346 | int eng = 0; |
---|
347 | if ( size(#)>0 ) |
---|
348 | { |
---|
349 | if ( typeof(#[1]) == "int" ) |
---|
350 | { |
---|
351 | eng = int(#[1]); |
---|
352 | } |
---|
353 | } |
---|
354 | // returns a list with a ring and an ideal LD in it |
---|
355 | int ppl = printlevel-voice+2; |
---|
356 | // printf("plevel :%s, voice: %s",printlevel,voice); |
---|
357 | def save = basering; |
---|
358 | int N = nvars(basering); |
---|
359 | int Nnew = 2*(N+2); |
---|
360 | int i,j; |
---|
361 | string s; |
---|
362 | list RL = ringlist(basering); |
---|
363 | list L, Lord; |
---|
364 | list tmp; |
---|
365 | intvec iv; |
---|
366 | L[1] = RL[1]; // char |
---|
367 | L[4] = RL[4]; // char, minpoly |
---|
368 | // check whether vars have admissible names |
---|
369 | list Name = RL[2]; |
---|
370 | list RName; |
---|
371 | RName[1] = "u"; |
---|
372 | RName[2] = "v"; |
---|
373 | RName[3] = "t"; |
---|
374 | RName[4] = "Dt"; |
---|
375 | for(i=1;i<=N;i++) |
---|
376 | { |
---|
377 | for(j=1; j<=size(RName);j++) |
---|
378 | { |
---|
379 | if (Name[i] == RName[j]) |
---|
380 | { |
---|
381 | "Variable names should not include u,v,t,Dt"; |
---|
382 | return(0); |
---|
383 | } |
---|
384 | } |
---|
385 | } |
---|
386 | // now, create the names for new vars |
---|
387 | tmp[1] = "u"; |
---|
388 | tmp[2] = "v"; |
---|
389 | list UName = tmp; |
---|
390 | list DName; |
---|
391 | for(i=1;i<=N;i++) |
---|
392 | { |
---|
393 | DName[i] = "D"+Name[i]; // concat |
---|
394 | } |
---|
395 | tmp = 0; |
---|
396 | tmp[1] = "t"; |
---|
397 | tmp[2] = "Dt"; |
---|
398 | list NName = UName + tmp + Name + DName; |
---|
399 | L[2] = NName; |
---|
400 | tmp = 0; |
---|
401 | // Name, Dname will be used further |
---|
402 | kill UName; |
---|
403 | kill NName; |
---|
404 | // block ord (a(1,1),dp); |
---|
405 | tmp[1] = "a"; // string |
---|
406 | iv = 1,1; |
---|
407 | tmp[2] = iv; //intvec |
---|
408 | Lord[1] = tmp; |
---|
409 | // continue with dp 1,1,1,1... |
---|
410 | tmp[1] = "dp"; // string |
---|
411 | s = "iv="; |
---|
412 | for(i=1;i<=Nnew;i++) |
---|
413 | { |
---|
414 | s = s+"1,"; |
---|
415 | } |
---|
416 | s[size(s)]= ";"; |
---|
417 | execute(s); |
---|
418 | tmp[2] = iv; |
---|
419 | Lord[2] = tmp; |
---|
420 | tmp[1] = "C"; |
---|
421 | iv = 0; |
---|
422 | tmp[2] = iv; |
---|
423 | Lord[3] = tmp; |
---|
424 | tmp = 0; |
---|
425 | L[3] = Lord; |
---|
426 | // we are done with the list |
---|
427 | def @R = ring(L); |
---|
428 | setring @R; |
---|
429 | matrix @D[Nnew][Nnew]; |
---|
430 | @D[3,4]=1; |
---|
431 | for(i=1; i<=N; i++) |
---|
432 | { |
---|
433 | @D[4+i,N+4+i]=1; |
---|
434 | } |
---|
435 | // @D[N+3,2*(N+2)]=1; old t,Dt stuff |
---|
436 | // L[5] = matrix(UpOneMatrix(Nnew)); |
---|
437 | // L[6] = @D; |
---|
438 | ncalgebra(1,@D); |
---|
439 | dbprint(ppl,"// -1-1- the ring @R(u,v,t,Dt,_x,_Dx) is ready"); |
---|
440 | dbprint(ppl-1, @R); |
---|
441 | // create the ideal I |
---|
442 | poly F = imap(save,F); |
---|
443 | ideal I = u*F-t,u*v-1; |
---|
444 | poly p; |
---|
445 | for(i=1; i<=N; i++) |
---|
446 | { |
---|
447 | p = u*Dt; // u*Dt |
---|
448 | p = diff(F,var(4+i))*p; |
---|
449 | I = I, var(N+4+i) + p; |
---|
450 | } |
---|
451 | // -------- the ideal I is ready ---------- |
---|
452 | dbprint(ppl,"// -1-2- starting the elimination of u,v in @R"); |
---|
453 | dbprint(ppl-1, I); |
---|
454 | ideal J = engine(I,eng); |
---|
455 | ideal K = nselect(J,1,2); |
---|
456 | dbprint(ppl,"// -1-3- u,v are eliminated"); |
---|
457 | dbprint(ppl-1, K); // K is without u,v |
---|
458 | setring save; |
---|
459 | // ------------ new ring @R2 ------------------ |
---|
460 | // without u,v and with the elim.ord for t,Dt |
---|
461 | // tensored with the K[s] |
---|
462 | // keep: N, i,j,s, tmp, RL |
---|
463 | Nnew = 2*N+2+1; |
---|
464 | // list RL = ringlist(save); // is defined earlier |
---|
465 | L = 0; // kill L; |
---|
466 | kill Lord, tmp, iv, RName; |
---|
467 | list Lord, tmp; |
---|
468 | intvec iv; |
---|
469 | L[1] = RL[1]; L[4] = RL[4]; // char, minpoly |
---|
470 | // check whether vars have admissible names -> done earlier |
---|
471 | // list Name = RL[2]; |
---|
472 | list RName; |
---|
473 | RName[1] = "t"; |
---|
474 | RName[2] = "Dt"; |
---|
475 | // now, create the names for new var (here, s only) |
---|
476 | tmp[1] = "s"; |
---|
477 | // DName is defined earlier |
---|
478 | list NName = RName + Name + DName + tmp; |
---|
479 | L[2] = NName; |
---|
480 | tmp = 0; |
---|
481 | // block ord (a(1,1),dp); |
---|
482 | tmp[1] = "a"; iv = 1,1; tmp[2] = iv; //intvec |
---|
483 | Lord[1] = tmp; |
---|
484 | // continue with a(1,1,1,1)... |
---|
485 | tmp[1] = "dp"; s = "iv="; |
---|
486 | for(i=1; i<= Nnew; i++) |
---|
487 | { |
---|
488 | s = s+"1,"; |
---|
489 | } |
---|
490 | s[size(s)]= ";"; execute(s); |
---|
491 | kill NName; |
---|
492 | tmp[2] = iv; |
---|
493 | Lord[2] = tmp; |
---|
494 | // extra block for s |
---|
495 | // tmp[1] = "dp"; iv = 1; |
---|
496 | // s[size(s)]= ","; s = s + "1,1,1;"; execute(s); tmp[2] = iv; |
---|
497 | // Lord[3] = tmp; |
---|
498 | kill s; |
---|
499 | tmp[1] = "C"; iv = 0; tmp[2] = iv; |
---|
500 | Lord[3] = tmp; tmp = 0; |
---|
501 | L[3] = Lord; |
---|
502 | // we are done with the list. Now, add a Plural part |
---|
503 | def @R2 = ring(L); |
---|
504 | setring @R2; |
---|
505 | matrix @D[Nnew][Nnew]; |
---|
506 | @D[1,2] = 1; |
---|
507 | for(i=1; i<=N; i++) |
---|
508 | { |
---|
509 | @D[2+i,2+N+i] = 1; |
---|
510 | } |
---|
511 | ncalgebra(1,@D); |
---|
512 | dbprint(ppl,"// -2-1- the ring @R2(t,Dt,_x,_Dx,s) is ready"); |
---|
513 | dbprint(ppl-1, @R2); |
---|
514 | ideal MM = maxideal(1); |
---|
515 | MM = 0,0,MM; |
---|
516 | map R01 = @R, MM; |
---|
517 | ideal K = R01(K); |
---|
518 | // ideal K = imap(@R,K); // names of vars are important! |
---|
519 | poly G = t*Dt+s+1; // s is a variable here |
---|
520 | K = NF(K,std(G)),G; |
---|
521 | // -------- the ideal K_(@R2) is ready ---------- |
---|
522 | dbprint(ppl,"// -2-2- starting the elimination of t,Dt in @R2"); |
---|
523 | dbprint(ppl-1, K); |
---|
524 | ideal M = engine(K,eng); |
---|
525 | ideal K2 = nselect(M,1,2); |
---|
526 | dbprint(ppl,"// -2-3- t,Dt are eliminated"); |
---|
527 | dbprint(ppl-1, K2); |
---|
528 | // dbprint(ppl-1+1," -2-4- std of K2"); |
---|
529 | // option(redSB); option(redTail); K2 = std(K2); |
---|
530 | // K2; // without t,Dt, and with s |
---|
531 | // -------- the ring @R3 ---------- |
---|
532 | // _x, _Dx, s; elim.ord for _x,_Dx. |
---|
533 | // keep: N, i,j,s, tmp, RL |
---|
534 | setring save; |
---|
535 | Nnew = 2*N+1; |
---|
536 | // list RL = ringlist(save); // is defined earlier |
---|
537 | // kill L; |
---|
538 | kill Lord, tmp, iv, RName; |
---|
539 | list Lord, tmp; |
---|
540 | intvec iv; |
---|
541 | L[1] = RL[1]; L[4] = RL[4]; // char, minpoly |
---|
542 | // check whether vars have admissible names -> done earlier |
---|
543 | // list Name = RL[2]; |
---|
544 | // now, create the names for new var (here, s only) |
---|
545 | tmp[1] = "s"; |
---|
546 | // DName is defined earlier |
---|
547 | list NName = Name + DName + tmp; |
---|
548 | L[2] = NName; |
---|
549 | tmp = 0; |
---|
550 | // block ord (a(1,1...),dp); |
---|
551 | string s = "iv="; |
---|
552 | for(i=1; i<=Nnew-1; i++) |
---|
553 | { |
---|
554 | s = s+"1,"; |
---|
555 | } |
---|
556 | s[size(s)]= ";"; |
---|
557 | execute(s); |
---|
558 | tmp[1] = "a"; // string |
---|
559 | tmp[2] = iv; //intvec |
---|
560 | Lord[1] = tmp; |
---|
561 | // continue with dp 1,1,1,1... |
---|
562 | tmp[1] = "dp"; // string |
---|
563 | s[size(s)]=","; s= s+"1;"; |
---|
564 | execute(s); |
---|
565 | kill s; |
---|
566 | kill NName; |
---|
567 | tmp[2] = iv; |
---|
568 | Lord[2] = tmp; |
---|
569 | tmp[1] = "C"; iv = 0; tmp[2] = iv; |
---|
570 | Lord[3] = tmp; tmp = 0; |
---|
571 | L[3] = Lord; |
---|
572 | // we are done with the list. Now add a Plural part |
---|
573 | def @R3 = ring(L); |
---|
574 | setring @R3; |
---|
575 | matrix @D[Nnew][Nnew]; |
---|
576 | for(i=1; i<=N; i++) |
---|
577 | { |
---|
578 | @D[i,N+i]=1; |
---|
579 | } |
---|
580 | ncalgebra(1,@D); |
---|
581 | dbprint(ppl,"// -3-1- the ring @R3(_x,_Dx,s) is ready"); |
---|
582 | dbprint(ppl-1, @R3); |
---|
583 | ideal MM = maxideal(1); |
---|
584 | MM = 0,0,MM; |
---|
585 | map R12 = @R2, MM; |
---|
586 | ideal K = R12(K2); |
---|
587 | poly F = imap(save,F); |
---|
588 | K = K,F; |
---|
589 | dbprint(ppl,"// -3-2- starting the elimination of _x,_Dx in @R3"); |
---|
590 | dbprint(ppl-1, K); |
---|
591 | ideal M = engine(K,eng); |
---|
592 | ideal K3 = nselect(M,1,Nnew-1); |
---|
593 | dbprint(ppl,"// -3-3- _x,_Dx are eliminated in @R3"); |
---|
594 | dbprint(ppl-1, K3); |
---|
595 | // the ring @R4 and the search for minimal negative int s |
---|
596 | ring @R4 = 0,(s),dp; |
---|
597 | dbprint(ppl,"// -4-1- the ring @R4 is ready"); |
---|
598 | ideal K4 = imap(@R3,K3); |
---|
599 | poly p = K4[1]; |
---|
600 | dbprint(ppl,"// -4-2- factorization"); |
---|
601 | ideal P = factorize(p,1); // without constants and multiplicities |
---|
602 | // "------ b-function factorizes into ----------"; P; |
---|
603 | int sP = minIntRoot(P, 1); |
---|
604 | dbprint(ppl,"// -4-3- minimal integer root found"); |
---|
605 | dbprint(ppl-1, sP); |
---|
606 | // convert factors to a list of their roots |
---|
607 | // assume all factors are linear |
---|
608 | ideal BS = normalize(P); |
---|
609 | BS = subst(BS,s,0); |
---|
610 | BS = -BS; |
---|
611 | // TODO: sort BS! |
---|
612 | // ------ substitute s found in the ideal ------ |
---|
613 | // ------- going back to @R2 and substitute -------- |
---|
614 | setring @R2; |
---|
615 | ideal K3 = subst(K2,s,sP); |
---|
616 | // create the ordinary Weyl algebra and put the result into it, |
---|
617 | // thus creating the ring @R5 |
---|
618 | // keep: N, i,j,s, tmp, RL |
---|
619 | setring save; |
---|
620 | Nnew = 2*N; |
---|
621 | // list RL = ringlist(save); // is defined earlier |
---|
622 | kill Lord, tmp, iv; |
---|
623 | L = 0; |
---|
624 | list Lord, tmp; |
---|
625 | intvec iv; |
---|
626 | L[1] = RL[1]; L[4] = RL[4]; // char, minpoly |
---|
627 | // check whether vars have admissible names -> done earlier |
---|
628 | // list Name = RL[2]; |
---|
629 | // DName is defined earlier |
---|
630 | list NName = Name + DName; |
---|
631 | L[2] = NName; |
---|
632 | // dp ordering; |
---|
633 | string s = "iv="; |
---|
634 | for(i=1;i<=Nnew;i++) |
---|
635 | { |
---|
636 | s = s+"1,"; |
---|
637 | } |
---|
638 | s[size(s)]= ";"; |
---|
639 | execute(s); |
---|
640 | tmp = 0; |
---|
641 | tmp[1] = "dp"; // string |
---|
642 | tmp[2] = iv; //intvec |
---|
643 | Lord[1] = tmp; |
---|
644 | kill s; |
---|
645 | tmp[1] = "C"; |
---|
646 | iv = 0; |
---|
647 | tmp[2] = iv; |
---|
648 | Lord[2] = tmp; |
---|
649 | tmp = 0; |
---|
650 | L[3] = Lord; |
---|
651 | // we are done with the list |
---|
652 | // Add: Plural part |
---|
653 | def @R5 = ring(L); |
---|
654 | setring @R5; |
---|
655 | matrix @D[Nnew][Nnew]; |
---|
656 | for(i=1; i<=N; i++) |
---|
657 | { |
---|
658 | @D[i,N+i]=1; |
---|
659 | } |
---|
660 | ncalgebra(1,@D); |
---|
661 | dbprint(ppl,"// -5-1- the ring @R5 is ready"); |
---|
662 | dbprint(ppl-1, @R5); |
---|
663 | ideal K5 = imap(@R2,K3); |
---|
664 | option(redSB); |
---|
665 | dbprint(ppl,"// -5-2- the final cosmetic std"); |
---|
666 | K5 = engine(K5,eng); // std does the job too |
---|
667 | // total cleanup |
---|
668 | kill @R; |
---|
669 | kill @R2; |
---|
670 | kill @R3; |
---|
671 | ideal BS = imap(@R4,BS); |
---|
672 | export BS; |
---|
673 | kill @R4; |
---|
674 | ideal LD = K5; |
---|
675 | export LD; |
---|
676 | return(@R5); |
---|
677 | } |
---|
678 | example |
---|
679 | { |
---|
680 | "EXAMPLE:"; echo = 2; |
---|
681 | ring r = 0,(x,y,z),Dp; |
---|
682 | poly F = x^2+y^3+z^5; |
---|
683 | printlevel = 0; |
---|
684 | def A = annfsOT(F); |
---|
685 | setring A; |
---|
686 | LD; |
---|
687 | print(matrix(BS)); |
---|
688 | } |
---|
689 | |
---|
690 | |
---|
691 | proc annfsgms(poly F, list #) |
---|
692 | "USAGE: annfsgms(f [,eng]); f a poly, eng an optional int |
---|
693 | ASSUME: f has an isolated critical point at 0 |
---|
694 | RETURN: ring |
---|
695 | PURPOSE: compute the D-module structure of basering[f^s] |
---|
696 | NOTE: activate this ring with the @code{setring} command. In this ring, |
---|
697 | @* - the ideal LD is the needed D-mod structure, |
---|
698 | @* - the ideal BS is the list of roots of a Bernstein polynomial of f. |
---|
699 | @* If eng <>0, @code{std} is used for Groebner basis computations, |
---|
700 | @* otherwise (and by default) @code{slimgb} is used. |
---|
701 | @* If printlevel=1, progress debug messages will be printed, |
---|
702 | @* if printlevel>=2, all the debug messages will be printed. |
---|
703 | EXAMPLE: example annfsgms; shows examples |
---|
704 | " |
---|
705 | { |
---|
706 | LIB "gmssing.lib"; |
---|
707 | int eng = 0; |
---|
708 | if ( size(#)>0 ) |
---|
709 | { |
---|
710 | if ( typeof(#[1]) == "int" ) |
---|
711 | { |
---|
712 | eng = int(#[1]); |
---|
713 | } |
---|
714 | } |
---|
715 | int ppl = printlevel-voice+2; |
---|
716 | // returns a ring with the ideal LD in it |
---|
717 | def save = basering; |
---|
718 | // compute the Bernstein poly from gmssing.lib |
---|
719 | list RL = ringlist(basering); |
---|
720 | // in the descr. of the ordering, replace "p" by "s" |
---|
721 | list NL = convloc(RL); |
---|
722 | // create a ring with the ordering, converted to local |
---|
723 | def @LR = ring(NL); |
---|
724 | setring @LR; |
---|
725 | poly F = imap(save, F); |
---|
726 | ideal B = bernstein(F)[1]; |
---|
727 | // since B may not contain (s+1) [following gmssing.lib] |
---|
728 | // add it! |
---|
729 | B = B,-1; |
---|
730 | B = simplify(B,2+4); // erase zero and repeated entries |
---|
731 | // find the minimal integer value |
---|
732 | int S = minIntRoot(B,0); |
---|
733 | dbprint(ppl,"// -0- minimal integer root found"); |
---|
734 | dbprint(ppl-1,S); |
---|
735 | setring save; |
---|
736 | int N = nvars(basering); |
---|
737 | int Nnew = 2*(N+2); |
---|
738 | int i,j; |
---|
739 | string s; |
---|
740 | // list RL = ringlist(basering); |
---|
741 | list L, Lord; |
---|
742 | list tmp; |
---|
743 | intvec iv; |
---|
744 | L[1] = RL[1]; // char |
---|
745 | L[4] = RL[4]; // char, minpoly |
---|
746 | // check whether vars have admissible names |
---|
747 | list Name = RL[2]; |
---|
748 | list RName; |
---|
749 | RName[1] = "u"; |
---|
750 | RName[2] = "v"; |
---|
751 | RName[3] = "t"; |
---|
752 | RName[4] = "Dt"; |
---|
753 | for(i=1;i<=N;i++) |
---|
754 | { |
---|
755 | for(j=1; j<=size(RName);j++) |
---|
756 | { |
---|
757 | if (Name[i] == RName[j]) |
---|
758 | { |
---|
759 | "Variable names should not include u,v,t,Dt"; |
---|
760 | return(0); |
---|
761 | } |
---|
762 | } |
---|
763 | } |
---|
764 | // now, create the names for new vars |
---|
765 | // tmp[1] = "u"; tmp[2] = "v"; tmp[3] = "t"; tmp[4] = "Dt"; |
---|
766 | list UName = RName; |
---|
767 | list DName; |
---|
768 | for(i=1;i<=N;i++) |
---|
769 | { |
---|
770 | DName[i] = "D"+Name[i]; // concat |
---|
771 | } |
---|
772 | list NName = UName + Name + DName; |
---|
773 | L[2] = NName; |
---|
774 | tmp = 0; |
---|
775 | // Name, Dname will be used further |
---|
776 | kill UName; |
---|
777 | kill NName; |
---|
778 | // block ord (a(1,1),dp); |
---|
779 | tmp[1] = "a"; // string |
---|
780 | iv = 1,1; |
---|
781 | tmp[2] = iv; //intvec |
---|
782 | Lord[1] = tmp; |
---|
783 | // continue with dp 1,1,1,1... |
---|
784 | tmp[1] = "dp"; // string |
---|
785 | s = "iv="; |
---|
786 | for(i=1; i<=Nnew; i++) // need really all vars! |
---|
787 | { |
---|
788 | s = s+"1,"; |
---|
789 | } |
---|
790 | s[size(s)]= ";"; |
---|
791 | execute(s); |
---|
792 | tmp[2] = iv; |
---|
793 | Lord[2] = tmp; |
---|
794 | tmp[1] = "C"; |
---|
795 | iv = 0; |
---|
796 | tmp[2] = iv; |
---|
797 | Lord[3] = tmp; |
---|
798 | tmp = 0; |
---|
799 | L[3] = Lord; |
---|
800 | // we are done with the list |
---|
801 | def @R = ring(L); |
---|
802 | setring @R; |
---|
803 | matrix @D[Nnew][Nnew]; |
---|
804 | @D[3,4] = 1; // t,Dt |
---|
805 | for(i=1; i<=N; i++) |
---|
806 | { |
---|
807 | @D[4+i,4+N+i]=1; |
---|
808 | } |
---|
809 | // L[5] = matrix(UpOneMatrix(Nnew)); |
---|
810 | // L[6] = @D; |
---|
811 | ncalgebra(1,@D); |
---|
812 | dbprint(ppl,"// -1-1- the ring @R is ready"); |
---|
813 | dbprint(ppl-1,@R); |
---|
814 | // create the ideal |
---|
815 | poly F = imap(save,F); |
---|
816 | ideal I = u*F-t,u*v-1; |
---|
817 | poly p; |
---|
818 | for(i=1; i<=N; i++) |
---|
819 | { |
---|
820 | p = u*Dt; // u*Dt |
---|
821 | p = diff(F,var(4+i))*p; |
---|
822 | I = I, var(N+4+i) + p; // Dx, Dy |
---|
823 | } |
---|
824 | // add the relations between t,Dt and s |
---|
825 | // I = I, t*Dt+1+S; |
---|
826 | // -------- the ideal I is ready ---------- |
---|
827 | dbprint(ppl,"// -1-2- starting the elimination of u,v in @R"); |
---|
828 | ideal J = engine(I,eng); |
---|
829 | ideal K = nselect(J,1,2); |
---|
830 | dbprint(ppl,"// -1-3- u,v are eliminated in @R"); |
---|
831 | dbprint(ppl-1,K); // without u,v: not yet our answer |
---|
832 | //----- create a ring with elim.ord for t,Dt ------- |
---|
833 | setring save; |
---|
834 | // ------------ new ring @R2 ------------------ |
---|
835 | // without u,v and with the elim.ord for t,Dt |
---|
836 | // keep: N, i,j,s, tmp, RL |
---|
837 | Nnew = 2*N+2; |
---|
838 | // list RL = ringlist(save); // is defined earlier |
---|
839 | kill Lord,tmp,iv, RName; |
---|
840 | L = 0; |
---|
841 | list Lord, tmp; |
---|
842 | intvec iv; |
---|
843 | L[1] = RL[1]; // char |
---|
844 | L[4] = RL[4]; // char, minpoly |
---|
845 | // check whether vars have admissible names -> done earlier |
---|
846 | // list Name = RL[2]; |
---|
847 | list RName; |
---|
848 | RName[1] = "t"; |
---|
849 | RName[2] = "Dt"; |
---|
850 | // DName is defined earlier |
---|
851 | list NName = RName + Name + DName; |
---|
852 | L[2] = NName; |
---|
853 | tmp = 0; |
---|
854 | // block ord (a(1,1),dp); |
---|
855 | tmp[1] = "a"; // string |
---|
856 | iv = 1,1; |
---|
857 | tmp[2] = iv; //intvec |
---|
858 | Lord[1] = tmp; |
---|
859 | // continue with dp 1,1,1,1... |
---|
860 | tmp[1] = "dp"; // string |
---|
861 | s = "iv="; |
---|
862 | for(i=1;i<=Nnew;i++) |
---|
863 | { |
---|
864 | s = s+"1,"; |
---|
865 | } |
---|
866 | s[size(s)]= ";"; |
---|
867 | execute(s); |
---|
868 | kill s; |
---|
869 | kill NName; |
---|
870 | tmp[2] = iv; |
---|
871 | Lord[2] = tmp; |
---|
872 | tmp[1] = "C"; |
---|
873 | iv = 0; |
---|
874 | tmp[2] = iv; |
---|
875 | Lord[3] = tmp; |
---|
876 | tmp = 0; |
---|
877 | L[3] = Lord; |
---|
878 | // we are done with the list |
---|
879 | // Add: Plural part |
---|
880 | def @R2 = ring(L); |
---|
881 | setring @R2; |
---|
882 | matrix @D[Nnew][Nnew]; |
---|
883 | @D[1,2]=1; |
---|
884 | for(i=1; i<=N; i++) |
---|
885 | { |
---|
886 | @D[2+i,2+N+i]=1; |
---|
887 | } |
---|
888 | ncalgebra(1,@D); |
---|
889 | dbprint(ppl,"// -2-1- the ring @R2 is ready"); |
---|
890 | dbprint(ppl-1,@R2); |
---|
891 | ideal MM = maxideal(1); |
---|
892 | MM = 0,0,MM; |
---|
893 | map R01 = @R, MM; |
---|
894 | ideal K2 = R01(K); |
---|
895 | // add the relations between t,Dt and s |
---|
896 | // K2 = K2, t*Dt+1+S; |
---|
897 | poly G = t*Dt+S+1; |
---|
898 | K2 = NF(K2,std(G)),G; |
---|
899 | dbprint(ppl,"// -2-2- starting elimination for t,Dt in @R2"); |
---|
900 | ideal J = engine(K2,eng); |
---|
901 | ideal K = nselect(J,1,2); |
---|
902 | dbprint(ppl,"// -2-3- t,Dt are eliminated"); |
---|
903 | dbprint(ppl-1,K); |
---|
904 | // "------- produce a final result --------"; |
---|
905 | // ----- create the ordinary Weyl algebra and put |
---|
906 | // ----- the result into it |
---|
907 | // ------ create the ring @R5 |
---|
908 | // keep: N, i,j,s, tmp, RL |
---|
909 | setring save; |
---|
910 | Nnew = 2*N; |
---|
911 | // list RL = ringlist(save); // is defined earlier |
---|
912 | kill Lord, tmp, iv; |
---|
913 | // kill L; |
---|
914 | L = 0; |
---|
915 | list Lord, tmp; |
---|
916 | intvec iv; |
---|
917 | L[1] = RL[1]; // char |
---|
918 | L[4] = RL[4]; // char, minpoly |
---|
919 | // check whether vars have admissible names -> done earlier |
---|
920 | // list Name = RL[2]; |
---|
921 | // DName is defined earlier |
---|
922 | list NName = Name + DName; |
---|
923 | L[2] = NName; |
---|
924 | // dp ordering; |
---|
925 | string s = "iv="; |
---|
926 | for(i=1;i<=2*N;i++) |
---|
927 | { |
---|
928 | s = s+"1,"; |
---|
929 | } |
---|
930 | s[size(s)]= ";"; |
---|
931 | execute(s); |
---|
932 | tmp = 0; |
---|
933 | tmp[1] = "dp"; // string |
---|
934 | tmp[2] = iv; //intvec |
---|
935 | Lord[1] = tmp; |
---|
936 | kill s; |
---|
937 | tmp[1] = "C"; |
---|
938 | iv = 0; |
---|
939 | tmp[2] = iv; |
---|
940 | Lord[2] = tmp; |
---|
941 | tmp = 0; |
---|
942 | L[3] = Lord; |
---|
943 | // we are done with the list |
---|
944 | // Add: Plural part |
---|
945 | def @R5 = ring(L); |
---|
946 | setring @R5; |
---|
947 | matrix @D[Nnew][Nnew]; |
---|
948 | for(i=1; i<=N; i++) |
---|
949 | { |
---|
950 | @D[i,N+i]=1; |
---|
951 | } |
---|
952 | ncalgebra(1,@D); |
---|
953 | dbprint(ppl,"// -3-1- the ring @R5 is ready"); |
---|
954 | dbprint(ppl-1,@R5); |
---|
955 | ideal K5 = imap(@R2,K); |
---|
956 | option(redSB); |
---|
957 | dbprint(ppl,"// -3-2- the final cosmetic std"); |
---|
958 | K5 = engine(K5,eng); // std does the job too |
---|
959 | // total cleanup |
---|
960 | kill @R; |
---|
961 | kill @R2; |
---|
962 | ideal LD = K5; |
---|
963 | ideal BS = imap(@LR,B); |
---|
964 | kill @LR; |
---|
965 | export BS; |
---|
966 | export LD; |
---|
967 | return(@R5); |
---|
968 | } |
---|
969 | example |
---|
970 | { |
---|
971 | "EXAMPLE:"; echo = 2; |
---|
972 | ring r = 0,(x,y,z),Dp; |
---|
973 | poly F = x^2+y^3+z^5; |
---|
974 | def A = annfsgms(F); |
---|
975 | setring A; |
---|
976 | LD; |
---|
977 | print(matrix(BS)); |
---|
978 | } |
---|
979 | |
---|
980 | proc convloc(list @NL) |
---|
981 | "USAGE: convloc(L); L a list |
---|
982 | RETURN: list |
---|
983 | PURPOSE: convert a ringlist L into another ringlist, |
---|
984 | where all the 'p' orderings are replaced with the 's' orderings. |
---|
985 | ASSUME: L is a result of a ringlist command |
---|
986 | EXAMPLE: example minIntRoot; shows examples |
---|
987 | " |
---|
988 | { |
---|
989 | list NL = @NL; |
---|
990 | // given a ringlist, returns a new ringlist, |
---|
991 | // where all the p-orderings are replaced with s-ord's |
---|
992 | int i,j,k,found; |
---|
993 | int nblocks = size(NL[3]); |
---|
994 | for(i=1; i<=nblocks; i++) |
---|
995 | { |
---|
996 | for(j=1; j<=size(NL[3][i]); j++) |
---|
997 | { |
---|
998 | if (typeof(NL[3][i][j]) == "string" ) |
---|
999 | { |
---|
1000 | found = 0; |
---|
1001 | for (k=1; k<=size(NL[3][i][j]); k++) |
---|
1002 | { |
---|
1003 | if (NL[3][i][j][k]=="p") |
---|
1004 | { |
---|
1005 | NL[3][i][j][k]="s"; |
---|
1006 | found = 1; |
---|
1007 | // printf("replaced at %s,%s,%s",i,j,k); |
---|
1008 | } |
---|
1009 | } |
---|
1010 | } |
---|
1011 | } |
---|
1012 | } |
---|
1013 | return(NL); |
---|
1014 | } |
---|
1015 | example |
---|
1016 | { |
---|
1017 | "EXAMPLE:"; echo = 2; |
---|
1018 | ring r = 0,(x,y,z),(Dp(2),dp(1)); |
---|
1019 | list L = ringlist(r); |
---|
1020 | list N = convloc(L); |
---|
1021 | def rs = ring(N); |
---|
1022 | setring rs; |
---|
1023 | rs; |
---|
1024 | } |
---|
1025 | |
---|
1026 | |
---|
1027 | proc minIntRoot(ideal P, int fact) |
---|
1028 | "USAGE: minIntRoot(P, fact); P an ideal, fact an int |
---|
1029 | RETURN: int |
---|
1030 | PURPOSE: minimal integer root of a maximal ideal P |
---|
1031 | NOTE: if fact==1, P is the result of some 'factorize' call, |
---|
1032 | @* else P is treated as the result of bernstein::gmssing.lib |
---|
1033 | @* in both cases without constants and multiplicities |
---|
1034 | EXAMPLE: example minIntRoot; shows examples |
---|
1035 | " |
---|
1036 | { |
---|
1037 | // ideal P = factorize(p,1); |
---|
1038 | // or ideal P = bernstein(F)[1]; |
---|
1039 | intvec vP; |
---|
1040 | number nP; |
---|
1041 | int i; |
---|
1042 | if ( fact ) |
---|
1043 | { |
---|
1044 | // the result comes from "factorize" |
---|
1045 | P = normalize(P); // now leadcoef = 1 |
---|
1046 | P = lead(P)-P; |
---|
1047 | // nP = leadcoef(P[i]-lead(P[i])); // for 1 var only, extract the coeff |
---|
1048 | } |
---|
1049 | else |
---|
1050 | { |
---|
1051 | // bernstein deletes -1 usually |
---|
1052 | P = P, -1; |
---|
1053 | } |
---|
1054 | // for both situations: |
---|
1055 | // now we have an ideal of fractions of type "number" |
---|
1056 | int sP = size(P); |
---|
1057 | for(i=1; i<=sP; i++) |
---|
1058 | { |
---|
1059 | nP = leadcoef(P[i]); |
---|
1060 | if ( (nP - int(nP)) == 0 ) |
---|
1061 | { |
---|
1062 | vP = vP,int(nP); |
---|
1063 | } |
---|
1064 | } |
---|
1065 | if ( size(vP)>=2 ) |
---|
1066 | { |
---|
1067 | vP = vP[2..size(vP)]; |
---|
1068 | } |
---|
1069 | sP = -Max(-vP); |
---|
1070 | if (sP == 0) |
---|
1071 | { |
---|
1072 | "Warning: zero root!"; |
---|
1073 | } |
---|
1074 | return(sP); |
---|
1075 | } |
---|
1076 | example |
---|
1077 | { |
---|
1078 | "EXAMPLE:"; echo = 2; |
---|
1079 | ring r = 0,(x,y),ds; |
---|
1080 | poly f1 = x*y*(x+y); |
---|
1081 | ideal I1 = bernstein(f1)[1]; |
---|
1082 | minIntRoot(I1,0); |
---|
1083 | poly f2 = x2-y3; |
---|
1084 | ideal I2 = bernstein(f2)[1]; |
---|
1085 | minIntRoot(I2,0); |
---|
1086 | // now we illustrate the behaviour of factorize |
---|
1087 | ring r2 = 0,x,dp; |
---|
1088 | poly f3 = 9*(x+2/3)*(x+1)*(x+4/3); // b-poly of f1=x*y*(x+y) |
---|
1089 | ideal I3 = factorize(f3,1); |
---|
1090 | minIntRoot(I3,1); |
---|
1091 | // and a more interesting situation |
---|
1092 | ring s = 0,(x,y,z),ds; |
---|
1093 | poly f = x3 + y3 + z3; |
---|
1094 | ideal I = bernstein(f)[1]; |
---|
1095 | minIntRoot(I,0); |
---|
1096 | } |
---|
1097 | |
---|
1098 | proc isHolonomic(def M) |
---|
1099 | "USAGE: isHolonomic(M); M an ideal/module/matrix |
---|
1100 | RETURN: int, 1 if M is holonomic and 0 otherwise |
---|
1101 | PURPOSE: check the modules for the property of holonomy |
---|
1102 | NOTE: M is holonomic, if 2*dim(M) = dim(R), where R is a |
---|
1103 | ground ring |
---|
1104 | EXAMPLE: example isHolonomic; shows examples |
---|
1105 | " |
---|
1106 | { |
---|
1107 | if ( (typeof(M) != "ideal") && (typeof(M) != "module") && (typeof(M) != "matrix") ) |
---|
1108 | { |
---|
1109 | // print(typeof(M)); |
---|
1110 | ERROR("an argument of type ideal/module/matrix expected"); |
---|
1111 | } |
---|
1112 | if (attrib(M,"isSB")!=1) |
---|
1113 | { |
---|
1114 | M = std(M); |
---|
1115 | } |
---|
1116 | int dimR = gkdim(std(0)); |
---|
1117 | int dimM = gkdim(M); |
---|
1118 | return( (dimR==2*dimM) ); |
---|
1119 | } |
---|
1120 | example |
---|
1121 | { |
---|
1122 | "EXAMPLE:"; echo = 2; |
---|
1123 | ring R = 0,(x,y),dp; |
---|
1124 | poly F = x*y*(x+y); |
---|
1125 | def A = annfsBM(F,0); |
---|
1126 | setring A; |
---|
1127 | LD; |
---|
1128 | isHolonomic(LD); |
---|
1129 | ideal I = std(LD[1]); |
---|
1130 | isHolonomic(I); |
---|
1131 | } |
---|
1132 | |
---|
1133 | proc reiffen(int p, int q) |
---|
1134 | "USAGE: reiffen(p, q); int p, int q |
---|
1135 | RETURN: ring |
---|
1136 | PURPOSE: set up the polynomial, describing a Reiffen curve |
---|
1137 | NOTE: activate this ring with the @code{setring} command and find the curve as a polynomial RC |
---|
1138 | @* a Reiffen curve is defined as F = x^p + y^q + xy^{q-1}, q >= p+1 >= 5 |
---|
1139 | ASSUME: q >= p+1 >= 5. Otherwise an error message is returned |
---|
1140 | EXAMPLE: example reiffen; shows examples |
---|
1141 | " |
---|
1142 | { |
---|
1143 | // a Reiffen curve is defined as |
---|
1144 | // F = x^p + y^q +x*y^{q-1}, q \geq p+1 \geq 5 |
---|
1145 | |
---|
1146 | if ( (p<4) || (q<5) || (q-p<1) ) |
---|
1147 | { |
---|
1148 | "Some of conditions p>=4, q>=5 or q>=p+1 is not satisfied!"; |
---|
1149 | return(0); |
---|
1150 | } |
---|
1151 | ring @r = 0,(x,y),ds; |
---|
1152 | poly RC = y^q +x^p + x*y^(q-1); |
---|
1153 | export RC; |
---|
1154 | return(@r); |
---|
1155 | } |
---|
1156 | example |
---|
1157 | { |
---|
1158 | "EXAMPLE:"; echo = 2; |
---|
1159 | def r = reiffen(4,5); |
---|
1160 | setring r; |
---|
1161 | RC; |
---|
1162 | } |
---|
1163 | |
---|
1164 | proc arrange(int p) |
---|
1165 | "USAGE: arrange(p); int p |
---|
1166 | RETURN: ring |
---|
1167 | PURPOSE: set up the polynomial, describing a generic hyperplane arrangement |
---|
1168 | NOTE: must be executed in a ring |
---|
1169 | ASSUME: basering is present |
---|
1170 | EXAMPLE: example arrange; shows examples |
---|
1171 | " |
---|
1172 | { |
---|
1173 | int UseBasering = 0 ; |
---|
1174 | if (p<2) |
---|
1175 | { |
---|
1176 | ERROR("p>=2 is required!"); |
---|
1177 | } |
---|
1178 | if ( nameof(basering)!="basering" ) |
---|
1179 | { |
---|
1180 | if (p > nvars(basering)) |
---|
1181 | { |
---|
1182 | ERROR("too big p"); |
---|
1183 | } |
---|
1184 | else |
---|
1185 | { |
---|
1186 | def @r = basering; |
---|
1187 | UseBasering = 1; |
---|
1188 | } |
---|
1189 | } |
---|
1190 | else |
---|
1191 | { |
---|
1192 | // no basering found |
---|
1193 | ERROR("no basering found!"); |
---|
1194 | // ring @r = 0,(x(1..p)),dp; |
---|
1195 | } |
---|
1196 | int i,j,sI; |
---|
1197 | poly q = 1; |
---|
1198 | list ar; |
---|
1199 | ideal tmp; |
---|
1200 | tmp = ideal(var(1)); |
---|
1201 | ar[1] = tmp; |
---|
1202 | for (i = 2; i<=p; i++) |
---|
1203 | { |
---|
1204 | // add i-nary sums to the product |
---|
1205 | ar = indAR(ar,i); |
---|
1206 | } |
---|
1207 | for (i = 1; i<=p; i++) |
---|
1208 | { |
---|
1209 | tmp = ar[i]; sI = size(tmp); |
---|
1210 | for (j = 1; j<=sI; j++) |
---|
1211 | { |
---|
1212 | q = q*tmp[j]; |
---|
1213 | } |
---|
1214 | } |
---|
1215 | if (UseBasering) |
---|
1216 | { |
---|
1217 | return(q); |
---|
1218 | } |
---|
1219 | // poly AR = q; export AR; |
---|
1220 | // return(@r); |
---|
1221 | } |
---|
1222 | example |
---|
1223 | { |
---|
1224 | "EXAMPLE:"; echo = 2; |
---|
1225 | ring X = 0,(x,y,z,t),dp; |
---|
1226 | poly q = arrange(3); |
---|
1227 | factorize(q,1); |
---|
1228 | } |
---|
1229 | |
---|
1230 | static proc indAR(list L, int n) |
---|
1231 | "USAGE: indAR(L,n); list L, int n |
---|
1232 | RETURN: list |
---|
1233 | PURPOSE: computes arrangement inductively, using L and varn(n) as the |
---|
1234 | next variable |
---|
1235 | ASSUME: L has a structure of an arrangement |
---|
1236 | EXAMPLE: example indAR; shows examples |
---|
1237 | " |
---|
1238 | { |
---|
1239 | if ( (n<2) || (n>nvars(basering)) ) |
---|
1240 | { |
---|
1241 | ERROR("incorrect n"); |
---|
1242 | } |
---|
1243 | int sl = size(L); |
---|
1244 | list K; |
---|
1245 | ideal tmp; |
---|
1246 | poly @t = L[sl][1] + var(n); //1 elt |
---|
1247 | K[sl+1] = ideal(@t); |
---|
1248 | tmp = L[1]+var(n); |
---|
1249 | K[1] = tmp; tmp = 0; |
---|
1250 | int i,j,sI; |
---|
1251 | ideal I; |
---|
1252 | for(i=sl; i>=2; i--) |
---|
1253 | { |
---|
1254 | I = L[i-1]; sI = size(I); |
---|
1255 | for(j=1; j<=sI; j++) |
---|
1256 | { |
---|
1257 | I[j] = I[j] + var(n); |
---|
1258 | } |
---|
1259 | tmp = L[i],I; |
---|
1260 | K[i] = tmp; |
---|
1261 | I = 0; tmp = 0; |
---|
1262 | } |
---|
1263 | kill I; kill tmp; |
---|
1264 | return(K); |
---|
1265 | } |
---|
1266 | example |
---|
1267 | { |
---|
1268 | "EXAMPLE:"; echo = 2; |
---|
1269 | ring r = 0,(x,y,z,t,v),dp; |
---|
1270 | list L; |
---|
1271 | L[1] = ideal(x); |
---|
1272 | list K = indAR(L,2); |
---|
1273 | K; |
---|
1274 | list M = indAR(K,3); |
---|
1275 | M; |
---|
1276 | M = indAR(M,4); |
---|
1277 | M; |
---|
1278 | } |
---|
1279 | |
---|
1280 | static proc exCheckGenericity() |
---|
1281 | { |
---|
1282 | LIB "control.lib"; |
---|
1283 | ring r = (0,a,b,c),x,dp; |
---|
1284 | poly p = (x-a)*(x-b)*(x-c); |
---|
1285 | def A = annfsBM(p); |
---|
1286 | setring A; |
---|
1287 | ideal J = slimgb(LD); |
---|
1288 | matrix T = lift(LD,J); |
---|
1289 | T = normalize(T); |
---|
1290 | genericity(T); |
---|
1291 | // Ann =x^3*Dx+3*x^2*t*Dt+(-a-b-c)*x^2*Dx+(-2*a-2*b-2*c)*x*t*Dt+3*x^2+(a*b+a*c+b*c)*x*Dx+(a*b+a*c+b*c)*t*Dt+(-2*a-2*b-2*c)*x+(-a*b*c)*Dx+(a*b+a*c+b*c) |
---|
1292 | // genericity: g = a2-ab-ac+b2-bc+c2 =0 |
---|
1293 | // g = (a -(b+c)/2)^2 + (3/4)*(b-c)^2; |
---|
1294 | // g ==0 <=> a=b=c |
---|
1295 | // indeed, Ann = (x-a)^2*(x*Dx+3*t*Dt+(-a)*Dx+3) |
---|
1296 | // -------------------------------------------- |
---|
1297 | // BUT a direct computation shows |
---|
1298 | // when a=b=c, |
---|
1299 | // Ann = x*Dx+3*t*Dt+(-a)*Dx+3 |
---|
1300 | } |
---|
1301 | |
---|
1302 | static proc exOT_17() |
---|
1303 | { |
---|
1304 | // Oaku-Takayama, p.208 |
---|
1305 | ring R = 0,(x,y),dp; |
---|
1306 | poly F = x^3-y^2; // x^2+x*y+y^2; |
---|
1307 | option(prot); |
---|
1308 | option(mem); |
---|
1309 | // option(redSB); |
---|
1310 | def A = annfsOT(F,0); |
---|
1311 | setring A; |
---|
1312 | LD; |
---|
1313 | LIB "gkdim.lib"; |
---|
1314 | gkdim(LD); // a holonomic check |
---|
1315 | // poly F = x^3-y^2; // = x^7 - y^5; // x^3-y^4; // x^5 - y^4; |
---|
1316 | } |
---|
1317 | |
---|
1318 | static proc exOT_16() |
---|
1319 | { |
---|
1320 | // Oaku-Takayama, p.208 |
---|
1321 | ring R = 0,(x),dp; |
---|
1322 | poly F = x*(1-x); |
---|
1323 | option(prot); |
---|
1324 | option(mem); |
---|
1325 | // option(redSB); |
---|
1326 | def A = annfsOT(F,0); |
---|
1327 | setring A; |
---|
1328 | LD; |
---|
1329 | LIB "gkdim.lib"; |
---|
1330 | gkdim(LD); // a holonomic check |
---|
1331 | // poly F = x^3-y^2; // = x^7 - y^5; // x^3-y^4; // x^5 - y^4; |
---|
1332 | } |
---|
1333 | |
---|
1334 | static proc ex_bcheck() |
---|
1335 | { |
---|
1336 | ring R = 0,(x,y),dp; |
---|
1337 | poly F = x*y*(x+y); |
---|
1338 | option(prot); |
---|
1339 | option(mem); |
---|
1340 | int eng = 0; |
---|
1341 | // option(redSB); |
---|
1342 | def A = annfsOT(F,eng); |
---|
1343 | // def A = annfsgms(F,0); |
---|
1344 | setring A; |
---|
1345 | LD; |
---|
1346 | } |
---|
1347 | |
---|
1348 | static proc ex_bcheck2() |
---|
1349 | { |
---|
1350 | ring R = 0,(x,y),dp; |
---|
1351 | poly F = x*y*(x+y); |
---|
1352 | int eng = 0; |
---|
1353 | def A = annfsOT(F,eng); |
---|
1354 | // def A = annfsgms(F,0); |
---|
1355 | setring A; |
---|
1356 | LD; |
---|
1357 | } |
---|