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