# source:git/Singular/LIB/dmodvar.lib@f4490f5

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