source: git/Singular/LIB/dmodapp.lib @ 6fe9a5

spielwiese
Last change on this file since 6fe9a5 was 6fe9a5, checked in by Viktor Levandovskyy <levandov@…>, 12 years ago
*levandov: bugfixes for bfun_lib and dmodapp_lib by Daniel Andres git-svn-id: file:///usr/local/Singular/svn/trunk@14018 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 90.4 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Noncommutative";
4info="
5LIBRARY: dmodapp.lib     Applications of algebraic D-modules
6AUTHORS: Viktor Levandovskyy,  levandov@math.rwth-aachen.de
7@*       Daniel Andres,   daniel.andres@math.rwth-aachen.de
8
9Support: DFG Graduiertenkolleg 1632 'Experimentelle und konstruktive Algebra'
10
11OVERVIEW:
12Let K be a field of characteristic 0, R = K[x1,...,xN] and
13D be the Weyl algebra in variables x1,...,xN,d1,...,dN.
14In this library there are the following procedures for algebraic D-modules:
15
16@* - given a cyclic representation D/I of a holonomic module and a polynomial
17 F in R, it is proved that the localization of D/I with respect to the mult.
18 closed set of all powers of F is a holonomic D-module. Thus we aim to compute
19 its cyclic representaion D/L for an ideal L in D. The procedures for the
20 localization are DLoc, SDLoc and DLoc0.
21
22@* - annihilator in D of a given polynomial F from R as well as
23 of a given rational function G/F from Quot(R). These can be computed via
24 procedures annPoly resp. annRat.
25
26@* - Groebner bases with respect to weights (according to (SST), given an
27 arbitrary integer vector containing weights for variables, one computes the
28 homogenization of a given ideal relative to this vector, then one computes a
29 Groebner basis and returns the dehomogenization of the result), initial
30 forms and initial ideals in Weyl algebras with respect to a given weight
31 vector can be computed with GBWeight, inForm, initialMalgrange and
32 initialIdealW.
33
34@* - restriction and integration of a holonomic module D/I. Suppose I
35 annihilates a function F(x1,...,xn). Our aim is to compute an ideal J
36 directly from I, which annihilates
37@*   - F(0,...,0,xk,...,xn) in case of restriction or
38@*   - the integral of F with respect to x1,...,xm in case of integration.
39 The corresponding procedures are restrictionModule, restrictionIdeal,
40 integralModule and integralIdeal.
41
42@* - characteristic varieties defined by ideals in Weyl algebras can be computed
43 with charVariety and charInfo.
44
45@* - appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebras,
46 which annihilate corresponding Appel hypergeometric functions.
47
48
49References:
50@* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric
51         Differential Equations', Springer, 2000
52@* (OTW) Oaku, Takayama, Walther 'A Localization Algorithm for D-modules',
53         Journal of Symbolic Computation, 2000
54@* (OT)  Oaku, Takayama 'Algorithms for D-modules',
55         Journal of Pure and Applied Algebra, 1998
56
57
58PROCEDURES:
59
60annPoly(f);  computes annihilator of a polynomial f in the corr. Weyl algebra
61annRat(f,g); computes annihilator of rational function f/g in corr. Weyl algebra
62DLoc(I,f);   computes presentation of localization of D/I wrt symbolic power f^s
63SDLoc(I,f);  computes generic presentation of the localization of D/I wrt f^s
64DLoc0(I,f);  computes presentation of localization of D/I wrt f^s based on SDLoc
65
66GBWeight(I,u,v[,a,b]);       computes Groebner basis of I wrt a weight vector
67initialMalgrange(f[,s,t,v]); computes Groebner basis of initial Malgrange ideal
68initialIdealW(I,u,v[,s,t]);  computes initial ideal of  wrt a given weight
69inForm(f,w);                 computes initial form of poly/ideal wrt a weight
70
71restrictionIdeal(I,w[,eng,m,G]);  computes restriction ideal of I wrt w
72restrictionModule(I,w[,eng,m,G]); computes restriction module of I wrt w
73integralIdeal(I,w[,eng,m,G]);     computes integral ideal of I wrt w
74integralModule(I,w[,eng,m,G]);    computes integral module of I wrt w
75deRhamCohom(f[,eng,m]);           computes basis of n-th de Rham cohom. group
76deRhamCohomIdeal(I[,w,eng,m,G]);  computes basis of n-th de Rham cohom. group
77
78charVariety(I); computes characteristic variety of the ideal I
79charInfo(I);    computes char. variety, singular locus and primary decomp.
80isFsat(I,F);    checks whether the ideal I is F-saturated
81
82
83
84appelF1();     creates an ideal annihilating Appel F1 function
85appelF2();     creates an ideal annihilating Appel F2 function
86appelF4();     creates an ideal annihilating Appel F4 function
87
88fourier(I[,v]);        applies Fourier automorphism to ideal
89inverseFourier(I[,v]); applies inverse Fourier automorphism to ideal
90
91bFactor(F);    computes the roots of irreducible factors of an univariate poly
92intRoots(L);   dismisses non-integer roots from list in bFactor format
93poly2list(f);  decomposes the polynomial f into a list of terms and exponents
94fl2poly(L,s);  reconstructs a monic univariate polynomial from its factorization
95
96insertGenerator(id,p[,k]); inserts an element into an ideal/module
97deleteGenerator(id,k);     deletes the k-th element from an ideal/module
98
99engine(I,i);   computes a Groebner basis with the algorithm specified by i
100isInt(n);      checks whether number n is actually an int
101sortIntvec(v); sorts intvec
102
103KEYWORDS: D-module; annihilator of polynomial; annihilator of rational function;
104D-localization; localization of D-module; D-restriction; restriction of
105D-module; D-integration; integration of D-module; characteristic variety;
106Appel function; Appel hypergeometric function
107
108SEE ALSO: bfun_lib, dmod_lib, dmodvar_lib, gmssing_lib
109";
110
111/*
112  Changelog
113  21.09.10 by DA:
114  - restructured library for better readability
115  - new / improved procs:
116  - toolbox: isInt, intRoots, sortIntvec
117  - GB wrt weights: GBWeight, initialIdealW rewritten using GBWeight
118  - restriction/integration: restrictionX, integralX where X in {Module, Ideal},
119  fourier, inverseFourier, deRhamCohom, deRhamCohomIdeal
120  - characteristic variety: charVariety, charInfo
121  - added keywords for features
122  - reformated help strings and examples in existing procs
123  - added SUPPORT in header
124  - added reference (OT)
125
126  04.10.10 by DA:
127  - incorporated suggestions by Oleksandr Motsak, among other:
128  - bugfixes for fl2poly, sortIntvec, annPoly, GBWeight
129  - enhanced functionality for deleteGenerator, inForm
130
131  11.10.10 by DA:
132  - procedure bFactor now sorts the roots using new static procedure sortNumberIdeal
133
134  17.03.11 by DA:
135  - bugfixes for inForm with polynomial input, typo in restrictionIdealEngine
136*/
137
138
139LIB "bfun.lib";     // for pIntersect etc
140LIB "dmod.lib";     // for SannfsBM etc
141LIB "general.lib";  // for sort
142LIB "gkdim.lib";
143LIB "nctools.lib";  // for isWeyl etc
144LIB "poly.lib";
145LIB "primdec.lib";  // for primdecGKZ
146LIB "qhmoduli.lib"; // for Max
147LIB "sing.lib";     // for slocus
148
149
150///////////////////////////////////////////////////////////////////////////////
151// testing for consistency of the library:
152proc testdmodapp()
153{
154  example annPoly;
155  example annRat;
156  example DLoc;
157  example SDLoc;
158  example DLoc0;
159  example GBWeight;
160  example initialMalgrange;
161  example initialIdealW;
162  example inForm;
163  example restrictionIdeal;
164  example restrictionModule;
165  example integralIdeal;
166  example integralModule;
167  example deRhamCohom;
168  example deRhamCohomIdeal;
169  example charVariety;
170  example charInfo;
171  example isFsat;
172  example appelF1;
173  example appelF2;
174  example appelF4;
175  example fourier;
176  example inverseFourier;
177  example bFactor;
178  example intRoots;
179  example poly2list;
180  example fl2poly;
181  example insertGenerator;
182  example deleteGenerator;
183  example engine;
184  example isInt;
185  example sortIntvec;
186}
187
188
189// general assumptions ////////////////////////////////////////////////////////
190
191static proc dmodappAssumeViolation()
192{
193  // char K <> 0 or qring
194  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
195  {
196    ERROR("Basering is inappropriate: characteristic>0 or qring present");
197  }
198  return();
199}
200
201static proc dmodappMoreAssumeViolation()
202{
203  // char K <> 0, qring
204  dmodappAssumeViolation();
205  // no Weyl algebra
206  if (isWeyl() == 0)
207  {
208    ERROR("Basering is not a Weyl algebra");
209  }
210  // wrong sequence of vars
211  int i,n;
212  n = nvars(basering)/2;
213  for (i=1; i<=n; i++)
214  {
215    if (bracket(var(i+n),var(i))<>1)
216    {
217      ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i)));
218    }
219  }
220  return();
221}
222
223static proc safeVarName (string s, string cv)
224// assumes 's' to be a valid variable name
225// returns valid var name string @@..@s
226{
227  string S;
228  if (cv == "v")  { S = "," + "," + varstr(basering)  + ","; }
229  if (cv == "c")  { S = "," + "," + charstr(basering) + ","; }
230  if (cv == "cv") { S = "," + charstr(basering) + "," + varstr(basering) + ","; }
231  s = "," + s + ",";
232  while (find(S,s) <> 0)
233  {
234    s[1] = "@";
235    s = "," + s;
236  }
237  s = s[2..size(s)-1];
238  return(s)
239    }
240
241static proc intLike (def i)
242{
243  string str = typeof(i);
244  if (str == "int" || str == "number" || str == "bigint")
245  {
246    return(1);
247  }
248  else
249  {
250    return(0);
251  }
252}
253
254
255// toolbox ////////////////////////////////////////////////////////////////////
256
257proc engine(def I, int i)
258  "USAGE:  engine(I,i);  I  ideal/module/matrix, i an int
259RETURN:  the same type as I
260PURPOSE: compute the Groebner basis of I with the algorithm, chosen via i
261NOTE:    By default and if i=0, slimgb is used; otherwise std does the job.
262EXAMPLE: example engine; shows examples
263"
264{
265  /* std - slimgb mix */
266  def J;
267  //  ideal J;
268  if (i==0)
269  {
270    J = slimgb(I);
271  }
272  else
273  {
274    // without options -> strange! (ringlist?)
275    intvec v = option(get);
276    option(redSB);
277    option(redTail);
278    J = std(I);
279    option(set, v);
280  }
281  return(J);
282}
283example
284{
285  "EXAMPLE:"; echo = 2;
286  ring r = 0,(x,y),Dp;
287  ideal I  = y*(x3-y2),x*(x3-y2);
288  engine(I,0); // uses slimgb
289  engine(I,1); // uses std
290}
291
292proc poly2list (poly f)
293  "USAGE:  poly2list(f); f a poly
294RETURN:  list of exponents and corresponding terms of f
295PURPOSE: converts a poly to a list of pairs consisting of intvecs (1st entry)
296@*       and polys (2nd entry), where the i-th pair contains the exponent of the
297@*       i-th term of f and the i-th term (with coefficient) itself.
298EXAMPLE: example poly2list; shows examples
299"
300{
301  list l;
302  int i = 1;
303  if (f == 0) // just for the zero polynomial
304  {
305    l[1] = list(leadexp(f), lead(f));
306  }
307  else
308  {
309    l[size(f)] = list(); // memory pre-allocation
310    while (f != 0)
311    {
312      l[i] = list(leadexp(f), lead(f));
313      f = f - lead(f);
314      i++;
315    }
316  }
317  return(l);
318}
319example
320{
321  "EXAMPLE:"; echo = 2;
322  ring r = 0,x,dp;
323  poly F = x;
324  poly2list(F);
325  ring r2 = 0,(x,y),dp;
326  poly F = x2y+5xy2;
327  poly2list(F);
328  poly2list(0);
329}
330
331proc fl2poly(list L, string s)
332  "USAGE:  fl2poly(L,s); L a list, s a string
333RETURN:  poly
334PURPOSE: reconstruct a monic polynomial in one variable from its factorization
335ASSUME:  s is a string with the name of some variable and
336@*       L is supposed to consist of two entries:
337@*        - L[1] of the type ideal with the roots of a polynomial
338@*        - L[2] of the type intvec with the multiplicities of corr. roots
339EXAMPLE: example fl2poly; shows examples
340"
341{
342  if (varNum(s)==0)
343  {
344    ERROR(s+ " is no variable in the basering");
345  }
346  poly x = var(varNum(s));
347  poly P = 1;
348  ideal RR = L[1];
349  int sl = ncols(RR);
350  intvec IV = L[2];
351  if (sl <> nrows(IV))
352  {
353    ERROR("number of roots doesn't match number of multiplicites");
354  }
355  for(int i=1; i<=sl; i++)
356  {
357    if (IV[i] > 0)
358    {
359      P = P*((x-RR[i])^IV[i]);
360    }
361    else
362    {
363      printf("Ignored the root with incorrect multiplicity %s",string(IV[i]));
364    }
365  }
366  return(P);
367}
368example
369{
370  "EXAMPLE:"; echo = 2;
371  ring r = 0,(x,y,z,s),Dp;
372  ideal I = -1,-4/3,0,-5/3,-2;
373  intvec mI = 2,1,2,1,1;
374  list BS = I,mI;
375  poly p = fl2poly(BS,"s");
376  p;
377  factorize(p,2);
378}
379
380proc insertGenerator (list #)
381  "USAGE:  insertGenerator(id,p[,k]);
382@*       id an ideal/module, p a poly/vector, k an optional int
383RETURN:  of the same type as id
384PURPOSE: inserts p into id at k-th position and returns the enlarged object
385NOTE:    If k is given, p is inserted at position k, otherwise (and by default),
386@*       p is inserted at the beginning (k=1).
387SEE ALSO: deleteGenerator
388EXAMPLE: example insertGenerator; shows examples
389"
390{
391  if (size(#) < 2)
392  {
393    ERROR("insertGenerator has to be called with at least 2 arguments (ideal/module,poly/vector)");
394  }
395  string inp1 = typeof(#[1]);
396  if (inp1 == "ideal" || inp1 == "module")
397  {
398    def id = #[1];
399  }
400  else { ERROR("first argument has to be of type ideal or module"); }
401  string inp2 = typeof(#[2]);
402  if (inp2 == "poly" || inp2 == "vector")
403  {
404    def f = #[2];
405  }
406  else { ERROR("second argument has to be of type poly/vector"); }
407  if (inp1 == "ideal" && inp2 == "vector")
408  {
409    ERROR("second argument has to be a polynomial if first argument is an ideal");
410  }
411  // don't check module/poly combination due to auto-conversion
412  //   if (inp1 == "module" && inp2 == "poly")
413  //   {
414  //     ERROR("second argument has to be a vector if first argument is a module");
415  //   }
416  int n = ncols(id);
417  int k = 1; // default
418  if (size(#)>=3)
419  {
420    if (intLike(#[3]))
421    {
422      k = int(#[3]);
423      if (k<=0)
424      {
425        ERROR("third argument has to be positive");
426      }
427    }
428    else { ERROR("third argument has to be of type int"); }
429  }
430  execute(inp1 +" J;");
431  if (k == 1) { J = f,id; }
432  else
433  {
434    if (k>n)
435    {
436      J = id;
437      J[k] = f;
438    }
439    else // 1<k<=n
440    {
441      J[n+1] = id[n]; // preinit
442      J[1..k-1] = id[1..k-1];
443      J[k] = f;
444      J[k+1..n+1] = id[k..n];
445    }
446  }
447  return(J);
448}
449example
450{
451  "EXAMPLE:"; echo = 2;
452  ring r = 0,(x,y,z),dp;
453  ideal I = x^2,z^4;
454  insertGenerator(I,y^3);
455  insertGenerator(I,y^3,2);
456  module M = I*gen(3);
457  insertGenerator(M,[x^3,y^2,z],2);
458  insertGenerator(M,x+y+z,4);
459}
460
461proc deleteGenerator (def id, int k)
462  "USAGE:   deleteGenerator(id,k);  id an ideal/module, k an int
463RETURN:   of the same type as id
464PURPOSE:  deletes the k-th generator from the first argument and returns
465@*        the altered object
466SEE ALSO: insertGenerator
467EXAMPLE:  example deleteGenerator; shows examples
468"
469{
470  string inp1 = typeof(id);
471  if (inp1 <> "ideal" && inp1 <> "module")
472  {
473    ERROR("first argument has to be of type ideal or module");
474  }
475  execute(inp1 +" J;");
476  int n = ncols(id);
477  if (n == 1 && k == 1) { return(J); }
478  if (k<=0 || k>n)
479  {
480    ERROR("second argument has to be in the range 1,...,"+string(n));
481  }
482  J[n-1] = 0; // preinit
483  if (k == 1) { J = id[2..n]; }
484  else
485  {
486    if (k == n) { J = id[1..n-1]; }
487    else
488    {
489      J[1..k-1] = id[1..k-1];
490      J[k..n-1] = id[k+1..n];
491    }
492  }
493  return(J);
494}
495example
496{
497  "EXAMPLE:"; echo = 2;
498  ring r = 0,(x,y,z),dp;
499  ideal I = x^2,y^3,z^4;
500  deleteGenerator(I,2);
501  module M = [x,y,z],[x2,y2,z2],[x3,y3,z3];
502  print(deleteGenerator(M,2));
503  M = M[1];
504  deleteGenerator(M,1);
505}
506
507static proc sortNumberIdeal (ideal I)
508// sorts ideal of constant polys (ie numbers), returns according permutation
509{
510  int i;
511  int nI = ncols(I);
512  intvec dI;
513  for (i=nI; i>0; i--)
514  {
515    dI[i] = int(denominator(leadcoef(I[i])));
516  }
517  int lcmI = lcm(dI);
518  for (i=nI; i>0; i--)
519  {
520    dI[i] = int(lcmI*leadcoef(I[i]));
521  }
522  intvec perm = sortIntvec(dI)[2];
523  return(perm);
524}
525example
526{
527  "EXAMPLE:"; echo = 2;
528  ring r = 0,s,dp;
529  ideal I = -9/20,-11/20,-23/20,-19/20,-1,-13/10,-27/20,-13/20,-21/20,-17/20,
530    -11/10,-9/10,-7/10; // roots of BS poly of reiffen(4,5)
531  intvec v = sortNumberIdeal(I); v;
532  I[v];
533}
534
535proc bFactor (poly F)
536  "USAGE:  bFactor(f);  f poly
537RETURN:  list of ideal and intvec and possibly a string
538PURPOSE: tries to compute the roots of a univariate poly f
539NOTE:    The output list consists of two or three entries:
540@*       roots of f as an ideal, their multiplicities as intvec, and,
541@*       if present, a third one being the product of all irreducible factors
542@*       of degree greater than one, given as string.
543@*       If f is the zero polynomial or if f has no roots in the ground field,
544@*       this is encoded as root 0 with multiplicity 0.
545DISPLAY: If printlevel=1, progress debug messages will be printed,
546@*       if printlevel>=2, all the debug messages will be printed.
547EXAMPLE: example bFactor; shows examples
548"
549{
550  int ppl = printlevel - voice +2;
551  def save = basering;
552  ideal LI = variables(F);
553  list L;
554  int i = size(LI);
555  if (i>1) { ERROR("poly has to be univariate")}
556  if (i == 0)
557  {
558    dbprint(ppl,"// poly is constant");
559    L = list(ideal(0),intvec(0),string(F));
560    return(L);
561  }
562  poly v = LI[1];
563  L = ringlist(save); L = L[1..4];
564  L[2] = list(string(v));
565  L[3] = list(list("dp",intvec(1)),list("C",intvec(0)));
566  def @S = ring(L);
567  setring @S;
568  poly ir = 1; poly F = imap(save,F);
569  list L = factorize(F);
570  ideal I = L[1];
571  intvec m = L[2];
572  ideal II; intvec mm;
573  for (i=2; i<=ncols(I); i++)
574  {
575    if (deg(I[i]) > 1)
576    {
577      ir = ir * I[i]^m[i];
578    }
579    else
580    {
581      II[size(II)+1] = I[i];
582      mm[size(II)] = m[i];
583    }
584  }
585  II = normalize(II);
586  II = subst(II,var(1),0);
587  II = -II;
588  intvec perm = sortNumberIdeal(II);
589  II = II[perm];
590  mm = mm[perm];
591  if (size(II)>0)
592  {
593    dbprint(ppl,"// found roots");
594    dbprint(ppl-1,string(II));
595  }
596  else
597  {
598    dbprint(ppl,"// found no roots");
599  }
600  L = list(II,mm);
601  if (ir <> 1)
602  {
603    dbprint(ppl,"// found irreducible factors");
604    ir = cleardenom(ir);
605    dbprint(ppl-1,string(ir));
606    L = L + list(string(ir));
607  }
608  else
609  {
610    dbprint(ppl,"// no irreducible factors found");
611  }
612  setring save;
613  L = imap(@S,L);
614  return(L);
615}
616example
617{
618  "EXAMPLE:"; echo = 2;
619  ring r = 0,(x,y),dp;
620  bFactor((x^2-1)^2);
621  bFactor((x^2+1)^2);
622  bFactor((y^2+1/2)*(y+9)*(y-7));
623  bFactor(1);
624  bFactor(0);
625}
626
627proc isInt (number n)
628  "USAGE:  isInt(n); n a number
629RETURN:  int, 1 if n is an integer or 0 otherwise
630PURPOSE: check whether given object of type number is actually an int
631NOTE:    Parameters are treated as integers.
632EXAMPLE: example isInt; shows an example
633"
634{
635  number d = denominator(n);
636  if (d<>1)
637  {
638    return(0);
639  }
640  else
641  {
642    return(1);
643  }
644}
645example
646{
647  "EXAMPLE:"; echo = 2;
648  ring r = 0,x,dp;
649  number n = 4/3;
650  isInt(n);
651  n = 11;
652  isInt(n);
653}
654
655proc intRoots (list l)
656  "USAGE:  isInt(L); L a list
657RETURN:  list
658PURPOSE: extracts integer roots from a list given in @code{bFactor} format
659ASSUME:  The input list must be given in the format of @code{bFactor}.
660NOTE:    Parameters are treated as integers.
661SEE ALSO: bFactor
662EXAMPLE: example intRoots; shows an example
663"
664{
665  int wronginput;
666  int sl = size(l);
667  if (sl>0)
668  {
669    if (typeof(l[1])<>"ideal"){wronginput = 1;}
670    if (sl>1)
671    {
672      if (typeof(l[2])<>"intvec"){wronginput = 1;}
673      if (sl>2)
674      {
675        if (typeof(l[3])<>"string"){wronginput = 1;}
676        if (sl>3){wronginput = 1;}
677      }
678    }
679  }
680  if (sl<2){wronginput = 1;}
681  if (wronginput)
682  {
683    ERROR("Given list has wrong format.");
684  }
685  int i,j;
686  ideal l1 = l[1];
687  int n = ncols(l1);
688  j = 1;
689  ideal I;
690  intvec v;
691  for (i=1; i<=n; i++)
692  {
693    if (size(l1[j])>1) // poly not number
694    {
695      ERROR("Ideal in list has wrong format.");
696    }
697    if (isInt(leadcoef(l1[i])))
698    {
699      I[j] = l1[i];
700      v[j] = l[2][i];
701      j++;
702    }
703  }
704  return(list(I,v));
705}
706example
707{
708  "EXAMPLE:"; echo = 2;
709  ring r = 0,x,dp;
710  list L = bFactor((x-4/3)*(x+3)^2*(x-5)^4); L;
711  intRoots(L);
712}
713
714proc sortIntvec (intvec v)
715  "USAGE:  sortIntvec(v); v an intvec
716RETURN:  list of two intvecs
717PURPOSE: sorts an intvec
718NOTE:    In the output list L, the first entry consists of the entries of v
719@*       satisfying L[1][i] >= L[1][i+1]. The second entry is a permutation
720@*       such that v[L[2]] = L[1].
721@*       Unlike in the procedure @code{sort}, zeros are not dismissed.
722SEE ALSO: sort
723EXAMPLE: example sortIntvec; shows examples
724"
725{
726  int i;
727  intvec vpos,vzero,vneg,vv,sortv,permv;
728  list l;
729  for (i=1; i<=nrows(v); i++)
730  {
731    if (v[i]>0)
732    {
733      vpos = vpos,i;
734    }
735    else
736    {
737      if (v[i]==0)
738      {
739        vzero = vzero,i;
740      }
741      else // v[i]<0
742      {
743        vneg = vneg,i;
744      }
745    }
746  }
747  if (size(vpos)>1)
748  {
749    vpos = vpos[2..size(vpos)];
750    vv = v[vpos];
751    l = sort(vv);
752    vv = l[1];
753    vpos = vpos[l[2]];
754    sortv = vv[size(vv)..1];
755    permv = vpos[size(vv)..1];
756  }
757  if (size(vzero)>1)
758  {
759    vzero = vzero[2..size(vzero)];
760    permv = permv,vzero;
761    sortv = sortv,0:size(vzero);
762  }
763  if (size(vneg)>1)
764  {
765    vneg = vneg[2..size(vneg)];
766    vv = -v[vneg];
767    l = sort(vv);
768    vv = -l[1];
769    vneg = vneg[l[2]];
770    sortv = sortv,vv;
771    permv = permv,vneg;
772  }
773  if (permv[1]==0)
774  {
775    sortv = sortv[2..size(sortv)];
776    permv = permv[2..size(permv)];
777  }
778  return(list(sortv,permv));
779}
780example
781{
782  "EXAMPLE:"; echo = 2;
783  ring r = 0,x,dp;
784  intvec v = -1,0,1,-2,0,2;
785  list L = sortIntvec(v); L;
786  v[L[2]];
787  v = -3,0;
788  sortIntvec(v);
789  v = 0,-3;
790  sortIntvec(v);
791}
792
793
794// F-saturation ///////////////////////////////////////////////////////////////
795
796proc isFsat(ideal I, poly F)
797  "USAGE:  isFsat(I, F);  I an ideal, F a poly
798RETURN:  int, 1  if I is F-saturated and 0 otherwise
799PURPOSE: checks whether the ideal I is F-saturated
800NOTE:    We check indeed that Ker(D--> F--> D/I) is 0, where D is the basering.
801EXAMPLE: example isFsat; shows examples
802"
803{
804  /* checks whether I is F-saturated, that is Ke  (D -F-> D/I) is 0 */
805  /* works in any algebra */
806  /*  for simplicity : later check attrib */
807  /* returns 1 if I is F-sat */
808  if (attrib(I,"isSB")!=1)
809  {
810    I = groebner(I);
811  }
812  matrix @M = matrix(I);
813  matrix @F[1][1] = F;
814  def S = modulo(module(@F),module(@M));
815  S = NF(S,I);
816  S = groebner(S);
817  return( (gkdim(S) == -1) );
818}
819example
820{
821  "EXAMPLE:"; echo = 2;
822  ring r = 0,(x,y),dp;
823  poly G = x*(x-y)*y;
824  def A = annfs(G);
825  setring A;
826  poly F = x3-y2;
827  isFsat(LD,F);
828  ideal J = LD*F;
829  isFsat(J,F);
830}
831
832
833// annihilators ///////////////////////////////////////////////////////////////
834
835proc annRat(poly g, poly f)
836  "USAGE:   annRat(g,f);  f, g polynomials
837RETURN:   ring (a Weyl algebra) containing an ideal 'LD'
838PURPOSE:  compute the annihilator of the rational function g/f in the
839@*        corresponding Weyl algebra
840ASSUME:   basering is commutative and over a field of characteristic 0
841NOTE:     Activate the output ring with the @code{setring} command.
842@*        In the output ring, the ideal 'LD' (in Groebner basis) is the
843@*        annihilator of g/f.
844@*        The algorithm uses the computation of Ann(f^{-1}) via D-modules,
845@*        see (SST).
846DISPLAY:  If printlevel =1, progress debug messages will be printed,
847@*        if printlevel>=2, all the debug messages will be printed.
848SEE ALSO: annPoly
849EXAMPLE:  example annRat; shows examples
850"
851{
852  // assumption check
853  dmodappAssumeViolation();
854  if (!isCommutative())
855  {
856    ERROR("Basering must be commutative.");
857  }
858  // assumptions: f is not a constant
859  if (f==0) { ERROR("the denominator f cannot be zero"); }
860  if ((leadexp(f) == 0) && (size(f) < 2))
861  {
862    // f = const, so use annPoly
863    g = g/f;
864    def @R = annPoly(g);
865    return(@R);
866  }
867  // computes the annihilator of g/f
868  def save = basering;
869  int ppl = printlevel-voice+2;
870  dbprint(ppl,"// -1-[annRat] computing the ann f^s");
871  def  @R1 = SannfsBM(f);
872  setring @R1;
873  poly f = imap(save,f);
874  int i,mir;
875  int isr = 0; // checkRoot1(LD,f,1); // roots are negative, have to enter positive int
876  if (!isr)
877  {
878    // -1 is not the root
879    // find the m.i.r iteratively
880    mir = 0;
881    for(i=nvars(save)+1; i>=1; i--)
882    {
883      isr =  checkRoot1(LD,f,i);
884      if (isr) { mir =-i; break; }
885    }
886    if (mir ==0)
887    {
888      ERROR("No integer root found! Aborting computations, inform the authors!");
889    }
890    // now mir == i is m.i.r.
891  }
892  else
893  {
894    // -1 is the m.i.r
895    mir = -1;
896  }
897  dbprint(ppl,"// -2-[annRat] the minimal integer root is ");
898  dbprint(ppl-1, mir);
899  // use annfspecial
900  dbprint(ppl,"// -3-[annRat] running annfspecial ");
901  ideal AF = annfspecial(LD,f,mir,-1); // ann f^{-1}
902  //  LD = subst(LD,s,j);
903  //  LD = engine(LD,0);
904  // modify the ring: throw s away
905  // output ring comes from SannfsBM
906  list U = ringlist(@R1);
907  list tmp; // variables
908  for(i=1; i<=size(U[2])-1; i++)
909  {
910    tmp[i] = U[2][i];
911  }
912  U[2] = tmp;
913  tmp = 0;
914  tmp[1] = U[3][1]; // x,Dx block
915  tmp[2] = U[3][3]; // module block
916  U[3] = tmp;
917  tmp = 0;
918  tmp = U[1],U[2],U[3],U[4];
919  def @R2 = ring(tmp);
920  setring @R2;
921  // now supply with Weyl algebra relations
922  int N = nvars(@R2)/2;
923  matrix @D[2*N][2*N];
924  for(i=1; i<=N; i++)
925  {
926    @D[i,N+i]=1;
927  }
928  def @R3 = nc_algebra(1,@D);
929  setring @R3;
930  dbprint(ppl,"// - -[annRat] ring without s is ready:");
931  dbprint(ppl-1,@R3);
932  poly g = imap(save,g);
933  matrix G[1][1] = g;
934  matrix LL = matrix(imap(@R1,AF));
935  kill @R1;   kill @R2;
936  dbprint(ppl,"// -4-[annRat] running modulo");
937  ideal LD  = modulo(G,LL);
938  dbprint(ppl,"// -4-[annRat] running GB on the final result");
939  LD  = engine(LD,0);
940  export LD;
941  return(@R3);
942}
943example
944{
945  "EXAMPLE:"; echo = 2;
946  ring r = 0,(x,y),dp;
947  poly g = 2*x*y;  poly f = x^2 - y^3;
948  def B = annRat(g,f);
949  setring B;
950  LD;
951  // Now, compare with the output of Macaulay2:
952  ideal tst = 3*x*Dx + 2*y*Dy + 1, y^3*Dy^2 - x^2*Dy^2 + 6*y^2*Dy + 6*y,
953    9*y^2*Dx^2*Dy-4*y*Dy^3+27*y*Dx^2+2*Dy^2, 9*y^3*Dx^2-4*y^2*Dy^2+10*y*Dy -10;
954  option(redSB); option(redTail);
955  LD = groebner(LD);
956  tst = groebner(tst);
957  print(matrix(NF(LD,tst)));  print(matrix(NF(tst,LD)));
958  // So, these two answers are the same
959}
960
961proc annPoly(poly f)
962  "USAGE:   annPoly(f);  f a poly
963RETURN:   ring (a Weyl algebra) containing an ideal 'LD'
964PURPOSE:  compute the complete annihilator ideal of f in the corresponding
965@*        Weyl algebra
966ASSUME:   basering is commutative and over a field of characteristic 0
967NOTE:     Activate the output ring with the @code{setring} command.
968@*        In the output ring, the ideal 'LD' (in Groebner basis) is the
969@*        annihilator.
970DISPLAY:  If printlevel =1, progress debug messages will be printed,
971@*        if printlevel>=2, all the debug messages will be printed.
972SEE ALSO: annRat
973EXAMPLE:  example annPoly; shows examples
974"
975{
976  // assumption check
977  dmodappAssumeViolation();
978  if (!isCommutative())
979  {
980    ERROR("Basering must be commutative.");
981  }
982  // computes a system of linear PDEs with polynomial coeffs for f
983  def save = basering;
984  list L = ringlist(save);
985  list Name = L[2];
986  int N = nvars(save);
987  int i;
988  for (i=1; i<=N; i++)
989  {
990    Name[N+i] = safeVarName("D"+Name[i],"cv"); // concat
991  }
992  L[2] = Name;
993  def @R = ring(L);
994  setring @R;
995  def @@R = Weyl();
996  setring @@R;
997  kill @R;
998  matrix M[1][N];
999  for (i=1; i<=N; i++)
1000  {
1001    M[1,i] = var(N+i);
1002  }
1003  matrix F[1][1] = imap(save,f);
1004  def I = modulo(module(F),module(M));
1005  ideal LD = I;
1006  LD = groebner(LD);
1007  export LD;
1008  return(@@R);
1009}
1010example
1011{
1012  "EXAMPLE:"; echo = 2;
1013  ring r = 0,(x,y,z),dp;
1014  poly f = x^2*z - y^3;
1015  def A = annPoly(f);
1016  setring A;    // A is the 3rd Weyl algebra in 6 variables
1017  LD;           // the Groebner basis of annihilator
1018  gkdim(LD);    // must be 3 = 6/2, since A/LD is holonomic module
1019  NF(Dy^4, LD); // must be 0 since Dy^4 clearly annihilates f
1020  poly f = imap(r,f);
1021  NF(LD*f,std(ideal(Dx,Dy,Dz))); // must be zero if LD indeed annihilates f
1022}
1023
1024
1025
1026// localizations //////////////////////////////////////////////////////////////
1027
1028proc DLoc(ideal I, poly F)
1029  "USAGE:  DLoc(I, f);  I an ideal, f a poly
1030RETURN:  list of ideal and list
1031ASSUME:  the basering is a Weyl algebra
1032PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s
1033NOTE:    In the output list L,
1034@*       - L[1] is an ideal (given as Groebner basis), the presentation of the
1035@*       localization,
1036@*       - L[2] is a list containing roots with multiplicities of Bernstein
1037@*       polynomial of (D/I)_f.
1038DISPLAY: If printlevel =1, progress debug messages will be printed,
1039@*       if printlevel>=2, all the debug messages will be printed.
1040EXAMPLE: example DLoc; shows examples
1041"
1042{
1043  /* runs SDLoc and DLoc0 */
1044  /* assume: run from Weyl algebra */
1045  dmodappAssumeViolation();
1046  if (!isWeyl())
1047  {
1048    ERROR("Basering is not a Weyl algebra");
1049  }
1050  int old_printlevel = printlevel;
1051  printlevel=printlevel+1;
1052  def @R = basering;
1053  def @R2 = SDLoc(I,F);
1054  setring @R2;
1055  poly F = imap(@R,F);
1056  def @R3 = DLoc0(LD,F);
1057  setring @R3;
1058  ideal bs = BS[1];
1059  intvec m = BS[2];
1060  setring @R;
1061  ideal LD0 = imap(@R3,LD0);
1062  ideal bs = imap(@R3,bs);
1063  list BS; BS[1] = bs; BS[2] = m;
1064  kill @R3;
1065  printlevel = old_printlevel;
1066  return(list(LD0,BS));
1067}
1068example;
1069{
1070  "EXAMPLE:"; echo = 2;
1071  ring r = 0,(x,y,Dx,Dy),dp;
1072  def R = Weyl();    setring R; // Weyl algebra in variables x,y,Dx,Dy
1073  poly F = x2-y3;
1074  ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
1075  // I is not holonomic, since its dimension is not 4/2=2
1076  gkdim(I);
1077  list L = DLoc(I, x2-y3);
1078  L[1]; // localized module (R/I)_f is isomorphic to R/LD0
1079  L[2]; // description of b-function for localization
1080}
1081
1082proc DLoc0(ideal I, poly F)
1083  "USAGE:  DLoc0(I, f);  I an ideal, f a poly
1084RETURN:  ring (a Weyl algebra) containing an ideal 'LD0' and a list 'BS'
1085PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s,
1086@*       where D is a Weyl Algebra, based on the output of procedure SDLoc
1087ASSUME:  the basering is similar to the output ring of SDLoc procedure
1088NOTE:    activate the output ring with the @code{setring} command. In this ring,
1089@*       - the ideal LD0 (given as Groebner basis) is the presentation of the
1090@*       localization,
1091@*       - the list BS contains roots and multiplicities of Bernstein
1092@*       polynomial of (D/I)_f.
1093DISPLAY: If printlevel =1, progress debug messages will be printed,
1094@*       if printlevel>=2, all the debug messages will be printed.
1095EXAMPLE: example DLoc0; shows examples
1096"
1097{
1098  dmodappAssumeViolation();
1099  /* assume: to be run in the output ring of SDLoc */
1100  /* doing: add F, eliminate vars*Dvars, factorize BS */
1101  /* analogue to annfs0 */
1102  def @R2 = basering;
1103  // we're in D_n[s], where the elim ord for s is set
1104  ideal J = NF(I,std(F));
1105  // make leadcoeffs positive
1106  int i;
1107  for (i=1; i<= ncols(J); i++)
1108  {
1109    if (leadcoef(J[i]) <0 )
1110    {
1111      J[i] = -J[i];
1112    }
1113  }
1114  J = J,F;
1115  ideal M = groebner(J);
1116  int Nnew = nvars(@R2);
1117  ideal K2 = nselect(M,1..Nnew-1);
1118  int ppl = printlevel-voice+2;
1119  dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering");
1120  dbprint(ppl-1, K2);
1121  // the ring @R3 and the search for minimal negative int s
1122  ring @R3 = 0,s,dp;
1123  dbprint(ppl,"// -2-1- the ring @R3 = K[s] is ready");
1124  ideal K3 = imap(@R2,K2);
1125  poly p = K3[1];
1126  dbprint(ppl,"// -2-2- attempt the factorization");
1127  list PP = factorize(p);          //with constants and multiplicities
1128  ideal bs; intvec m;             //the Bernstein polynomial is monic, so
1129  // we are not interested in constants
1130  for (i=2; i<= size(PP[1]); i++)  //we delete P[1][1] and P[2][1]
1131  {
1132    bs[i-1] = PP[1][i];
1133    m[i-1]  = PP[2][i];
1134  }
1135  ideal bbs; int srat=0; int HasRatRoots = 0;
1136  int sP;
1137  for (i=1; i<= size(bs); i++)
1138  {
1139    if (deg(bs[i]) == 1)
1140    {
1141      bbs = bbs,bs[i];
1142    }
1143  }
1144  if (size(bbs)==0)
1145  {
1146    dbprint(ppl-1,"// -2-3- factorization: no rational roots");
1147    //    HasRatRoots = 0;
1148    HasRatRoots = 1; // s0 = -1 then
1149    sP = -1;
1150    // todo: return ideal with no subst and a b-function unfactorized
1151  }
1152  else
1153  {
1154    // exist rational roots
1155    dbprint(ppl-1,"// -2-3- factorization: rational roots found");
1156    HasRatRoots = 1;
1157    //    dbprint(ppl-1,bbs);
1158    bbs = bbs[2..ncols(bbs)];
1159    ideal P = bbs;
1160    dbprint(ppl-1,P);
1161    srat = size(bs) - size(bbs);
1162    // define minIntRoot on linear factors or find out that it doesn't exist
1163    intvec vP;
1164    number nP;
1165    P = normalize(P); // now leadcoef = 1
1166    P = ideal(matrix(lead(P))-matrix(P));
1167    sP = size(P);
1168    int cnt = 0;
1169    for (i=1; i<=sP; i++)
1170    {
1171      nP = leadcoef(P[i]);
1172      if ( (nP - int(nP)) == 0 )
1173      {
1174        cnt++;
1175        vP[cnt] = int(nP);
1176      }
1177    }
1178    //     if ( size(vP)>=2 )
1179    //     {
1180    //       vP = vP[2..size(vP)];
1181    //     }
1182    if ( size(vP)==0 )
1183    {
1184      // no roots!
1185      dbprint(ppl,"// -2-4- no integer root, setting s0 = -1");
1186      sP = -1;
1187      //      HasRatRoots = 0; // older stuff, here we do substitution
1188      HasRatRoots = 1;
1189    }
1190    else
1191    {
1192      HasRatRoots = 1;
1193      sP = -Max(-vP);
1194      dbprint(ppl,"// -2-4- minimal integer root found");
1195      dbprint(ppl-1, sP);
1196      //    int sP = minIntRoot(bbs,1);
1197      //       P =  normalize(P);
1198      //       bs = -subst(bs,s,0);
1199      if (sP >=0)
1200      {
1201        dbprint(ppl,"// -2-5- nonnegative root, setting s0 = -1");
1202        sP = -1;
1203      }
1204      else
1205      {
1206        dbprint(ppl,"// -2-5- the root is negative");
1207      }
1208    }
1209  }
1210
1211  if (HasRatRoots)
1212  {
1213    setring @R2;
1214    K2 = subst(I,s,sP);
1215    // IF min int root exists ->
1216    // create the ordinary Weyl algebra and put the result into it,
1217    // thus creating the ring @R5
1218    // ELSE : return the same ring with new objects
1219    // keep: N, i,j,s, tmp, RL
1220    Nnew = Nnew - 1; // former 2*N;
1221    // list RL = ringlist(save);  // is defined earlier
1222    //  kill Lord, tmp, iv;
1223    list L = 0;
1224    list Lord, tmp;
1225    intvec iv;
1226    list RL = ringlist(basering);
1227    L[1] = RL[1];
1228    L[4] = RL[4];  //char, minpoly
1229    // check whether vars have admissible names -> done earlier
1230    // list Name = RL[2]M
1231    // DName is defined earlier
1232    list NName; // = RL[2]; // skip the last var 's'
1233    for (i=1; i<=Nnew; i++)
1234    {
1235      NName[i] =  RL[2][i];
1236    }
1237    L[2] = NName;
1238    // dp ordering;
1239    string s = "iv=";
1240    for (i=1; i<=Nnew; i++)
1241    {
1242      s = s+"1,";
1243    }
1244    s[size(s)] = ";";
1245    execute(s);
1246    tmp     = 0;
1247    tmp[1]  = "dp";  // string
1248    tmp[2]  = iv;  // intvec
1249    Lord[1] = tmp;
1250    kill s;
1251    tmp[1]  = "C";
1252    iv = 0;
1253    tmp[2]  = iv;
1254    Lord[2] = tmp;
1255    tmp     = 0;
1256    L[3]    = Lord;
1257    // we are done with the list
1258    // Add: Plural part
1259    def @R4@ = ring(L);
1260    setring @R4@;
1261    int N = Nnew/2;
1262    matrix @D[Nnew][Nnew];
1263    for (i=1; i<=N; i++)
1264    {
1265      @D[i,N+i]=1;
1266    }
1267    def @R4 = nc_algebra(1,@D);
1268    setring @R4;
1269    kill @R4@;
1270    dbprint(ppl,"// -3-1- the ring @R4 is ready");
1271    dbprint(ppl-1, @R4);
1272    ideal K4 = imap(@R2,K2);
1273    intvec vopt = option(get);
1274    option(redSB);
1275    dbprint(ppl,"// -3-2- the final cosmetic std");
1276    K4 = groebner(K4);  // std does the job too
1277    option(set,vopt);
1278    // total cleanup
1279    setring @R2;
1280    ideal bs = imap(@R3,bs);
1281    bs = -normalize(bs); // "-" for getting correct coeffs!
1282    bs = subst(bs,s,0);
1283    kill @R3;
1284    setring @R4;
1285    ideal bs = imap(@R2,bs); // only rationals are the entries
1286    list BS; BS[1] = bs; BS[2] = m;
1287    export BS;
1288    //    list LBS = imap(@R3,LBS);
1289    //    list BS; BS[1] = sbs; BS[2] = m;
1290    //    BS;
1291    //    export BS;
1292    ideal LD0 = K4;
1293    export LD0;
1294    return(@R4);
1295  }
1296  else
1297  {
1298    /* SHOULD NEVER GET THERE */
1299    /* no rational/integer roots */
1300    /* return objects in the copy of current ring */
1301    setring @R2;
1302    ideal LD0 = I;
1303    poly BS = normalize(K2[1]);
1304    export LD0;
1305    export BS;
1306    return(@R2);
1307  }
1308}
1309example;
1310{
1311  "EXAMPLE:"; echo = 2;
1312  ring r = 0,(x,y,Dx,Dy),dp;
1313  def R = Weyl();    setring R; // Weyl algebra in variables x,y,Dx,Dy
1314  poly F = x2-y3;
1315  ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
1316  // moreover I is not holonomic, since its dimension is not 2 = 4/2
1317  gkdim(I); // 3
1318  def W = SDLoc(I,F);  setring W; // creates ideal LD in W = R[s]
1319  def U = DLoc0(LD, x2-y3);  setring U; // compute in R
1320  LD0; // Groebner basis of the presentation of localization
1321  BS; // description of b-function for localization
1322}
1323
1324proc SDLoc(ideal I, poly F)
1325  "USAGE:  SDLoc(I, f);  I an ideal, f a poly
1326RETURN:  ring (basering extended by a new variable) containing an ideal 'LD'
1327PURPOSE: compute a generic presentation of the localization of D/I w.r.t. f^s
1328ASSUME:  the basering D is a Weyl algebra over a field of characteristic 0
1329NOTE:    Activate this ring with the @code{setring} command. In this ring,
1330@*       the ideal LD (given as Groebner basis) is the presentation of the
1331@*       localization.
1332DISPLAY: If printlevel =1, progress debug messages will be printed,
1333@*       if printlevel>=2, all the debug messages will be printed.
1334EXAMPLE: example SDLoc; shows examples
1335"
1336{
1337  /* analogue to Sannfs */
1338  /* printlevel >=4 gives debug info */
1339  /* assume: we're in the Weyl algebra D  in x1,x2,...,d1,d2,... */
1340
1341  dmodappAssumeViolation();
1342  if (!isWeyl())
1343  {
1344    ERROR("Basering is not a Weyl algebra");
1345  }
1346  def save = basering;
1347  /* 1. create D <t, dt, s > as in LOT */
1348  /* ordering: eliminate t,dt */
1349  int ppl = printlevel-voice+2;
1350  int N = nvars(save); N = N div 2;
1351  int Nnew = 2*N + 3; // t,Dt,s
1352  int i,j;
1353  string s;
1354  list RL = ringlist(save);
1355  list L, Lord;
1356  list tmp;
1357  intvec iv;
1358  L[1] = RL[1]; // char
1359  L[4] = RL[4]; // char, minpoly
1360  // check whether vars have admissible names
1361  list Name  = RL[2];
1362  list RName;
1363  RName[1] = "@t";
1364  RName[2] = "@Dt";
1365  RName[3] = "@s";
1366  for(i=1;i<=N;i++)
1367  {
1368    for(j=1; j<=size(RName);j++)
1369    {
1370      if (Name[i] == RName[j])
1371      {
1372        ERROR("Variable names should not include @t,@Dt,@s");
1373      }
1374    }
1375  }
1376  // now, create the names for new vars
1377  tmp    =  0;
1378  tmp[1] = "@t";
1379  tmp[2] = "@Dt";
1380  list SName ; SName[1] = "@s";
1381  list NName = tmp + Name + SName;
1382  L[2]   = NName;
1383  tmp    = 0;
1384  kill NName;
1385  // block ord (a(1,1),dp);
1386  tmp[1]  = "a"; // string
1387  iv      = 1,1;
1388  tmp[2]  = iv; //intvec
1389  Lord[1] = tmp;
1390  // continue with dp 1,1,1,1...
1391  tmp[1]  = "dp"; // string
1392  s       = "iv=";
1393  for(i=1;i<=Nnew;i++)
1394  {
1395    s = s+"1,";
1396  }
1397  s[size(s)]= ";";
1398  execute(s);
1399  tmp[2]    = iv;
1400  Lord[2]   = tmp;
1401  tmp[1]    = "C";
1402  iv        = 0;
1403  tmp[2]    = iv;
1404  Lord[3]   = tmp;
1405  tmp       = 0;
1406  L[3]      = Lord;
1407  // we are done with the list
1408  def @R@ = ring(L);
1409  setring @R@;
1410  matrix @D[Nnew][Nnew];
1411  @D[1,2]=1;
1412  for(i=1; i<=N; i++)
1413  {
1414    @D[2+i,N+2+i]=1;
1415  }
1416  // ADD [s,t]=-t, [s,Dt]=Dt
1417  @D[1,Nnew] = -var(1);
1418  @D[2,Nnew] = var(2);
1419  def @R = nc_algebra(1,@D);
1420  setring @R;
1421  kill @R@;
1422  dbprint(ppl,"// -1-1- the ring @R(@t,@Dt,_x,_Dx,@s) is ready");
1423  dbprint(ppl-1, @R);
1424  poly  F = imap(save,F);
1425  ideal I = imap(save,I);
1426  dbprint(ppl-1, "the ideal after map:");
1427  dbprint(ppl-1, I);
1428  poly p = 0;
1429  for(i=1; i<=N; i++)
1430  {
1431    p = diff(F,var(2+i))*@Dt + var(2+N+i);
1432    dbprint(ppl-1, p);
1433    I = subst(I,var(2+N+i),p);
1434    dbprint(ppl-1, var(2+N+i));
1435    p = 0;
1436  }
1437  I = I, @t - F;
1438  // t*Dt + s +1 reduced with t-f gives f*Dt + s
1439  I = I, F*var(2) + var(Nnew); // @s
1440  // -------- the ideal I is ready ----------
1441  dbprint(ppl,"// -1-2- starting the elimination of @t,@Dt in @R");
1442  dbprint(ppl-1, I);
1443  //  ideal J = engine(I,eng);
1444  ideal J = groebner(I);
1445  dbprint(ppl-1,"// -1-2-1- result of the  elimination of @t,@Dt in @R");
1446  dbprint(ppl-1, J);;
1447  ideal K = nselect(J,1..2);
1448  dbprint(ppl,"// -1-3- @t,@Dt are eliminated");
1449  dbprint(ppl-1, K);  // K is without t, Dt
1450  K = groebner(K);  // std does the job too
1451  // now, we must change the ordering
1452  // and create a ring without t, Dt
1453  setring save;
1454  // ----------- the ring @R3 ------------
1455  // _x, _Dx,s;  elim.ord for _x,_Dx.
1456  // keep: N, i,j,s, tmp, RL
1457  Nnew = 2*N+1;
1458  kill Lord, tmp, iv, RName;
1459  list Lord, tmp;
1460  intvec iv;
1461  L[1] = RL[1];
1462  L[4] = RL[4];  // char, minpoly
1463  // check whether vars hava admissible names -> done earlier
1464  // now, create the names for new var
1465  tmp[1] = "s";
1466  list NName = Name + tmp;
1467  L[2] = NName;
1468  tmp = 0;
1469  // block ord (dp(N),dp);
1470  // string s is already defined
1471  s = "iv=";
1472  for (i=1; i<=Nnew-1; i++)
1473  {
1474    s = s+"1,";
1475  }
1476  s[size(s)]=";";
1477  execute(s);
1478  tmp[1] = "dp";  // string
1479  tmp[2] = iv;   // intvec
1480  Lord[1] = tmp;
1481  // continue with dp 1,1,1,1...
1482  tmp[1] = "dp";  // string
1483  s[size(s)] = ",";
1484  s = s+"1;";
1485  execute(s);
1486  kill s;
1487  kill NName;
1488  tmp[2]      = iv;
1489  Lord[2]     = tmp;
1490  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
1491  Lord[3]     = tmp;  tmp = 0;
1492  L[3]        = Lord;
1493  // we are done with the list. Now add a Plural part
1494  def @R2@ = ring(L);
1495  setring @R2@;
1496  matrix @D[Nnew][Nnew];
1497  for (i=1; i<=N; i++)
1498  {
1499    @D[i,N+i]=1;
1500  }
1501  def @R2 = nc_algebra(1,@D);
1502  setring @R2;
1503  kill @R2@;
1504  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
1505  dbprint(ppl-1, @R2);
1506  ideal MM = maxideal(1);
1507  MM = 0,s,MM;
1508  map R01 = @R, MM;
1509  ideal K = R01(K);
1510  // total cleanup
1511  ideal LD = K;
1512  // make leadcoeffs positive
1513  for (i=1; i<= ncols(LD); i++)
1514  {
1515    if (leadcoef(LD[i]) <0 )
1516    {
1517      LD[i] = -LD[i];
1518    }
1519  }
1520  export LD;
1521  kill @R;
1522  return(@R2);
1523}
1524example;
1525{
1526  "EXAMPLE:"; echo = 2;
1527  ring r = 0,(x,y,Dx,Dy),dp;
1528  def R = Weyl(); // Weyl algebra on the variables x,y,Dx,Dy
1529  setring R;
1530  poly F = x2-y3;
1531  ideal I = Dx*F, Dy*F;
1532  // note, that I is not holonomic, since it's dimension is not 2
1533  gkdim(I);  // 3, while dim R = 4
1534  def W = SDLoc(I,F);
1535  setring W; // = R[s], where s is a new variable
1536  LD;        // Groebner basis of s-parametric presentation
1537}
1538
1539
1540// Groebner basis wrt weights and initial ideal business //////////////////////
1541
1542proc GBWeight (ideal I, intvec u, intvec v, list #)
1543  "USAGE:  GBWeight(I,u,v [,s,t,w]);
1544@*       I ideal, u,v intvecs, s,t optional ints, w an optional intvec
1545RETURN:  ideal, Groebner basis of I w.r.t. the weights u and v
1546ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
1547@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
1548@*       holds, i.e. the sequence of variables is given by
1549@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
1550@*       belonging to x(i).
1551PURPOSE: computes a Groebner basis with respect to given weights
1552NOTE:    The weights u and v are understood as weight vectors for x(i) and D(i),
1553@*       respectively. According to (SST), one computes the homogenization of a
1554@*       given ideal relative to (u,v), then one computes a Groebner basis and
1555@*       returns the dehomogenization of the result.
1556@*       If s<>0, @code{std} is used for Groebner basis computations,
1557@*       otherwise, and by default, @code{slimgb} is used.
1558@*       If t<>0, a matrix ordering is used for Groebner basis computations,
1559@*       otherwise, and by default, a block ordering is used.
1560@*       If w is given and consists of exactly 2*n strictly positive entries,
1561@*       w is used for constructing the weighted homogenized Weyl algebra,
1562@*       see Noro (2002). Otherwise, and by default, the homogenization weight
1563@*       (1,...,1) is used.
1564DISPLAY: If printlevel=1, progress debug messages will be printed,
1565@*       if printlevel>=2, all the debug messages will be printed.
1566EXAMPLE: example GBWeight; shows examples
1567"
1568{
1569  dmodappMoreAssumeViolation();
1570  int ppl = printlevel - voice +2;
1571  def save = basering;
1572  int n = nvars(save)/2;
1573  int whichengine = 0;           // default
1574  int methodord   = 0;           // default
1575  intvec homogweights = 1:(2*n); // default
1576  if (size(#)>0)
1577  {
1578    if (intLike(#[1]))
1579    {
1580      whichengine = int(#[1]);
1581    }
1582    if (size(#)>1)
1583    {
1584      if (intLike(#[2]))
1585      {
1586        methodord = int(#[2]);
1587      }
1588      if (size(#)>2)
1589      {
1590        if (typeof(#[3])=="intvec")
1591        {
1592          if (size(#[3])==2*n && allPositive(#[3])==1)
1593          {
1594            homogweights = #[3];
1595          }
1596          else
1597          {
1598            print("// Homogenization weight vector must consist of positive entries and be");
1599            print("// of size " + string(n) + ". Using weight (1,...,1).");
1600          }
1601        }
1602      }
1603    }
1604  }
1605  // 1. create  homogenized Weyl algebra
1606  // 1.1 create ordering
1607  int i;
1608  list RL = ringlist(save);
1609  int N = 2*n+1;
1610  intvec uv = u,v,0;
1611  homogweights = homogweights,1;
1612  list Lord = list(list("a",homogweights));
1613  list C0 = list("C",intvec(0));
1614  if (methodord == 0) // default: blockordering
1615  {
1616    Lord[5] = C0;
1617    Lord[4] = list("lp",intvec(1));
1618    Lord[3] = list("dp",intvec(1:(N-1)));
1619    Lord[2] = list("a",uv);
1620  }
1621  else                // M() ordering
1622  {
1623    intmat @Ord[N][N];
1624    @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
1625    for (i=1; i<=N-2; i++)
1626    {
1627      @Ord[2+i,N - i] = -1;
1628    }
1629    dbprint(ppl-1,"// the ordering matrix:",@Ord);
1630    Lord[2] = list("M",intvec(@Ord));
1631    Lord[3] = C0;
1632  }
1633  // 1.2 the homog var
1634  list Lvar = RL[2]; Lvar[N] = safeVarName("h","cv");
1635  // 1.3 create commutative ring
1636  list L@@Dh = RL; L@@Dh = L@@Dh[1..4];
1637  L@@Dh[2] = Lvar; L@@Dh[3] = Lord;
1638  def @Dh = ring(L@@Dh); kill L@@Dh;
1639  setring @Dh;
1640  // 1.4 create non-commutative relations
1641  matrix @relD[N][N];
1642  for (i=1; i<=n; i++)
1643  {
1644    @relD[i,n+i] = var(N)^(homogweights[i]+homogweights[n+i]);
1645  }
1646  def Dh = nc_algebra(1,@relD);
1647  setring Dh; kill @Dh;
1648  dbprint(ppl-1,"// computing in ring",Dh);
1649  // 2. Compute the initial ideal
1650  ideal I = imap(save,I);
1651  I = homog(I,var(N));
1652  // 2.1 the hard part: Groebner basis computation
1653  dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine));
1654  I = engine(I, whichengine);
1655  dbprint(ppl, "// finished Groebner basis computation");
1656  dbprint(ppl-1, "// ideal before dehomogenization is " +string(I));
1657  I = subst(I,var(N),1); // dehomogenization
1658  setring save;
1659  I = imap(Dh,I); kill Dh;
1660  return(I);
1661}
1662example
1663{
1664  "EXAMPLE:"; echo = 2;
1665  ring r = 0,(x,y,Dx,Dy),dp;
1666  def D2 = Weyl();
1667  setring D2;
1668  ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6;
1669  intvec u = -2,-3;
1670  intvec v = -u;
1671  GBWeight(I,u,v);
1672  ideal J = std(I);
1673  GBWeight(J,u,v); // same as above
1674  u = 0,1;
1675  GBWeight(I,u,v);
1676}
1677
1678proc inForm (def I, intvec w)
1679"USAGE:  inForm(I,w);  I ideal or poly, w intvec
1680RETURN:  ideal, generated by initial forms of generators of I w.r.t. w, or
1681@*       poly, initial form of input poly w.r.t. w
1682PURPOSE: computes the initial form of an ideal or a poly w.r.t. the weight w
1683NOTE:    The size of the weight vector must be equal to the number of variables
1684@*       of the basering.
1685EXAMPLE: example inForm; shows examples
1686"
1687{
1688  string inp1 = typeof(I);
1689  if ((inp1 <> "ideal") && (inp1 <> "poly"))
1690  {
1691    ERROR("first argument has to be an ideal or a poly");
1692  }
1693  if (size(w) != nvars(basering))
1694  {
1695    ERROR("weight vector has wrong dimension");
1696  }
1697  ideal II = I;
1698  if (simplify(II,2) == 0)
1699  {
1700    return(I);
1701  }
1702  int j,i;
1703  bigint s,m;
1704  list l;
1705  poly g;
1706  ideal J;
1707  for (j=1; j<=ncols(II); j++)
1708  {
1709    l = poly2list(II[j]);
1710    m = scalarProd(w,l[1][1]);
1711    g = l[1][2];
1712    for (i=2; i<=size(l); i++)
1713    {
1714      s = scalarProd(w,l[i][1]);
1715      if (s == m)
1716      {
1717        g = g + l[i][2];
1718      }
1719      else
1720      {
1721        if (s > m)
1722        {
1723          m = s;
1724          g = l[i][2];
1725        }
1726      }
1727    }
1728    J[j] = g;
1729  }
1730  if (inp1 == "ideal")
1731  {
1732    return(J);
1733  }
1734  else
1735  {
1736    return(J[1]);
1737  }
1738}
1739example
1740{
1741  "EXAMPLE:"; echo = 2;
1742  ring r = 0,(x,y,Dx,Dy),dp;
1743  def D = Weyl(); setring D;
1744  poly F = 3*x^2*Dy+2*y*Dx;
1745  poly G = 2*x*Dx+3*y*Dy+6;
1746  ideal I = F,G;
1747  intvec w1 = -1,-1,1,1;
1748  intvec w2 = -1,-2,1,2;
1749  intvec w3 = -2,-3,2,3;
1750  inForm(I,w1);
1751  inForm(I,w2);
1752  inForm(I,w3);
1753  inForm(F,w1);
1754}
1755
1756proc initialIdealW(ideal I, intvec u, intvec v, list #)
1757  "USAGE:  initialIdealW(I,u,v [,s,t,w]);
1758@*       I ideal, u,v intvecs, s,t optional ints, w an optional intvec
1759RETURN:  ideal, GB of initial ideal of the input ideal wrt the weights u and v
1760ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and for all
1761@*       1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the
1762@*       sequence of variables is given by x(1),...,x(n),D(1),...,D(n),
1763@*       where D(i) is the differential operator belonging to x(i).
1764PURPOSE: computes the initial ideal with respect to given weights.
1765NOTE:    u and v are understood as weight vectors for x(1..n) and D(1..n)
1766@*       respectively.
1767@*       If s<>0, @code{std} is used for Groebner basis computations,
1768@*       otherwise, and by default, @code{slimgb} is used.
1769@*       If t<>0, a matrix ordering is used for Groebner basis computations,
1770@*       otherwise, and by default, a block ordering is used.
1771@*       If w is given and consists of exactly 2*n strictly positive entries,
1772@*       w is used as homogenization weight.
1773@*       Otherwise, and by default, the homogenization weight (1,...,1) is used.
1774DISPLAY: If printlevel=1, progress debug messages will be printed,
1775@*       if printlevel>=2, all the debug messages will be printed.
1776EXAMPLE: example initialIdealW; shows examples
1777"
1778{
1779  // assumption check in GBWeight
1780  int ppl = printlevel - voice + 2;
1781  printlevel = printlevel + 1;
1782  I = GBWeight(I,u,v,#);
1783  printlevel = printlevel - 1;
1784  intvec uv = u,v;
1785  I = inForm(I,uv);
1786  int eng;
1787  if (size(#)>0)
1788  {
1789    if(typeof(#[1])=="int" || typeof(#[1])=="number")
1790    {
1791      eng = int(#[1]);
1792    }
1793  }
1794  dbprint(ppl,"// starting cosmetic Groebner basis computation");
1795  I = engine(I,eng);
1796  dbprint(ppl,"// finished cosmetic Groebner basis computation");
1797  return(I);
1798}
1799example
1800{
1801  "EXAMPLE:"; echo = 2;
1802  ring r = 0,(x,y,Dx,Dy),dp;
1803  def D2 = Weyl();
1804  setring D2;
1805  ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6;
1806  intvec u = -2,-3;
1807  intvec v = -u;
1808  initialIdealW(I,u,v);
1809  ideal J = std(I);
1810  initialIdealW(J,u,v); // same as above
1811  u = 0,1;
1812  initialIdealW(I,u,v);
1813}
1814
1815proc initialMalgrange (poly f,list #)
1816  "USAGE:  initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec
1817RETURN:  ring, Weyl algebra induced by basering, extended by two new vars t,Dt
1818PURPOSE: computes the initial Malgrange ideal of a given polynomial w.r.t. the
1819@*       weight vector (-1,0...,0,1,0,...,0) such that the weight of t is -1
1820@*       and the weight of Dt is 1.
1821ASSUME:  The basering is commutative and over a field of characteristic 0.
1822NOTE:    Activate the output ring with the @code{setring} command.
1823@*       The returned ring contains the ideal 'inF', being the initial ideal
1824@*       of the Malgrange ideal of f.
1825@*       Varnames of the basering should not include t and Dt.
1826@*       If a<>0, @code{std} is used for Groebner basis computations,
1827@*       otherwise, and by default, @code{slimgb} is used.
1828@*       If b<>0, a matrix ordering is used for Groebner basis computations,
1829@*       otherwise, and by default, a block ordering is used.
1830@*       If a positive weight vector v is given, the weight
1831@*       (d,v[1],...,v[n],1,d+1-v[1],...,d+1-v[n]) is used for homogenization
1832@*       computations, where d denotes the weighted degree of f.
1833@*       Otherwise and by default, v is set to (1,...,1). See Noro (2002).
1834DISPLAY: If printlevel=1, progress debug messages will be printed,
1835@*       if printlevel>=2, all the debug messages will be printed.
1836EXAMPLE: example initialMalgrange; shows examples
1837"
1838{
1839  dmodappAssumeViolation();
1840  if (!isCommutative())
1841  {
1842    ERROR("Basering must be commutative.");
1843  }
1844  int ppl = printlevel - voice + 2;
1845  def save = basering;
1846  int n = nvars(save);
1847  int i;
1848  int whichengine = 0; // default
1849  int methodord   = 0; // default
1850  intvec u0 = 1:n;     // default
1851  if (size(#)>0)
1852  {
1853    if (intLike(#[1]))
1854    {
1855      whichengine = int(#[1]);
1856    }
1857    if (size(#)>1)
1858    {
1859      if (intLike(#[2]))
1860      {
1861        methodord = int(#[2]);
1862      }
1863      if (size(#)>2)
1864      {
1865        if ((typeof(#[3])=="intvec") && (size(#[3])==n) && (allPositive(#[3])==1))
1866        {
1867          u0 = #[3];
1868        }
1869      }
1870    }
1871  }
1872  list RL = ringlist(save);
1873  RL = RL[1..4]; // if basering is commutative nc_algebra
1874  list C0 = list("C",intvec(0));
1875  // 1. create homogenization weights
1876  // 1.1. get the weighted degree of f
1877  list Lord = list(list("wp",u0),C0);
1878  list L@r = RL;
1879  L@r[3] = Lord;
1880  def r = ring(L@r); kill L@r,Lord;
1881  setring r;
1882  poly f = imap(save,f);
1883  int d = deg(f);
1884  setring save; kill r;
1885  // 1.2 the homogenization weights
1886  intvec homogweights = d;
1887  homogweights[n+2] = 1;
1888  for (i=1; i<=n; i++)
1889  {
1890    homogweights[i+1]   = u0[i];
1891    homogweights[n+2+i] = d+1-u0[i];
1892  }
1893  // 2. create extended Weyl algebra
1894  int N = 2*n+2;
1895  // 2.1 create names for vars
1896  string vart  = safeVarName("t","cv");
1897  string varDt = safeVarName("D"+vart,"cv");
1898  while (varDt <> "D"+vart)
1899  {
1900    vart  = safeVarName("@"+vart,"cv");
1901    varDt = safeVarName("D"+vart,"cv");
1902  }
1903  list Lvar;
1904  Lvar[1] = vart; Lvar[n+2] = varDt;
1905  for (i=1; i<=n; i++)
1906  {
1907    Lvar[i+1]   = string(var(i));
1908    Lvar[i+n+2] = safeVarName("D" + string(var(i)),"cv");
1909  }
1910  //  2.2 create ordering
1911  list Lord = list(list("dp",intvec(1:N)),C0);
1912  // 2.3 create the (n+1)-th Weyl algebra
1913  list L@D = RL; L@D[2] = Lvar; L@D[3] = Lord;
1914  def @D = ring(L@D); kill L@D;
1915  setring @D;
1916  def D = Weyl();
1917  setring D; kill @D;
1918  dbprint(ppl,"// the (n+1)-th Weyl algebra :" ,D);
1919  // 3. compute the initial ideal
1920  // 3.1 create the Malgrange ideal
1921  poly f = imap(save,f);
1922  ideal I = var(1)-f;
1923  for (i=1; i<=n; i++)
1924  {
1925    I = I, var(n+2+i)+diff(f,var(i+1))*var(n+2);
1926  }
1927  // I = engine(I,whichengine); // todo efficient to compute GB wrt dp first?
1928  // 3.2 computie the initial ideal
1929  intvec w = 1,0:n;
1930  printlevel = printlevel + 1;
1931  I = initialIdealW(I,-w,w,whichengine,methodord,homogweights);
1932  printlevel = printlevel - 1;
1933  ideal inF = I; attrib(inF,"isSB",1);
1934  export(inF);
1935  setring save;
1936  return(D);
1937}
1938example
1939{
1940  "EXAMPLE:"; echo = 2;
1941  ring r = 0,(x,y),dp;
1942  poly f = x^2+y^3+x*y^2;
1943  def D = initialMalgrange(f);
1944  setring D;
1945  inF;
1946  setring r;
1947  intvec v = 3,2;
1948  def D2 = initialMalgrange(f,1,1,v);
1949  setring D2;
1950  inF;
1951}
1952
1953
1954// restriction and integration ////////////////////////////////////////////////
1955
1956static proc restrictionModuleEngine (ideal I, intvec w, list #)
1957// returns list L with 2 entries of type ideal
1958// L[1]=basis of free module, L[2]=generating system of submodule
1959// #=eng,m,G; eng=engine; m=min int root of bfctIdeal(I,w); G=GB of I wrt (-w,w)
1960{
1961  dmodappMoreAssumeViolation();
1962  if (!isHolonomic(I))
1963  {
1964    ERROR("Given ideal is not holonomic");
1965  }
1966  int l0,l0set,Gset;
1967  ideal G;
1968  int whichengine = 0;         // default
1969  if (size(#)>0)
1970  {
1971    if (intLike(#[1]))
1972    {
1973      whichengine = int(#[1]);
1974    }
1975    if (size(#)>1)
1976    {
1977      if (intLike(#[2]))
1978      {
1979        l0 = int(#[2]);
1980        l0set = 1;
1981      }
1982      if (size(#)>2)
1983      {
1984        if (typeof(#[3])=="ideal")
1985        {
1986          G = #[3];
1987          Gset = 1;
1988        }
1989      }
1990    }
1991  }
1992  int ppl = printlevel;
1993  int i,j,k;
1994  int n = nvars(basering)/2;
1995  if (w == 0:size(w))
1996  {
1997    ERROR("weight vector must not be zero");
1998  }
1999  if (size(w)<>n)
2000  {
2001    ERROR("weight vector must have exactly " + string(n) + " entries");
2002  }
2003  for (i=1; i<=n; i++)
2004  {
2005    if (w[i]<0)
2006    {
2007      ERROR("weight vector must not have negative entries");
2008    }
2009  }
2010  intvec ww = -w,w;
2011  if (!Gset)
2012  {
2013    G = GBWeight(I,-w,w,whichengine);
2014    dbprint(ppl,"// found GB wrt weight " +string(w));
2015    dbprint(ppl-1,"// " + string(G));
2016  }
2017  if (!l0set)
2018  {
2019    ideal inG = inForm(G,ww);
2020    inG = engine(inG,whichengine);
2021    poly s = 0;
2022    for (i=1; i<=n; i++)
2023    {
2024      s = s + w[i]*var(i)*var(i+n);
2025    }
2026    vector v = pIntersect(s,inG);
2027    list L = bFactor(vec2poly(v));
2028    dbprint(ppl,"// found b-function of given ideal wrt given weight");
2029    dbprint(ppl-1,"// roots: "+string(L[1]));
2030    dbprint(ppl-1,"// multiplicities: "+string(L[2]));
2031    kill inG,v,s;
2032    L = intRoots(L);           // integral roots of b-function
2033    if (L[2]==0:size(L[2]))       // no integral roots
2034    {
2035      return(list(ideal(0),ideal(0)));
2036    }
2037    intvec v;
2038    for (i=1; i<=ncols(L[1]); i++)
2039    {
2040      v[i] = int(L[1][i]);
2041    }
2042    l0 = Max(v);
2043    dbprint(ppl,"// maximal integral root is " +string(l0));
2044    kill L,v;
2045  }
2046  if (l0 < 0)                  // maximal integral root is < 0
2047  {
2048    return(list(ideal(0),ideal(0)));
2049  }
2050  intvec m;
2051  for (i=ncols(G); i>0; i--)
2052  {
2053    m[i] = deg(G[i],ww);
2054  }
2055  dbprint(ppl,"// weighted degree of generators of GB is " +string(m));
2056  def save = basering;
2057  list RL = ringlist(save);
2058  RL = RL[1..4];
2059  list Lvar;
2060  j = 1;
2061  intvec neww;
2062  for (i=1; i<=n; i++)
2063  {
2064    if (w[i]>0)
2065    {
2066      Lvar[j] = string(var(i+n));
2067      neww[j] = w[i];
2068      j++;
2069    }
2070  }
2071  list Lord;
2072  Lord[1] = list("dp",intvec(1:n));
2073  Lord[2] = list("C", intvec(0));
2074  RL[2] = Lvar;
2075  RL[3] = Lord;
2076  def r = ring(RL);
2077  kill Lvar, Lord, RL;
2078  setring r;
2079  ideal B;
2080  list Blist;
2081  intvec mm = l0,-m+l0;
2082  for (i=0; i<=Max(mm); i++)
2083  {
2084    B = weightKB(std(0),i,neww);
2085    Blist[i+1] = B;
2086  }
2087  setring save;
2088  list Blist = imap(r,Blist);
2089  ideal ff = maxideal(1);
2090  for (i=1; i<=n; i++)
2091  {
2092    if (w[i]<>0)
2093    {
2094      ff[i] = 0;
2095    }
2096  }
2097  map f = save,ff;
2098  ideal B,M;
2099  poly p;
2100  for (i=1; i<=size(G); i++)
2101  {
2102    for (j=1; j<=l0-m[i]+1; j++)
2103    {
2104      B = Blist[j];
2105      for (k=1; k<=ncols(B); k++)
2106      {
2107        p = B[k]*G[i];
2108        p = f(p);
2109        M[size(M)+1] = p;
2110      }
2111    }
2112  }
2113  ideal Bl0 = Blist[1..(l0+1)];
2114  dbprint(ppl,"// found basis of free module");
2115  dbprint(ppl-1,"// " + string(Bl0));
2116  dbprint(ppl,"// found generators of submodule");
2117  dbprint(ppl-1,"// " + string(M));
2118  return(list(Bl0,M));
2119}
2120
2121static proc restrictionModuleOutput (ideal B, ideal N, intvec w, int eng, string str)
2122// returns ring, which contains module "str"
2123{
2124  int n = nvars(basering)/2;
2125  int i,j;
2126  def save = basering;
2127  // 1: create new ring
2128  list RL = ringlist(save);
2129  RL = RL[1..4];
2130  list V = RL[2];
2131  poly prodvar = 1;
2132  int zeropresent;
2133  j = 0;
2134  for (i=1; i<=n; i++)
2135  {
2136    if (w[i]<>0)
2137    {
2138      V = delete(V,i-j);
2139      V = delete(V,i-2*j-1+n);
2140      j = j+1;
2141      prodvar = prodvar*var(i)*var(i+n);
2142    }
2143    else
2144    {
2145      zeropresent = 1;
2146    }
2147  }
2148  if (!zeropresent) // restrict/integrate all vars, return input ring
2149  {
2150    def newR = save;
2151  }
2152  else
2153  {
2154    RL[2] = V;
2155    V = list();
2156    V[1] = list("C", intvec(0));
2157    V[2] = list("dp",intvec(1:(2*size(ideal(w)))));
2158    RL[3] = V;
2159    def @D = ring(RL);
2160    setring @D;
2161    def newR = Weyl();
2162    setring save;
2163    kill @D;
2164  }
2165  // 2. get coker representation of module
2166  module M = coeffs(N,B,prodvar);
2167  if (zeropresent)
2168  {
2169    setring newR;
2170    module M = imap(save,M);
2171  }
2172  M = engine(M,eng);
2173  M = prune(M);
2174  M = engine(M,eng);
2175  execute("module " + str + " = M;");
2176  execute("export(" + str + ");");
2177  setring save;
2178  return(newR);
2179}
2180
2181proc restrictionModule (ideal I, intvec w, list #)
2182  "USAGE:  restrictionModule(I,w,[,eng,m,G]);
2183@*       I ideal, w intvec, eng and m optional ints, G optional ideal
2184RETURN:  ring (a Weyl algebra) containing a module 'resMod'
2185ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
2186@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
2187@*       holds, i.e. the sequence of variables is given by
2188@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
2189@*       belonging to x(i).
2190@*       Further, assume that I is holonomic and that w is n-dimensional with
2191@*       non-negative entries.
2192PURPOSE: computes the restriction module of a holonomic ideal to the subspace
2193@*       defined by the variables corresponding to the non-zero entries of the
2194@*       given intvec
2195NOTE:    The output ring is the Weyl algebra defined by the zero entries of w.
2196@*       It contains a module 'resMod' being the restriction module of I wrt w.
2197@*       If there are no zero entries, the input ring is returned.
2198@*       If eng<>0, @code{std} is used for Groebner basis computations,
2199@*       otherwise, and by default, @code{slimgb} is used.
2200@*       The minimal integer root of the b-function of I wrt the weight (-w,w)
2201@*       can be specified via the optional argument m.
2202@*       The optional argument G is used for specifying a Groebner Basis of I
2203@*       wrt the weight (-w,w), that is, the initial form of G generates the
2204@*       initial ideal of I wrt the weight (-w,w).
2205@*       Further note, that the assumptions on m and G (if given) are not
2206@*       checked.
2207DISPLAY: If printlevel=1, progress debug messages will be printed,
2208@*       if printlevel>=2, all the debug messages will be printed.
2209EXAMPLE: example restrictionModule; shows examples
2210"
2211{
2212  list L = restrictionModuleEngine(I,w,#);
2213  int eng;
2214  if (size(#)>0)
2215  {
2216    if (intLike(#[1]))
2217    {
2218      eng = int(#[1]);
2219    }
2220  }
2221  def newR = restrictionModuleOutput(L[1],L[2],w,eng,"resMod");
2222  return(newR);
2223}
2224example
2225{
2226  "EXAMPLE:"; echo = 2;
2227  ring r = 0,(a,x,b,Da,Dx,Db),dp;
2228  def D3 = Weyl();
2229  setring D3;
2230  ideal I = a*Db-Dx+2*Da, x*Db-Da, x*Da+a*Da+b*Db+1,
2231    x*Dx-2*x*Da-a*Da, b*Db^2+Dx*Da-Da^2+Db,
2232    a*Dx*Da+2*x*Da^2+a*Da^2+b*Dx*Db+Dx+2*Da;
2233  intvec w = 1,0,0;
2234  def rm = restrictionModule(I,w);
2235  setring rm; rm;
2236  print(resMod);
2237}
2238
2239static proc restrictionIdealEngine (ideal I, intvec w, string cf, list #)
2240{
2241  int eng;
2242  if (size(#)>0)
2243  {
2244    if(intLike(#[1]))
2245    {
2246      eng = int(#[1]);
2247    }
2248  }
2249  def save = basering;
2250  if (cf == "restriction")
2251  {
2252    def newR = restrictionModule(I,w,#);
2253    setring newR;
2254    matrix M = resMod;
2255    kill resMod;
2256  }
2257  if (cf == "integral")
2258  {
2259    def newR = integralModule(I,w,#);
2260    setring newR;
2261    matrix M = intMod;
2262    kill intMod;
2263  }
2264  int i,r,c;
2265  r = nrows(M);
2266  c = ncols(M);
2267  ideal J;
2268  if (r == 1) // nothing to do
2269  {
2270    J = M;
2271  }
2272  else
2273  {
2274    matrix zm[r-1][1]; // zero matrix
2275    matrix v[r-1][1];
2276    for (i=1; i<=c; i++)
2277    {
2278      if (M[1,i]<>0)
2279      {
2280        v = M[2..r,i];
2281        if (v == zm)
2282        {
2283          J[size(J)+1] = M[1,i];
2284        }
2285      }
2286    }
2287  }
2288  J = engine(J,eng);
2289  if (cf == "restriction")
2290  {
2291    ideal resIdeal = J;
2292    export(resIdeal);
2293  }
2294  if (cf == "integral")
2295  {
2296    ideal intIdeal = J;
2297    export(intIdeal);
2298  }
2299  setring save;
2300  return(newR);
2301}
2302
2303proc restrictionIdeal (ideal I, intvec w, list #)
2304  "USAGE:  restrictionIdeal(I,w,[,eng,m,G]);
2305@*       I ideal, w intvec, eng and m optional ints, G optional ideal
2306RETURN:  ring (a Weyl algebra) containing an ideal 'resIdeal'
2307ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
2308@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
2309@*       holds, i.e. the sequence of variables is given by
2310@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
2311@*       belonging to x(i).
2312@*       Further, assume that I is holonomic and that w is n-dimensional with
2313@*       non-negative entries.
2314PURPOSE: computes the restriction ideal of a holonomic ideal to the subspace
2315@*       defined by the variables corresponding to the non-zero entries of the
2316@*       given intvec
2317NOTE:    The output ring is the Weyl algebra defined by the zero entries of w.
2318@*       It contains an ideal 'resIdeal' being the restriction ideal of I wrt w.
2319@*       If there are no zero entries, the input ring is returned.
2320@*       If eng<>0, @code{std} is used for Groebner basis computations,
2321@*       otherwise, and by default, @code{slimgb} is used.
2322@*       The minimal integer root of the b-function of I wrt the weight (-w,w)
2323@*       can be specified via the optional argument m.
2324@*       The optional argument G is used for specifying a Groebner basis of I
2325@*       wrt the weight (-w,w), that is, the initial form of G generates the
2326@*       initial ideal of I wrt the weight (-w,w).
2327@*       Further note, that the assumptions on m and G (if given) are not
2328@*       checked.
2329DISPLAY: If printlevel=1, progress debug messages will be printed,
2330@*       if printlevel>=2, all the debug messages will be printed.
2331EXAMPLE: example restrictionIdeal; shows examples
2332"
2333{
2334  def rm = restrictionIdealEngine(I,w,"restriction",#);
2335  return(rm);
2336}
2337example
2338{
2339  "EXAMPLE:"; echo = 2;
2340  ring r = 0,(a,x,b,Da,Dx,Db),dp;
2341  def D3 = Weyl();
2342  setring D3;
2343  ideal I = a*Db-Dx+2*Da,
2344    x*Db-Da,
2345    x*Da+a*Da+b*Db+1,
2346    x*Dx-2*x*Da-a*Da,
2347    b*Db^2+Dx*Da-Da^2+Db,
2348    a*Dx*Da+2*x*Da^2+a*Da^2+b*Dx*Db+Dx+2*Da;
2349  intvec w = 1,0,0;
2350  def D2 = restrictionIdeal(I,w);
2351  setring D2; D2;
2352  resIdeal;
2353}
2354
2355proc fourier (ideal I, list #)
2356  "USAGE:  fourier(I[,v]); I an ideal, v an optional intvec
2357RETURN:  ideal
2358PURPOSE: computes the Fourier transform of an ideal in a Weyl algebra
2359ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
2360@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
2361@*       holds, i.e. the sequence of variables is given by
2362@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
2363@*       belonging to x(i).
2364NOTE:    The Fourier automorphism is defined by mapping x(i) to -D(i) and
2365@*       D(i) to x(i).
2366@*       If v is an intvec with entries ranging from 1 to n, the Fourier
2367@*       transform of I restricted to the variables given by v is computed.
2368SEE ALSO: inverseFourier
2369EXAMPLE: example fourier; shows examples
2370"
2371{
2372  dmodappMoreAssumeViolation();
2373  intvec v;
2374  if (size(#)>0)
2375  {
2376    if(typeof(#[1])=="intvec")
2377    {
2378      v = #[1];
2379    }
2380  }
2381  int n = nvars(basering)/2;
2382  int i;
2383  if(v <> 0:size(v))
2384  {
2385    v = sortIntvec(v)[1];
2386    for (i=1; i<size(v); i++)
2387    {
2388      if (v[i] == v[i+1])
2389      {
2390        ERROR("No double entries allowed in intvec");
2391      }
2392    }
2393  }
2394  else
2395  {
2396    v = 1..n;
2397  }
2398  ideal m = maxideal(1);
2399  for (i=1; i<=size(v); i++)
2400  {
2401    if (v[i]<0 || v[i]>n)
2402    {
2403      ERROR("Entries of intvec must range from 1 to "+string(n));
2404    }
2405    m[v[i]] = -var(v[i]+n);
2406    m[v[i]+n] = var(v[i]);
2407  }
2408  map F = basering,m;
2409  ideal FI = F(I);
2410  return(FI);
2411}
2412example
2413{
2414  "EXAMPLE:"; echo = 2;
2415  ring r = 0,(x,y,Dx,Dy),dp;
2416  def D2 = Weyl();
2417  setring D2;
2418  ideal I = x*Dx+2*y*Dy+2, x^2*Dx+y*Dx+2*x;
2419  intvec v = 2;
2420  fourier(I,v);
2421  fourier(I);
2422}
2423
2424proc inverseFourier (ideal I, list #)
2425  "USAGE:  inverseFourier(I[,v]); I an ideal, v an optional intvec
2426RETURN:  ideal
2427PURPOSE: computes the inverse Fourier transform of an ideal in a Weyl algebra
2428ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
2429@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
2430@*       holds, i.e. the sequence of variables is given by
2431@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
2432@*       belonging to x(i).
2433NOTE:    The Fourier automorphism is defined by mapping x(i) to -D(i) and
2434@*       D(i) to x(i).
2435@*       If v is an intvec with entries ranging from 1 to n, the inverse Fourier
2436@*       transform of I restricted to the variables given by v is computed.
2437SEE ALSO: fourier
2438EXAMPLE: example inverseFourier; shows examples
2439"
2440{
2441  dmodappMoreAssumeViolation();
2442  intvec v;
2443  if (size(#)>0)
2444  {
2445    if(typeof(#[1])=="intvec")
2446    {
2447      v = #[1];
2448    }
2449  }
2450  int n = nvars(basering)/2;
2451  int i;
2452  if(v <> 0:size(v))
2453  {
2454    v = sortIntvec(v)[1];
2455    for (i=1; i<size(v); i++)
2456    {
2457      if (v[i] == v[i+1])
2458      {
2459        ERROR("No double entries allowed in intvec");
2460      }
2461    }
2462  }
2463  else
2464  {
2465    v = 1..n;
2466  }
2467  ideal m = maxideal(1);
2468  for (i=1; i<=size(v); i++)
2469  {
2470    if (v[i]<0 || v[i]>n)
2471    {
2472      ERROR("Entries of intvec must range between 1 and "+string(n));
2473    }
2474    m[v[i]] = var(v[i]+n);
2475    m[v[i]+n] = -var(v[i]);
2476  }
2477  map F = basering,m;
2478  ideal FI = F(I);
2479  return(FI);
2480}
2481example
2482{
2483  "EXAMPLE:"; echo = 2;
2484  ring r = 0,(x,y,Dx,Dy),dp;
2485  def D2 = Weyl();
2486  setring D2;
2487  ideal I = x*Dx+2*y*Dy+2, x^2*Dx+y*Dx+2*x;
2488  intvec v = 2;
2489  ideal FI = fourier(I);
2490  inverseFourier(FI);
2491}
2492
2493proc integralModule (ideal I, intvec w, list #)
2494  "USAGE:  integralModule(I,w,[,eng,m,G]);
2495@*       I ideal, w intvec, eng and m optional ints, G optional ideal
2496RETURN:  ring (a Weyl algebra) containing a module 'intMod'
2497ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
2498@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
2499@*       holds, i.e. the sequence of variables is given by
2500@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
2501@*       belonging to x(i).
2502@*       Further, assume that I is holonomic and that w is n-dimensional with
2503@*       non-negative entries.
2504PURPOSE: computes the integral module of a holonomic ideal w.r.t. the subspace
2505@*       defined by the variables corresponding to the non-zero entries of the
2506@*       given intvec
2507NOTE:    The output ring is the Weyl algebra defined by the zero entries of w.
2508@*       It contains a module 'intMod' being the integral module of I wrt w.
2509@*       If there are no zero entries, the input ring is returned.
2510@*       If eng<>0, @code{std} is used for Groebner basis computations,
2511@*       otherwise, and by default, @code{slimgb} is used.
2512@*       Let F(I) denote the Fourier transform of I w.r.t. w.
2513@*       The minimal integer root of the b-function of F(I) w.r.t. the weight
2514@*       (-w,w) can be specified via the optional argument m.
2515@*       The optional argument G is used for specifying a Groebner Basis of F(I)
2516@*       wrt the weight (-w,w), that is, the initial form of G generates the
2517@*       initial ideal of F(I) w.r.t. the weight (-w,w).
2518@*       Further note, that the assumptions on m and G (if given) are not
2519@*       checked.
2520DISPLAY: If printlevel=1, progress debug messages will be printed,
2521@*       if printlevel>=2, all the debug messages will be printed.
2522EXAMPLE: example integralModule; shows examples
2523"
2524{
2525  int l0,l0set,Gset;
2526  ideal G;
2527  int whichengine = 0;         // default
2528  if (size(#)>0)
2529  {
2530    if (intLike(#[1]))
2531    {
2532      whichengine = int(#[1]);
2533    }
2534    if (size(#)>1)
2535    {
2536      if (intLike(#[2]))
2537      {
2538        l0 = int(#[2]);
2539        l0set = 1;
2540      }
2541      if (size(#)>2)
2542      {
2543        if (typeof(#[3])=="ideal")
2544        {
2545          G = #[3];
2546          Gset = 1;
2547        }
2548      }
2549    }
2550  }
2551  int ppl = printlevel;
2552  int i;
2553  int n = nvars(basering)/2;
2554  intvec v;
2555  for (i=1; i<=n; i++)
2556  {
2557    if (w[i]>0)
2558    {
2559      if (v == 0:size(v))
2560      {
2561        v[1] = i;
2562      }
2563      else
2564      {
2565        v[size(v)+1] = i;
2566      }
2567    }
2568  }
2569  ideal FI = fourier(I,v);
2570  dbprint(ppl,"// computed Fourier transform of given ideal");
2571  dbprint(ppl-1,"// " + string(FI));
2572  list L;
2573  if (l0set)
2574  {
2575    if (Gset) // l0 and G given
2576    {
2577      L = restrictionModuleEngine(FI,w,whichengine,l0,G);
2578    }
2579    else      // l0 given, G not
2580    {
2581      L = restrictionModuleEngine(FI,w,whichengine,l0);
2582    }
2583  }
2584  else        // nothing given
2585  {
2586    L = restrictionModuleEngine(FI,w,whichengine);
2587  }
2588  ideal B,N;
2589  B = inverseFourier(L[1],v);
2590  N = inverseFourier(L[2],v);
2591  def newR = restrictionModuleOutput(B,N,w,whichengine,"intMod");
2592  return(newR);
2593}
2594example
2595{
2596  "EXAMPLE:"; echo = 2;
2597  ring r = 0,(x,b,Dx,Db),dp;
2598  def D2 = Weyl();
2599  setring D2;
2600  ideal I = x*Dx+2*b*Db+2, x^2*Dx+b*Dx+2*x;
2601  intvec w = 1,0;
2602  def im = integralModule(I,w);
2603  setring im; im;
2604  print(intMod);
2605}
2606
2607proc integralIdeal (ideal I, intvec w, list #)
2608  "USAGE:  integralIdeal(I,w,[,eng,m,G]);
2609@*       I ideal, w intvec, eng and m optional ints, G optional ideal
2610RETURN:  ring (a Weyl algebra) containing an ideal 'intIdeal'
2611ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
2612@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
2613@*       holds, i.e. the sequence of variables is given by
2614@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
2615@*       belonging to x(i).
2616@*       Further, assume that I is holonomic and that w is n-dimensional with
2617@*       non-negative entries.
2618PURPOSE: computes the integral ideal of a holonomic ideal w.r.t. the subspace
2619@*       defined by the variables corresponding to the non-zero entries of the
2620@*       given intvec.
2621NOTE:    The output ring is the Weyl algebra defined by the zero entries of w.
2622@*       It contains ideal 'intIdeal' being the integral ideal of I w.r.t. w.
2623@*       If there are no zero entries, the input ring is returned.
2624@*       If eng<>0, @code{std} is used for Groebner basis computations,
2625@*       otherwise, and by default, @code{slimgb} is used.
2626@*       The minimal integer root of the b-function of I wrt the weight (-w,w)
2627@*       can be specified via the optional argument m.
2628@*       The optional argument G is used for specifying a Groebner basis of I
2629@*       wrt the weight (-w,w), that is, the initial form of G generates the
2630@*       initial ideal of I wrt the weight (-w,w).
2631@*       Further note, that the assumptions on m and G (if given) are not
2632@*       checked.
2633DISPLAY: If printlevel=1, progress debug messages will be printed,
2634@*       if printlevel>=2, all the debug messages will be printed.
2635EXAMPLE: example integralIdeal; shows examples
2636"
2637{
2638  def im = restrictionIdealEngine(I,w,"integral",#);
2639  return(im);
2640}
2641example
2642{
2643  "EXAMPLE:"; echo = 2;
2644  ring r = 0,(x,b,Dx,Db),dp;
2645  def D2 = Weyl();
2646  setring D2;
2647  ideal I = x*Dx+2*b*Db+2, x^2*Dx+b*Dx+2*x;
2648  intvec w = 1,0;
2649  def D1 = integralIdeal(I,w);
2650  setring D1; D1;
2651  intIdeal;
2652}
2653
2654proc deRhamCohomIdeal (ideal I, list #)
2655  "USAGE:  deRhamCohomIdeal (I[,w,eng,k,G]);
2656@*       I ideal, w optional intvec, eng and k optional ints, G optional ideal
2657RETURN:  ideal
2658ASSUME:  The basering is the n-th Weyl algebra D over a field of characteristic
2659@*       zero and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
2660@*       holds, i.e. the sequence of variables is given by
2661@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
2662@*       belonging to x(i).
2663@*       Further, assume that I is of special kind, namely let f in K[x] and
2664@*       consider the module K[x,1/f]f^m, where m is smaller than or equal to
2665@*       the minimal integer root of the Bernstein-Sato polynomial of f.
2666@*       Since this module is known to be a holonomic D-module, it has a cyclic
2667@*       presentation D/I.
2668PURPOSE: computes a basis of the n-th de Rham cohomology group of the complement
2669@*       of the hypersurface defined by f
2670NOTE:    If I does not satisfy the assumptions described above, the result might
2671@*       have no meaning. Note that I can be computed with @code{annfs}.
2672@*       If w is an intvec with exactly n strictly positive entries, w is used
2673@*       in the computation. Otherwise, and by default, w is set to (1,...,1).
2674@*       If eng<>0, @code{std} is used for Groebner basis computations,
2675@*       otherwise, and by default, @code{slimgb} is used.
2676@*       Let F(I) denote the Fourier transform of I wrt w.
2677@*       An integer smaller than or equal to the minimal integer root of the
2678@*       b-function of F(I) wrt the weight (-w,w) can be specified via the
2679@*       optional argument k.
2680@*       The optional argument G is used for specifying a Groebner Basis of F(I)
2681@*       wrt the weight (-w,w), that is, the initial form of G generates the
2682@*       initial ideal of F(I) wrt the weight (-w,w).
2683@*       Further note, that the assumptions on I, k and G (if given) are not
2684@*       checked.
2685THEORY:  (SST) pp. 232-235
2686DISPLAY: If printlevel=1, progress debug messages will be printed,
2687@*       if printlevel>=2, all the debug messages will be printed.
2688SEE ALSO: deRhamCohom
2689EXAMPLE: example deRhamCohomIdeal; shows examples
2690"
2691{
2692  intvec w = 1:(nvars(basering)/2);
2693  int l0,l0set,Gset;
2694  ideal G;
2695  int whichengine = 0;         // default
2696  if (size(#)>0)
2697  {
2698    if (typeof(#[1])=="intvec")
2699    {
2700      if (allPositive(#[1])==1)
2701      {
2702        w = #[1];
2703      }
2704      else
2705      {
2706        print("// Entries of intvec must be strictly positive");
2707        print("// Using weight " + string(w));
2708      }
2709      if (size(#)>1)
2710      {
2711        if (intLike(#[2]))
2712        {
2713          whichengine = int(#[2]);
2714        }
2715        if (size(#)>2)
2716        {
2717          if (intLike(#[3]))
2718          {
2719            l0 = int(#[3]);
2720            l0set = 1;
2721          }
2722          if (size(#)>3)
2723          {
2724            if (typeof(#[4])=="ideal")
2725            {
2726              G = #[4];
2727              Gset = 1;
2728            }
2729          }
2730        }
2731      }
2732    }
2733  }
2734  int ppl = printlevel;
2735  int i,j;
2736  int n = nvars(basering)/2;
2737  intvec v;
2738  for (i=1; i<=n; i++)
2739  {
2740    if (w[i]>0)
2741    {
2742      if (v == 0:size(v))
2743      {
2744        v[1] = i;
2745      }
2746      else
2747      {
2748        v[size(v)+1] = i;
2749      }
2750    }
2751  }
2752  ideal FI = fourier(I,v);
2753  dbprint(ppl,"// computed Fourier transform of given ideal");
2754  dbprint(ppl-1,"// " + string(FI));
2755  list L;
2756  if (l0set)
2757  {
2758    if (Gset) // l0 and G given
2759    {
2760      L = restrictionModuleEngine(FI,w,whichengine,l0,G);
2761    }
2762    else      // l0 given, G not
2763    {
2764      L = restrictionModuleEngine(FI,w,whichengine,l0);
2765    }
2766  }
2767  else        // nothing given
2768  {
2769    L = restrictionModuleEngine(FI,w,whichengine);
2770  }
2771  ideal B,N;
2772  B = inverseFourier(L[1],v);
2773  N = inverseFourier(L[2],v);
2774  dbprint(ppl,"// computed integral module of given ideal");
2775  dbprint(ppl-1,"// " + string(B));
2776  dbprint(ppl-1,"// " + string(N));
2777  ideal DR;
2778  poly p;
2779  poly Dt = 1;
2780  for (i=1; i<=n; i++)
2781  {
2782    Dt = Dt*var(i+n);
2783  }
2784  N = simplify(N,2+8);
2785  printlevel = printlevel-1;
2786  N = linReduceIdeal(N);
2787  N = simplify(N,2+8);
2788  for (i=1; i<=size(B); i++)
2789  {
2790    p = linReduce(B[i],N);
2791    if (p<>0)
2792    {
2793      DR[size(DR)+1] = B[i]*Dt;
2794      j=1;
2795      while (p<N[j])
2796      {
2797        j++;
2798      }
2799      N = insertGenerator(N,p,j+1);
2800    }
2801  }
2802  printlevel = printlevel + 1;
2803  return(DR);
2804}
2805example
2806{
2807  "EXAMPLE:"; echo = 2;
2808  ring r = 0,(x,y,z),dp;
2809  poly F = x^3+y^3+z^3;
2810  bfctAnn(F);            // Bernstein-Sato poly of F has minimal integer root -2
2811  def W = annRat(1,F^2); // so we compute the annihilator of 1/F^2
2812  setring W; W;          // Weyl algebra, contains LD = Ann(1/F^2)
2813  LD;                    // K[x,y,z,1/F]F^(-2) is isomorphic to W/LD as W-module
2814  deRhamCohomIdeal(LD);  // we see that the K-dim is 2
2815}
2816
2817proc deRhamCohom (poly f, list #)
2818  "USAGE:  deRhamCohom(f[,eng,m]);  f poly, eng and m optional ints
2819RETURN:  ring (a Weyl Algebra) containing an ideal 'DR'
2820ASSUME:  Basering is a commutative and over a field of characteristic 0.
2821PURPOSE: computes a basis of the n-th de Rham cohomology group of the complement
2822@*       of the hypersurface defined by f, where n denotes the number of
2823@*       variables of the basering
2824NOTE:    The output ring is the n-th Weyl algebra. It contains an ideal 'DR'
2825@*       being a basis of the n-th de Rham cohomology group of the complement of
2826@*       the hypersurface defined by f.
2827@*       If eng<>0, @code{std} is used for Groebner basis computations,
2828@*       otherwise, and by default, @code{slimgb} is used.
2829@*       If m is given, it is assumed to be less than or equal to the minimal
2830@*       integer root of the Bernstein-Sato polynomial of f. This assumption is
2831@*       not checked. If not specified, m is set to the minimal integer root of
2832@*       the Bernstein-Sato polynomial of f.
2833THEORY:  (SST) pp. 232-235
2834DISPLAY: If printlevel=1, progress debug messages will be printed,
2835@*       if printlevel>=2, all the debug messages will be printed.
2836SEE ALSO: deRhamCohomIdeal
2837EXAMPLE: example deRhamCohom; shows example
2838"
2839{
2840  int ppl = printlevel - voice + 2;
2841  int eng,l0,l0given;
2842  if (size(#)>0)
2843  {
2844    if(intLike(#[1]))
2845    {
2846      eng = int(#[1]);
2847    }
2848    if (size(#)>1)
2849    {
2850      if(intLike(#[2]))
2851      {
2852        l0 = int(#[2]);
2853        l0given = 1;
2854      }
2855    }
2856  }
2857  if (!isCommutative())
2858  {
2859    ERROR("Basering must be commutative.");
2860  }
2861  int i;
2862  def save = basering;
2863  int n = nvars(save);
2864  dbprint(ppl,"// Computing s-parametric annihilator Ann(f^s)...");
2865  def A = Sannfs(f);
2866  setring A;
2867  dbprint(ppl,"// ...done");
2868  dbprint(ppl-1,"//    Got: " + string(LD));
2869  poly f = imap(save,f);
2870  if (!l0given)
2871  {
2872    dbprint(ppl,"// Computing b-function of given poly...");
2873    ideal LDf = LD,f;
2874    LDf = engine(LDf,eng);
2875    vector v = pIntersect(var(2*n+1),LDf);   // BS poly of f
2876    list BS = bFactor(vec2poly(v));
2877    dbprint(ppl,"// ...done");
2878    dbprint(ppl-1,"// roots: " + string(BS[1]));
2879    dbprint(ppl-1,"// multiplicities: " + string(BS[2]));
2880    BS = intRoots(BS);
2881    intvec iv;
2882    for (i=1; i<=ncols(BS[1]); i++)
2883    {
2884      iv[i] = int(BS[1][i]);
2885    }
2886    l0 = Min(iv);
2887    kill v,iv,BS,LDf;
2888  }
2889  dbprint(ppl,"// Computing Ann(f^" + string(l0) + ")...");
2890  LD = annfspecial(LD,f,l0,l0); // Ann(f^l0)
2891  // create new ring without s
2892  list RL = ringlist(A);
2893  RL = RL[1..4];
2894  list Lt = RL[2];
2895  Lt = delete(Lt,2*n+1);
2896  RL[2] = Lt;
2897  Lt = RL[3];
2898  Lt = delete(Lt,2);
2899  RL[3] = Lt;
2900  def @B = ring(RL);
2901  setring @B;
2902  def B = Weyl();
2903  setring B;
2904  kill @B;
2905  ideal LD = imap(A,LD);
2906  LD = engine(LD,eng);
2907  dbprint(ppl,"// ...done");
2908  dbprint(ppl-1,"//    Got: " + string(LD));
2909  kill A;
2910  intvec w = 1:n;
2911  ideal DR = deRhamCohomIdeal(LD,w,eng);
2912  export(DR);
2913  setring save;
2914  return(B);
2915}
2916example
2917{
2918  "EXAMPLE:"; echo = 2;
2919  ring r = 0,(x,y,z),dp;
2920  poly f = x^3+y^3+z^3;
2921  def A = deRhamCohom(f); // we see that the K-dim is 2
2922  setring A;
2923  DR;
2924}
2925
2926// Appel hypergeometric functions /////////////////////////////////////////////
2927
2928proc appelF1()
2929  "USAGE:  appelF1();
2930RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel1'
2931PURPOSE: defines the ideal in a parametric Weyl algebra,
2932@*       which annihilates Appel F1 hypergeometric function
2933NOTE:    The output ring is a parametric Weyl algebra. It contains an ideal
2934@*       'IAappel1' annihilating Appel F1 hypergeometric function.
2935@*       See (SST) p. 48.
2936EXAMPLE: example appelF1; shows example
2937"
2938{
2939  // Appel F1, d = b', SST p.48
2940  ring @r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
2941  def @S = Weyl();
2942  setring @S;
2943  ideal IAppel1 =
2944    (x*Dx)*(x*Dx+y*Dy+c-1) - x*(x*Dx+y*Dy+a)*(x*Dx+b),
2945    (y*Dy)*(x*Dx+y*Dy+c-1) - y*(x*Dx+y*Dy+a)*(y*Dy+d),
2946    (x-y)*Dx*Dy - d*Dx + b*Dy;
2947  export IAppel1;
2948  kill @r;
2949  return(@S);
2950}
2951example
2952{
2953  "EXAMPLE:"; echo = 2;
2954  def A = appelF1();
2955  setring A;
2956  IAppel1;
2957}
2958
2959proc appelF2()
2960  "USAGE:  appelF2();
2961RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel2'
2962PURPOSE: defines the ideal in a parametric Weyl algebra,
2963@*       which annihilates Appel F2 hypergeometric function
2964NOTE:    The output ring is a parametric Weyl algebra. It contains an ideal
2965@*       'IAappel2' annihilating Appel F2 hypergeometric function.
2966@*       See (SST) p. 85.
2967EXAMPLE: example appelF2; shows example
2968"
2969{
2970  // Appel F2, c = b', SST p.85
2971  ring @r = (0,a,b,c),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
2972  def @S = Weyl();
2973  setring @S;
2974  ideal IAppel2 =
2975    (x*Dx)^2 - x*(x*Dx+y*Dy+a)*(x*Dx+b),
2976    (y*Dy)^2 - y*(x*Dx+y*Dy+a)*(y*Dy+c);
2977  export IAppel2;
2978  kill @r;
2979  return(@S);
2980}
2981example
2982{
2983  "EXAMPLE:"; echo = 2;
2984  def A = appelF2();
2985  setring A;
2986  IAppel2;
2987}
2988
2989proc appelF4()
2990  "USAGE:  appelF4();
2991RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel4'
2992PURPOSE: defines the ideal in a parametric Weyl algebra,
2993@*       which annihilates Appel F4 hypergeometric function
2994NOTE:    The output ring is a parametric Weyl algebra. It contains an ideal
2995@*       'IAappel4' annihilating Appel F4 hypergeometric function.
2996@*       See (SST) p. 39.
2997EXAMPLE: example appelF4; shows example
2998"
2999{
3000  // Appel F4, d = c', SST, p. 39
3001  ring @r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
3002  def @S = Weyl();
3003  setring @S;
3004  ideal IAppel4 =
3005    Dx*(x*Dx+c-1) - (x*Dx+y*Dy+a)*(x*Dx+y*Dy+b),
3006    Dy*(y*Dy+d-1) - (x*Dx+y*Dy+a)*(x*Dx+y*Dy+b);
3007  export IAppel4;
3008  kill @r;
3009  return(@S);
3010}
3011example
3012{
3013  "EXAMPLE:"; echo = 2;
3014  def A = appelF4();
3015  setring A;
3016  IAppel4;
3017}
3018
3019
3020// characteric variety ////////////////////////////////////////////////////////
3021
3022proc charVariety(ideal I, list #)
3023  "USAGE:  charVariety(I [,eng]); I an ideal, eng an optional int
3024RETURN:  ring (commutative) containing an ideal 'charVar'
3025PURPOSE: computes an ideal whose zero set is the characteristic variety of I in
3026@*       the sense of D-module theory
3027ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
3028@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
3029@*       holds, i.e. the sequence of variables is given by
3030@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
3031@*       belonging to x(i).
3032NOTE:    The output ring is commutative. It contains an ideal 'charVar'.
3033@*       If eng<>0, @code{std} is used for Groebner basis computations,
3034@*       otherwise, and by default, @code{slimgb} is used.
3035DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
3036@*       if @code{printlevel}>=2, all the debug messages will be printed.
3037SEE ALSO: charInfo
3038EXAMPLE: example charVariety; shows examples
3039"
3040{
3041  // assumption check is done in GBWeight
3042  int eng;
3043  if (size(#)>0)
3044  {
3045    if (intLike(#[1]))
3046    {
3047      eng = int(#[1]);
3048    }
3049  }
3050  int ppl = printlevel - voice + 2;
3051  def save = basering;
3052  int n = nvars(save)/2;
3053  intvec u = 0:n;
3054  intvec v = 1:n;
3055  dbprint(ppl,"// Computing Groebner basis wrt weights...");
3056  I = GBWeight(I,u,v,eng);
3057  dbprint(ppl,"// ...finished");
3058  dbprint(ppl-1,"//    Got: " + string(I));
3059  list RL = ringlist(save);
3060  RL = RL[1..4];
3061  def newR = ring(RL);
3062  setring newR;
3063  ideal charVar = imap(save,I);
3064  intvec uv = u,v;
3065  charVar = inForm(charVar,uv);
3066  charVar = groebner(charVar);
3067  export(charVar);
3068  setring save;
3069  return(newR);
3070}
3071example
3072{
3073  "EXAMPLE:"; echo = 2;
3074  ring r = 0,(x,y),Dp;
3075  poly F = x3-y2;
3076  printlevel = 0;
3077  def A  = annfs(F);
3078  setring A;      // Weyl algebra
3079  LD;             // the annihilator of F
3080  def CA = charVariety(LD);
3081  setring CA; CA; // commutative ring
3082  charVar;
3083  dim(charVar);   // hence I is holonomic
3084}
3085
3086proc charInfo(ideal I)
3087  "USAGE:  charInfo(I);  I an ideal
3088RETURN:  ring (commut.) containing ideals 'charVar','singLoc' and list 'primDec'
3089PURPOSE: computes characteristic variety of I (in the sense of D-module theory),
3090@*       its singular locus and primary decomposition
3091ASSUME:  The basering is the n-th Weyl algebra over a field of characteristic 0
3092@*       and for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
3093@*       holds, i.e. the sequence of variables is given by
3094@*       x(1),...,x(n),D(1),...,D(n), where D(i) is the differential operator
3095@*       belonging to x(i).
3096NOTE:    In the output ring, which is commutative:
3097@*       - the ideal 'charVar' is the characteristic variety char(I),
3098@*       - the ideal 'SingLoc' is the singular locus of char(I),
3099@*       - the list 'primDec' is the primary decomposition of char(I).
3100DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
3101@*       if @code{printlevel}>=2, all the debug messages will be printed.
3102EXAMPLE: example charInfo; shows examples
3103"
3104{
3105  int ppl = printlevel - voice + 2;
3106  def save = basering;
3107  dbprint(ppl,"// computing characteristic variety...");
3108  def A = charVariety(I);
3109  setring A;
3110  dbprint(ppl,"// ...done");
3111  dbprint(ppl-1,"//    Got: " + string(charVar));
3112  dbprint(ppl,"// computing singular locus...");
3113  ideal singLoc = slocus(charVar);
3114  singLoc = groebner(singLoc);
3115  dbprint(ppl,"// ...done");
3116  dbprint(ppl-1,"//    Got: " + string(singLoc));
3117  dbprint(ppl,"// computing primary decomposition...");
3118  list primDec = primdecGTZ(charVar);
3119  dbprint(ppl,"// ...done");
3120  //export(charVar,singLoc,primDec);
3121  export(singLoc,primDec);
3122  setring save;
3123  return(A);
3124}
3125example
3126{
3127  "EXAMPLE:"; echo = 2;
3128  ring r = 0,(x,y),Dp;
3129  poly F = x3-y2;
3130  printlevel = 0;
3131  def A  = annfs(F);
3132  setring A;      // Weyl algebra
3133  LD;             // the annihilator of F
3134  def CA = charInfo(LD);
3135  setring CA; CA; // commutative ring
3136  charVar;        // characteristic variety
3137  singLoc;        // singular locus
3138  primDec;        // primary decomposition
3139}
3140
3141
3142// examples ///////////////////////////////////////////////////////////////////
3143
3144/*
3145  static proc exCusp()
3146  {
3147  "EXAMPLE:"; echo = 2;
3148  ring r = 0,(x,y,Dx,Dy),dp;
3149  def R = Weyl();   setring R;
3150  poly F = x2-y3;
3151  ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
3152  def W = SDLoc(I,F);
3153  setring W;
3154  LD;
3155  def U = DLoc0(LD,x2-y3);
3156  setring U;
3157  LD0;
3158  BS;
3159  // the same with DLoc:
3160  setring R;
3161  DLoc(I,F);
3162  }
3163
3164  static proc exWalther1()
3165  {
3166  // p.18 Rem 3.10
3167  ring r = 0,(x,Dx),dp;
3168  def R = nc_algebra(1,1);
3169  setring R;
3170  poly F = x;
3171  ideal I = x*Dx+1;
3172  def W = SDLoc(I,F);
3173  setring W;
3174  LD;
3175  ideal J = LD, x;
3176  eliminate(J,x*Dx); // must be [1]=s // agree!
3177  // the same result with Dloc0:
3178  def U = DLoc0(LD,x);
3179  setring U;
3180  LD0;
3181  BS;
3182  }
3183
3184  static proc exWalther2()
3185  {
3186  // p.19 Rem 3.10 cont'd
3187  ring r = 0,(x,Dx),dp;
3188  def R = nc_algebra(1,1);
3189  setring R;
3190  poly F = x;
3191  ideal I = (x*Dx)^2+1;
3192  def W = SDLoc(I,F);
3193  setring W;
3194  LD;
3195  ideal J = LD, x;
3196  eliminate(J,x*Dx); // must be [1]=s^2+2*s+2 // agree!
3197  // the same result with Dloc0:
3198  def U = DLoc0(LD,x);
3199  setring U;
3200  LD0;
3201  BS;
3202  // almost the same with DLoc
3203  setring R;
3204  DLoc(I,F);
3205  }
3206
3207  static proc exWalther3()
3208  {
3209  // can check with annFs too :-)
3210  // p.21 Ex 3.15
3211  LIB "nctools.lib";
3212  ring r = 0,(x,y,z,w,Dx,Dy,Dz,Dw),dp;
3213  def R = Weyl();
3214  setring R;
3215  poly F = x2+y2+z2+w2;
3216  ideal I = Dx,Dy,Dz,Dw;
3217  def W = SDLoc(I,F);
3218  setring W;
3219  LD;
3220  ideal J = LD, x2+y2+z2+w2;
3221  eliminate(J,x*y*z*w*Dx*Dy*Dz*Dw); // must be [1]=s^2+3*s+2 // agree
3222  ring r2 =  0,(x,y,z,w),dp;
3223  poly F = x2+y2+z2+w2;
3224  def Z = annfs(F);
3225  setring Z;
3226  LD;
3227  BS;
3228  // the same result with Dloc0:
3229  setring W;
3230  def U = DLoc0(LD,x2+y2+z2+w2);
3231  setring U;
3232  LD0;  BS;
3233  // the same result with DLoc:
3234  setring R;
3235  DLoc(I,F);
3236  }
3237
3238  static proc ex_annRat()
3239  {
3240  // more complicated example for annRat
3241  ring r = 0,(x,y,z),dp;
3242  poly f = x3+y3+z3; // mir = -2
3243  poly g = x*y*z;
3244  def A = annRat(g,f);
3245  setring A;
3246  }
3247*/
Note: See TracBrowser for help on using the repository browser.