One Hat Cyber Team

View File Name : paramBpoly.m2

-- Copyright 1999-2002 by Anton Leykin and Harrison Tsai

-----------------------------------------------------------------------------
-- paramBpoly attempts to compute the b-polys for a poly with parameters
-- IN>  f: RingElement (an element in a Weyl algebra with a polynomial ring for a coefficientRing)
-- IN>  file: File (e.g., stdio)
-- OUT> list of all possible Bernstein-Sato polynomials (also write the tree of b-polys into file)
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-- keys for HashTables that are used locally
-----------------------------------------------------------------------------
protect children
protect ringOpts
protect tempI
protect tempN
protect tempT
protect tempBF
protect poly
protect inv
protect iList
protect tempV'
protect parentNODE
protect Itype
protect tempV''
protect level
protect coeff
protect cSet
 
-----------------------------------------------------------------------------
-- "lucky" prime & QQflag (says when to use Q)
-----------------------------------------------------------------------------
DBIGPRIME = 32749;
QQflag := false;

-----------------------------------------------------------------------------
-- next comes a ton of small functions needed further...
-----------------------------------------------------------------------------

-- get gens from the paramGB output
gensGB := GB -> matrix{GB#0}
-- Matrix to List
toL := u -> (flatten entries (u))
 
-- QtoZ (transforms Q-polys into Z-polys)
QtoZ := (p, Z) -> (
     c := product(listForm p / (u -> denominator u#1));
     sum(listForm p, (u -> promote(numerator (c*u#1), 
		    coefficientRing Z)*Z_(u#0)))
     );

-- QtoZmodP (transforms Q-polys into Z/P-polys)
QtoZmodP := (p, Z) -> (
     K2 := coefficientRing Z;
     ZmodP := coefficientRing K2; 
     sum(listForm p, (u -> (
		    c := product(listForm u#1 / (v -> denominator v#1));
		    (QtoZ(c*u#1, K2) * (1_K2//c))*Z_(u#0)
		    )))
     );

-- ZtoQ (transforms Z-polys into Q-polys)
ZtoQ := (p, Q) -> (
	       sum(listForm p, (u -> promote(u#1, coefficientRing Q)*Q_(u#0)))
	       );
-- (Z/p to Q)
ZmodP2Q := c -> (
     if class c =!= ZZ then c = substitute(c,ZZ);
     F1 := DBIGPRIME//100;
     F2 := 100;
     
     b := DBIGPRIME;
     i := 1;
     while (abs(i * c % b)) > F1 and i < F2 do i = i+1;
     j := 1;
     while (abs(j * (b-c) % b)) > F1 and j < F2 do j = j+1;
     if j < i then ((j * (b-c) % b) / (-j)) else ((i * c % b) / i) 
     );

-- (applies ZmodP2Q to coeffs of polynomial)
ZmodP2Qpoly := f -> (
     if f == 0 then f else(
     	  A := QQ(monoid [options ring f]);   
     	  sum(listForm f / (u -> (
			 ZmodP2Q(substitute(u#1, ZZ))
			 * A_(u#0))))
     	  )
     );

promoteIdeal := (I, R) -> ( 
     if I == 0 then ideal 0_R
     else ideal mingens ideal ((toL gens I) / (u->promote(u, R)))     
     );	  

liftIdeal := (I, R) -> ( 
     K := (ideal ((toL gens I) / (u->lift(u, R)))); 
     ideal mingens (K + promoteIdeal(ideal presentation ring I, R)) 
     );	  

-- liftList
liftList := (u, A) -> (u / (v -> lift(v, A)) );

--least common multiplier
lcm := l -> (
     if #l == 0 then null
     else ( 
	  I := ideal l#0;
	  i := 1;
	  while i<#l do (
	       I = intersect(I, ideal l#i); i = i + 1
	       );
	  I
	  )
     );

----------------------------------------------------------------------------
-- END of "ton of small functions"
---------------------------------------------------------------------------


-------------------------------------------------------------------           
-- calculateAss : (I) -> (list)
-------------------------------------------------------------------
-- I = ideal
-- list = the list of minimal primes associated to I
calculateAss := (I) -> (
     pInfo(2,"calculating minimal primes associated to I = " | 
	  toString I | "...");  
     R := ring I;
     local r;
     if isQuotientRing(R) then (	  
	  J := ideal presentation R;
	  S := ring J;
	  J = J + liftIdeal(I, S); 
	  r = (calculateAss J) / (u -> promoteIdeal(u, R)); 
	  ) 
     else (
	  local l;
	  if QQflag then (
	       --for QQ
	       Z := ZZ(monoid [(entries vars R)#0]);
	       
	       I' := toL(gens I) / (u->QtoZ(u,Z));
	       I' = ideal I';
	       decomp := minimalPrimes(I');
	       l = decomp / (u -> (
			 Ide := ideal (toL(gens u) / (u->ZtoQ(u,R)));
		    	 Ide
		    ));
	       )
	  else ( 
	       -- for ZZ/DBIGPRIME
	       l = minimalPrimes(I);
	       );
	  r = select(l, u -> u != (ideal 1_R));
	  );
     r
     )

factorPolyDropExponents := f -> (
     exp := drop(factor f, -1);
     exp = apply(exp, u -> u#0);
     toString exp
     );


----------------------------------------------------------------------------
-- Data structures needed for b-1.m2
----------------------------------------------------------------------------

-- tree of varieties:
--
-- NODE := (	
--     	    	I: ideal, 
--     	    	Itype: integer (COMPUTED, NOTCOMPUTED 
--                     or the reference for a NODE n such that n#tempI <= I),
--     	    	level: integer,
--     	   	bf: b-polynomial, 
--     	    	inv: elements inverted when computing bf,
--     	    	children: list of NODEs below this NODE,   
--     	    	parentNODE: parent, 
--     	    	t: time used by b-poly algorithm
--     	   ); 

-- Itypes:
COMPUTED := -1;
NOTCOMPUTED := -2;

createVTree := f -> (
     R := ring f;
     new HashTable from {
	  coeff => coefficientRing R,
	  ringOpts => options R,
	  poly => f,
	  tempN => {} 
	  }  
     ); 

saveVTree = method()
saveVTree(HashTable, String) := (V, filename) -> ( 
     filename << toString V#coeff << endl << toString V#ringOpts << endl
     << toString V#poly << endl << toString V#tempN	<< endl )

deleteDollars = method()
deleteDollars(String) := (s) -> (
     cs := characters s;
     ret := "";
     scan(cs, c->(if c!="$" then ret=ret|c));
     ret 
     ) 
loadVTree = method()
loadVTree(String) := filename -> (
     f := lines get filename;
     f = apply(f, deleteDollars);
     K := value f#0;
     opts := value f#1;
     A := K(monoid [opts]);
     new HashTable from {
	  coeff => K,
	  ringOpts => opts,
	  poly => value f#2,
	  tempN => value f#3 
	  }
     )

search4ABigGuy := (V, I, pNODE) -> (
     dad := V#tempN#pNODE;
     grandpa := dad#parentNODE;
     pInfo(3, {"grandpa: ideal =", I, " parent = ", pNODE});
     if grandpa < 0 then NOTCOMPUTED
     else (
	  pInfo(3, "grandpa children: " | 
	       toString ((V#tempN#grandpa)#children / (i->
			 (i, isSubset((V#tempN#i)#tempI, I)))));
	  temp := select(1, (V#tempN#grandpa)#children, 
--	       i-> (i < parentNODE) and isSubset((V#tempN#i)#tempI, I) ); 
	       i-> (i < pNODE) and isSubset((V#tempN#i)#tempI, I) ); 
	  if #temp > 0 then temp#0
	  else search4ABigGuy(V, I, grandpa)
	  ) 
     );

add2VTree := (V, II, p) -> (
     pInfo(3, "add2VTree"|toString {V,II,p});
     new HashTable from {
	  coeff => V#coeff,
	  ringOpts => V#ringOpts,
	  poly => V#poly,
	  tempN  => append(V#tempN, new HashTable from {
	       	    tempI => II,
		    Itype => (
			 if p < 0 then NOTCOMPUTED 
		    	 else search4ABigGuy(V, II, p)
			 ),
		    level => if p < 0 then 0 else (V#tempN#p)#level + 1,
		    parentNODE => p
		    })  
	  } 
     );

modifyNODE := (V, num, l) -> (
     ---------(V, num, {bf, inv, children, t})
     kidNumbers := l#2 / (J -> (
	       V = add2VTree(V, J, num);
	       #(V#tempN) - 1     -- the number of the last element
	       )); 
     new HashTable from {
	  coeff => V#coeff,
	  ringOpts => V#ringOpts,
	  poly => V#poly,
	  tempN => take(V#tempN, num) | {new HashTable from {
	       	    tempI => (V#tempN#num)#tempI,
		    Itype => COMPUTED,
		    level => (V#tempN#num)#level,
		    tempBF => l#0,
		    inv => l#1,
		    children => kidNumbers,
		    parentNODE => (V#tempN#num)#parentNODE,
		    tempT => l#3
	       	    }} | drop(V#tempN, num + 1)   
	  } 
     );

modifyCSET := (V, num, l) -> (
     old := V#tempN#num;
     new HashTable from {
	  coeff => V#coeff,
	  ringOpts => V#ringOpts,
	  poly => V#poly,
	  tempN  => take(V#tempN, num) | {new HashTable from {
	       	    tempI => old#tempI,
		    Itype => old#Itype,
		    level => old#level,
		    tempBF => old#tempBF,
		    inv => l,
		    children => old#children,
		    parentNODE => old#parentNODE,
		    tempT => old#tempT
	       	    }} | drop(V#tempN, num + 1)   
	  } 
     );


cutCrapOut := l -> (
     i := 0;
     while i < #l do (
	  l' := drop(l, {i,i});
	  if any(l', j -> isSubset(j, l#i)) then (
		    << "throwing away " << l#i << endl;
		    l = l';
		    )
	  else i = i + 1;
     	  );
     l
     );

refineVTree := (V, num) -> (
     temp := V#tempN#num;
     if temp#Itype < 0 and #temp#children != 0 then (
	  scan(temp#children, u -> (
		    V = refineVTree(V, u);
	       	    temp2 := V#tempN#u;
		    if temp2#Itype < 0 and temp2#tempBF == temp#tempBF then (
		    	 << "patching node #" << num << endl;
			 p := position(temp#inv, v-> v == temp2#tempI);
			 if p===null then error "VTree is damaged!";
			 V = modifyCSET(V, num, drop(temp#inv, {p,p}) | temp2#inv);
			 temp = V#tempN#num;
			 );
	       ));
     	  V = modifyCSET(V, num, cutCrapOut(temp#inv));   
	  ); 
     V
     );

appendCSet := (c, l) -> (
     i := 0;
     isDifferent := true;
     while isDifferent do (
	  if (c#i)#tempV' == l#tempV' then isDifferent = false
	  else i = i + 1; 
	  );
     if isDifferent then take(c,i) | {l} | drop(c,i)
     else take(c,i) | {new HashTable from{
		    tempV'' => cutCrapOut (l#tempV'' | (c#i)#tempV''),
		    tempV' => l#tempV'
		    }} | drop(c,i+1)
     );

chooseMinList := l -> (
     select(l, u->(
	       all(l, v->(
			 (not isSubset(v,u) )
			 or v == u
			 ))
	       ))
     );
chooseMinListList := (l1, l2) -> (
     chooseMinList select(l1, u->any(l2, v-> isSubset(v,u)))
     );
 
VTree2bSet2 := V -> (
     --Itype = symbol Itype;
     if any(V#tempN, u -> (u#Itype == NOTCOMPUTED)) then
     error "VTree is not complete"; 
     
     -- for each b(s) \in B(n,d)
     -- make the list of ideals for which b(s) is a generic b-function 
     r := {};
     scan(V#tempN, node->( 
	       if node#Itype == COMPUTED then 
	       (
		    t := position(r, u ->u#tempBF == node#tempBF);
		    r = (
			 if t =!= null 
		    	 then take(r, t) | { 
			      new HashTable from {
			      	   tempBF => node#tempBF,
			      	   iList => (r#t)#iList | {
					if node#tempI == 0 then ideal 0_(V#coeff)
			      		else node#tempI  
			      		}
				   }
			      } | drop (r, t+1)
		    	 else r | { 
			      new HashTable from {
				   tempBF => node#tempBF,
			      	   iList => {
					if node#tempI == 0 then ideal 0_(V#coeff)
			      		else node#tempI  
					} 
				   }
			      }
			 ); 
     	       	    )
	       ));      
     i := 0;
     while i < #r do (
	  l1 := (r#i)#iList;
	  l2 := flatten(drop(r, {i,i}) / ( u -> u#iList));
	  c := { new HashTable from {
	       	    tempV' => {"not empty"},
		    tempV'' => {ideal 0_(V#coeff)} 
		    }};
	  while #((last c)#tempV') > 0 do (
	       -- choose min ideals in arg1 containing 
	       -- at least one ideal from arg2 
	       temp := chooseMinListList(l1, (last c)#tempV'');
	       c = c | { new HashTable from {
			 tempV' => temp,
			 tempV''=> chooseMinListList(l2, temp)
			 }}; 
	       );
     	  r = take(r,i) | {  new HashTable from {
			      	   tempBF => (r#i)#tempBF,
			      	   iList => (r#i)#iList,
				   cSet => take(c,{1,#c-2}) 
				   } 
			      } | drop(r, i+1); 
     	  i = i + 1; 
	  );
     r
     ); 

isNonEmpty := h -> ( #h#tempV'' == 0 or not isSubset( intersect(h#tempV''), h#tempV') );


------------------------------------
-- TeX output
------------------------------------     
Sequence2String := s -> (
          if # s === 1 then concatenate("(",toString s#0,")")
          else concatenate("(",between(",",toString \ s),")")
          )


VSet2tex := (l, TeXfile) -> (
     j := 0;
     while j < #l do(
	  vv := l#j;
	  TeXfile << "$V" 
	  << Sequence2String (first entries gens vv / ZmodP2Qpoly)
	  << (if j < #l - 1 then " \\cup $ " else "$");
	  j = j + 1;
	  );  
     );

factorBFunctionZmodP = method()
factorBFunctionZmodP RingElement := Product => f -> (
     R := ring f;
     f = select(factor f, u->first degree u#0 > 0);
     f = select(f, u->first degree u#0 > 0);
     S := symbol S;
     QR := QQ[S];
     result := apply(f, u->(
	       if first degree u#0 != 1 then error "internal error: incorrect b-function";
	       coeff := listForm u#0 / (v->v#1);
	       pInfo(666, {"coeff = ", coeff});
     	       Power(QR_0 + (if #coeff> 1 then ZmodP2Q(coeff#1//coeff#0) else 0), u#1)
	       ));
     result
     );-- end factorBFunctionZmodP

BSet2tex := (s, TeXfile) -> (
     TeXfile << "\\begin{itemize}\n"; 
     s / (u -> ( 
	       TeXfile << "\\item $b(s)=";
	       TeXfile << (
		    -- b-function
		    S := symbol S; 
		    R := (ZZ/DBIGPRIME)[S];
		    toString factorBFunctionZmodP sum(u#tempBF, 
			v->(if class v#1 === ZZ 
			     -- hashtable was loaded from file
			     then v#1
			     -- original hashtable
			     else substitute(v#1,ZZ)) 
			* R_(v#0))
	       ) << "$ corresponds to" << endl;
	       i := 0;
	       while i < #u#cSet do( 
		    VSet2tex((u#cSet#i)#tempV', TeXfile);
	       	    if #(u#cSet#i)#tempV'' > 0 then TeXfile << "$\\setminus$";
	       	    VSet2tex((u#cSet#i)#tempV'', TeXfile);
	       	    if i < #u#cSet - 1 then TeXfile << "$\\cup$";
		    i = i + 1;
		    );
		TeXfile << endl;	    
	       ));
     TeXfile << "\\end{itemize}" << endl;
     );

------------------------------------------------------------------------------
-- ***************************************************************************
-- globalBFunction (2 versions) adapted for b-1
-- ***************************************************************************
------------------------------------------------------------------------------
getChangeMat := g->g#1
getInv := GB -> (
     l := flatten entries getChangeMat GB;
     l = l / (u -> (entries transpose (coefficients u)#1)#0);
     l = (flatten l) / (u -> lift(u, coefficientRing ring GB#0#0));
     select(l / (u -> denominator u), u -> not isUnit u)
     ); 
-- a local function used by gbW and inW  

homGBparam := (I, w) -> (
     W := ring I;
     createHomWeylAlgebra W;
     HW := W.HomWeylAlgebra;
     dpairs := W.monoid.Options.WeylAlgebra;
     -- Do some sanity checking
     if dpairs === {} 
     then error "expected a Weyl algebra";
     if any(dpairs, v -> class v =!= List)
     then error "expected non-homogenized Weyl algebra";
     if #w =!= numgens W
     then error ("expected weight vector of length " | numgens W);
     -- Make the new weight vector
     wts := prepend(-1,w);

     -- Homogenize I
     I1 := W.WAtoHWA I;
     homogenize(gens I1, HW_0, wts)
     );

-- computes the GB of ideal "I" with respect to weight vector "w"  
gbWparam = method()
gbWparam(Ideal, List) := (I, w) -> (
     I2 := homGBparam(I, w);
     -- Do the computation *************************************************
     ideal compress (ring I).HWAtoWA gens paramGB I2
     );-- end gbW

-- computes in_w(I). The result is a WA ideal (as opposed to 
-- an ideal in the associated commutative ring) 
inWparam = method()
inWparam(Ideal,List) := (I, w) -> (
     I2 := homGBparam(I, w);
     GB := paramGB (I2);    
     I3 := leadTerm(1, gensGB GB);
     (ideal compress (ring I).HWAtoWA I3, getInv(GB)) 
     );

local globalBFunctionParam;
local globalBFunctionParam2;
----------------------------------------------------------------------------
-- This function calls 2 versions of globalBFunction with parameters
BSWITCH := 1
setBSwitch = method()
setBSwitch ZZ := n -> (
     if n!=1 and n!=2 then error "invalid BSWITCH";
     BSWITCH = n;     
     )
CALLglobalBFunctionParam = method()
CALLglobalBFunctionParam(Thing) := a -> (
     if BSWITCH == 1 then globalBFunctionParam a
     else globalBFunctionParam2 a 
     );

---------------------------------------------------------------------------
 
makeMonic := f -> ( (1 / (leadCoefficient f)) * f );
ZZtoQQ := u -> (
     c := substitute(lift(u, ZZ/DBIGPRIME), ZZ);
     b := DBIGPRIME;
     i := 1;
     while (abs(i * c % b)) > 100 and i < 100 do i = i+1;
     ((i * c % b) / i) 
     );

ZZmakeQQ := f -> (
     Q := QQ(monoid [first entries vars ring f]);
     try sum(listForm f / (u-> ZZtoQQ(u#1) * Q_(u#0)))
     else f
     );

makeQQ := f -> (
     s := symbol s;
     R := QQ[s];
     if not member(QQ, R.baseRings) then
     error "QQ is not a base ring of R";
     sum(listForm f, u -> lift(u#1, QQ) * s^(sum u#0))
     );

---------------------------------------------------------------------
-- 1st version of globalB with parameters 
---------------------------------------------------------------------
bFunctionParam := (I, w) -> (
     if not (ring I).?ThetaRing then
     createThetaRing (ring I);

     if not (ring I).ThetaRing.?IntRing then
     createIntRing (ring I).ThetaRing;
     
     createDpairs (ring I);

     W := ring I;
     T := W.ThetaRing;
     TI := T.IntRing;
     dpV := W.dpairVars;
     dpI := W.dpairInds;
     
     -- sanity check
     if (#(dpI#2) != 0) then
     error "expected no central variables in Weyl algebra";
     if (#w != ((numgens W) // 2)) then
     error "expected weight vector of length " | ((numgens W) // 2);
     
     w = apply(numgens W, i -> (
	       p := position(dpI#1, u -> u == i); 
	       if p =!= null  then w#p
	       else (
		    p = position(dpI#0, u -> u == i);
		    -w#p
		    )  
	       ));
     -- compute in_(-w,w) (I) 
     temp := inWparam(I, w); 
     inv := temp#1;
     inI := temp#0;
     n := #(dpI#0);
     eulerOp := sum(n, i -> w_(dpI#1#i)*(dpV#0#i)*(dpV#1#i));
     	  
     (     
	  -- NON-GENERIC
	  dpI' := {select(dpI#0, i -> w#i != 0), 
	       select(dpI#1, i -> w#i != 0)};
	  dpI'' := {select(dpI#0, i -> w#i == 0), 
	       select(dpI#1, i -> w#i == 0)};
	  -- want: eliminate all u_i, v_i as well as x_i, dx_i of weight 0
	  u := symbol u;
	  v := symbol v;	  
	  UV := (coefficientRing W)(monoid [ (dpI'#0) / (i -> u_i), 
	       (dpI'#1) / (i -> v_i), 
	       (dpI''#0) / (i -> W_i), (dpI''#1) / (i -> W_i),
	       (dpI'#0) / (i -> W_i), (dpI'#1) / (i -> W_i),
	       WeylAlgebra => W.monoid.Options.WeylAlgebra,
	       MonomialSize => 16,
	       MonomialOrder => Eliminate (2 * #dpI#0) ]);
	  WtoUV := map(UV, W, matrix { apply(numgens W, i -> (
			      if member(i, dpI'#0) then 
			      (UV_(u_i) * substitute(W_i, UV))
			      else if member(i, dpI'#1) 
			      then (substitute(W_i, UV) * UV_(v_i))
			      else substitute(W_i, UV)
			      )) 
		    });
	  
	  GB :=  paramGB ((WtoUV inI) 
	       	    + ideal apply(#dpI'#0, 
			 i -> (UV_(u_(dpI'#0#i)) * UV_(v_(dpI'#1#i)) - 1)
		    	 ));
	  inv = inv | getInv(GB);
	  intGB := gensGB GB;
	  intIdeal := ideal substitute(selectInSubring(1, intGB), W);
	  );
     	  -- compute J = intIdeal \cap K[\theta]
     	  genJ := (flatten entries gens intIdeal) / W.WtoT;
	  elimIdeal := (T.RtoIR ideal genJ)
     	  + ideal(TI_(numgens TI - 1) - T.RtoIR W.WtoT eulerOp);
     	  
	  GB = paramGB(elimIdeal);
	  inv = inv | getInv(GB);
	  elimIdealGB := gensGB GB;
	  -- take the generator of J and cook up the b-function
	  bfcn := (mingens ideal selectInSubring(1,elimIdealGB))_(0,0);
	  s := symbol s;
	  TItoS := map((coefficientRing TI)[s], TI, 
	       matrix{toList(numgens T:0) | {s}});
	  bfcn = TItoS bfcn;
	  --(if QQflag 
	  --     then makeQQ makeMonic bfcn
	  --     else ZZmakeQQ makeMonic bfcn, 
	  --     
	  --     inv)
	  (bfcn, inv)     	
     );-- end of bFunctionParam

globalBFunctionParam = f -> (
     W := ring f;
     createDpairs W;
     dpI := W.dpairInds;
     
     -- sanity check
     if (#(W.dpairInds#2) != 0) then
     error "expected no central variables in Weyl algebra";
     if any(listForm f, m -> any(dpI#1, i -> m#0#i != 0)) then
     error "expected no differentials in the polynomial";
     
     t := symbol t;
     dt := symbol dt;     
     WT := (coefficientRing W)(monoid [ t, dt, (entries vars W)#0,
	       MonomialSize => 16,
	       WeylAlgebra => W.monoid.Options.WeylAlgebra | {t => dt}]);
     t = WT_t;
     dt= WT_dt;
     w := {1} | toList (((numgens W) // 2):0);
     f' := substitute(f,WT);
     If := ideal ({t - f'} 
     	  | (dpI#1 / (i->(
	       	    	 DX := WT_(i+2);
	       	    	 (DX * f' - f' * DX) * dt + DX
	       	    	 )))
	  );
     temp := bFunctionParam(If, w);
     bfunc := temp#0;
     inv := temp#1;
     s := (ring bfunc)_0;
     (makeMonic substitute(bfunc, { s => -s - 1 }), inv)
     );--end of globalbFunctionParam

-------------------------------------------------------------
-- 2nd version of globalB with parameters 
--------------------------------------------------------------
AnnFsParam := f -> (
     pInfo(1, "computing AnnFsParam... ");
     W := ring f;
     K := last (coefficientRing W).baseRings;
     np := numgens K; 
     ord := W.monoid.Options.MonomialOrder;
     -- number of parameters
     n := numgens W;
              	       
     createDpairs W;
     dpI := W.dpairInds;
     dpV := W.dpairVars;
     I := ideal dpV#1;
     
     t := symbol t;
     dt := symbol dt;
     WAopts := W.monoid.Options.WeylAlgebra | {t => dt};
     WT := (coefficientRing W)(monoid [ t, dt, (entries vars W)#0, 
	  WeylAlgebra => WAopts,
	  MonomialOrder => Eliminate 2 ]);
     u := symbol u;
     v := symbol v;
     WTUV := (coefficientRing W)(monoid [ u, v, t, dt, (entries vars W)#0,
	  WeylAlgebra => WAopts,
	  MonomialOrder => Eliminate 2 ]);
     WtoWTUV := map(WTUV, W, (vars WTUV)_{4..n+3});
     -- twist generators of I into generators of KI
     f' := substitute(f,WTUV);
     twistList := apply( toList(3..n+3), 
	  i -> WTUV_i + (WTUV_i*f' - f'*WTUV_i)*dt);
     twistMap := map(WTUV, WTUV, matrix{{u,v,t-f'}|twistList});
     tempKI := twistMap(ideal t + WtoWTUV I);
     wts := {1,-1,1,-1} | toList(n:0);
     KI := ideal homogenize(gens tempKI, u, wts);
     
     g := (entries gens KI)#0 | { u * v - 1 };
     GB := paramGB(ideal g);
     inv := getInv GB;
     preGens := flatten entries substitute(
	  selectInSubring(1, gensGB GB), WT);
     s := symbol s;
     WS := (coefficientRing W)(monoid [(entries vars W)#0, s,
	  WeylAlgebra => W.monoid.Options.WeylAlgebra]);
     WTtoWS := g -> (
	  e := first exponents leadMonomial g;
	  if e#0 > e#1 then g = dt^(e#0-e#1) * g
	  else g = t^(e#1-e#0) * g;
	  g' := 0_WS;
	  while (d := first exponents leadMonomial g; d#0 * d#1 != 0) do(
	       c := leadCoefficient g;
	       g' = g' + c * (-s-1)^(d#1) * WS_(drop(d, 2) | {0}); -- >%-0	
	       g = g - c * (t*dt)^(d#1) * WT_({0,0} | drop(d, 2));
	       ); 
	  g' + substitute(g, WS)
	  );
     pInfo(666, {"AnnFSparam = ", (preGens / WTtoWS)});
     (ideal (preGens / WTtoWS), inv) 
     )

globalBFunctionParam2 = f -> (
     W := ring f;
     
     temp := AnnFsParam f;
     AnnI := temp#0;
     inv := temp#1;
     
     Ws := ring AnnI;
     ns := numgens Ws;
     
     
     elimWs := (coefficientRing Ws)(monoid [(entries vars Ws)#0,
	  WeylAlgebra => Ws.monoid.Options.WeylAlgebra,
	  MonomialOrder => Eliminate (ns-1)]);
     ff := substitute(f,elimWs);
     elimAnnI := substitute(AnnI, elimWs);
     H := (gens elimAnnI) | matrix{{ff}};
     
     gbH := paramGB(H);
     inv = inv | getInv gbH; 
     
     bpolys := selectInSubring(1, gensGB gbH);
     if (bpolys == 0) then error "module not specializable";
     if (rank source bpolys > 1) then error "ideal principal but not
     realized as such.  Need better implementation";
     bpoly := bpolys_(0,0);
     
     Ks := (coefficientRing W)(monoid [Ws_(ns-1)]);
     bpoly = substitute(bpoly, Ks);
     (bpoly, inv)
     )

---------------------------------------------
-- takes care of a node in a VTree
---------------------------------------------
takeCareOf = method()
takeCareOf(HashTable, ZZ) := (V, num) -> (
     pInfo(1, "computing node #" | toString num | "...");
     node := V#tempN#num;
     K := frac(V#coeff / node#tempI);
     A := K(monoid [ V#ringOpts ]); 
     n := numgens A;
     f := sum(listForm V#poly / (u->promote(u#1,K) * A_(u#0)));
     
     timeSpent := timing (bf := CALLglobalBFunctionParam f);
     inv := bf#1;
     bf = bf#0;
          
     assocPrimes := select(
     	  if inv != {} 
     	  then calculateAss liftIdeal(lcm(inv), V#coeff) 
     	  else {},
	  I ->  (
	       TK := frac(V#coeff / I);
     	       TA := TK(monoid [ V#ringOpts ]); 
	       sum(listForm V#poly / (u->promote(u#1,TK) * TA_(u#0))) !=0
	       )
     	  );
     modifyNODE(V, num, {apply(listForm bf, u->(u#0,substitute(u#1,ZZ))), 
	       assocPrimes, assocPrimes, timeSpent#0})
     );

----------
-- MAIN
----------

paramBpoly = method(Options => {"ground field" => 32749}) -- 0 stays for QQ
paramBpoly(RingElement, File) := List => o -> (f, file) -> (
     QQflag = (o#"ground field" == 0);
     if not QQflag then (
	  if isPrime o#"ground field" 
     	  then DBIGPRIME = o#"ground field"
     	  else error "need a prime";
	  )
     else error "algorithm is implemented over finite field so far"; 
     R := ring f;
     K := coefficientRing R;
     if coefficientRing K =!= QQ then 
     error "base ring = QQ expected";
	  
     Z := ((ZZ/DBIGPRIME)(monoid [(entries vars K)#0]))(monoid [
	  (entries vars R)#0, WeylAlgebra => R.monoid.Options.WeylAlgebra]);
     f = QtoZmodP(f,Z);  
     V := createVTree(f);
     V = add2VTree(V, ideal 0_(coefficientRing Z), -1);     
     local num;
     while (num = position(V#tempN, u -> u#Itype == NOTCOMPUTED)) =!= null 
     do ( 
     	  V = takeCareOf(V, num);
     	  );
     -- transform into bSet
     bs := VTree2bSet2 V;
     BSet2tex(bs,file);
     ret := apply(bs, u->(
	       s := symbol s;
	       RZ := (ZZ/DBIGPRIME)[s];
	       factorBFunctionZmodP(
		    sum(u#tempBF, v->
			 (if class v#1 === ZZ 
			     -- hashtable was loaded from file
			     then v#1
			     -- original hashtable
			     else substitute(v#1,ZZ)) 
			* RZ_(v#0))) 
	       ));
     ret
     );


--VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
-- ex-"paramGB.m2" contents
--VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV  

--------------------------------------------------------------
-- Computes parametric GB
--------------------------------------------------------------

--some helpful funcs
--ring(List) := g -> ring g#0
--gens (List) := g -> matrix{g#0}

myLCM := l -> (
     if #l == 0 then null
     else ( 
	  I := ideal l#0;
	  i := 1;
	  while i<#l do (
	       I = intersect(I, ideal l#i); i = i + 1
	       );
	  first first entries mingens I
	  )
     );

in2 := f -> if f!=0 then leadTerm f else 0_(ring f) 

isdivisible := (f,g) -> (
     if (f==0) or (g==0) 
     then false 
     else(
     	  a := listForm f;
       	  b := listForm g;	
       	  all(a#0#0-b#0#0, u->(u>=0))
     	  )
     )

---------------------------------------------------------------
-- paramGB
-- input: I: ideal
--        ChangeMatrix=>false: option
-- output: { Groebner basis generators for I [, changematrix] }

paramGB = I -> (
    if class I === Ideal then I = gens I;
    R := ring I;
    WAflag := R.monoid.Options.WeylAlgebra =!= {};
    if WAflag then (
	 createDpairs R;
	 xVar := R.dpairInds#0;
	 dxVar := R.dpairInds#1;
	 );
    K := last (coefficientRing R).baseRings;
    isQuotFlag := isQuotientRing(K);  
    L := if isQuotFlag then ring presentation K else K;
    varP := first entries vars L;
    np := #varP;
    ord := R.monoid.Options.MonomialOrder;
    newMonomialOrder := ord -> (
	 append(ord, GRevLex=>np));
    Rgens := (options R).Variables;
    Lgens := (options L).Variables;
    
    -- make a new ring with parameters:
    -- The only point of using "xxx", "mmm" was the fact 
    -- that M2 crashes if the old variables are used
    xxx := symbol xxx;
    mmm := symbol mmm;
    NewR := (coefficientRing L) [
	 xxx_0..xxx_(numgens R - 1),
	 mmm_0..mmm_(np-1),  
	 WeylAlgebra => 
	 (if WAflag 
	      then apply(#xVar, i -> (xxx_(xVar#i) => xxx_(dxVar#i)))
	      else {}
	      ),
	 MonomialOrder => newMonomialOrder R.monoid.Options.MonomialOrder
	 ];
    L2NewR := map( NewR, L, toList(0..np-1) / (i->mmm_i) );
    Zero := if isQuotFlag then L2NewR ideal presentation K
    else ideal 0_NewR;
    R2NewR := f -> (
	 l := listForm f;
	 mlcm := myLCM (l / (u -> lift(denominator u#1, L)));
	 l = l / (u -> ( u#0, lift(numerator u#1, L) * (
			mlcm // lift(denominator u#1, L)) ));
	 sum(l, (u -> (L2NewR u#1) * NewR_( toList u#0 | toList(np : 0) )))
	 );
    J := first entries I / R2NewR;
    g := first entries gens gb (ideal J + Zero);
    
    -- go back to the original ring
    NewR2R := f -> (
	 l := listForm f;
	 l = l / (u -> 
	      (drop(toList u#0, -np), u#1 * L_(drop(toList u#0, numgens R))));
	 sum(l, (u -> promote(u#1, K) * R_(u#0)))    
	 );
    g = g / NewR2R;      	  
    
    -- clean up
    i := 0;
    while i < #g do (
	 if g#i == 0_R
	 then (
	      g = drop(g, {i,i})
	      )
	 else i = i + 1; 
	 );
    i = 0;
    while i < #g do (
	 if any(drop(g, {i,i}), u->isdivisible(in2 g#i, in2 u)) 
	 then g = drop(g, {i,i})
	 else i = i + 1; 
	 );
    
    -- construct the matrix of invertibles
    inv := matrix { g / (f -> (
		   (1_K / (leadCoefficient f))*1_R
		   ))};
    
    -- make every generator monic
    g = g / ( u -> makeMonic u);
    
    {g, inv}
    );


TEST ///
A =  (QQ [a,b,c,d]) [x, y, Dx, Dy, WeylAlgebra => {x=>Dx, y=>Dy}]
Dtrace 1
bf = paramBpoly(a*x^2,stdio)
assert(listForm value first bf == {({2},1/1), ({1},3/2), ({0},1/2)})
///