多項式の計算で乗算を増やさない方法

数値計算なんかをプログラムで書いてると、多項式の値を求めたいことがよくある。
たとえば次の多項式

$f(x) = x^2+2x+3$

これを計算するコードは、たとえば次のようになる。

// C/C++
f = x*x + 2*x + 3;

もう少し一般的に、 $n$ 次の多項式

$f(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0$

について、 $x$ と係数 $a_i$ の配列から $f(x)$ を計算する関数を考える。

double poly1(double x, double a[], int n)
{
	double f = 0.0;
	for (int i = 0; i <= n; i++)
	{
		double g = a[i];
		for (int j = 0; j < i; j++)
			g *= x;
		f += g;
	}
	return f;
}

これでもいいのだが、次数の高い多項式を計算しようとすると、この関数の計算量は非常に大きくなる。
理由は簡単で、各 $x^i$ を直接乗算で求めているからだ。
たとえばこの関数で100次の多項式を計算しようとすると、加算100回、乗算5050回ということになる。この関数の乗算の計算回数は $O(n^2)$ である。
また、乗算を何度も繰り返すことになるので、浮動小数点数を使った計算では誤差が問題になる。

というわけで、多項式の計算から乗算の回数を減らしたい。

pow関数を使う

普通のプログラミング言語なら数学ライブラリが付属しているはずで、
その中にある（と思われる） $a^b$ を求めるpowを使えば、計算の回数は減らすことができる。

#include <math.h>

double poly2(double x, double a[], int n)
{
	double f = 0.0;
	for (int i = 0; i <= n; i++)
		f += a[i] * pow(x, (double)i);
	return f;
}

ちなみにpowは、普通次のように実装されている。powが使えなければこの式を代用すればいい。

double pow(double a, double b)
{
	return exp(b * log(a));
}

ただし、powで計算する方法は浮動小数点数が使えない環境や、誤差が出ては困るような状況（整数での計算など）では使うことができない。

高次から順に計算する（ホーナー法）

他に使われる方法としては、計算の順序を変える方法がある。
まず、最初の式を少し変形してみる。

$x^2+2x+3 = ((1)x+2)x+3$

$x$ を順番に括り出しただけで特別なことではないが、式から高次の $x$ が消えている。
一般の多項式に当てはめるなら、

$a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0 = ( (\cdots ((a_n)x + a_{n-1})x + \cdots)x + a_1)x + a_0$

となる。これを一番内側の括弧から順に計算すればいい。

double poly3(double x, double a[], int n)
{
	double f = a[n];
	for (int i = n-1; i >= 0; i--)
		f = f*x + a[i];
	return f;
}

このアルゴリズムはホーナー法（Horner法）と呼ばれているもので、高次の $x$ をうまくまとめることで乗算を減らす方法である。
関数をながめると、乗算が加算と同じ回数まで減少していることが分かる。たとえば100次の多項式なら5050回の乗算が100回になるわけで、速度の差は歴然である。
またpoly2と比べると、特殊な関数を使用していない分、こちらの方がいくらか高速に動くし、整数で計算した場合も誤差は出ない。

ただし、ホーナー法でも乗算の繰り返しを行っているので、浮動小数点数で計算する場合、最初のpoly1ほどではないにしろ誤差は大きくなる。
浮動小数点数で高速な数学関数が使える環境なら、powを使った実装を選ぶほうがいいかもしれない。