# source:git/Singular/LIB/dmod.lib@437f97

spielwiese
Last change on this file since 437f97 was 437f97, checked in by Viktor Levandovskyy <levandov@…>, 17 years ago
*levandov: cosmetic changes to names and descroption git-svn-id: file:///usr/local/Singular/svn/trunk@9113 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100755`
File size: 32.4 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id: dmod.lib,v 1.3 2006-05-08 14:43:21 levandov Exp \$";
3category="Noncommutative";
4info="
5LIBRARY: dmod.lib     Algorithms for algebraic D-modules
6AUTHORS: Viktor Levandovskyy,     levandov@mathematik.uni-kl.de
7@*         Jorge Martin Morales,      jorge@unizar.es
8
9THEORY: 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
12that 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
17PROCEDURES:
18annfsOT(F[,eng]);    compute Ann F^s for a poly F with the algorithm of Oaku-Takayama
19annfsBM(F[,eng]);    compute Ann F^s for a poly F with the algorithm of Briancon-Maisonobe
20reiffen(p,q);      create the polynomial, describing a Reiffen curve
21arrange(p);        create the polynomial, describing a generic hyperplane arrangement
22isHolonomic(M);    check whether a module is holonomic
23convloc(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
31LIB "nctools.lib";
32LIB "elim.lib";
33LIB "qhmoduli.lib"; // for Max
34LIB "gkdim.lib";
35
36proc 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
54proc annfsBM(poly F, list #)
55"USAGE:  annfsBM(f [,eng]);  f a poly, eng an optional int
56RETURN:  ring
57PURPOSE: compute the D-module structure of basering[f^s], according
58to the algorithm by Briancon and Maisonobe
59NOTE:    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.
66EXAMPLE: 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
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}
319example
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
331proc annfsOT(poly F, list #)
332"USAGE:  annfsOT(f [,eng]);  f a poly, eng an optional int
333RETURN:  ring
334PURPOSE: compute the D-module structure of basering[f^s], according
335to the algorithm by Oaku and Takayama
336NOTE:    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.
343EXAMPLE: 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
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}
678example
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
691proc annfsgms(poly F, list #)
692"USAGE:  annfsgms(f [,eng]);  f a poly, eng an optional int
693ASSUME:  f has an isolated critical point at 0
694RETURN:  ring
695PURPOSE: compute the D-module structure of basering[f^s]
696NOTE:    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.
703EXAMPLE: 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]
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
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
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}
969example
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
980proc convloc(list @NL)
981"USAGE:  convloc(L); L a list
982RETURN:  list
983PURPOSE: convert a ringlist L into another ringlist,
984where all the 'p' orderings are replaced with the 's' orderings.
985ASSUME: L is a result of a ringlist command
986EXAMPLE: 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}
1015example
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
1027proc minIntRoot(ideal P, int fact)
1028"USAGE:  minIntRoot(P, fact); P an ideal, fact an int
1029RETURN:  int
1030PURPOSE: minimal integer root of a maximal ideal P
1031NOTE:    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
1034EXAMPLE: 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
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  {
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}
1076example
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
1098proc isHolonomic(def M)
1099"USAGE:  isHolonomic(M); M an ideal/module/matrix
1100RETURN:  int, 1 if M is holonomic and 0 otherwise
1101PURPOSE: check the modules for the property of holonomy
1102NOTE:    M is holonomic, if 2*dim(M) = dim(R), where R is a
1103ground ring
1104EXAMPLE: 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}
1120example
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
1133proc reiffen(int p, int q)
1134"USAGE:  reiffen(p, q);  int p, int q
1135RETURN:  ring
1136PURPOSE: set up the polynomial, describing a Reiffen curve
1137NOTE:    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
1139ASSUME: q >= p+1 >= 5. Otherwise an error message is returned
1140EXAMPLE: 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}
1156example
1157
1158  "EXAMPLE:"; echo = 2;
1159  def r = reiffen(4,5);
1160  setring r;
1161  RC;
1162}
1163
1164proc arrange(int p)
1165"USAGE:  arrange(p);  int p
1166RETURN:  ring
1167PURPOSE: set up the polynomial, describing a generic hyperplane arrangement
1168NOTE:   must be executed in a ring
1169ASSUME: basering is present
1170EXAMPLE: 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}
1222example
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
1230static proc indAR(list L, int n)
1231"USAGE:  indAR(L,n);  list L, int n
1232RETURN:  list
1233PURPOSE: computes arrangement inductively, using L and varn(n) as the
1234next variable
1235ASSUME: L has a structure of an arrangement
1236EXAMPLE: 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}
1266example
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
1280static 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
1302static 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
1318static 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
1334static 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
1348static 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}
Note: See TracBrowser for help on using the repository browser.