source: git/Singular/LIB/dmodvar.lib @ 6fe9a5

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