2016/12/28(水)中点はどうやって計算するべきか

今回は非常に単純な話で、IEEE754のdoubleで表現された数a,bが与えられたとき、その中点(a+b)/2はどう計算するのが最善でしょうか、という問題です。もちろん誤差は入りますが、出来ればnearestにしたい、すなわち、真の(a+b)/2をIEEE754のルールで丸めたものを計算することを目標にします。IEEE754の内部表現は2進数ですから、2で割る操作は指数部を1減らすだけなので普通は誤差は入りません。しかし、極端な場合を考えると案外難しいのです。

(a+b)*0.5

まずは普通に足して2で割る方法。「2で割る」は「0.5を掛ける」にした方がコンパイラに優しいでしょう。ごく普通の方法で安定度抜群ですが、例えば
a=21023
b=21023
の場合、a+bがオーバーフローしてしまう問題があります。

a*0.5+b*0.5

オーバーフローを避けて先に2で割る方法です。計算量は少し増えてしまいますが。こちらは、アンダーフロー領域で問題になります。例えば、
a=2-1074
b=2-1074
のとき、a*0.5、b*0.5はともにアンダーフローで0になってしまい、結果も0になってしまいます。

a+(b-a)*0.5

これもよく見る式です。IEEE754では関係ありませんが、内部演算が10進数の場合にこうやって計算することが推奨されているのを見かけたことがあります。しかし、
a=21023
b=-21023
でオーバーフローしてしまいます。

a+b*0.5-a*0.5

上の式でオーバーフローを避けたものです。今度こそ良さそうですが、
a=5*2-1074
b=7*2-1074
としてみましょう。このとき、IEEE754の偶数丸めルールの関係で、
a*0.5=2*2-1074
b*0.5=4*2-1074
となります。よって計算結果は7*2-1074となってしまい、真の6*2-1074になりません。

a+b*0.5-a*0.5 (ただし上向き丸めで計算)

INTLAB version 9で採用している計算方法です。a*0.5とb*0.5に入る丸めの方向が逆向きになることが無いので、先の例だとうまく計算できます。しかし残念ながら必ずnearestにはならないようで、Mさんがちゃんと計算できない例を発見しました。
a=(253-1)*2-1073
b=7*2-1074
とすると、
(nearest)=2-1021+2-1073
(計算値)=2-1021+2-1072
と値がずれてしまいます。a*0.5は無誤差ですが、b*0.5と加算で両方上向き丸めになった結果、真値(2-1021+2-1073+2-1075)を丸めたものとずれてしまいました。

結論?

kvライブラリの区間の中点を計算するmid関数では、aとbの絶対値が両方共1より大きいときはa*0.5+b*0.5で、それ以外のときは(a+b)*0.5で計算しています。このようにifで場合分けをしていいなら話は簡単なのですが、matlabの場合はベクトル単位で一括処理しないとパフォーマンスが出ないのでなるべくifを使いたくないという事情があります。if文禁止で完璧な中点の計算アルゴリズムはまだ見つかっていません。誰か考えて!

2016/12/10(土)kv-0.4.38

kvライブラリを0.4.38に更新しました。

今回は、区間演算にmidrad、max、minを追加しました。また、dka.hppの中にあるvdka (精度保証付きDKA法) で、解が重複する場合の処理をきちんと行なうようにしました。

midradとmax/minの追加は、どちらも同僚のS先生の要請によるものです。少し詳しく説明します。

midrad

midradは、区間 I に対して、
midrad(I, m, r);
とすると中心の近似mと、半径rを同時に計算します。これは、中心m、半径rの区間で I を包含する、すなわち、上端下端型の区間から中心半径型の区間への変換に使われることを意図しています。単純に
m = mid(I);
r = rad(I);
とすれば良さそうなものですが、中心mに丸め誤差が入るので、これだと I を包含しない可能性があるのです。例えば、ε=2-52として、I=[1,1+3ε]とします。このとき、midとradを別々に計算すると、
m = 1+2ε
r = 1.5ε
となりますが、[m-r,m+r]は I を含まないことが分かります。midradを使うと、ちゃんと
m = 1+2ε
r = 2ε
になります。ミスが起こりやすいポイントなので気をつけて下さい。

max, min

max, minは区間同士のmax, minです。max(I1, I2)は、x1∈I1、x2∈I2としたときのmax(x1,x2)の取り得る範囲の区間を返します。minも同様です。計算は、
max([a,b],[c,d]) = [max(a,c),max(b,d)]
min([a,b],[c,d]) = min(a,c),max(b.d)]
のように行います。

この関数は、後述するようにVisual C++と少し相性が悪いこともあって実装していなかったのですが、実装しないと思わぬバグに遭遇する可能性があると同僚のS先生に指摘され、実装することにしました。例えば、max, minを実装していないversionのkvで、
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>

typedef kv::interval<double> itv;

using namespace std;

int main()
{
    cout << max(itv(2., 4.), itv(3., 5.)) << endl;
    cout << min(itv(3., 5.), itv(2., 4.)) << endl;
}
のようなプログラムを動かすと、予想に反して
[2,4]
[3,5]
という結果が帰って来ます。これは、std::max, std::minが呼び出されてしまっており、恐らくmaxは
template <class T> 
T max(const T& a, const T& b)
{
    return (a < b) ? b : a;
}
のような実装になっていて、[2,4]<[3,5]は偽なので[2,4]が返されてしまっていると推測されます。よく考えれば分かるとは言え、これは事故を誘発しやすいので、ちゃんと区間用のmaxとminを実装しました。kv-0.4.38で上のプログラムを動かすと、ちゃんと
[3,5]
[2,4]
となります。

なお、Visual C++はwindows.hというヘッダファイルを持っており、これをincludeするとマクロでmaxとminを定義してしまうという大変極悪な仕様になっています。これがincludeされていると、max,minを定義している箇所がマクロ置換されてしまいinterval.hppがコンパイル出来ません。windows.hとkvをどうしても同時に使いたければ、
#define NOMINMAX
#include <windows.h>
のようにwindows.hをincludeする前にNOMINMAXを定義して下さい。

2016/12/03(土)double-double演算が異常に高精度になる!?

いわゆるdouble-doubleによる4倍精度演算が不思議な挙動を示した例を見つけたので、メモしておきます。

いわゆる普通の電卓で、適当な数(例えば100)を入れて、平方根のボタンを何回か押して、次に二乗(多くの電卓で[×][=]という操作)を同じ回数だけ行います。このとき、その回数がある程度以上多いと、丸め誤差でちゃんと100に戻ってきません。100円ショップに売っていた8桁のごく普通の電卓で試してみたところ、
100 → (10回平方根) → 1.0045072 → (10回二乗) → 99.9806
と、誤差が観測されました。更に回数を増やしてみると、
100 → (20回平方根) → 1.0000042 → (20回二乗) → 81.635475
100 → (25回平方根) → 1 → (25回二乗) → 1
のようになりました。平方根を取った値は徐々に1に近づき、1+εのεを保持する桁数が徐々に小さくなっていって、誤差がひどくなっていく様子がよく分かります。

もちろん、普通のPCで倍精度(double)を用いても同じことで、誤差が入ります。
#include <iostream>
#include <cmath>

int main()
{
	int i;
	double x;

	std::cout.precision(17);

	x = 100;
	
	for (i=0; i<10; i++) {
		x = sqrt(x);
	}

	for (i=0; i<10; i++) {
		x = x * x;
	}

	std::cout << x << "\n";
}
を実行すると、
100.00000000000637
のように誤差が入りました。kvライブラリを使って区間演算にしてみます。
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>

typedef kv::interval<double> itv;

int main()
{
	int i;
	itv x;

	std::cout.precision(17);

	x = 100;
	
	for (i=0; i<10; i++) {
		x = sqrt(x);
	}

	for (i=0; i<10; i++) {
		x = x * x;
	}

	std::cout << x << "\n";
}
すると、
[99.99999999997776,100.00000000004346]
のように100を含む区間が得られます。

さて、ここからが本題です。このdoubleで区間演算をしたときの区間幅を、平方根と二乗の回数を変化させながらプロットしてみます。
sqrt1.png

60回手前で計算が破綻していることが分かります。どこまで行けるかは仮数部の長さで決まる筈なので、mpfrを使って仮数部長を変化させて比較してみます。
sqrt2.png

すると、mpfrの仮数部を53bit(doubleと同じ)にした場合はdoubleとぴったり同じ、mpfrの仮数部長を長くすればそれだけ破綻までの回数が大きくなっています。ここまでは予想通り。ここで、このグラフにdd(double-double演算による擬似4倍精度、仮数部は106bit相当)を追加してみましょう。
sqrt3.png

なんだか異常なグラフが得られました。仮数部106bit相当なのでmpfr106と同じ動きをすると思いきや、最初は同じ挙動を示すものの途中から全く精度の悪化が見られず、mpfr212をも凌ぐ精度を叩きだしています。そんな馬鹿なとmpfrの超高精度を追加してみます。
sqrt4.png

すると、mpfrの仮数部1100bitで、ようやくddに勝つことができました。この現象は、
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

typedef kv::interval< kv::dd > itv;

int main()
{
	int i;
	itv x;

	std::cout.precision(17);

	x = 100;
	
	for (i=0; i<10; i++) {
		x = sqrt(x);
	}

	for (i=0; i<10; i++) {
		x = x * x;
	}

	std::cout << x << "\n";
}
のようなプログラムで区間幅を観察すれば容易に確かめられます。

さて、なぜこんなことが起きたのでしょうか。最初はバグを疑ったのですが、バグではありませんでした。後日この現象の原因を追記しようと思いますので、少し考えてみて下さい。

解答 (12月4日追記)

一日経ったので理由を簡単に説明します。

まず、ddとmpfr106で、平方根を40回行った場合の値を見比べてみます。表示は40桁にしました。
[1.000000000004188377884927590880100368969,1.000000000004188377884927590880168161704] (dd)
[1.000000000004188377884927590880118857896,1.000000000004188377884927590880168161704] (mpfr106)
ここではほぼ違いは見られません。上限と下限で一致している桁数は32桁で、4倍精度としては普通です。次に、平方根を80回にしてみます。
[1.000000000000000000000003809307495356531,1.000000000000000000000003809307495356571] (dd)
[1.000000000000000000000003809307449647879,1.000000000000000000000003809307498951686] (mpfr106)
こちらは顕著な違いが見られます。mpfrの方は32桁で変わっていませんが、ddの方は38桁も一致していることが分かります。

これは、ddの内部表現の特殊性によるものです。ddは、簡単にいうと仮数部の上位53bitと下位53bitを分割して格納するようなフォーマットです。よって、上位と下位の指数部は53ずれるのが普通なのですが(正確には下位の符号が利用できるので54ずれる)、上位と下位の間に0(上位と下位の符号が違う場合は1)が連続するような場合、それを省略することができ、上位と下位の指数部のずれが大きくなって仮数部長が大きくなることがあります。この場合も、40回のときの値の内部表現を見ると、
1.000000000004188427382700865564-4.9497773274684e-17
ですが、80回のときは
1+3.8093074953565e-24
で上位と下位が大きく離れています。たまたま収束先が1(=doubleで正確に表現可能)で、1+εのεの精密な表現力が問われるような問題だったので、ddが異常に高精度になったという仕掛けでした。

極めてレアな現象でいつもこういうことが起きるわけではないですが、偶然出会ったので記事に残しておきたくなったのでした。

2016/11/20(日)kv-0.4.37

kvライブラリを久しぶりに0.4.37にアップデートしました。

今回は、double-double (dd) 関連のアップデートです。ddのsqrtが(多分)精度が上がって速くなっています。また、sqrtに無限大を入れた時にNaNになってしまっていたバグを修正しました。

そして、重大なバグ修正を含んでいます。0.4.36までは、ddを内部に持つ区間演算 interval<dd> の除算において、ある特定の条件のときに丸めの向きを間違うバグがあり、精度保証されていなかった可能性があります。interval<dd> を使って何らかの精度保証を行っている方は、速やかに0.4.37にアップデートをお願いします。近似計算としてddを使っている場合は問題ありません。

また、ddに関してはそれなりに利用者がいるにもかかわらずきちんとした形でアルゴリズムを記載していませんでした。今回、
を書きましたので、興味のある方は是非お読み下さい。

2016/10/31(月)bash on Ubuntu on Windowsを試してみる

Windows10の夏の大型アップデート(Anniversary Update)で搭載された、Windows Subsystem for Linux (bash on Ubuntu on Windows)を使ってみたので、記録を残しておきます。

cygwinやmsys2など、Windows上でunixツールを使うためのものは以前からいろいろありますが、Microsoft本家が出してきたこいつは、「ubuntuのバイナリがそのまま動く」という点が今までと違います。使うための条件は、
  • Windows10の64bit版であること
  • Windows10のバージョン1607以降であること
です。バージョンは、「スタート→歯車アイコン→システム→バージョン情報」で確認できます。
version1607.png

8月のリリース以降、少しずつ時間をずらしながらWindows Updateを降らせていたようですが、そろそろほとんどのWindows10が1607になった頃ではないでしょうか。

インストール

インストール方法の情報は検索すればたくさん出てきますが、次のような手順です。
  • 「スタートを右クリック→プログラムと機能→Windowsの機能の有効化または無効化」で、「Windows Subsysyem for Linux (Beta)」をチェックする。再起動。
    wsl.png

  • 「スタート→歯車アイコン→更新とセキュリティ→開発者向け」で、「開発者モードを使う」を、「サイドロードアプリ」から「開発者モード」に変更。再起動。
    developermode.png

  • 「スタートを右クリック→コマンドプロンプト」でコマンドプロンプトを起動し、bashとタイプ。"y"でダウンロードとインストールが始まります。数分かかります。ユーザIDとパスワードを聞かれてインストール完了。
    install.png

  • 次回以降は、スタートメニューに「Bash on Ubuntu on Windows」が登録されているのでそこから起動できます。
    startmenu.png

アンインストール

いろいろ試しにパッケージを入れたりしてシステムが壊れることもあるかと思いますが、コマンドプロンプトで
lxrun /uninstall /full
と入れるときれいさっぱり削除できます。上の「bashとタイプ」のところからやりなおすことが出来ます。
uninstall.png

簡単な使い方など

  • 中身はubuntu 14.04です。いつものubuntuの作法通り、
    sudo apt update
    sudo apt upgrade
    
    で最新に更新しておきましょう。
  • これで、端末内で完結するような作業は大体できます。最初は最低限しかインストールされていないので、いつものubuntuの作法で必要なパッケージをインストールしましょう。とりあえず
    sudo apt install build-essential
    
    でCコンパイラなど最低限の開発環境を入れることをお勧めします。
  • ubuntu側からは、windows側のファイルが例えばCドライブなら
    /mnt/c/
    
    以下に見えます。逆にwindows側から見てubuntuのファイルシステムは
    c:\Users\(windowsユーザ名)\AppData\Local\lxss\ 
    
    にあり、このフォルダがubuntu側の"/"に対応しています。隠しファイルになっているので普通の状態では見えないかも知れません。恐らく、ubuntu側からwindowsのファイルシステムを操作することはOK、windows側からubuntuのファイルシステムを操作するのはNG、だと思われます。
プログラミングして楽しむだけならこれで十分と思われます。が、次のような問題点があり、生活の全てをこの環境で行なうのは難しそうです。
  • X windowを使うソフトウェアが動かない。gnuplotくらい使いたい!
  • (文字幅を正しく扱えないせいか)日本語が頻繁に文字化けする。
  • そもそも日本語が入力できない。

X環境を作る

買ったままのwindowsでソフトウェア追加無しにunix環境が使えるのがbash on ubuntu on windowsのメリットですが、windows用のX serverを入れてしまえば使えるソフトウェアが一気に増え、また日本語の問題も解決できる可能性があります。そこで、上記問題点を解決すべく、X serverを入れてみます。

X serverは有料、無料いろいろあると思いますが、いろいろ検索してみると無料では
の2つがよく使われているようです。今回はXmingの方を使ってみました。Xming X Serverによると、Public Domain Releasesという無料版と、Website Releaseという新しいが寄付が必要な版があるようです。今回はPublic Domain Releasesの方を使いました。
  • Xming 6.9.0.31
  • Xming-fonts 7.7.0.10
をダウンロードしてインストールしました。途中sshクライアントを入れるか聞かれますが、自分は不要だったので「Don't install SSH client」としました。

スタートメニューから「Xming」を起動します。(「XLaunch」の方だといろいろオプションを設定してから起動します。今回はデフォルトのままで十分。) すると、右下に
xming.png

のようなアイコンが出ます。これで、bash on ubuntu on windowsからの描画命令を受け止める準備が出来たことになります。動作チェックします。
sudo apt install x11-apps
sudo apt install x11-utils
sudo apt install x11-xserver-utils
といくつかのx11の基本アプリをインストールして、
DISPLAY=localhost:0.0 xeyes &
とxeyesを起動してみます。
xeyes.png

のようにマウスカーソルを追う目玉が表示されたら、正常に動作しています。

gnuplotを試してみましょう。
sudo apt install gnuplot-x11
としてインストールし、
DISPLAY=localhost:0.0 gnuplot
として起動します。
gnuplot.png

のように、ちゃんと動作しました。毎回「DISPLAY=localhost:0.0」とするのが面倒なら、
export DISPLAY=localhost:0.0
とすると、端末を閉じるまで有効になります。

驚くべきことに、日本語フォントを追加してやると、少しエラーが出るもののfirefoxを動かすことができます。
sudo apt install fonts-ipafont
sudo apt install firefox
firefox.png

日本語環境を整える

ここまでちゃんと動作するとなると、日本語の読み書きがまともにできないのが惜しくなってきます。そこで、少し頑張って環境を整えてみました。いろいろ試行錯誤した結果ではありますが、もっといい方法もありそうなので情報が欲しいところです。以下、自分が試した方法を書きます。X serverがwindows側にインストールしてあって、また上で書いたように、
sudo apt install fonts-ipafont
で日本語フォントを追加してあるものとします。

まず、端末は、windows側を捨ててXの方で動かすことにします。いろいろ試しましたが、lxterminalが良さそうでした。
sudo apt install lxterminal
でインストールし、
DISPLAY=localhost:0.0 lxterminal &
で起動します。起動後に「編集→設定」でフォントをMonoSpace 10からMonospace 15くらいにしてあげると見やすい感じになりました。
lxterminal2.png

windows側の端末を使わずにこちらを使うことにします。こちらからだと「DISPLAY=locahost:0.0」をいちいち打たなくてよくなります。windows側の端末を閉じると全部落ちてしまうので、アイコンにでもしておきましょう。

かな漢字変換のシステムを入れます。いろいろ試しましたが、uim-anthyが何とか動作しました。
sudo apt install uim uim-xim uim-anthy
のようにインストールします。そして、windows側のbashターミナルで、
DISPLAY=localhost:0.0 UIM_CANDWIN_PROG=uim-candwin-gtk uim-xim &
のようにかな漢字変換サーバを起動し、lxterminalの起動は
DISPLAY=localhost:0.0 XMODIFIERS="@im=uim" GTK_IM_MODULE=uim QT_IM_MODULE=uim lxterminal &
とします。これで、「半角/全角」キーで日本語入力ができるようになりました。
uimanthy.png

terminal起動時の設定が長いですが、このterminalから起動したものにはこの設定が伝わるので、ここからいろいろ起動することにすれば楽です。

その他もろもろ

他にもいろいろ入れてみました。

TeX環境。
sudo apt install texlive-lang-cjk
で簡単に入ります(ちょっと時間がかかります)。ついでに
sudo apt install evince
でpdfビューアも。
tex.png

java。
sudo apt install default-jdk
あちこちでjavaは動かないという記述を見かけましたが、普通に動いているように見えます。

自分はvimで十分ですが、もう少し普通のエディタを使いたいなら、
sudo apt install geany
あたりはいかがでしょうか。
geany.png

ま、windows側のお気に入りのエディタを使えば済むことではありますが。

vmwareなどの仮想化ソフトを使うよりずっと軽いのが嬉しいです。windowsとの分業がしやすいのも大きな利点かと。次はsshdなどサーバ系のソフトをいろいろ試してみたいと思います。

追記

上で書いたかな漢字変換サーバとlxterminalの起動を自動化するなら、例えばhome directoryの.bashrcの末尾に
if [ $SHLVL -eq 1 ]; then
  if DISPLAY=localhost:0.0 xset q > /dev/null 2>&1 ; then
    DISPLAY=localhost:0.0 UIM_CANDWIN_PROG=uim-candwin-gtk uim-xim &
    DISPLAY=localhost:0.0 XMODIFIERS="@im=uim" GTK_IM_MODULE=uim QT_IM_MODULE=uim lxterminal &
  fi
fi
のように書けばいいでしょう。最初に起動されたbashで、なおかつXが利用可能なら、かな漢字変換サーバとlxterminalを起動します。

次のwindows10の大型アップデートでubuntu 16.04になるとか日本語入力も普通に出来るようになるとか噂が聞こえてくるので、そのときにはここに書いたことの大半は無意味になってしまうかもしれません。

この記事は、次のページを参考に書きました。
貴重な情報を公開して下さった皆様に感謝します。
OK キャンセル 確認 その他