Ticket #55: bfun.lib.diff2

File bfun.lib.diff2, 38.0 KB (added by Oleksandr , 15 years ago)

Inspecting my changes/sugestions vs latest author's version

Line 
12c2
2< version="$Id: bfct.lib,v 1.7 2008/12/09 16:50:21 levandov Exp $";
3---
4> version="$Id: bfun.lib,v 1.1 2009/02/12 20:25:22 levandov Exp $";
55c5
6< LIBRARY: bfct.lib     M. Noro's Algorithm for Bernstein-Sato polynomial
7---
8> LIBRARY: bfun.lib     Algorithms for b-functions and Bernstein-Sato polynomial
910,12c10,12
10< @*      one is interested in the global Bernstein-Sato polynomial b(s) in K[s],
11< @*      defined to be the monic polynomial, satisfying a functional identity
12< @*      L * F^{s+1} = b(s) F^s,  for some operator L in D[s], where D is the
13---
14> @*      one is interested in the global b-Function (also known as Bernstein-Sato
15> @*      polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying a functional
16> @*      identity L * F^{s+1} = b(s) F^s,   for some operator L in D[s]. Here, D stands for an
1715,30c15,20
18< @*      = global Bernstein-Sato polynomial b(s) in K[s]
19< @*      = the list of its roots, which are known to be rational
20< @*      = roots multiplicities.
21<
22< // TODO: IN THE ABOVE IS NOT CLEAR: WHO IS L? ARE b(s) AND L UNIQUE?
23<
24< // TODO: MAYBE ADD MORE THEORY?
25<
26< // TODO: GLOBAL ASSUMPTIONS? E.G: WHAT ABOUT PRIME CHAR. OF K? BASERING IS COMM?
27<
28< // TODO: EXPLAIN WHY IS THIS IMPORTANT?
29<
30< // TODO: ADD APPLICATION REFERENCES?
31<
32< // TODO: ADD THEORY / ALGORITHM REFERENCES
33<
34---
35> @*      - Bernstein-Sato polynomial b(s) in K[s],
36> @*      - the list of its roots, which are known to be rational
37> @*      - the multiplicities of the roots.
38> @*      References: Saito, Sturmfels, Takayama: Groebner Deformations of Hypergeometric
39> @*      Differential Equations (2000), Noro: An Efficient Modular Algorithm for Computing
40> @*      the Global b-function, (2002).
4135,51c25,34
42< // TODO: WHAT IS THE DIFFERENCE BETWEEN THE FOLLOWING PROCS?
43<
44< bfct(f[,s,t,v]);            compute the global Bernstein-Sato poly. of poly. f
45<
46< // TODO: BAD NAMES: RENAME! (SEE 3.9.1 Procedures in a library, 5^th rule)
47<
48< bfctsyz(f[,r,s,t,u,v]);     compute the global Bernstein-Sato poly. of poly. f
49< bfctann(f[,s]);             compute the global Bernstein-Sato poly. of poly. f
50< bfctonestep(f[,s,t]);       compute the global Bernstein-Sato poly. of poly. f
51<
52< // TODO: WHAT IS THE GLOBAL B-FUNCTION OF AN IDEAL (WRT. WEIGHT)
53<
54< bfctideal(I,w[,s,t]);       compute the global b-function of ideal I wrt weight
55< pintersect(f,I);            compute the intersection of ideals <f> and I
56< pintersectsyz(f,I[,p,s,t]); compute the intersection of ideals <f> and I
57< linreduce(f,I[,s]);         reduce a poly by linear reductions only
58< ncsolve(I[,s]);             compute a linear dependency of the elements of I
59---
60> bfct(f[,s,t,v]);            computes the Bernstein-Sato polynomial of poly f
61> bfctSyz(f[,r,s,t,u,v]);     computes the Bernstein-Sato polynomial of poly f
62> bfctAnn(f[,s]);             computes the Bernstein-Sato polynomial of poly f
63> bfctOneGB(f[,s,t]);         computes the Bernstein-Sato polynomial of poly f
64> bfctIdeal(I,w[,s,t]);       computes the global b-function of ideal I w.r.t. the weight w
65> pIntersect(f,I);            intersection of the ideal I with the subalgebra K[f] for a poly f
66> pIntersectSyz(f,I[,p,s,t]); intersection of the ideal I with the subalgebra K[f] for a poly f
67> linReduce(f,I[,s]);         reduces a poly by linear reductions only
68> linReduceIdeal(I,[s,t]);    reduces generators of ideal by linear reductions only
69> linSyzSolve(I[,s]);         computes a linear dependency of the elements of ideal I
7055,57c38
71< // TODO: BAD NAME: MAYBE 'allpositive'.
72< ispositive(v);   check whether all entries of an intvec are positive
73< isin(l,i);       check whether an element is a member of a list
74---
75> allPositive(v);  checks whether all entries of an intvec are positive
7667c48
77< LIB "dmodapp.lib";  // for initialideal etc
78---
79> LIB "dmodapp.lib";  // for initialIdealW etc
8071,72c52,53
81< // The following procedure is for testing the procedures of the library:
82< proc testbfctlib ()
83---
84> // testing for consistency of the library:
85> proc testbfunlib ()
8676,77c57
87<   example ispositive;
88<   example isin;
89---
90>   example allPositive;
9184c64
92<   example bfctonestep;
93---
94>   example bfctOneGB;
9589c69,70
96<   example ncsolve;
97---
98>   example linReduceIdeal;
99>   example linSyzSolve;
10092,93c73
101<
102< //--------------- auxiliary procedures ---------------------------------------------------------
103---
104> //--------------- auxiliary procedures ----------------------------------------
10598a79
106> ASSUME: basering is a Weyl algebra
107150d130
108<   LIB "bfct.lib";
109160,166c140,142
110< ///////////////////////////////////////////////////////////////////////////////
111< // TODO: BAD NAME: INTVEC CANNOT BE POSITIVE... MAYBE 'allPositive'.
112< // (see 3.9.1 Procedures in a library, 5^th rule)
113<
114< proc ispositive (intvec v)
115< "USAGE:  ispositive(v);  v an intvec
116< RETURN:  int, 1 if all entries of v are positive, or 0 otherwise
117---
118> proc allPositive (intvec v)
119> "USAGE:  allPositive(v);  v an intvec
120> RETURN:  int, 1 if all components of v are positive, or 0 otherwise
121168c144
122< EXAMPLE: example ispositive; shows an example
123---
124> EXAMPLE: example allPositive; shows an example
125186c162
126<   ispositive(v);
127---
128>   allPositive(v);
129188c164
130<   ispositive(w);
131---
132>   allPositive(w);
133191,199c167,169
134<
135< ///////////////////////////////////////////////////////////////////////////////
136< // TODO: IMPLEMENTATION IS OVERSIMPLIFIED! WITHOUT FURTHER ASSUMPTIONS ON THE
137< // LIST/ARGUMENT STRUCTURE IT MAY BE CONSIDERED BUGGY!
138< // TODO: BAD NAME: MAYBE 'findFirst' (see 3.9.1 Procedures in a library, 5^th rule)
139<
140< proc isin (list l, def i)
141< "USAGE:  isin(l,i);  l a list, i an argument of any type
142< RETURN:  int, the position of the first appearance of i in l, or 0 if not there
143---
144> static proc findFirst (list l, i)
145> "USAGE:  findFirst(l,i);  l a list, i an argument of any type
146> RETURN:  int, the position of the first appearance of i in l, or 0 if i is not a member of l
147201c171
148< EXAMPLE: example isin; shows an example
149---
150> EXAMPLE: example findFirst; shows an example
151214a185,188
152> //   findFirst(list(1, 2, list(1)),2);       // seems to be a bit buggy,
153> //   findFirst(list(1, 2, list(1)),3);       // but works for the purposes
154> //   findFirst(list(1, 2, list(1)),list(1)); // of this library
155> //   findFirst(l,list(2));
156220,226c194,195
157<   isin(l,4);
158<   isin(l,2);
159<
160<   isin(list(1, 2, list(1)),2); // OK, BUT:
161<   isin(list(1, 2, list(1)),3); // TODO: FIX THIS BUG!
162<   isin(list(1, 2, list(1)),list(1)); // TODO: FIX THIS BUG!
163<   isin(l,list(2)); // TODO: FIX THIS BUG!
164---
165>   findFirst(l,4);
166>   findFirst(l,2);
167229,231d197
168<
169< // TODO: BAD NAME: (see 3.9.1 Procedures in a library, 5^th rule)
170<
171235,236c201,202
172< PURPOSE: compute the scalar product of two intvecs
173< NOTE:    the arguments must have the same size
174---
175> PURPOSE: computes the scalar product of two intvecs
176> ASSUME:  the arguments are of the same size
177240c206
178<   int i; int sp = 0;
179---
180>   int i; int sp;
181264,275c229,236
182<
183< // TODO: BAD NAME: (see 3.9.1 Procedures in a library, 5^th rule)
184<
185< proc linreduce(poly f, ideal I, list #)
186< "USAGE:  linreduce(f, I [,s,t]);  f a poly, I an ideal, s,t optional ints
187< RETURN:  poly/list, linear reductum (over field) of f by elements from I
188< PURPOSE: reduce a poly only by linear reductions (no monomial multiplications)
189< NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient
190< @*       vector of the used reductions is returned, otherwise (and by default)
191< @*       only reduced poly is returned.
192< // TODO: WHAT IS THE OUTPUT FORMAT IF S <> 0?
193< @*       If t==0 (and by default) all monomials are reduced (if possible),
194---
195> proc linReduceIdeal(ideal I, list #)
196> "USAGE:  linReduceIdeal(I [,s,t,u]); I an ideal, s,t,u optional ints
197> RETURN:  ideal or list, linear reductum (over field) of I by its elements
198> PURPOSE: reduce a list of polys only by linear reductions (no monomial multiplications)
199> NOTE:    If s<>0, a list consisting of the reduced ideal and the coefficient
200> @*       vectors of the used reductions given as module is returned.
201> @*       Otherwise (and by default), only the reduced ideal is returned.
202> @*       If t<>0 (and by default) all monomials are reduced (if possible),
203277,278c238,243
204< // TODO: WRITE DOWN WHAT DOES THE ABOVE MEAN? 
205< EXAMPLE: example linreduce; shows examples
206---
207> @*       If u<>0 (and by default), the ideal is first sorted in increasing order.
208> @*       If u is set to 0 and the given ideal is not sorted in the way described,
209> @*       the result might not be as expected.
210> DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
211> @*       if printlevel>=2, all the debug messages will be printed.
212> EXAMPLE: example linReduceIdeal; shows examples
213280a246,248
214>   // #[1] = remembercoeffs
215>   // #[2] = redtail
216>   // #[3] = sortideal
217283c251,252
218<   int redlm          = 0; // default
219---
220>   int redtail        = 1; // default
221>   int sortideal      = 1; // default
222294c263,270
223<         redlm = #[2];
224---
225>         redtail = #[2];
226>       }
227>       if (size(#)>2)
228>       {
229>         if (typeof(#[3])=="int" || typeof(#[3])=="number")
230>         {
231>           sortideal = #[3];
232>         }
233298,309d273
234<   // TODO: SEEMS A BIT OVERCOMPLICATED:
235<   // = ARRAY OF NORMALIZED TERMS OF F (monomf) HAS TO BE UPDATED AFTER EVERY
236<   //   REDUCTION => NO SPEED UP (JUST MORE WORK)
237<   //   MOREOVER YOO TRY TO REDUCE JUST ANY TERM IN F INSTEAD OF
238<   //   REDUCING THE LEADING TERM AS FAR AS POSSIBLE AND THEN REDUCE THE
239<   //   REMAINING TAIL (THIS IS THE USUAL APPROACH E.G. IN REDUCED NORMAL FORM)!
240<   //   THIS WAY YOU CAN ALSO OMIT TWO INDEPENDENT CASES...
241<   // SPEED-UPS:
242<   // = YOU CAN USE GRADING BY DEGREE TO SEED UP SEARCHES
243<   // = YOU COULD PREREDUCE I IN ORDER TO OMIT CIRCULAR REDUCTIOS
244<   //   E.G. linreduce(X, IDEAL(X-Y, Y-Z, Z-...)
245<   int i,j,k;
246311c275,287
247<   ideal lmI,lcI;
248---
249>   int sZeros = sI - size(I);
250>   int i,j;
251>   ideal J,lmJ,ordJ;
252>   list lJ = sort(I);
253>   module M; // for the coefficients
254>   // step 1: prepare, e.g. sort I
255>   if (sortideal <> 0)
256>   {
257>     if (sZeros > 0) // I contains zeros
258>     {
259>       if (remembercoeffs <> 0)
260>       {
261>         j = 0;
262314,315c290,295
263<     lmI[i] = leadmonom(I[i]);
264<     lcI[i] = leadcoef(I[i]);
265---
266>           if (I[i] == 0)
267>           {
268>             j++;
269>             J[j] = 0;
270>             ordJ[j] = -1;
271>             M[j] = gen(i);
272317,320c297
273<   vector v;
274<   poly c;
275<   int reduction = 1;
276<   if (redlm == 0)
277---
278>           else
279322,323c299,303
280<     ideal monomf;
281<     for (k=1; k<=size(f); k++)
282---
283>             M[i+sZeros-j] = gen(lJ[2][i-j]+j);
284>           }
285>         }
286>       }
287>       else
288325c305,326
289<       monomf[k] = normalize(f[k]);
290---
291>         for (i=1; i<=sZeros; i++)
292>         {
293>           J[i] = 0;
294>           ordJ[i] = -1;
295>         }
296>       }
297>       I = J,lJ[1];
298>     }
299>     else // I contains no zeros
300>     {
301>       I = lJ[1];
302>       if (remembercoeffs <> 0)
303>       {
304>         for (i=1; i<=size(lJ[1]); i++) { M[i+sZeros] = gen(lJ[2][i]); }
305>       }
306>     }
307>   }
308>   else // assume I is already sorted
309>   {
310>     if (remembercoeffs <> 0)
311>     {
312>       for (i=1; i<=ncols(I); i++) { M[i] = gen(i);  }
313326a328,347
314>   } 
315>   dbprint(ppl-1,"initially sorted ideal:", I);
316>   if (remembercoeffs <> 0) { dbprint(ppl-1," used permutations:", M); }
317>   // step 2: reduce leading monomials by linear reductions
318>   poly lm,c,redpoly,maxlmJ;
319>   J[sZeros+1] = I[sZeros+1];
320>   lm = I[sZeros+1];
321>   maxlmJ = leadmonom(lm);
322>   lmJ[sZeros+1] = maxlmJ;
323>   int ordlm,reduction;
324>   ordJ[sZeros+1] = ord(lm);
325>   vector v;
326>   int noRedPast;
327>   for (i=sZeros+2; i<=sI; i++)
328>   {
329>     redpoly = I[i];
330>     lm = leadmonom(redpoly);
331>     ordlm = ord(lm);
332>     if (remembercoeffs <> 0) { v = M[i]; }
333>     reduction = 1;
334328a350
335>       noRedPast = i;
336330c352
337<       for (i=sI; i>=1; i--)
338---
339>       for (j=sZeros+1; j<noRedPast; j++)
340332,333c354,356
341<         dbprint(ppl,"testing ideal entry:",i);
342<         for (j=1; j<=size(f); j++)
343---
344>         if (lm == 0) { break; } // nothing more to reduce
345>         if (lm > maxlmJ) { break; } //lm is bigger than maximal monomial to reduce with
346>         if (ordlm == ordJ[j])
347335c358
348<           if (monomf[j] == lmI[i])
349---
350>           if (lm == lmJ[j])
351337,341c360,374
352<             c = leadcoef(f[j])/lcI[i];
353<             f = f - c*I[i];
354<             dbprint(ppl,"reducing poly to ",f);
355<             monomf = 0;
356<             for (k=1; k<=size(f); k++)
357---
358>             dbprint(ppl-1,"reducing " + string(redpoly));
359>             dbprint(ppl-1," with " + string(J[j]));
360>             c = leadcoef(redpoly)/leadcoef(J[j]);
361>             redpoly = redpoly - c*J[j];
362>             dbprint(ppl-1," to " + string(redpoly));
363>             lm = leadmonom(redpoly);
364>             ordlm = ord(lm);
365>             if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; }
366>             noRedPast = j;
367>             reduction = 1;
368>           }
369>         }
370>       }
371>     }
372>     for (j=sZeros+1; j<i; j++)
373343c376
374<               monomf[k] = normalize(f[k]);
375---
376>       if (redpoly < J[j]) { break; }
377345c378,380
378<             reduction = 1;
379---
380>     J = insertGenerator(J,redpoly,j);
381>     lmJ = insertGenerator(lmJ,lm,j);
382>     ordJ = insertGenerator(ordJ,poly(ordlm),j);
383348c383,385
384<               v = v - c * gen(i);
385---
386>       v = M[i];
387>       M = deleteGenerator(M,i);
388>       M = insertGenerator(M,v,j);
389350c387
390<             break;
391---
392>     maxlmJ = lmJ[i];
393351a389,412
394>   // step 3: reduce tails by linear reductions as well
395>   if (redtail <> 0)
396>   {
397>     dbprint(ppl,"finished reducing leading monomials");
398>     dbprint(ppl-1,string(J));
399>     if (remembercoeffs <> 0) { dbprint(ppl-1,"used reductions:" + string(M)); }
400>     int k;
401>     for (i=sZeros+1; i<=sI; i++)
402>     {
403>       lm = lmJ[i];
404>       for (j=i+1; j<=sI; j++)
405>       {
406>         for (k=2; k<=size(J[j]); k++) // run over all terms in J[j]
407>         {
408>           if (ordJ[i] == ord(J[j][k]))
409>           {   
410>             if (lm == normalize(J[j][k]))
411>             {
412>               c = leadcoef(J[j][k])/leadcoef(J[i]);
413>               dbprint(ppl-1,"reducing " + string(J[j]));
414>               dbprint(ppl-1," with " + string(J[i]));
415>               J[j] = J[j] - c*J[i];
416>               dbprint(ppl-1," to " + string(J[j]));
417>               if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; }
418353c414,422
419<         if (reduction == 1)
420---
421>           }
422>         }
423>       }
424>     }
425>   }
426>   if (remembercoeffs <> 0) { return(list(J,M)); }
427>   else { return(J); }
428> }
429> example
430355c424,467
431<           break;
432---
433>   "EXAMPLE:"; echo = 2;
434>   ring r = 0,(x,y),dp;
435>   ideal I = 3,x+9,y4+5x,2y4+7x+2;
436>   linReduceIdeal(I);     // reduces tails
437>   linReduceIdeal(I,0,0); // no reductions of tails
438>   list l = linReduceIdeal(I,1); // reduces tails and shows reductions used
439>   l;
440>   module M = I;
441>   l[1] - ideal(M*l[2]);
442> }
443>
444> proc linReduce(poly f, ideal I, list #)
445> "USAGE:  linReduce(f, I [,s,t,u]); f a poly, I an ideal, s,t,u optional ints
446> RETURN:  poly or list, linear reductum (over field) of f by elements from I
447> PURPOSE: reduce a poly only by linear reductions (no monomial multiplications)
448> NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient
449> @*       vector of the used reductions is returned, otherwise (and by default)
450> @*       only reduced poly is returned.
451> @*       If t<>0 (and by default) all monomials are reduced (if possible),
452> @*       otherwise, only leading monomials are reduced.
453> @*       If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. instead
454> @*       of the given ideal, the output of @code{linReduceIdeal} is used. If u is
455> @*       set to 0 and the given ideal does not equal the output of
456> @*       @code{linReduceIdeal}, the result might not be as expected.
457> DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
458> @*       if printlevel>=2, all the debug messages will be printed.
459> EXAMPLE: example linReduce; shows examples
460> "
461> {
462>   int ppl = printlevel - voice + 2;
463>   int remembercoeffs = 0; // default
464>   int redtail        = 1; // default
465>   int prepareideal   = 1; // default
466>   if (size(#)>0)
467>   {
468>     if (typeof(#[1])=="int" || typeof(#[1])=="number")
469>     {
470>       remembercoeffs = #[1];
471>     }
472>     if (size(#)>1)
473>     {
474>       if (typeof(#[2])=="int" || typeof(#[2])=="number")
475>       {
476>         redtail = #[2];
477356a469,473
478>       if (size(#)>2)
479>       {
480>         if (typeof(#[3])=="int" || typeof(#[3])=="number")
481>         {
482>           prepareideal = #[3];
483360c477,484
484<   else // reduce only leading monomials
485---
486>   }
487>   int i,j,k;
488>   int sI = ncols(I);
489>   // pre-reduce I:
490>   module M;
491>   if (prepareideal <> 0)
492>   {
493>     if (remembercoeffs <> 0)
494361a486,510
495>       list sortedI = linReduceIdeal(I,1,redtail);
496>       I = sortedI[1];
497>       M = sortedI[2];
498>       dbprint(ppl-1,"prepared ideal:",I);
499>       dbprint(ppl-1," with operations:",M);
500>     }
501>     else { I = linReduceIdeal(I,0,redtail); }
502>   }
503>   else
504>   {
505>     if (remembercoeffs <> 0)
506>     {
507>       for (i=1; i<=sI; i++) { M[i] = gen(i); }
508>     }
509>   }
510>   ideal lmI,lcI,ordI;
511>   for (i=1; i<=sI; i++)
512>   {
513>     lmI[i] = leadmonom(I[i]);
514>     lcI[i] = leadcoef(I[i]);
515>     ordI[i] = ord(lmI[i]);
516>   }
517>   vector v;
518>   poly c;
519>   // === reduce leading monomials ===
520363c512,513
521<     while (reduction == 1) // while there was a reduction
522---
523>   int ordf = ord(lm);
524>   for (i=sI; i>=1; i--)  // I is sorted: smallest lm's on top
525365,366c515,516
526<       reduction = 0;
527<       for (i=sI;i>=1;i--)
528---
529>     if (lm == 0) { break; }
530>     else
531368c518,520
532<         if (lm <> 0 && lm == lmI[i])
533---
534>       if (ordf == ordI[i])
535>       {
536>         if (lm == lmI[i])
537373,374c525,532
538<           reduction = 1;
539<           if (remembercoeffs <> 0)
540---
541>           ordf = ord(lm);
542>           if (remembercoeffs <> 0) { v = v - c * M[i]; }
543>         }
544>       }
545>     }
546>   }
547>   // === reduce tails as well ===
548>   if (redtail <> 0)
549376c534,550
550<             v = v - c * gen(i);
551---
552>     dbprint(ppl,"finished reducing leading monomials");
553>     dbprint(ppl-1,string(f));
554>     if (remembercoeffs <> 0) { dbprint(ppl-1,"used reductions:" + string(v)); }
555>     for (i=1; i<=sI; i++)
556>     {
557>       dbprint(ppl,"testing ideal entry:",i);
558>       for (j=1; j<=size(f); j++)
559>       {
560>         if (ord(f[j]) == ordI[i])
561>         {
562>           if (normalize(f[j]) == lmI[i])
563>           {
564>             c = leadcoef(f[j])/lcI[i];
565>             f = f - c*I[i];
566>             dbprint(ppl-1,"reducing with " + string(I[i]));
567>             dbprint(ppl-1," to " + string(f));
568>             if (remembercoeffs <> 0) { v = v - c * M[i]; }
569387,390c561
570<   else
571<   {
572<     return(f);
573<   }
574---
575>   else { return(f); }
576397d567
577<   module M = module(I);
578400,408c570,572
579<   linreduce(g,I);
580<   linreduce(g,I,0,1);
581<   list l = linreduce(g,I,1);
582<   l;
583<   l[1] - (module(M[1..nrows(l[2])]) * l[2])[1][1] - (g);
584<   linreduce(f,I);
585<   l = linreduce(f,I,1);
586<   l;
587<   l[1] - (module(M[1..nrows(l[2])]) * l[2])[1][1] - (f);
588---
589>   linReduce(g,I);     // reduces tails
590>   linReduce(g,I,0,0); // no reductions of tails
591>   linReduce(f,I,1);   // reduces tails and shows reductions used
592410,411c574,575
593<   I = ideal(x3 - y3, y3 - x2, x3 - y2, x2 - y, y2-x); M = I;
594<   l = linreduce(f, I, 1);
595---
596>   I = x3-y3, y3-x2,x3-y2,x2-y,y2-x;
597>   list l = linReduce(f,I,1);
598413,414c577,578
599<   l[1] - (module(M[1..nrows(l[2])]) * l[2])[1][1] - (f);
600<
601---
602>   module M = I;
603>   f - (l[1]-(M*l[2])[1,1]);
604417,426c581,583
605<
606<
607< // TODO: WHAT ARE THE ASSUMPTIONS ON THE BASERING: COMM? NON-COMM?
608<
609< // TODO: BAD NAME: (see 3.9.1 Procedures in a library, 5^th rule)
610<
611<
612< proc ncsolve (ideal I, list #)
613< "USAGE:  ncsolve(I[,s]);  I an ideal, s an optional int
614< RETURN:  coefficient vector of a linear combination of 0 in the elements of I
615---
616> proc linSyzSolve (ideal I, list #)
617> "USAGE:  linSyzSolve(I[,s]);  I an ideal, s an optional int
618> RETURN:  vector, coefficient vector of a linear combination of 0 in the elements of I
619428,431c585,587
620<          if such one exists
621< NOTE:    The optional integer s determines the engine, that computes the GB:
622< // TODO: NOTE INCONSISTENCE WITH YOUR OTHER LIBS... (engine, smith, jacobson)
623< @*       0 means 'slimgb', otherwise 'std'.
624---
625> @*       if such one exists
626> NOTE:    If s<>0, @code{std} is used for Groebner basis computations,
627> @*       otherwise, @code{slimgb} is used.
628435c591
629< EXAMPLE: example ncsolve; shows examples
630---
631> EXAMPLE: example linSyzSolve; shows examples
632471c627,628
633<     def @B = save + @aA; setring @B;
634---
635>     def @B = save + @aA;
636>     setring @B;
637490c645
638<     dbprint(ppl, "ncsolve: starting Groebner basis computation with engine:", whichengine);
639---
640>     dbprint(ppl, "linSyzSolve: starting Groebner basis computation with engine:", whichengine);
641495c650
642<       dbprint(ppl+1, "no solutions by ncsolve");
643---
644>       dbprint(ppl+1, "no solutions by linSyzSolve");
645525c680
646<   ncsolve(I);
647---
648>   linSyzSolve(I);
649527,533c682
650<   ncsolve(J);
651<   // the following is just a bug check:
652<   ring r = (0,@a(1)),(@a(2)),dp;
653<   ideal I = @a(2),2@a(2);
654<   ncsolve(I);
655<   ideal J = @a(2),@a(2)2;
656<   ncsolve(J);
657---
658>   linSyzSolve(J);
659536,542c685,689
660< // TODO: WHAT ARE THE ASSUMPTIONS ON THE BASERING: COMM? NON-COMM?
661<
662< proc pintersect (poly s, ideal I)
663< "USAGE:  pintersect(f, I);  f a poly, I an ideal
664< RETURN:  vector, coefficient vector of the monic generator of the intersection
665<          of ideals <f> and I
666< PURPOSE: compute the intersection of the ideal I with the ideal <f>
667---
668> proc pIntersect (poly s, ideal I, list #)
669> "USAGE:  pIntersect(f, I [,s]);  f a poly, I an ideal, s an optional int
670> RETURN:  vector, coefficient vector of the monic polynomial
671> PURPOSE: compute the intersection of ideal I with the subalgebra K[f]
672> ASSUME:  I is given as Groebner basis.
673544c691,692
674< @*       I should be given as standard basis.
675---
676> @*       If s>0 is given, it is searched for the generator of the intersection
677> @*       only up to degree s. Otherwise (and by default), no bound is assumed.
678550a699,702
679>   if (attrib(I,"isSB") <> 1)
680>   {
681>     print("WARNING: The input has no SB attribute!");
682>     print(" Treating it as if it were a Groebner basis and proceeding...");
683552c704,713
684<   int ppl = printlevel-voice+2;
685---
686>   }
687>   int bound = 0; // default
688>   if (size(#)>0)
689>   {
690>     if (typeof(#[1])=="int" || typeof(#[1])=="number")
691>     {
692>       bound = #[1];
693>     }
694>   }
695>   int ppl = printlevel-voice+1;
696556c717
697<     if (simplify(I,1) == ideal(1))
698---
699>     if (simplify(I,2)[1] == 1)
700571,574c732
701<     if (i == ncols(I)+1)
702<     {
703<       break;
704<     }
705---
706>     if (i == ncols(I)+1) { break; }
707580,583c738
708<         if (degI[i][j] <> 0)
709<         {
710<           break;
711<         }
712---
713>         if (degI[i][j] <> 0) { break; }
714593a749,753
715>   if (bound>0 && bound<degbound) // given bound is too small
716>   {
717>     print("// Try a bound of at least " + string(degbound));
718>     return(vector(0));
719>   } 
720597a758
721>     print("// Intersection is zero");
722609a771
723>     if (bound>0 && i>bound) { return(vector(0)); }
724647,648c809
725<   // we obtain the coefficients of the generator of the intersection by the used
726<   // reductions:
727---
728>   // we obtain the coefficients of the generator of the intersection by the used reductions:
729685a847
730>   pIntersect(t*Dt,inF,1);
731688,692d849
732<
733< // TODO: WHAT ARE THE ASSUMPTIONS ON THE BASERING: COMM? NON-COMM?
734< // OVER A FIELD OF CHAR 0? WHAT ABOUT PRIME CHAR?
735< // TODO: BAD NAME: (see 3.9.1 Procedures in a library, 5^th rule)
736< 
737694,705c851,859
738< "USAGE:  pintersectsyz(f, I [,p,s,t]);  f poly, I ideal, opt p,t ints
739< RETURN:  vector, coefficient vector of the monic generator of the intersection
740< @*       of ideals <f> and I
741< PURPOSE: compute the intersection of an ideal I with a principal ideal <f>
742< NOTE:    If the intersection is zero, this proc might not terminate.
743< @*       I should be given as standard basis.
744< @*       If p>0 is given (must be prime), this proc computes the generator of
745< @*       the intersection in char p first and then only searches for a generator
746< @*       of the obtained degree in the basering. Otherwise, it searches for all
747< @*       degrees by computing syzygies.
748< // TODO: OVERCOMPLICATED: DO YOU NEED BOTH s AND t? ONE DOES KNOW p BEFOREHAND!?
749< @*       If s<>0, @code{std} is used for GB computations in char 0,
750---
751> "USAGE:  pIntersectSyz(f, I [,p,s,t]);  f a poly, I an ideal, p, t optial ints, p a prime number
752> RETURN:  vector, coefficient vector of the monic polynomial
753> PURPOSE: compute the intersection of an ideal I with the subalgebra K[f]
754> ASSUME:  I is given as Groebner basis.
755> NOTE:    If the intersection is zero, this procedure might not terminate.
756> @*       If p>0 is given, this proc computes the generator of the intersection in char p first
757> @*       and then only searches for a generator of the obtained degree in the basering.
758> @*       Otherwise, it searched for all degrees by computing syzygies.
759> @*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
760707c861
761< @*       If t<>0 and by default, @code{std} is used for GB computations in char >0,
762---
763> @*       If t<>0 and by default, @code{std} is used for Groebner basis computations in char >0,
764714,715d867
765<   // TODO: ADD ASSUMPTION CHECKS: I - GB, P - ZERO OR PRIME...
766<
767717a870,873
768>   if (attrib(I,"isSB") <> 1)
769>   {
770>     print("WARNING: The input has no SB attribute!");
771>     print(" Treating it as if it were a Groebner basis and proceeding...");
772718a875
773>   }
774768c925
775<     poly s = phi(s);  // TODO: WHAT ABOUT "imap(save, s);" INSTEAD?
776---
777>     poly s = phi(s);
778793c950
779<     dbprint(ppl,"ncsolve starts with: "+string(matrix(NI)));
780---
781>     dbprint(ppl,"linSyzSolve starts with: "+string(matrix(NI)));
782798c955
783<       v = ncsolve(NI,modengine);
784---
785>       v = linSyzSolve(NI,modengine);
786803c960
787<         v = ncsolve(NI,whichengine);
788---
789>         v = linSyzSolve(NI,whichengine);
790817c974
791<       v = ncsolve(NI,whichengine);
792---
793>       v = linSyzSolve(NI,whichengine);
794820c977
795<     dbprint(ppl,"ncsolve ready  with: "+string(MM));
796---
797>     dbprint(ppl,"linSyzSolve ready  with: "+string(MM));
798822c979
799<     //  "ncsolve ready with"; print(v);
800---
801>     //  "linSyzSolve ready with"; print(v);
802834c991
803<         dbprint(ppl,"ncsolve: bad solution!");
804---
805>         dbprint(ppl,"linSyzSolve: bad solution!");
806838c995
807<         dbprint(ppl,"ncsolve: got solution!");
808---
809>         dbprint(ppl,"linSyzSolve: got solution!");
810894c1050
811<   if ( (ringvar > nvars(save)) || (ringvar < 1) )
812---
813>   if (ringvar > nvars(save))
814898,899c1054
815<
816<   poly p = 0;
817---
818>   poly p;
819914,917d1068
820<   iv = 0;
821<   v = 0;
822<   vec2poly(v,2);
823<   vec2poly(iv); 
824920,921d1070
825< // TODO: ADD INPUT/OUTPUT FORMAT EXPLANATION!
826< // TODO: BAD NAME: (see 3.9.1 Procedures in a library, 5^th rule)
827970,976d1118
828<
829< // TODO: ADD HELP AS IF THIS IS A FULL-SCALE GLOBAL PROC!
830< // TODO: INCLUDE ASSUMPTIONS AND INPUT/OUTPUT FORMATS!!!
831< // TODO: ADD ALGORITHMS REFERENCES
832< // TODO: ADD STEP-BY-STEP COMMENTS/EXPLANATIONS!
833< // TODO: BAD NAME (see 3.9.1 Procedures in a library, 5^th rule)
834<
835982a1125,1128
836>   if (size(variables(f)) == 0) // f is constant
837>   {
838>     return(list(ideal(0),intvec(0)));
839>   } 
840993c1139
841<     def D = SannfsBFCT(f,whichengine); // TODO: WHAT IS THIS? WHERE IS THIS FROM???setring
842---
843>     def D = SannfsBFCT(f,whichengine);
8441012c1158
845<         if (isin(usedprimes,q)==0) // if q was not already used
846---
847>         if (findFirst(usedprimes,q)==0) // if q was not already used
8481043,1053d1188
849<
850<
851< // TODO: ASSUMPTIONS?
852< // TODO: ALGORITHM REFERENCE?
853< // TODO: MAYBE MORE EXPLICITE THEORY?
854< // TODO: WHAT IS A 'Malgrange ideal'?
855<
856< // TODO: DO YOU REALLY COMPUTE THE global Bernstein-Sato polynomial
857< // according to the algorithm by Masayuki Noro?
858< // OR ONLY ROOTS?
859<
8601057,1063c1192,1199
861< PURPOSE: computes the roots (as an ideal, together with their multiplicities
862< @*       as an intvec) of the global Bernstein-Sato polynomial bs(s) for the
863< @*       hypersurface, defined by f, computed according to the algorithm by
864< @*       Masayuki Noro.
865< BACKGROUND: In this proc, the initial Malgrange ideal is computed and than
866< @*       a system of linear equations is solved by linear reductions.
867< NOTE:    If s<>0, @code{std} is used for GB computations,
868---
869> PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
870> @*       for the hypersurface defined by f.
871> ASSUME:  The basering is a commutative polynomial ring in char 0.
872> BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm
873> @*       by Masayuki Noro and then a system of linear equations is solved by linear reductions.
874> NOTE:    In the output list, the ideal contains all the roots
875> @*       and the intvec their multiplicities.
876> @*       If s<>0, @code{std} is used for GB computations,
8771067,1068c1203,1204
878< @*       If v is a positive weight vector, v is used for homogenization
879< @*       computations, otherwise and by default, no weights are used.
880---
881> @*       If v is a positive weight vector, v is used for homogenization computations,
882> @*       otherwise and by default, no weights are used.
8831098c1234
884<         if (typeof(#[3])=="intvec" && size(#[3])==n && ispositive(#[3])==1)
885---
886>         if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
8871118,1131d1253
888<
889< // TODO: WRONG NAME : RENAME TO 'bfctSyz' (see 3.9.1 Procedures in a library, 5^th rule)
890<
891< // TODO AS ABOVE:
892< // TODO: ASSUMPTIONS?
893< // TODO: ALGORITHM REFERENCE?
894< // TODO: MAYBE MORE EXPLICITE THEORY?
895< // TODO: WHAT IS A 'Malgrange ideal'?
896<
897< // TODO: DO YOU REALLY COMPUTE THE global Bernstein-Sato polynomial
898< // according to the algorithm by Masayuki Noro?
899< // OR ONLY ROOTS?
900<
901<
9021133c1255
903< "USAGE:  bfctsyz(f [,r,s,t,u,v]);  f a poly, optional  r,s,t,u ints, v an intvec
904---
905> "USAGE:  bfctSyz(f [,r,s,t,u,v]);  f a poly, r,s,t,u optional ints, v an optional intvec
9061135,1141c1257,1264
907< PURPOSE: computes the roots (as an ideal, together with their multiplicities
908< @*       as an intvec) of the global Bernstein-Sato polynomial bs(s) for the
909< @*       hypersurface, defined by f, computed according to the algorithm by
910< @*       Masayuki Noro.
911< BACKGROUND: In this proc, the initial Malgrange ideal is computed and than
912< @*       a system of linear equations is solved by computing syzygies.
913< NOTE:    If r<>0, @code{std} is used for GB computations in characteristic 0,
914---
915> PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
916> @*       for the hypersurface defined by f
917> ASSUME:  The basering is a commutative polynomial ring in char 0.
918> BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm
919> @*       by Masayuki Noro and then a system of linear equations is solved by computing syzygies.
920> NOTE:    In the output list, the ideal contains all the roots
921> @*       and the intvec their multiplicities.
922> @*       If r<>0, @code{std} is used for GB computations in characteristic 0,
9231153c1276
924< EXAMPLE: example bfct; shows examples
925---
926> EXAMPLE: example bfctSyz; shows examples
9271197c1320
928<             if (typeof(#[5])=="intvec" && size(#[5])==n && ispositive(#[5])==1)
929---
930>             if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1)
9311219,1225d1341
932<
933< // TODO: WRONG NAME : RENAME TO 'bfctIdeal' (see 3.9.1 Procedures in a library, 5^th rule)
934< // TODO: MORE ASSUMPTIONS?
935< // TODO: ALGORITHM REFERENCE?
936< // TODO: EXAPLAIN WHAT IS A b-function of I w.r.t. a weight vector! IS IT UNIQUE?
937< // TODO: MAYBE MORE EXPLICITE THEORY?
938<
9391229,1235c1345,1352
940< ASSUME:  Assume, I is an ideal in the basering, which is the n-th Weyl algebra,
941< @*       where the sequence of the variables is x(1),...,x(n),D(1),...,D(n).
942< PURPOSE: computes the roots (as an ideal, together with their multiplicities as
943< @*       an intvec) of the global b-function of I wrt the weight vector (-w,w),
944< @*       computed according to the algorithm by Masayuki Noro.
945< BACKGROUND: // TODO: EXPLAIN HERE!
946< NOTE:    If s<>0, @code{std} is used for GB computations in characteristic 0,
947---
948> PURPOSE: computes the roots of the global b-function of I wrt the weight vector (-w,w).
949> ASSUME:  The basering is the n-th Weyl algebra in char 0, where the sequence of
950> @*       the variables is x(1),...,x(n),D(1),...,D(n).
951> BACKGROUND:  In this proc, the initial ideal of I is computed according to the algorithm by
952> @*       Masayuki Noro and then a system of linear equations is solved by linear reductions.
953> NOTE:    In the output list, the ideal contains all the roots
954> @*       and the intvec their multiplicities.
955> @*       If s<>0, @code{std} is used for GB computations in characteristic 0,
9561264c1381
957<   ideal J = initialideal(I,-w,w,whichengine,methodord);
958---
959>   ideal J = initialIdealW(I,-w,w,whichengine,methodord);
9601280c1397
961<   ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6;
962---
963>   ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6; I = std(I);
9641289,1297c1406,1407
965< // TODO: WRONG NAME : RENAME TO 'bfctOneStep' (see 3.9.1 Procedures in a library, 5^th rule)
966<
967< // TODO: ASSUMPTIONS?
968< // TODO: ALGORITHM REFERENCE?
969< // TODO: MAYBE MORE EXPLICITE THEORY?
970<
971<
972< proc bfctonestep (poly f,list #)
973< "USAGE:  bfctonestep(f [,s,t]);  f a poly, s,t optional ints
974---
975> proc bfctOneGB (poly f,list #)
976> "USAGE:  bfctOneGB(f [,s,t]);  f a poly, s,t optional ints
9771299,1303c1409,1416
978< PURPOSE: computes the roots (as an ideal, together with their multiplicities
979< @*       as an intvec) of the global Bernstein-Sato polynomial bs(s) for the
980< @*       hypersurface, defined by f, using only one GB computation
981< BACKGROUND: // TODO: EXPLAIN HERE!
982< NOTE:    If s<>0, @code{std} is used for the GB computation, otherwise,
983---
984> PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the
985> @*       hypersurface defined by f, using only one GB computation
986> ASSUME:  The basering is a commutative polynomial ring in char 0.
987> BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the
988> @*       algorithm by Masayuki Noro and combined with an elimination ordering.
989> NOTE:    In the output list, the ideal contains all the roots
990> @*       and the intvec their multiplicities.
991> @*       If s<>0, @code{std} is used for the GB computation, otherwise,
9921309c1422
993< EXAMPLE: example bfctonestep; shows examples
994---
995> EXAMPLE: example bfctOneGB; shows examples
9961314a1428
997>   int noofvars = 2*n+4;
9981332c1446,1482
999<   def DDh = initialidealengine("bfctonestep", whichengine, methodord, f);
1000---
1001>   intvec uv;
1002>   uv[n+3] = 1;
1003>   ring r = 0,(x(1..n)),dp;
1004>   poly f = fetch(save,f);
1005>   uv[1] = -1; uv[noofvars] = 0;
1006>   //  for the ordering
1007>   intvec @a; @a = 1:noofvars; @a[2] = 2;
1008>   intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0;
1009>   if (methodord == 0) // default: block ordering
1010>   {
1011>     ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(noofvars-1),lp(1));
1012>   }
1013>   else // M() ordering
1014>   {
1015>     intmat @Ord[noofvars][noofvars];
1016>     @Ord[1,1..noofvars] = uv;
1017>     @Ord[2,1..noofvars] = 1:(noofvars-1);
1018>     for (i=1; i<=noofvars-2; i++)
1019>     {
1020>       @Ord[2+i,noofvars - i] = -1;
1021>     }
1022>     dbprint(ppl,"weights for ordering:",transpose(@a));
1023>     dbprint(ppl,"the ordering matrix:",@Ord);
1024>     ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord));
1025>   }
1026>   dbprint(ppl,"the ring Dh:",Dh);
1027>   // non-commutative relations
1028>   matrix @relD[noofvars][noofvars];
1029>   @relD[1,2] = t*h^2;// s*t = t*s+t*h^2
1030>   @relD[2,n+3] = Dt*h^2;// Dt*s = s*Dt+h^2
1031>   @relD[1,n+3] = h^2;
1032>   for (i=1; i<=n; i++)
1033>   {
1034>     @relD[i+2,n+3+i] = h^2;
1035>   }
1036>   dbprint(ppl,"nc relations:",@relD);
1037>   def DDh = nc_algebra(1,@relD);
10381334,1336c1484,1490
1039<   dbprint(ppl, "the initial ideal:", string(matrix(inF)));
1040<   intvec tonselect = 1;
1041<   for (i=3; i<=2*n+4; i++)
1042---
1043>   dbprint(ppl,"computing in ring",DDh);
1044>   ideal I;
1045>   poly f = imap(r,f);
1046>   kill r;
1047>   f = homog(f,h);
1048>   I = t - f, t*Dt - s;  // defining the Malgrange ideal
1049>   for (i=1; i<=n; i++)
10501338c1492
1051<     tonselect = tonselect,i;
1052---
1053>     I = I, D(i)+diff(f,x(i))*Dt;
10541340,1342c1494,1506
1055<   inF = nselect(inF,tonselect);
1056<   dbprint(ppl, "generators containing only s:", string(matrix(inF)));
1057<   inF = engine(inF, whichengine); // is now a principal ideal;
1058---
1059>   dbprint(ppl, "starting Groebner basis computation with engine:", whichengine);
1060>   I = engine(I, whichengine);
1061>   dbprint(ppl, "finished Groebner basis computation");
1062>   dbprint(ppl, "I before dehomogenization is" ,I);
1063>   I = subst(I,h,1); // dehomogenization
1064>   dbprint(ppl, "I after dehomogenization is" ,I);
1065>   I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv
1066>   dbprint(ppl, "the initial ideal:", string(matrix(I)));
1067>   intvec tonselect = 1;
1068>   for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; }
1069>   I = nselect(I,tonselect);
1070>   dbprint(ppl, "generators containing only s:", string(matrix(I)));
1071>   I = engine(I, whichengine); // is now a principal ideal;
10721346,1347c1510,1511
1073<   ideal inF = @m(inF);
1074<   poly p = inF[1];
1075---
1076>   ideal I = @m(I);
1077>   poly p = I[1];
10781356,1357c1520,1521
1079<   bfctonestep(f);
1080<   bfctonestep(f,1,1);
1081---
1082>   bfctOneGB(f);
1083>   bfctOneGB(f,1,1);
10841360,1367d1523
1085<
1086< // TODO: WRONG NAME : RENAME TO 'bfctAnn'
1087<
1088< // TODO: ASSUMPTIONS?
1089< // TODO: ALGORITHM REFERENCE?
1090< // TODO: MAYBE MORE EXPLICITE THEORY?
1091<
1092<
10931371,1375c1527,1534
1094< PURPOSE: computes the roots (as an ideal, together with their multiplicities
1095< @*       as an intvec) of the global Bernstein-Sato polynomial bs(s) for the
1096< @*       hypersurface, defined by f
1097< BACKGROUND: In this proc, ann(f^s) is computed. // TODO: EXPLAIN MORE!
1098< NOTE:    If r<>0, @code{std} is used for GB computations,
1099---
1100> PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
1101> @*       for the hypersurface defined by f
1102> ASSUME:  The basering is a commutative polynomial ring in char 0.
1103> BACKGROUND: In this proc, ann(f^s) is computed and then a system of linear
1104> @*       equations is solved by linear reductions.
1105> NOTE:    In the output list, the ideal contains all the roots
1106> @*       and the intvec their multiplicities.
1107> @*       If r<>0, @code{std} is used for GB computations,
11081403c1562,1563
1109< static proc hardexamples () // TODO: HOW LONG DOES IT TAKE ON SOME COMPUTER?
1110---
1111> /*
1112> //static proc hardexamples ()
11131425c1585
1114<
1115---
1116> */