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))

      

+3


source to share


3 answers


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 the E2 = (b0, ..., bm)

    corresponding set of E1 + 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 the i

    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 the E2 = (b0, ..., bm)

    corresponding set of E1 * 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 the i

    th element is calculated using

    sigma 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.

+1


source


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)

+3


source


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))

      

+2


source







All Articles