source: git/Singular/LIB/dmodapp.lib @ 0af200

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