# Multivariate remainder of polynomial in respect to a set of polynomials

I would like to have a really fast routine that computes the so called Normal Form of a multivariate polynomial f in respect to a set of other multivariate polynomials W. I was not able to find a built in function to do this in Mathematica, so I wrote the following code:

NormalForm[f_, W_List] :=
su = Length[W];
r = 0;
p = f;
Do[
, {qq,1, su}];
While[(! (p === 0)),
i = 1;
dividing = True;
While[i < su + 1 && dividing,
i = i + 1;
,
p = (p - Sum[leadp/leadW[i] (W[[i]])[[qq]], {qq, 1, Length[W[[i]]]}]);
dividing = False;
];
];
If[dividing, r = r + leadp; p = p - leadp;]
];
r
];


This implementation is basically a direct Mathematica interpretation of corresponding pseudo code. To find the leading term in each polynomial I use the following LT[p_] function (n is the number of different variables):

LT[p_] := Block[{},
frst = (p - (p /. Subscript[x, 1] -> 0) // Expand) /.Subscript[x, y_] -> Subscript[x, y] a;
If[frst === 0, frst = p /. Subscript[x, y_] -> Subscript[x, y] a;];
scnd = Coefficient[frst, a, Exponent[frst, a]];
MonomialList[scnd, Table[Subscript[x, ind], {ind, 2, n}]][[1]]
];


It might look a bit weird to you because the variable $x_1$ is special. It can only appear as $x_1^0=1$ or $x_1^1=x_1$, and has highest weight in the monomial order. The remaining $x_i$ are ordered in grevlex (a is used to realize this).

The above code works conceptually. However, whenever the set W grows to above 20 or so polynomials (with a couple dozen monomials each) and the integer coefficients of some of the polynomials become thousands of digits long, the code becomes too slow for my practical application. Therefore, I wonder if there is any way to speed it up, which I have missed? A speed increase by a factor of 10 could make my life so much easier. Thanks for any suggestion!

• Check the documentation for PolynomialReduce. It is intended for finding normal forms with respect to a given set of polynomials and term order. – Daniel Lichtblau Nov 20 '15 at 23:29
• Thank you! Turns out, PolynomialReduce is about ten times faster than my code, so hopefully this solves the efficiency issue. I wish they would name functions more intuitively, though! – Kagaratsch Nov 21 '15 at 0:23
• "NormalForm" is used in the Groebner basis literature, but as naming goes it really is less descriptive than "PolynomialReduce". – Daniel Lichtblau Nov 21 '15 at 3:39