# source:git/Singular/LIB/dmodapp.lib@66d68c

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