0
$\begingroup$

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] := 
  Module[{su, r, p, leadW, leadWM, qq, i, dividing, leadp, leadpM},
   su = Length[W];
   r = 0;
   p = f;
   Do[
    leadW[qq] = LT[W[[qq]]]; 
    leadWM[qq] = leadW[qq]/(leadW[qq] /. Subscript[x, _] -> 1);
   , {qq,1, su}];
   While[(! (p === 0)),
    i = 1;
    dividing = True;
    leadp = LT[p];
    leadpM = leadp/(leadp /. Subscript[x, _] -> 1);
    While[i < su + 1 && dividing,
     If[! (FreeQ[Denominator[leadpM/leadWM[i]], Subscript[x, _]]),
       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!

$\endgroup$
  • 2
    $\begingroup$ Check the documentation for PolynomialReduce. It is intended for finding normal forms with respect to a given set of polynomials and term order. $\endgroup$ – Daniel Lichtblau Nov 20 '15 at 23:29
  • $\begingroup$ 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! $\endgroup$ – Kagaratsch Nov 21 '15 at 0:23
  • 2
    $\begingroup$ "NormalForm" is used in the Groebner basis literature, but as naming goes it really is less descriptive than "PolynomialReduce". $\endgroup$ – Daniel Lichtblau Nov 21 '15 at 3:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.