# source:git/Singular/LIB/dmodvar.lib@6a07eb

spielwiese
Last change on this file since 6a07eb was 6a07eb, checked in by Viktor Levandovskyy <levandov@…>, 13 years ago
• Property mode set to 100755
File size: 22.9 KB
Line
1////////////////////////////////////////////////////////////////////////
2version="\$Id\$";
3category="Noncommutative";
4info="
5LIBRARY: dmodvar.lib     Algebraic D-modules for varieties
6
7AUTHORS: Daniel Andres,          daniel.andres@math.rwth-aachen.de
8       Viktor Levandovskyy,    levandov@math.rwth-aachen.de
9       Jorge Martin-Morales,   jorge@unizar.es
10
11THEORY: Let K be a field of characteristic 0. Given a polynomial ring R = K[x_1,...,x_n] and
12@* a set of polynomial f_1,..., f_r in R, define F = f_1 * ... * f_r and F^s:=f_1^s_1*...*f_r^s_r
13@* for symbolic discrete (that is shiftable) variables s_1,..., s_r.
14@* The module R[1/F]*F^s has a structure of a D<S>-module, where D<S> := D(R) tensored with S over K, where
15@*   - D(R) is an n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>
16@*   - S is the universal enveloping algebra of gl_r, generated by s_{ij}, where s_{ii}=s_i.
17@* One is interested in the following data:
18@*   - the left ideal Ann F^s in D<S>, usually denoted by LD in the output
19@*   - global Bernstein polynomial in one variable s = s_1 + ...+ s_r, denoted by bs,
20@*   - its minimal integer root s0, the list of all roots of bs, which are known
21@*     to be negative rational numbers, with their multiplicities, which is denoted by BS
22@*   - an r-tuple of operators in D<S>, denoted by PS, such that the functional equality
23@*     sum(k=1 to k=r) P_k*f_k*F^s = bs*F^s holds in R[1/F]*F^s.
24
25REFERENCES:
26  (BMS06) Budur, Mustata, Saito: Bernstein-Sato polynomials of arbitrary varieties (2006).
27  (ALM09) Andres, Levandovskyy, Martin-Morales : Principal Intersection and Bernstein-Sato Polynomial of an Affine Variety (2009).
28
29MAIN PROCEDURES:
30bfctVarIn(F[,L]);     compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using initial ideal approach
31bfctVarAnn(F[,L]);  compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using Sannfs approach
32SannfsVar(F[,O,e]); compute the annihilator of F^s in the ring D<S>
33
34AUXILIARY PROCEDURES:
35makeIF(F[,ORD]);    create the Malgrange ideal, associated with F = F[1],..,F[P]
36
38
39KEYWORDS: D-module; D-module structure; Bernstein-Sato polynomial for variety; global Bernstein-Sato polynomial for variety;
40Weyl algebra; parametric annihilator for variety; Budur-Mustata-Saito approach; initial ideal approach;
41";
42
43// Static procs:
44//coDim(I);           compute the codimension of the leading ideal of I
45// dmodvarAssumeViolation()
46// ORDstr2list (ORD, NN)
47// smallGenCoDim(I,k)
48
49
50LIB "bfun.lib";    // for pIntersect
51LIB "dmodapp.lib"; // for isCommutative etc.
52
53
54///////////////////////////////////////////////////////////////////////////////
55
56// testing for consistency of the library:
57proc testdmodvarlib ()
58{
59  "AUXILIARY PROCEDURES:";
60  example makeIF;
61  "MAIN PROCEDURES:";
62  example bfctVarIn;
63  example bfctVarAnn;
64  example SannfsVar;
65}
66
67//   example coDim;
68
69///////////////////////////////////////////////////////////////////////////////
70
71static proc dmodvarAssumeViolation()
72{
73  // returns Boolean : yes/no [for assume violation]
74  // char K = 0
75  // no qring
76  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
77  {
78    //    "ERROR: no qring is allowed";
79    return(1);
80  }
81  return(0);
82}
83
84// da: in smallGenCoDim(), rewritten using mstd business
85static proc coDim (ideal I)
86"USAGE:  coDim (I);  I an ideal
87RETURN:  int
88PURPOSE: computes the codimension of the ideal generated by the leading monomials
89@*       of the given generators of the ideal. This is also the codimension of
90@*       the ideal if it is represented by a standard basis.
91NOTE:    The codimension of an ideal I means the number of variables minus the
92@*       Krull dimension of the basering modulo I.
93EXAMPLE: example SannfsVar; shows examples
94"
95{
96  int n = nvars(basering);
97  int d = dim(I); // to insert: check whether I is in GB
98  return(n-d);
99}
100example
101{
102  "EXAMPLE:"; echo = 2;
103  ring R = 0,(x,y,z),Dp;
104  ideal I = x^2+y^3, z;
105  coDim(std(I));
106}
107
108static proc ORDstr2list (string ORD, int NN)
109{
110  /* convert an ordering defined in NN variables in the  */
111  /* string form into the same ordering in the list form */
112  string st;
113  st = "ring @Z = 0,z(1.." + string(NN) + "),";
114  st = st + ORD + ";";
115  execute(st); kill st;
116  list L = ringlist(@Z)[3];
117  kill @Z;
118  return(L);
119}
120
121proc SannfsVar (ideal F, list #)
122"USAGE:  SannfsVar(F [,ORD,eng]);  F an ideal, ORD an optional string, eng an optional int
123RETURN:  ring
124PURPOSE: compute the D<S>-module structure of D<S>*f^s where f = F[1]*..*F[P]
125and D<S> is the Weyl algebra D tensored with K<S>=U(gl_P), according to the
126generalized algorithm by Briancon and Maisonobe for affine varieties.
127NOTE:    activate this ring with the @code{setring} command.
128@*       In the ring D<S>, the ideal LD is the needed D<S>-module structure.
129@*       The value of ORD must be an elimination ordering in D<Dt,S> for Dt
130@*       written in the string form, otherwise the result may have no meaning.
131@*       By default ORD = '(a(1..(P)..1),a(1..(P+P^2)..1),dp)'.
132@*       If eng<>0, @code{std} is used for Groebner basis computations,
133@*       otherwise, and by default @code{slimgb} is used.
134@*       If printlevel=1, progress debug messages will be printed,
135@*       if printlevel>=2, all the debug messages will be printed.
136EXAMPLE: example SannfsVar; shows examples
137"
138{
139  if (dmodvarAssumeViolation())
140  {
141    ERROR("Basering is inappropriate: characteristic>0 or qring present");
142  }
143  if (!isCommutative())
144  {
145    ERROR("Basering must be commutative");
146  }
147  def save = basering;
148  int N = nvars(basering);
149  int P = ncols(F);  //ncols better than size, since F[i] could be zero
150  // P is needed for default ORD
151  int i,j,k,l;
152  // st = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";
153  string st = "(a(" + string(1:P);
154  st = st + "),a(" + string(1:(P+P^2));
155  st = st + "),dp)";
156  // default values
157  string ORD = st;
158  int eng = 0;
159  if ( size(#)>0 )
160  {
161    if ( typeof(#[1]) == "string" )
162    {
163      ORD = string(#[1]);
164      // second arg
165      if (size(#)>1)
166      {
167        // exists 2nd arg
168        if ( typeof(#[2]) == "int" )
169        {
170          // the case: given ORD, given engine
171          eng = int(#[2]);
172        }
173        else
174        {
175          eng = 0;
176        }
177      }
178      else
179      {
180        // no second arg
181        eng = 0;
182      }
183    }
184    else
185    {
186      if ( typeof(#[1]) == "int" )
187      {
188        // the case: default ORD, engine given
189        eng = int(#[1]);
190        // ORD = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";  //is already set
191      }
192      else
193      {
194        // incorr. 1st arg
195        ORD = string(st);
196      }
197    }
198  }
199  // size(#)=0, i.e. there is no elimination ordering and no engine given
200  // eng = 0; ORD = "(a(1..(P)..1),a(1..(P^2)..1),dp)";  //are already set
201  int ppl = printlevel-voice+2;
202  // returns a list with a ring and an ideal LD in it
203  // save, N, P and the indices are already defined
204  int Nnew = 2*N+P+P^2;
205  list RL = ringlist(basering);
206  list L;
207  L[1] = RL[1];  //char
208  L[4] = RL[4];  //char, minpoly
209  // check whether vars have admissible names
210  list Name  = RL[2];
211  list RName;
212  // (i,j) <--> (i-1)*p+j
213  for (i=1; i<=P; i++)
214  {
215    RName[i] = "Dt("+string(i)+")";
216    for (j=1; j<=P; j++)
217    {
218      st = "s("+string(i)+")("+string(j)+")";
219      RName[P+(i-1)*P+j] = st;
220    }
221  }
222  for(i=1; i<=N; i++)
223  {
224    for(j=1; j<=size(RName); j++)
225    {
226      if (Name[i] == RName[j])
227      {
228        ERROR("Variable names should not include Dt(i), s(i)(j)");
229      }
230    }
231  }
232  // now, create the names for new vars
233  list DName;
234  for(i=1; i<=N; i++)
235  {
236    DName[i] = "D"+Name[i];  //concat
237  }
238  list NName = RName + Name + DName;
239  L[2] = NName;
240  // Name, Dname will be used further
241  kill NName;
242  //block ord (a(1..(P)..1),a(1..(P+P^2)..1),dp);
243  //export Nnew;
244  L[3] = ORDstr2list(ORD,Nnew);
245  // we are done with the list
246  def @R@ = ring(L);
247  setring @R@;
248  matrix @D[Nnew][Nnew];
249  // kronecker(i,j) equals (i==j)
250  // (i,j) <--> (i-1)*p+j
251  for (i=1; i<=P; i++)
252  {
253    for (j=1; j<=P; j++)
254    {
255      for (k=1; k<=P; k++)
256      {
257        //[sij,Dtk] = djk*Dti
258        @D[k,P+(i-1)*P+j] = (j==k)*Dt(i);
259        for (l=1; l<=P; l++)
260        {
261          if ( (i-k)*P < l-j )
262          {
263            //[sij,skl] = djk*sil - dil*skj
264            @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*s(i)(l) + (i==l)*s(k)(j);
265          }
266        }
267      }
268    }
269  }
270  for (i=1; i<=N; i++)
271  {
272    //[Dx,x]=1
273    @D[P+P^2+i,P+P^2+N+i] = 1;
274  }
275  def @R = nc_algebra(1,@D);
276  setring @R;
277  //@R@ will be used further
278  dbprint(ppl,"// -1-1- the ring @R(_Dt,_s,_x,_Dx) is ready");
279  dbprint(ppl-1, @R);
280  // create the ideal I
281  // (i,j) <--> (i-1)*p+j
282  ideal  F = imap(save,F);
283  ideal I;
284  for (i=1; i<=P; i++)
285  {
286    for (j=1; j<=P; j++)
287    {
288      I[(i-1)*P+j] = Dt(i)*F[j] + s(i)(j);
289    }
290  }
291  poly p,q;
292  for (i=1; i<=N; i++)
293  {
294    p=0;
295    for (j=1; j<=P; j++)
296    {
297      q = Dt(j);
298      q = q*diff(F[j],var(P+P^2+i));
299      p = p + q;
300    }
301    I = I, p + var(P+P^2+N+i);
302  }
303  // -------- the ideal I is ready ----------
304  dbprint(ppl,"// -1-2- starting the elimination of "+string(Dt(1..P))+" in @R");
305  dbprint(ppl-1, I);
306  ideal J = engine(I,eng);
307  ideal K = nselect(J,1..P);
308  kill I,J;
309  dbprint(ppl,"// -1-3- all Dt(i) are eliminated");
310  dbprint(ppl-1, K);  //K is without Dt(i)
311  // ----------- the ring @R2(_s,_x,_Dx) ------------
312  //come back to the ring save, recover L and remove all Dt(i)
313  //L[1],L[4] do not change
314  setring save;
315  list Lord, tmp;
316  // variables
317  tmp = L[2];
318  Lord = tmp[P+1..Nnew];
319  L[2] = Lord;
320  // ordering
321  // st = "(a(1..(P^2)..1),dp)";
322  st = "(a(" + string(1:P^2);
323  st = st + "),dp)";
324  tmp = ORDstr2list(st,Nnew-P);
325  L[3] = tmp;
326  def @R2@ = ring(L);
327  kill L;
328  // we are done with the list
329  setring @R2@;
330  matrix tmpM,LordM;
331  // non-commutative relations
332  intvec iv = P+1..Nnew;
333  tmpM = imap(@R@,@D);
334  kill @R@;
335  LordM = submat(tmpM,iv,iv);
336  matrix @D2 = LordM;
337  def @R2 = nc_algebra(1,@D2);
338  setring @R2;
339  kill @R2@;
340  dbprint(ppl,"// -2-1- the ring @R(_s,_x,_Dx) is ready");
341  dbprint(ppl-1, @R2);
342  ideal K = imap(@R,K);
343  kill @R;
344  dbprint(ppl,"// -2-2- starting cosmetic Groebner basis computation");
345  dbprint(ppl-1, K);
346  K = engine(K,eng);
347  dbprint(ppl,"// -2-3- the cosmetic Groebner basis has been computed");
348  dbprint(ppl-1,K);
349  ideal LD = K;
350  attrib(LD,"isSB",1);
351  export LD;
352  return(@R2);
353}
354example
355{
356  "EXAMPLE:"; echo = 2;
357  ring R = 0,(x,y),Dp;
358  ideal F = x^3, y^5;
359  //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
360  //eng = 0;
361  def A = SannfsVar(F);
362  setring A;
363  LD;
364}
365
366proc bfctVarAnn (ideal F, list #)
367"USAGE:  bfctVarAnn(F[,gid,eng]); F an ideal, gid,eng optional ints
368RETURN:  list of an ideal and an intvec
369PURPOSE: computes the roots of the Bernstein-Sato polynomial and their multiplicities
370@*       for an affine algebraic variety defined by F = F[1],..,F[r].
371ASSUME:  The basering is a commutative polynomial ring in char 0.
372BACKGROUND: In this proc, the annihilator of f^s in D[s] is computed and then a
373@*       system of linear equations is solved by linear reductions in order to
374@*       find the minimal polynomial of S = s(1)(1) + ... + s(P)(P)
375NOTE:    In the output list, the ideal contains all the roots and the intvec their multiplicities.
376@*       If gid<>0, the ideal is used as given. Otherwise, and by default, a
377@*       heuristically better suited generating set is used.
378@*       If eng<>0, @code{std} is used for GB computations,
379@*       otherwise, and by default, @code{slimgb} is used.
380DISPLAY: If printlevel=1, progress debug messages will be printed,
381@*       if printlevel=2, all the debug messages will be printed.
382COMPUTATIONAL REMARK: The time of computation can be very different depending
383@*       on the chosen generators of F, although the result is always the same.
384EXAMPLE: example bfctVarAnn; shows examples
385"
386{
387  if (dmodvarAssumeViolation())
388  {
389    ERROR("Basering is inappropriate: characteristic>0 or qring present");
390  }
391  if (!isCommutative())
392  {
393    ERROR("Basering must be commutative");
394  }
395  int gid = 0; // default
396  int eng = 0; // default
397  if (size(#)>0)
398  {
399    if (typeof(#[1])=="int" || typeof(#[1])=="number")
400    {
401      gid = int(#[1]);
402    }
403    if (size(#)>1)
404    {
405      if (typeof(#[2])=="int" || typeof(#[2])=="number")
406      {
407        eng = int(#[2]);
408      }
409    }
410  }
411  def save = basering;
412  int ppl = printlevel - voice + 2;
413  printlevel = printlevel+1;
414  list L = smallGenCoDim(F,gid);
415  F = L[1];
416  int cd = L[2];
417  kill L;
418  def @R2 = SannfsVar(F,eng);
419  printlevel = printlevel-1;
420  setring @R2;
421  // we are in D[s] and LD is a std of SannfsVar(F)
422  ideal F = imap(save,F);
423  ideal GF = std(F);
424  ideal J = NF(LD,GF);
425  J = J, F;
426  dbprint(ppl,"// -3-1- starting Groebner basis of ann F^s + F ");
427  dbprint(ppl-1,J);
428  ideal K = engine(J,eng);
429  dbprint(ppl,"// -3-2- finished Groebner basis of ann F^s + F ");
430  dbprint(ppl-1,K);
431  poly S;
432  int i;
433  for (i=1; i<=size(F); i++)
434  {
435    S = S + s(i)(i);
436  }
437  dbprint(ppl,"// -4-1- computing the minimal polynomial of S");
438  //dbprint(ppl-1,"S = "+string(S));
439  module M = pIntersect(S,K);
440  dbprint(ppl,"// -4-2- the minimal polynomial has been computed");
441  //dbprint(ppl-1,M);
442  ring @R3 = 0,s,dp;
443  dbprint(ppl,"// -5-1- the ring @R3(s) is ready");
444  dbprint(ppl-1,@R3);
445  ideal M = imap(@R2,M);
446  //kill @R2;
447  poly p;
448  for (i=1; i<=size(M); i++)
449  {
450    p = p + M[i]*s^(i-1);
451  }
452  dbprint(ppl,"// -5-2- factorization of the minimal polynomial");
453  list P = factorize(p);          //with constants and multiplicities
454  dbprint(ppl-1,P);               //the Bernstein polynomial is monic,
455  ideal bs; intvec m;             //so we are not interested in constants
456  for (i=2; i<=ncols(P[1]); i++)   //and that is why we delete P[1][1] and P[2][1]
457  {
458    bs[i-1] = P[1][i];
459    m[i-1]  = P[2][i];
460  }
461  // convert factors to a list of their roots and multiplicities
462  bs =  normalize(bs);
463  bs = -subst(bs,s,0);
464  setring save;
465//   ideal GF = imap(@R2,GF);
466//   attrib(GF,"isSB",1);
467  kill @R2;
468  dbprint(ppl,"// -5-3- codimension of the variety");
469//   int cd = coDim(GF);
470  dbprint(ppl-1,cd);
471  ideal bs = imap(@R3,bs);
472  dbprint(ppl,"// -5-4- shifting BS(s)=minpoly(s-codim+1)");
473  for (i=1; i<=ncols(bs); i++)
474  {
475    bs[i] = bs[i] + cd - 1;
476  }
477  kill @R3;
478  list BS = bs,m;
479  return(BS);
480}
481example
482{
483  "EXAMPLE:"; echo = 2;
484  ring R = 0,(x,y,z),Dp;
485  ideal F = x^2+y^3, z;
486  bfctVarAnn(F);
487}
488
489proc makeIF (ideal F, list #)
490"USAGE:  makeIF(F [,ORD]);  F an ideal, ORD an optional string
491RETURN:  ring
492PURPOSE: create the ideal by Malgrange associated with F = F[1],..,F[P].
493NOTE:    activate this ring with the @code{setring} command. In this ring,
494@*       - the ideal IF is the ideal by Malgrange corresponding to F.
495@*       The value of ORD must be an arbitrary ordering in K<_t,_x,_Dt,_Dx>
496@*       written in the string form. By default ORD = 'dp'.
497@*       If printlevel=1, progress debug messages will be printed,
498@*       if printlevel>=2, all the debug messages will be printed.
499EXAMPLE: example makeIF; shows examples
500"
501{
502  string ORD = "dp";
503  if ( size(#)>0 )
504  {
505    if ( typeof(#[1]) == "string" )
506    {
507      ORD = string(#[1]);
508    }
509  }
510  int ppl = printlevel-voice+2;
511  def save = basering;
512  int N = nvars(save);
513  int P = ncols(F);
514  int Nnew = 2*P+2*N;
515  int i,j;
516  string st;
517  list RL = ringlist(save);
518  list L,Lord;
519  list tmp;
520  intvec iv;
521  L[1] = RL[1];
522  L[4] = RL[4];
523  //check whether vars have admissible names
524  list Name = RL[2];
525  list TName, DTName;
526  for (i=1; i<=P; i++)
527  {
528    TName[i] = "t("+string(i)+")";
529    DTName[i] = "Dt("+string(i)+")";
530  }
531  for (i=1; i<=N; i++)
532  {
533    for (j=1; j<=P; j++)
534    {
535      if (Name[i] == TName[j])
536      {
537        ERROR("Variable names should not include t(i)");
538      }
539    }
540  }
541  //now, create the names for new vars
542  list DName;
543  for (i=1; i<=N; i++)
544  {
545    DName[i] = "D"+Name[i];  //concat
546  }
547  list NName = TName + Name + DTName + DName;
548  L[2]   = NName;
549  // Name, Dname will be used further
550  kill NName, TName, Name, DTName, DName;
551  // ORD already set, default ord dp;
552  L[3] = ORDstr2list(ORD,Nnew);
553  // we are done with the list
554  def @R@ = ring(L);
555  setring @R@;
556  matrix @D[Nnew][Nnew];
557  for (i=1; i<=N+P; i++)
558  {
559    @D[i,i+N+P]=1;
560  }
561  def @R = nc_algebra(1,@D);
562  setring @R;
563  kill @R@;
564  //dbprint(ppl,"// -1-1- the ring @R(_t,_x,_Dt,_Dx) is ready");
565  //dbprint(ppl-1, @R);
566  // create the ideal I
567  ideal  F = imap(save,F);
568  ideal I;
569  for (j=1; j<=P; j++)
570  {
571    I[j] = t(j) - F[j];
572  }
573  poly p,q;
574  for (i=1; i<=N; i++)
575  {
576    p=0;
577    for (j=1; j<=P; j++)
578    {
579      q = Dt(j);
580      q = diff(F[j],var(P+i))*q;
581      p = p + q;
582    }
583    I = I, p + var(2*P+N+i);
584  }
585  // -------- the ideal I is ready ----------
586  ideal IF = I;
587  export IF;
588  return(@R);
589}
590example
591{
592  "EXAMPLE:"; echo = 2;
593  ring R = 0,(x,y,z),Dp;
594  ideal I = x^2+y^3, z;
595  def W = makeIF(I);
596  setring W;
597  IF;
598}
599
600proc bfctVarIn (ideal I, list #)
601"USAGE:  bfctVarIn(I [,a,b,c]);  I an ideal, a,b,c optional ints
602RETURN:  list of ideal and intvec
603PURPOSE: computes the roots of the Bernstein-Sato polynomial and their
604@*       multiplicities for an affine algebraic variety defined by I.
605ASSUME:  The basering is commutative and of characteristic 0.
606@*       Varnames of the basering do not include t(1),...,t(r) and
607@*       Dt(1),...,Dt(r), where r is the number of entries of the input ideal.
608BACKGROUND:  In this proc, the initial ideal of the multivariate Malgrange ideal
609@*       defined by I is computed and then a system of linear equations is solved
610@*       by linear reductions following the ideas by Noro.
611NOTE:    In the output list, say L,
612@*       - L[1] of type ideal contains all the rational roots of a b-function,
613@*       - L[2] of type intvec contains the multiplicities of above roots,
614@*       - optional L[3] of type string is the part of b-function without
615@*        rational roots.
616@*       Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and
617@*       L[3] is 1 (for nonzero constant) or 0 (for zero b-function).
618@*       If a<>0, the ideal is used as given. Otherwise, and by default, a
619@*       heuristically better suited generating set is used to reduce computation
620@*       time.
621@*       If b<>0, @code{std} is used for GB computations in characteristic 0,
622@*       otherwise, and by default, @code{slimgb} is used.
623@*       If c<>0, a matrix ordering is used for GB computations, otherwise,
624@*       and by default, a block ordering is used.
625DISPLAY: If printlevel=1, progress debug messages will be printed,
626@*       if printlevel>=2, all the debug messages will be printed.
627EXAMPLE: example bfctVarIn; shows examples
628"
629{
630  if (dmodvarAssumeViolation())
631  {
632    ERROR("Basering is inappropriate: characteristic>0 or qring present");
633  }
634  if (!isCommutative())
635  {
636    ERROR("Basering must be commutative");
637  }
638  int ppl = printlevel - voice + 2;
639  int idealasgiven  = 0; // default
640  int whicheng      = 0; // default
641  int whichord      = 0; // default
642  if (size(#)>0)
643  {
644    if (typeof(#[1])=="int" || typeof(#[1])=="number")
645    {
646      idealasgiven = int(#[1]);
647    }
648    if (size(#)>1)
649    {
650      if (typeof(#[2])=="int" || typeof(#[2])=="number")
651      {
652        whicheng = int(#[2]);
653      }
654      if (size(#)>2)
655      {
656        if (typeof(#[3])=="int" || typeof(#[3])=="number")
657        {
658          whichord = int(#[3]);
659        }
660      }
661    }
662  }
663  def save = basering;
664  int i;
665  int n = nvars(basering);
666  // step 0: get small generating set
667  I = simplify(I,2);
668  list L = smallGenCoDim(I,idealasgiven);
669  I = L[1];
670  int c = L[2];
671  kill L;
672  // step 1: setting up the multivariate Malgrange ideal
673  int r = size(I);
674  def D = makeIF(I);
675  setring D;
676  dbprint(ppl-1,"// Computing in " + string(n+r) + "-th Weyl algebra:", D);
677  dbprint(ppl-1,"// The Malgrange ideal: ", IF);
678  // step 2: compute the b-function of the Malgrange ideal w.r.t. approriate weights
679  intvec w = 1:r;
680  w[r+n] = 0;
681  dbprint(ppl,"// Computing the b-function of the Malgrange ideal...");
682  list L = bfctIdeal(IF,w,whicheng,whichord);
683  dbprint(ppl,"// ... done.");
684  dbprint(ppl-1,"// The b-function: ",L);
685  // step 3: shift the result
686  ring S = 0,s,dp;
687  list L = imap(D,L);
688  kill D;
689  if (size(L)==2)
690  {
691    ideal B = L[1];
692    for (i=1; i<=ncols(B); i++)
693    {
694      B[i] = -B[i]+c-r-1;
695    }
696    L[1] = B;
697  }
698  else // should never get here: BS poly has non-rational roots
699  {
700    string str = L[3];
701    L = delete(L,3);
702    str = "poly @b = (" + str + ")*(" + string(fl2poly(L,"s")) + ");";
703    execute(str);
704    poly b = subst(@b,s,-s+c-r-1);
705    L = bFactor(b);
706  }
707  setring save;
708  list L = imap(S,L);
709  return(L);
710}
711example
712{
713  "EXAMPLE:"; echo = 2;
714  ring R = 0,(x,y,z),dp;
715  ideal F = x^2+y^3, z;
716  list L = bfctVarIn(F);
717  L;
718}
719
720static proc smallGenCoDim (ideal I, int Iasgiven)
721{
722  // call from K[x]
723  // returns list L
724  // L[1]=I or L[1]=smaller generating set of I
725  // L[2]=codimension(I)
726  int ppl = printlevel - voice + 3;
727  int n = nvars(basering);
728  int c;
729  if (attrib(I,"isSB") == 1)
730  {
731    c = n - dim(I);
732    if (!Iasgiven)
733    {
734      list L = mstd(I);
735    }
736  }
737  else
738  {
739    def save = basering;
740    list RL = ringlist(save);
741    list @ord;
742    @ord[1] = list("dp", intvec(1:n));
743    @ord[2] = list("C", intvec(0));
744    RL[3] = @ord;
745    kill @ord;
746    if (size(RL)>4) // commutative G-algebra present
747    {
748      RL = RL[1..4];
749    }
750    def R = ring(RL);
751    kill RL;
752    setring R;
753    ideal I = imap(save,I);
754    if (!Iasgiven)
755    {
756      list L = mstd(I);
757      c = n - dim(L[1]);
758      setring save;
759      list L = imap(R,L);
760    }
761    else
762    {
763      I = std(I);
764      c = n - dim(I);
765      setring save;
766    }
767    kill R;
768  }
769  if (!Iasgiven)
770  {
771    if (size(L[2]) < size(I))
772    {
773      I = L[2];
774      dbprint(ppl,"// Found smaller generating set of the given variety: ", I);
775    }
776    else
777    {
778      dbprint(ppl,"// Have not found smaller generating set of the given variety.");
779    }
780  }
781  dbprint(ppl-1,"// The codim of the given variety is " + string(c) + ".");
782  if (!defined(L))
783  {
784    list L;
785  }
786  L[1] = I;
787  L[2] = c;
788  return(L);
789}
790
791
792// Some more examples
793
794static proc TXcups()
795{
796  "EXAMPLE:"; echo = 2;
797  //TX tangent space of X=V(x^2+y^3)
798  ring R = 0,(x0,x1,y0,y1),Dp;
799  ideal F = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
800  printlevel = 0;
801  //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
802  //eng = 0;
803  def A = SannfsVar(F);
804  setring A;
805  LD;
806}
807
808static proc ex47()
809{
810  ring r7 = 0,(x0,x1,y0,y1),dp;
811  ideal I = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
812  bfctVarIn(I);
813  // second ex - too big
814  ideal J = x0^4+y0^5, 4*x0^3*x1+5*y0^4*y1;
815  bfctVarIn(J);
816}
817
818static proc ex48()
819{
820  ring r8 = 0,(x1,x2,x3),dp;
821  ideal I = x1^3-x2*x3, x2^2-x1*x3, x3^2-x1^2*x2;
822  bfctVarIn(I);
823}
824
825static proc ex49 ()
826{
827  ring r9 = 0,(z1,z2,z3,z4),dp;
828  ideal I = z3^2-z2*z4, z2^2*z3-z1*z4, z2^3-z1*z3;
829  bfctVarIn(I);
830}
831
832static proc ex410()
833{
834  LIB "toric.lib";
835  ring r = 0,(z(1..7)),dp;
836  intmat A[3][7];
837  A = 6,4,2,0,3,1,0,0,1,2,3,0,1,0,0,0,0,0,1,1,2;
838  ideal I = toric_ideal(A,"pt");
839  I = std(I);
840//ideal I = z(6)^2-z(3)*z(7), z(5)*z(6)-z(2)*z(7), z(5)^2-z(1)*z(7),
841//  z(4)*z(5)-z(3)*z(6), z(3)*z(5)-z(2)*z(6), z(2)*z(5)-z(1)*z(6),
842//  z(3)^2-z(2)*z(4), z(2)*z(3)-z(1)*z(4), z(2)^2-z(1)*z(3);
843  bfctVarIn(I,1); // no result yet
844}
Note: See TracBrowser for help on using the repository browser.