source: git/Singular/LIB/dmodapp.lib @ 33b509

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