Transform a parse tree of a polynomial to a parse tree of its evaluation according to Horner's scheme

Question

Could you please point me to an algorithm that takes a (binary) parse tree for evaluating a polynomial expression in a single variable and returns an equivalent parse tree that evaluates the polynomial according to Horner's rule.

The intended use case is in expression templates. The idea is that for a matrix x the parse tree obtained from

a + bx + c * x*x + d * x*x*x...

will get optimized into the corresponding parse tree of

a + x *( b + x( c + x*d))

Answer 1

You could use the following tranformation.

Assumption: the parse tree of the polynomial is in the order of increasing exponents -- if this assumption does not hold, the partial polynomes can be swapped around in the parse tree to make the assumption hold

Assumption: the parse tree holds exponential forms of the variable (e.g. x^2) instead of multiplicational forms (e.g. x*x), except for x^0 -- simple transformations can convert between the two in either direction

Assumption: the coefficients (if constant) in the polynom are non-zero -- this is to avoid having to deal with (a+0*x+c*x^2 -> a+x(cx) instead of a+cx^2)

Parse tree for a+b*x^1+c*x^2+d*x^3:

  .+..
 /    \
a   ...+....
   /        \
  *         .+..
 / \       /    \
b  ^      *      *
  / \    / \    / \
 x   1  c   ^  d   ^
           / \    / \
          x   2  x   3

Transformation to a+x(b+c*x^1+d*x^2)

  +
 / \
a   *
   / \
  x   +
     / \
    b   .+..
       /    \
      *      *
     / \    / \
    c   ^  d   ^
       / \    / \
      x   1  x   2

Transformation to a+x(b+x(c+d*x^1))

  +
 / \
a   *
   / \
  x   +
     / \
    b   *
       / \
      x   +
         / \
        c   *
           / \
          d   ^
             / \
            x   1

Then finally (a+x(b+x(c+d*x))):

  +
 / \
a   *
   / \
  x   +
     / \
    b   *
       / \
      x   +
         / \
        c   *
           / \
          d   x

The common transformation is:

 .            ->    .
  .           ->     .
   .          ->      .
    +         ->      .*..
   / \        ->     /    \
  *   N(k+1)  ->    ^      +
 / \          ->   / \    / \
ck  ^         ->  x   k  ck  N'(k+1)
   / \        -> 
  x   k       -> 

where N'(k+1) is the same subtree as N(k+1) with each exponent j replaced with j-k (if k is 1, replace the x^k subtree with x)

The algorithm is then applied again(*) on N'(k+1) until its root is * (instead of +), indicating that the final partial polynom is reached. If the final exponent is 1, replace the exponential part with x (x^1 -> x)

(*) here is the recursion part

Note: if you keep track of the exponent changes in the N(k+1) subtrees cumulatively, you can put the last two steps together to transform N(k+1) recursively in one go

Note: If you want to allow negative exponents, then either

a) extract the highest negative exponent as the first term:

a*x^-2 + b*x^-1 + c + d*x + e*x^2  ->  x^-2(a+b*x+c*x^2+d*x^3+d*x^4)

and apply the above transformation

or b) separate the positive and negative exponential parts of the equation and appply the above transformation on each separately (you will need to "flip" the operands for the negative-exponent side and replace multiplication with division):

a*x^-2 + b*x^-1 + c + d*x + e*x^2  ->  [a+x^-2 + b*x-1] + [c + d*x + e*x^2] ->
-> [(a/x + b)/x] + [c+x(d+ex)]

this approach seems more complicated than a)

Answer 2

You just have to apply the following rules until you can't apply them anymore.

((x*A)*B)     -> (x*(A*B))
((x*A)+(x*B)) -> (x*(A+B)))
(A+(n+B))     -> (n+(A+B))  if n is a number

where A and B are subtrees.

Here is an OCaml code to do this:

type poly = X | Num of int | Add of poly * poly | Mul of poly * poly

let rec horner = function
  | X -> X
  | Num n -> Num n
  | Add (X, X) -> Mul (X, Num 2)
  | Mul (X, p)
  | Mul (p, X) -> Mul (X, horner p)
  | Mul (p1, Mul (X, p2))
  | Mul (Mul (X, p1), p2) -> Mul (X, horner (Mul (p1, p2)))
  | Mul (p1, p2) -> Mul (horner p1, horner p2) (* This case should never be used *)
  | Add (Mul (X, p1), Mul (X, p2)) -> Mul (X, horner (Add (p1, p2)))
  | Add (X, Mul (X, p))
  | Add (Mul (X, p), X) -> Mul (X, Add (Num 1, horner p))
  | Add (Num n, p)
  | Add (p, Num n) -> Add (Num n, horner p)
  | Add (p1, Add (Num n, p2))
  | Add (Add (Num n, p1), p2) -> Add (Num n, horner (Add (p1, p2)))
  | Add (p1, p2) -> horner (Add (horner p1, horner p2))