Converting the syntax tree of a polynomial into a parse tree for its evaluation according to the Horner scheme
Could you please point out an algorithm that takes a (binary) parse tree to compute a polynomial expression in one variable and returns an equivalent parse tree that evaluates the polynomial according to Horner's rule.
The assumed use case is used in expression templates. The idea is that for a matrix the x
parse tree obtained from
a + bx + c * x*x + d * x*x*x...
will be optimized into the appropriate parse tree
a + x *( b + x( c + x*d))
source to share
You can get the monomial coefficients of a tree using a recursive function. Converting these coefficients and obtaining an expression according to Horner's Law would be simple .
I can give you a simple recursive function that calculates these values, although it might be more efficient.
Theoretical material
First, let's formulate the expressions. Expression E
:
E = a0 + a1x + a2x^2 + ... + anx^n
can be written as (n + 1) -the outfit:
(a0, a1, a2, ..., an)
Then we define two operations:
-
Addition: With two expressions
E1 = (a0, ..., an)
and theE2 = (b0, ..., bm)
corresponding set ofE1 + E2
is equal to:{(a0+b0, a1+b1, ..., am+bm, a(m+1), ..., an) (n > m) E1 + E2 = {(a0+b0, a1+b1, ..., an+bn, b(n+1), ..., bm) (n < m) {(a0+b0, a1+b1, ..., an+bn) (n = m)
that is, elements
max(n,m)+1
, and thei
th element is calculated using (with C-ish syntax):i<=n?ai:0 + i<=m?bi:0
-
Multiplication: for two expressions
E1 = (a0, ..., an)
and theE2 = (b0, ..., bm)
corresponding set ofE1 * E2
is:E1 * E2 = (a0*b0, a0*b1+a1*b0, a0*b2+a1*b1+a2*b0, ... , an*bm)
that is the elements
n+m+1
, and thei
th element is calculated usingsigma over {ar*bs | 0<=r<=n, 0<=s<=m, r+s=i}
So the recursive function is defined like this:
tuple get_monomial_coef(node)
if node == constant
return (node.value) // Note that this is a tuple
if node == variable
return (0, 1) // the expression is E = x
left_expr = get_monomial_coef(node.left)
right_expr = get_monomial_coef(node.right)
if node == +
return add(left_expr, right_expr)
if node == *
return mul(left_expr, right_expr)
Where
tuple add(tuple one, tuple other)
n = one.size
m = other.size
for i = 0 to max(n, m)
result[i] = i<=n?one[i]:0 + i<=m?other[i]:0
return result
tuple mul(tuple one, tuple other)
n = one.size
m = other.size
for i = 0 to n+m
result[i] = 0
for j=max(0,i-m) to min(i,n)
result[i] += one[j]*other[i-j]
return result
Note: in an implementation mul
, j
one should iterate from 0
to i
, while the following conditions should also be met:
j <= n (because of one[j])
i-j <= m (because of other[i-j]) ==> j >= i-m
Therefore, it j
can go from max(0,i-m)
and min(i,n)
(which is equal n
, since n <= i
)
C ++ implementation
Now that you have the pseudocode, the implementation doesn't have to be complicated. For a tuple type, simple enough std::vector
. Therefore:
vector<double> add(const vector<double> &one, const vector<double> &other)
{
size_t n = one.size() - 1;
size_t m = other.size() - 1;
vector<double> result((n>m?n:m) + 1);
for (size_t i = 0, size = result.size(); i < size; ++i)
result[i] = (i<=n?one[i]:0) + (i<=m?other[i]:0);
return result;
}
vector<double> mul(const vector<double> &one, const vector<double> &other)
{
size_t n = one.size() - 1;
size_t m = other.size() - 1;
vector<double> result(n + m + 1);
for (size_t i = 0, size = n + m + 1; i < size; ++i)
{
result[i] = 0.0;
for (size_t j = i>m?i-m:0; j <= n; ++j)
result[i] += one[j]*other[i-j];
}
return result;
}
vector<double> get_monomial_coef(const Node &node)
{
vector<double> result;
if (node.type == CONSTANT)
{
result.push_back(node.value);
return result;
}
if (node.type == VARIABLE)
{
result.push_back(0.0);
result.push_back(1); // be careful about floating point errors
// you might want to choose a better type than
// double for example a class that distinguishes
// between constants and variable coefficients
// and implement + and * for it
return result;
}
vector<double> left_expr = get_monomial_coef(node.left);
vector<double> right_expr = get_monomial_coef(node.right);
if (node.type == PLUS)
return add(left_expr, right_expr);
if (node.type == MULTIPLY)
return mul(left_expr, right_expr);
// unknown node.type
}
vector<double> get_monomial_coef(const Tree &tree)
{
return get_monomial_coef(tree.root);
}
Note: This code is untested. It may contain errors or insufficient error checking. Make sure you understand and implement this yourself, not copy.
Now you need to build an expression tree from the values this function gives.
source to share
You can use the following transformation.
Assumption: the parse tree of the polynomial is in ascending order of exponents - if this assumption fails, the partial polynomials can be swapped around in the parse tree to make an assumption
Assumption: The parse tree contains the exponential forms of the variable (for example, x^2
) instead of the multiplicative forms (for example, x*x
), except x^0
- simple transformations can transform between them in any direction
Assumption: the coefficients (if constant) in the polynomial are nonzero - this does not need to be dealt 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
Converting to a+x(b+c*x^1+d*x^2)
+
/ \
a *
/ \
x +
/ \
b .+..
/ \
* *
/ \ / \
c ^ d ^
/ \ / \
x 1 x 2
Converting 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
General conversion:
. -> .
. -> .
. -> .
+ -> .*..
/ \ -> / \
* N(k+1) -> ^ +
/ \ -> / \ / \
ck ^ -> x k ck N'(k+1)
/ \ ->
x k ->
where N'(k+1)
is the same subtree as and N(k+1)
, with each exponent j
replaced with j-k
(if k
equal to 1, replace subtree x^k
with x
)
The algorithm is then applied again (*) on N'(k+1)
until its root is *
(instead of +
), indicating that the final partial polynomial has been reached. If the end exponent is 1, replace the exponential part with x
( x^1
→ x
)
(*) - recurrent part
Note. If you are tracking changes in exponent across subtrees N(k+1)
in aggregate, you can put the last two steps together to transform N(k+1)
recursively in one pass
Note. If you want to allow negative metrics, then either
a) choose 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 transform
or b) separate the positive and negative exponential parts of the equation and apply the above transformation each separately (you need to "flip" the operands for the negative exponent side and replace multiplication by 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)
source to share
You just need to apply the following rules until you can no longer apply them.
((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 the OCaml code for that:
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))
source to share