# source:git/Singular/LIB/dmodvar.lib@0610f0e

jengelh-datetimespielwiese
Last change on this file since 0610f0e was 0610f0e, checked in by Hans Schoenemann <hannes@…>, 13 years ago
• Property mode set to `100755`
File size: 22.6 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
12      R = K[x_1,...,x_n] and a set of polynomial f_1,..., f_r in R, define
13      F = f_1 * ... * f_r and F^s:=f_1^s_1*...*f_r^s_r for symbolic discrete
14      (that is shiftable) variables s_1,..., s_r.
15      The module R[1/F]*F^s has a structure of a D<S>-module, where
16      D<S> := D(R) tensored with S over K, where
17      - 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>
18      - S is the universal enveloping algebra of gl_r, generated by s_{ij}, where s_{ii}=s_i.
19 One is interested in the following data:
20 - the left ideal Ann F^s in D<S>, usually denoted by LD in the output
21 - global Bernstein polynomial in one variable s = s_1 + ...+ s_r, denoted by bs,
22 - its minimal integer root s0, the list of all roots of bs, which are known
23   to be rational, with their multiplicities, which is denoted by BS
24 - an r-tuple of operators in D<S>, denoted by PS, such that the functional equality
25     sum(k=1 to k=r) P_k*f_k*F^s = bs*F^s holds in R[1/F]*F^s.
26
27REFERENCES:
28  (BMS06) Budur, Mustata, Saito: Bernstein-Sato polynomials of arbitrary varieties (2006).
29  (ALM09) Andres, Levandovskyy, Martin-Morales : Principal Intersection and Bernstein-Sato Polynomial of an Affine Variety (2009).
30
31MAIN PROCEDURES:
32bfctVarIn(F[,L]);     compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using initial ideal approach
33bfctVarAnn(F[,L]);  compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using Sannfs approach
34SannfsVar(F[,O,e]); compute the annihilator of F^s in the ring D<S>
35
36AUXILIARY PROCEDURES:
37makeIF(F[,ORD]);    create the Malgrange ideal, associated with F = F[1],..,F[P]
38";
39
40// Static procs:
41//coDim(I);           compute the codimension of the leading ideal of I
42// dmodvarAssumeViolation()
43// ORDstr2list (ORD, NN)
44// smallGenCoDim(I,k)
45
46
47LIB "bfun.lib";    // for pIntersect
48LIB "dmodapp.lib"; // for isCommutative etc.
49
50
51///////////////////////////////////////////////////////////////////////////////
52
53// testing for consistency of the library:
54proc testdmodvarlib ()
55{
56  "AUXILIARY PROCEDURES:";
57  example makeIF;
58  "MAIN PROCEDURES:";
59  example bfctVarIn;
60  example bfctVarAnn;
61  example SannfsVar;
62}
63
64//   example coDim;
65
66///////////////////////////////////////////////////////////////////////////////
67
68static proc dmodvarAssumeViolation()
69{
70  // returns Boolean : yes/no [for assume violation]
71  // char K = 0
72  // no qring
73  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
74  {
75    //    "ERROR: no qring is allowed";
76    return(1);
77  }
78  return(0);
79}
80
81// da: in smallGenCoDim(), rewritten using mstd business
82static proc coDim (ideal I)
83"USAGE:  coDim (I);  I an ideal
84RETURN:  int
85PURPOSE: computes the codimension of the ideal generated by the leading monomials
86@*       of the given generators of the ideal. This is also the codimension of
87@*       the ideal if it is represented by a standard basis.
88NOTE:    The codimension of an ideal I means the number of variables minus the
89@*       Krull dimension of the basering modulo I.
90EXAMPLE: example SannfsVar; shows examples
91"
92{
93  int n = nvars(basering);
94  int d = dim(I); // to insert: check whether I is in GB
95  return(n-d);
96}
97example
98{
99  "EXAMPLE:"; echo = 2;
100  ring R = 0,(x,y,z),Dp;
101  ideal I = x^2+y^3, z;
102  coDim(std(I));
103}
104
105static proc ORDstr2list (string ORD, int NN)
106{
107  /* convert an ordering defined in NN variables in the  */
108  /* string form into the same ordering in the list form */
109  string st;
110  st = "ring @Z = 0,z(1.." + string(NN) + "),";
111  st = st + ORD + ";";
112  execute(st); kill st;
113  list L = ringlist(@Z)[3];
114  kill @Z;
115  return(L);
116}
117
118proc SannfsVar (ideal F, list #)
119"USAGE:  SannfsVar(F [,ORD,eng]);  F an ideal, ORD an optional string, eng an optional int
120RETURN:  ring
121PURPOSE: compute the D<S>-module structure of D<S>*f^s where f = F[1]*..*F[P]
122and D<S> is the Weyl algebra D tensored with K<S>=U(gl_P), according to the
123generalized algorithm by Briancon and Maisonobe for affine varieties.
124NOTE:    activate this ring with the @code{setring} command.
125@*       In the ring D<S>, the ideal LD is the needed D<S>-module structure.
126@*       The value of ORD must be an elimination ordering in D<Dt,S> for Dt
127@*       written in the string form, otherwise the result may have no meaning.
128@*       By default ORD = '(a(1..(P)..1),a(1..(P+P^2)..1),dp)'.
129@*       If eng<>0, @code{std} is used for Groebner basis computations,
130@*       otherwise, and by default @code{slimgb} is used.
131@*       If printlevel=1, progress debug messages will be printed,
132@*       if printlevel>=2, all the debug messages will be printed.
133EXAMPLE: example SannfsVar; shows examples
134"
135{
136  if (dmodvarAssumeViolation())
137  {
138    ERROR("Basering is inappropriate: characteristic>0 or qring present");
139  }
140  if (!isCommutative())
141  {
142    ERROR("Basering must be commutative");
143  }
144  def save = basering;
145  int N = nvars(basering);
146  int P = ncols(F);  //ncols better than size, since F[i] could be zero
147  // P is needed for default ORD
148  int i,j,k,l;
149  // st = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";
150  string st = "(a(" + string(1:P);
151  st = st + "),a(" + string(1:(P+P^2));
152  st = st + "),dp)";
153  // default values
154  string ORD = st;
155  int eng = 0;
156  if ( size(#)>0 )
157  {
158    if ( typeof(#[1]) == "string" )
159    {
160      ORD = string(#[1]);
161      // second arg
162      if (size(#)>1)
163      {
164        // exists 2nd arg
165        if ( typeof(#[2]) == "int" )
166        {
167          // the case: given ORD, given engine
168          eng = int(#[2]);
169        }
170        else
171        {
172          eng = 0;
173        }
174      }
175      else
176      {
177        // no second arg
178        eng = 0;
179      }
180    }
181    else
182    {
183      if ( typeof(#[1]) == "int" )
184      {
185        // the case: default ORD, engine given
186        eng = int(#[1]);
187        // ORD = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";  //is already set
188      }
189      else
190      {
191        // incorr. 1st arg
192        ORD = string(st);
193      }
194    }
195  }
196  // size(#)=0, i.e. there is no elimination ordering and no engine given
197  // eng = 0; ORD = "(a(1..(P)..1),a(1..(P^2)..1),dp)";  //are already set
198  int ppl = printlevel-voice+2;
199  // returns a list with a ring and an ideal LD in it
200  // save, N, P and the indices are already defined
201  int Nnew = 2*N+P+P^2;
202  list RL = ringlist(basering);
203  list L;
204  L[1] = RL[1];  //char
205  L[4] = RL[4];  //char, minpoly
206  // check whether vars have admissible names
207  list Name  = RL[2];
208  list RName;
209  // (i,j) <--> (i-1)*p+j
210  for (i=1; i<=P; i++)
211  {
212    RName[i] = "Dt("+string(i)+")";
213    for (j=1; j<=P; j++)
214    {
215      st = "s("+string(i)+")("+string(j)+")";
216      RName[P+(i-1)*P+j] = st;
217    }
218  }
219  for(i=1; i<=N; i++)
220  {
221    for(j=1; j<=size(RName); j++)
222    {
223      if (Name[i] == RName[j])
224      {
225        ERROR("Variable names should not include Dt(i), s(i)(j)");
226      }
227    }
228  }
229  // now, create the names for new vars
230  list DName;
231  for(i=1; i<=N; i++)
232  {
233    DName[i] = "D"+Name[i];  //concat
234  }
235  list NName = RName + Name + DName;
236  L[2] = NName;
237  // Name, Dname will be used further
238  kill NName;
239  //block ord (a(1..(P)..1),a(1..(P+P^2)..1),dp);
240  //export Nnew;
241  L[3] = ORDstr2list(ORD,Nnew);
242  // we are done with the list
243  def @R@ = ring(L);
244  setring @R@;
245  matrix @D[Nnew][Nnew];
246  // kronecker(i,j) equals (i==j)
247  // (i,j) <--> (i-1)*p+j
248  for (i=1; i<=P; i++)
249  {
250    for (j=1; j<=P; j++)
251    {
252      for (k=1; k<=P; k++)
253      {
254        //[sij,Dtk] = djk*Dti
255        @D[k,P+(i-1)*P+j] = (j==k)*Dt(i);
256        for (l=1; l<=P; l++)
257        {
258          if ( (i-k)*P < l-j )
259          {
260            //[sij,skl] = djk*sil - dil*skj
261            @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*s(i)(l) + (i==l)*s(k)(j);
262          }
263        }
264      }
265    }
266  }
267  for (i=1; i<=N; i++)
268  {
269    //[Dx,x]=1
270    @D[P+P^2+i,P+P^2+N+i] = 1;
271  }
272  def @R = nc_algebra(1,@D);
273  setring @R;
274  //@R@ will be used further
275  dbprint(ppl,"// -1-1- the ring @R(_Dt,_s,_x,_Dx) is ready");
276  dbprint(ppl-1, @R);
277  // create the ideal I
278  // (i,j) <--> (i-1)*p+j
279  ideal  F = imap(save,F);
280  ideal I;
281  for (i=1; i<=P; i++)
282  {
283    for (j=1; j<=P; j++)
284    {
285      I[(i-1)*P+j] = Dt(i)*F[j] + s(i)(j);
286    }
287  }
288  poly p,q;
289  for (i=1; i<=N; i++)
290  {
291    p=0;
292    for (j=1; j<=P; j++)
293    {
294      q = Dt(j);
295      q = q*diff(F[j],var(P+P^2+i));
296      p = p + q;
297    }
298    I = I, p + var(P+P^2+N+i);
299  }
300  // -------- the ideal I is ready ----------
301  dbprint(ppl,"// -1-2- starting the elimination of "+string(Dt(1..P))+" in @R");
302  dbprint(ppl-1, I);
303  ideal J = engine(I,eng);
304  ideal K = nselect(J,1..P);
305  kill I,J;
306  dbprint(ppl,"// -1-3- all Dt(i) are eliminated");
307  dbprint(ppl-1, K);  //K is without Dt(i)
308  // ----------- the ring @R2(_s,_x,_Dx) ------------
309  //come back to the ring save, recover L and remove all Dt(i)
310  //L[1],L[4] do not change
311  setring save;
312  list Lord, tmp;
313  // variables
314  tmp = L[2];
315  Lord = tmp[P+1..Nnew];
316  L[2] = Lord;
317  // ordering
318  // st = "(a(1..(P^2)..1),dp)";
319  st = "(a(" + string(1:P^2);
320  st = st + "),dp)";
321  tmp = ORDstr2list(st,Nnew-P);
322  L[3] = tmp;
323  def @R2@ = ring(L);
324  kill L;
325  // we are done with the list
326  setring @R2@;
327  matrix tmpM,LordM;
328  // non-commutative relations
329  intvec iv = P+1..Nnew;
330  tmpM = imap(@R@,@D);
331  kill @R@;
332  LordM = submat(tmpM,iv,iv);
333  matrix @D2 = LordM;
334  def @R2 = nc_algebra(1,@D2);
335  setring @R2;
336  kill @R2@;
337  dbprint(ppl,"// -2-1- the ring @R(_s,_x,_Dx) is ready");
338  dbprint(ppl-1, @R2);
339  ideal K = imap(@R,K);
340  kill @R;
341  dbprint(ppl,"// -2-2- starting cosmetic Groebner basis computation");
342  dbprint(ppl-1, K);
343  K = engine(K,eng);
344  dbprint(ppl,"// -2-3- the cosmetic Groebner basis has been computed");
345  dbprint(ppl-1,K);
346  ideal LD = K;
347  attrib(LD,"isSB",1);
348  export LD;
349  return(@R2);
350}
351example
352{
353  "EXAMPLE:"; echo = 2;
354  ring R = 0,(x,y),Dp;
355  ideal F = x^3, y^5;
356  //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
357  //eng = 0;
358  def A = SannfsVar(F);
359  setring A;
360  LD;
361}
362
363proc bfctVarAnn (ideal F, list #)
364"USAGE:  bfctVarAnn(F[,gid,eng]);  F an ideal, gid,eng optional ints
365RETURN:  list of ideal and intvec
366PURPOSE: computes the roots of the Bernstein-Sato polynomial and their
367@*       multiplicities for an affine algebraic variety defined by
368@*        F = F[1],..,F[r].
369ASSUME:  The basering is a commutative polynomial ring in char 0.
370BACKGROUND: In this proc, the annihilator of f^s in D[s] is computed and then a
371@*       system of linear equations is solved by linear reductions in order to
372@*       find the minimal polynomial of S = s(1)(1) + ... + s(P)(P), following
373@*       the ideas by Noro.
374NOTE:    In the output list, the ideal contains all the roots and the invec their
375@*        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.