メモ‎ > ‎

Sollyaによる最良近似多項式の算出

概要

Sollya(http://sollya.gforge.inria.fr/)とは、浮動小数点演算ライブラリの作成を補助するツールであり、その機能の一つとして「与えられた式に対する最良多項式近似を求める」というものがある。

最良多項式近似とは、関数 f(x) を多項式 g(x) で近似するにあたり、差の最大値 maxx|f(x) - g(x)| が最小になるような多項式を用いるというものである。もっとも、「任意の x について」ということを考えては、多くの場合は差の最大値が無限大に発散するため、x の範囲をあらかじめ決めたうえで差の最大値を最小化する。一般に多項式近似を行う際には、x の範囲を狭くしたほうが精度を上げるのが容易になるとされる。

以下に例として、1/(1+x) を x∈[0, 1] に対する3次ならびに4次の最良近似多項式で表したものを示す(赤が 1/(1+x)、青がそれの近似)。x∈[0, 1] ではよく近似できているものの、それ以外では離れていくということがわかる。

Sollyaの導入

Sollyaを導入するには、依存するライブラリがいくつかあり、(おそらく)Windows環境ではそれらの導入が面倒である。

筆者はUbuntu環境で試したため、それを前提として説明する。
Sollyaは、Unix系OS向けのソフトウェアでよく見られる、./configure → make → make install でセットアップする形態を採っている。それにあたっては、Sollyaが必要とするライブラリが事前に導入されている必要がある。具体的には http://sollya.gforge.inria.fr/sollya-4.1/sollya.php にある通り、

  • GMP
  • MPFR
  • MPFI
  • fplll
  • libxml2
が必要となる。これはUbuntuの場合(Debianなど、APTが利用できる環境なら同様?)、以下のコマンドで導入できる。

sudo apt-get install libgmp-dev libmpfr-dev libmpfi-dev libfplll-dev libxml2-dev
# sudo
はもし必要なら

これができたら、あとはダウンロードしたパッケージに対して ./configure → make → make install の手順でインストールすることになる。

動かす

「sollya」コマンドでsollyaを起動すると、プロンプトが表示されるので、

remez(1/(1+x), 3, [0;1]);

と入力すると

0.99873734200703370255650542617590628500070645822908 + x * (-0.95079348766754853145575955686649968095374245776006 + x * (0.6862914918332974716537832986049998194156383639214 + x * (-0.23549800416574894019802374173850013846189590616134)))

と表示される。これが 1/(1+x) の x∈[0, 1] に対する3次の最良近似多項式である。また、

g=remez(1/(1+x), 3, [0;1]);
dirtyinfnorm(g-1/(1+x),[0;1]);

と入力すると、

1.26265896698753297463116454110043606729501037121705e-3

と表示される。これは3次の最良近似多項式 g ともとの式 1/(1+x) との無限ノルム(=差の絶対値の最大値)、すなわち近似による誤差の最大値である。誤差が高々0.00126程度であることがこれでわかった。

なお、これらのコマンドは直接入力するのではなく、ファイルに書いて「sollya ファイル名」として呼び出してもよい。

以下に例として、「1/(1+x) の x∈[0, 1] に対するn次の最良近似多項式(3≦n≦10)の近似誤差の最大値を求める」sollyaプログラムを示す。

for n from 3 to 10 do {
  // n次の最良近似多項式を求める
  g = remez(1/(1+x), n, [0;1]);
 
  // 近似誤差の最大値を求める
  d = dirtyinfnorm(g-1/(1+x),[0;1]);
  print("Degree of polynomial:", n, "Maximum error:", d);
};

実行結果は以下の通りであった。

Degree of polynomial: 3 Maximum error: 1.26265896698753297463116454110043606729501037121705e-3
Degree of polynomial: 4 Maximum error: 2.1663794570872775066869973372838816966726104338067e-4
Degree of polynomial: 5 Maximum error: 3.7169201823339135012268617208475286208581962591616e-5
Degree of polynomial: 6 Maximum error: 6.3772288046692996513460503060770401262076583237391e-6
Degree of polynomial: 7 Maximum error: 1.09415974052879634722835986716059189378349083647219e-6
Degree of polynomial: 8 Maximum error: 1.8772818942049454215926652193913745693705257612117e-7
Degree of polynomial: 9 Maximum error: 3.2209096266982243034298389062808592414787706250833e-8
Degree of polynomial: 10 Maximum error: 5.5262141508142946539095377863500835163586891226266e-9

利用可能な関数の一覧は Sollya User's Manual に掲載されている(英語)。

ą
Hiroyuki Hanada,
2015/07/02 2:03
ą
Hiroyuki Hanada,
2015/07/02 2:03
Comments