ホーナー法について - いろいろノート

たいしたことないネタでもちょいちょい記事を書いてみることにします。今日はホーナー法について。

多項式を計算する際、乗法の回数を削減することができます。上記リンク先によると、19世紀からある方法のようです。

やり方はとてもシンプルです。以下のような多項式があるとします。

$\displaystyle p(x) = a_0 + a_1x+a_2x^{2}+ \cdots +a_nx^{n}$

これを、次のように変形するだけです。

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

JavaScriptで実装すると、以下のような感じになります。

var horner = (coeffs = [], x) => coeffs.reduce((a, b) => a * x + b, 0);

もしくは、簡潔さより高速さということであれば、以下のような感じでしょうか。

var horner = (coeffs = [], x) => {
  var r = coeffs[0];
  for (let i = 1, l = coeffs.length; i < l; i++)
    r = r * x + coeffs[i];
  return r;
}

第一引数には、 $[a_n, a_{n-1},\cdots,a_0$ ]の順に係数の配列を渡します。

実際のところ、べき乗の計算は効率的なアルゴリズムが知られているため、素直に展開形を実装するのに比べて、ホーナー法でそこまでの高速化は期待できません。ただ、特に小数がある場合、計算精度の向上を望めます。

また、FMA (Fused Multiply-Add) をサポートしている環境では、更なる精度向上が期待できそうです。例えばJuliaはfma()関数を呼び出すことで、CPUが対応していればFMAを使うことができます。

function horner(coeffs, x)
  r = coeffs[1]
  for i in 1:length(coeffs)
    r = fma(r, x, coeffs[i])
  end
  r
end

FMAについては、以下のブログ記事が参考になりました。

ちなみに、JuliaについてはBase.Math中に@horner(x, p...)というマクロが予め定義されているようです。

以下のようにして呼び出すことができます。

julia> using Base.Math: @horner
julia> a0, a1, a2 = 0.123, 0.234, 0.345
julia> x = 0.999
julia> @horner x a0 a1 a2
0.701076345

とても簡単に使えるアルゴリズムなので、ぜひ活用したいです。

以下の本を参考に書きました。