2026/05/30(土)kv-0.4.60

kvライブラリを0.4.60にアップデートしました。

今回は、g++-13での仕様変更により0.4.56からg++-13以上で利用できなくなっていた、Intel 80bit浮動小数点数のサポートを復活させました。

以前は、Intel FPUに盲腸のように存在してる80bit浮動小数点数を_Float64xという型で記述して活用していましたが、g++-13から、
#include <iostream>

int main()
{
    _Float64x x;
    std::cin >> x;
}
みたいな簡単なプログラムさえコンパイルできなくなっていました。仕方ないのでg++-13以上ではこの機能を抹消していたのですが、long doubleを使えば問題なく動いていたので、今回ソースコード中の_Float64xを全部書き直しました。また将来状況が変わって振り回されるのは嫌なので、kv::fp80というaliasを作って全部それで書くことにして、状況が変わったらすぐにまた変更できるように備えました。

80bit(仮数部64bit)の浮動小数点数を使って区間演算をしたり、2つ束ねて仮数部128bitの疑似拡張4倍精度を使ったり、またその疑似拡張4倍精度を端点に持つ区間演算をしたり、結構面白いことができます。もちろんmpfrを使えばもっと高精度にすることもできますが、mpfrより速度低下がずっと小さいです。doubleやddじゃ微妙に精度が足りないが、mpfrでは遅すぎるというギリギリのケースで役に立ちことがあるかもしれません。

2026/03/16(月)kv-0.4.59

ものすごく久しぶりに、kvを0.4.59にアップデートしました。なんと前回からぴったり一年後でした。
いくつかバグが溜まってきて、そろそろ放出したくなりました。
  • double/ddで、DBL_MAX/3がNaNになってしまう問題の修正。
  • eig.hppのeig(近似固有値計算)で、稀に無限ループに陥ってしまうバグを修正 (Thanks to Codex CLI)。
  • 区間演算のsin,cosで、極値を含む判定を厳密に (Thanks to 岩波さん)。
  • exp(interval<dd>(1000.))がなぜか[-inf,inf]になってしまうバグの修正。
などです。

eigはずっと直したかったのですが、AI様に聞いたら一発で直してくれました。精度保証のプログラムをAIに書かせるのは怖いですが、近似計算の部分ならいいでしょう。

最後のバグは、dd区間演算でオーバーフローして端点が無限大になるようなケースでいろいろ不具合が見つかりました。区間が広がってしまうバグなので嘘は付いていないでしょうけど、dd区間演算でなぜか極端に区間幅が大きくなって断念してしまったようなことがあれば、もしかしたらこのせいかもしれません。

学生たちの頑張りで、affine arithmeticは大きく改善できることが分かり書き換えたいですが、大改造になりそうなのでそれは次回にします。

また、githubで、今までversion情報はcommit messageにしか入れていなかったのですが、"v0.4.59"みたいなtagも付加するようにしました。githubは何も分かっていないですがいくらかまともになった?

2025/04/06(日)非正規化数の計算は遅い?

はじめに

よく知られているように、IEEE754に従う浮動小数点数はフォーマットの一部に非正規化数(denormalized number, denormal number, subnormal number)という領域があります。倍精度では、絶対値が2-1074から2-1022までの0に近い領域で、非常に小さな数が表現できるのと引き換えに、この範囲では仮数部の長さが通常の53bitより短く精度が低下します。

ところで、我々の業界では、計算の途中に非正規化数が出てくると計算速度が極端に低下(50~100倍ほど)すると言われていて、いくつかの精度保証アルゴリズムは非正規化数がなるべく出現しないように設計されています。しかし、実際にそんなに遅くなるのかどうか、実測したことは無かったので、計測してみました。

使用プログラム

使ったのは次のようなプログラムです。加算、乗算、除算 (減算は加算と同じだろうから省略) について、引数に非正規化数が含まれている場合と含まれていない場合の計算速度の比を計測しています。10億回繰り返しています。volatileに代入したり、計算結果を表示したりして、最適化で計算本体が消されないように小細工しています。実際にやっている計算は
  • 正規化数 + 非正規化数 = 正規化数
  • 非正規化数 * 正規化数 = 非正規化数
  • 非正規化数 / 正規化数 = 非正規化数
です。もしかしたら計算結果が正規化数か否かでも変わるかもしれませんが、そこまでは調べていません。
#include <iostream>
#include <cmath>
#include <chrono>

double vc(double x)
{
	volatile double tmp = x;
	return tmp;
}

int main()
{
	std::chrono::system_clock::time_point ts0, ts1, ts2, ts3;
	int i, j;
        double a, b;
        double t0, t1, r0, r1, r2;
	

	a = vc(1);
	b = vc(std::pow(2., -1021));

	std::cout << "add normal\n";
        ts0 = std::chrono::system_clock::now();
	for (i=0; i<10000; i++) {
		for (j=0; j<100000; j++) {
			a += b;
		}
	}
        ts1 = std::chrono::system_clock::now();
        t0 = std::chrono::duration_cast<std::chrono::nanoseconds>(ts1-ts0).count()/1e9;
        std::cout << t0 << "\n";

	std::cout << a << "\r              \n"; // dummy

	a = vc(1);
	b = vc(std::pow(2., -1023));

	std::cout << "add denormal\n";
        ts2 = std::chrono::system_clock::now();
	for (i=0; i<10000; i++) {
		for (j=0; j<100000; j++) {
			a += b;
		}
	}
        ts3 = std::chrono::system_clock::now();

        t1 = std::chrono::duration_cast<std::chrono::nanoseconds>(ts3-ts2).count()/1e9;
        std::cout << t1 << "\n";

	r0 = t1 / t0;
	std::cout << "\nadd ratio\n";
        std::cout << "denormal / normal = " << r0 << "\n";

	std::cout << a << "\r              \n"; // dummy

	a = vc(std::pow(2., -1021));
	b = vc(1);

	std::cout << "mul normal\n";
        ts0 = std::chrono::system_clock::now();
	for (i=0; i<10000; i++) {
		for (j=0; j<100000; j++) {
			a *= b;
		}
	}
        ts1 = std::chrono::system_clock::now();
        t0 = std::chrono::duration_cast<std::chrono::nanoseconds>(ts1-ts0).count()/1e9;
        std::cout << t0 << "\n";

	std::cout << a << "\r              \n"; // dummy

	a = vc(std::pow(2., -1023));
	b = vc(1);

	std::cout << "mul denormal\n";
        ts2 = std::chrono::system_clock::now();
	for (i=0; i<10000; i++) {
		for (j=0; j<100000; j++) {
			a *= b;
		}
	}
        ts3 = std::chrono::system_clock::now();

        t1 = std::chrono::duration_cast<std::chrono::nanoseconds>(ts3-ts2).count()/1e9;
        std::cout << t1 << "\n";

	r1 = t1 / t0;
	std::cout << "\nmul ratio\n";
        std::cout << "denormal / normal = " << r1 << "\n";

	std::cout << a << "\r              \n"; // dummy

	a = vc(std::pow(2., -1021));
	b = vc(1);

	std::cout << "div normal\n";
        ts0 = std::chrono::system_clock::now();
	for (i=0; i<10000; i++) {
		for (j=0; j<100000; j++) {
			a /= b;
		}
	}
        ts1 = std::chrono::system_clock::now();
        t0 = std::chrono::duration_cast<std::chrono::nanoseconds>(ts1-ts0).count()/1e9;
        std::cout << t0 << "\n";

	std::cout << a << "\r              \n"; // dummy

	a = vc(std::pow(2., -1023));
	b = vc(1);

	std::cout << "div denormal\n";
        ts2 = std::chrono::system_clock::now();
	for (i=0; i<10000; i++) {
		for (j=0; j<100000; j++) {
			a /= b;
		}
	}
        ts3 = std::chrono::system_clock::now();

        t1 = std::chrono::duration_cast<std::chrono::nanoseconds>(ts3-ts2).count()/1e9;
        std::cout << t1 << "\n";

	r2 = t1 / t0;
	std::cout << "\ndiv ratio\n";
        std::cout << "denormal / normal = " << r2 << "\n";

	std::cout << a << "\r              \n"; // dummy

	std::cout << "ratio (add/mul/div): " << r0 << " " << r1 << " " << r2 << "\n";
}

実験結果

これをg++で-O3をつけてコンパイル、実行して、正規化数の場合の計算時間を1としたときの非正規化数の計算時間を、手元に転がっていたPCで片っ端から調べてみました。次の表にまとめます。
CPU加算乗算除算
Intel Core M-5Y10c (Broadwell)0.99933842.199514.3503
Intel Core i5-6500T (Skylake)0.92403433.542811.956
Intel Core i9-7900X (Skylake)0.96088633.876411.9609
Intel Pentium CPU 4415Y (Kaby Lake)0.99709934.407912.1162
Intel Core i5-8250U (Kaby Lake R)0.99220934.503712.1195
11th Gen Intel Core i7-1195G70.9998429.120210.577
Intel Celeron N4120 (Gemini Lake Refresh)57.132643.77513.7602
Intel N100 (Alder Lake-N)95.714171.82423.5843
AMD Ryzen 7 7700 (Zen4)0.9816351.553461.39528
AMD Ryzen 7 8840U (Zen4)0.990941.329991.40047
Apple M10.9648730.9992541.00117
これを見ると、IntelのCPUは加算は遅くならないが乗除算が遅い、ATOM系は加算も遅い、AMDやApple Siliconはほとんど遅くならない、ことが読み取れます。

おわりに

個人で所有しているCPUの種類には限界があってなかなか網羅的な調査は難しいですが、いろいろ検索して情報を集めてみた結果、どうやら次の表のような感じっぽいです。(o=ペナルティなし、x=ペナルティあり)
CPU加算乗算除算
Intel SandyBridgeより古いxxx
AMD Bulldozerxxx
Intel SandyBridge以降11世代までoxx
Intel ATOM系xxx
Intel KNLoox
AMD Zenooo
Apple Siliconooo
SandyBridgeより古いのとかBulldozerとかKNLとか、実機で試すのは大変ですが、いつかやってみたいものです。あるいはどなたか動かしてくれないかしら。

Intel 12世代以降はどうなってるの?

ところで、実験を見て気になるのはN100です。加算も遅いという散々な結果ですが、N100って、Intel 12世代以降のいわゆるEコアで出来ているはず。ということは、12世代以降のCPUで運悪くEコアに当たってしまうと加算ですら非正規化数のペナルティが発生する? 自分はIntelの12世代以降は所有していないので試せないのですが、気になります。

追記 (2025年4月7日)

未開封のCore Ultra 7 155UのPCがあったのに気づいたのでセットアップ試してみました。こいつはCore Ultra シリーズ1というやつで、Pコア、Eコア、LP(低電力)Eコアがそれぞれ2,8,2コアあり、PコアはHyper Threadingで倍になるので全部で14threadというなかなかややこしい構成になっています。/proc/cpuinfoを読み出してじっと睨んでると14の論理コアごとに微妙に違っていて、
processorcore idcache size
0812288 KB
1812288 KB
21212288 KB
31212288 KB
4012288 KB
5112288 KB
6212288 KB
7312288 KB
8412288 KB
9512288 KB
10612288 KB
11712288 KB
12322048 KB
13332048 KB
0-3がPコア、4-11がEコア、12-13がLP Eコアと推定できます。

これを使って、
$ taskset 0x00000001 ./a.out
$ taskset 0x00000010 ./a.out
$ taskset 0x00001000 ./a.out
のようにbitmaskで使用していいprocessorを指定して実行しました。結果は、
CPU加算乗算除算
Core Ultra 7 155U (P core)0.96551641.671411.1905
Core Ultra 7 155U (E core)94.278271.638723.4713
Core Ultra 7 155U (LP E core)95.088671.63223.4997
の通りでした。Eコアは全演算にペナルティがあって、Pコアは乗除算にペナルティがあります。

なお、るふぁ先生が調査して下さった結果を合わせると、Intelの12世代以降は、
CPU加算乗算除算
Intel 12-14世代 (Pコア)oxx
Intel 12-14世代 (Eコア)xxx
Intel Core Ultra 第1シリーズ (Pコア)oxx
Intel Core Ultra 第1シリーズ (Eコア)xxx
Intel Core Ultra 第2シリーズ (Pコア)oxx
Intel Core Ultra 第2シリーズ (Eコア)oxx
と、Core Ultra 第2シリーズでEコアが良くなったっぽいです。

2025/03/18(火)数学ソフトウェアとフリードキュメント

昨日は、数学ソフトウェアとフリードキュメント 37という集まりで、kvライブラリについて50分間、話をしてきました。数値計算の誤差がどのくらい怖いか、という「脅し」のパートが受けてたみたいです。kvライブラリのホームページに置いてあるスライドにちょっと追加しただけではありますが、スライドが欲しいという話があったのでここに置いておきます。

2025/03/15(土)kv-0.4.58

kvライブラリを0.4.58にアップデートしました。

今回も大した変更はありません。
  • ODE Solverで、精度保証に失敗した場合ごく稀にstep sizeを小さくしてのretryをせずにエラーで止まってしまうことがあったのを修正。
  • PSAのpow(x, y)で、yが整数でない定数だがpsa型でない場合の動作を改善。
  • vleq.hppの精度を改善 (高安先生の要請による)
くらいです。

数学ソフトウェアとフリードキュメント 37でkvについてしゃべるので、いい機会なのでここまで溜まっていたアップデートを吐き出したかったという理由もありました。
OK キャンセル 確認 その他