単一の変数で多項式を評価するために (バイナリ) 解析ツリーを取得し、ホーナーの規則に従って多項式を評価する同等の解析ツリーを返すアルゴリズムを教えてください。
意図された使用例は式テンプレートにあります。アイデアは、行列x
の場合、解析ツリーは次から取得されるということです
a + bx + c * x*x + d * x*x*x...
の対応する解析ツリーに最適化されます
a + x *( b + x( c + x*d))
単一の変数で多項式を評価するために (バイナリ) 解析ツリーを取得し、ホーナーの規則に従って多項式を評価する同等の解析ツリーを返すアルゴリズムを教えてください。
意図された使用例は式テンプレートにあります。アイデアは、行列x
の場合、解析ツリーは次から取得されるということです
a + bx + c * x*x + d * x*x*x...
の対応する解析ツリーに最適化されます
a + x *( b + x( c + x*d))
次の変換を使用できます。
仮定: 多項式の解析ツリーは、指数が増加する順序になっています。この仮定が成り立たない場合は、部分多項式を解析ツリー内で交換して、仮定を成り立たせることができます。
仮定: 解析ツリーはx^2
、乗法形式 (例 ) ではなくx*x
、変数の指数形式 (例 ) を保持x^0
します。
仮定: 多項式の係数 (定数の場合) は非ゼロです -- これは ( a+0*x+c*x^2
->a+x(cx)
の代わりにa+cx^2
)
の解析ツリーa+b*x^1+c*x^2+d*x^3
:
.+..
/ \
a ...+....
/ \
* .+..
/ \ / \
b ^ * *
/ \ / \ / \
x 1 c ^ d ^
/ \ / \
x 2 x 3
への変換a+x(b+c*x^1+d*x^2)
+
/ \
a *
/ \
x +
/ \
b .+..
/ \
* *
/ \ / \
c ^ d ^
/ \ / \
x 1 x 2
への変換a+x(b+x(c+d*x^1))
+
/ \
a *
/ \
x +
/ \
b *
/ \
x +
/ \
c *
/ \
d ^
/ \
x 1
そして最後に ( a+x(b+x(c+d*x))
):
+
/ \
a *
/ \
x +
/ \
b *
/ \
x +
/ \
c *
/ \
d x
一般的な変換は次のとおりです。
. -> .
. -> .
. -> .
+ -> .*..
/ \ -> / \
* N(k+1) -> ^ +
/ \ -> / \ / \
ck ^ -> x k ck N'(k+1)
/ \ ->
x k ->
ここで、各指数が に置き換えられN'(k+1)
たものと同じサブツリーです (が 1 の場合、サブツリーをに置き換えます)N(k+1)
j
j-k
k
x^k
x
N'(k+1)
アルゴリズムは、ルートが*
( ではなく) になるまで再び適用されます (*) +
。これは、最終的な部分多項式に到達したことを示します。最終指数が 1 の場合、指数部分をx
( x^1
-> x
)に置き換えます。
(*) 再帰部分はこちら
注: サブツリーの指数の変化を累積的に追跡している場合N(k+1)
は、最後の 2 つのステップをまとめてN(k+1)
再帰的に一度に変換できます。
注: 負の指数を許可する場合は、次のいずれかを行います。
a) 最初の項として最大の負の指数を抽出します。
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)
上記の変換を適用します
または b) 方程式の正と負の指数部分を分離し、それぞれに上記の変換を個別に適用します (負の指数側のオペランドを「反転」し、乗算を除算に置き換える必要があります)。
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)]
このアプローチは、a) よりも複雑に思えます。
次のルールを適用できなくなるまで適用する必要があります。
((x*A)*B) -> (x*(A*B))
((x*A)+(x*B)) -> (x*(A+B)))
(A+(n+B)) -> (n+(A+B)) if n is a number
とA
はB
サブツリーです。
これを行うための OCaml コードは次のとおりです。
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))
再帰関数で木の単項式係数を取得できます。これらの係数を変換し、ホーナーの法則に従って式を取得するのは簡単です。
もっと効率的な関数が存在するかもしれませんが、これらの値を計算する単純な再帰関数を提供できます。
まず、式を定式化しましょう。式E
:
E = a0 + a1x + a2x^2 + ... + anx^n
(n+1)-ary タプルとして記述できます。
(a0, a1, a2, ..., an)
次に、2 つの操作を定義します。
追加: 2 つの式E1 = (a0, ..., an)
andが与えられE2 = (b0, ..., bm)
た場合、対応する のタプルE1 + E2
は次のとおりです。
{(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)
つまり、max(n,m)+1
要素があり、i
th要素は (C っぽい構文で) によって計算されます。
i<=n?ai:0 + i<=m?bi:0
乗算: 2 つの式E1 = (a0, ..., an)
とが与えられた場合E2 = (b0, ..., bm)
、 の対応するタプルE1 * E2
は次のとおりです。
E1 * E2 = (a0*b0, a0*b1+a1*b0, a0*b2+a1*b1+a2*b0, ... , an*bm)
つまり、n+m+1
要素があり、i
th要素は次のように計算されます。
sigma over {ar*bs | 0<=r<=n, 0<=s<=m, r+s=i}
したがって、再帰関数は次のように定義されます。
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)
どこ
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
注: の実装ではmul
、からまでj
反復する必要がありますが、次の条件も保持する必要があります。0
i
j <= n (because of one[j])
i-j <= m (because of other[i-j]) ==> j >= i-m
したがって、およびj
から移動できます(以降に等しい)max(0,i-m)
min(i,n)
n
n <= i
疑似コードができたので、実装は難しくありません。タプル型の場合、シンプルstd::vector
で十分です。したがって:
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);
}
注: このコードはテストされていません。バグや不十分なエラー チェックが含まれている可能性があります。コピー&ペーストするのではなく、理解して自分で実装してください。
ここから先は、この関数が提供する値から式ツリーを構築するだけです。