2024/08/27(火)kv-0.4.57

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

今回は、
  • 以前にるふぁさんにご指摘頂いた、defint_autostepで使用する級数の次数を1にしたときに正しく計算できていなかったバグの修正。
  • ode_maffineに自動微分型を指定し、なおかつcallback関数を指定したとき、 callback関数が初期値に関する微分の情報をも渡してくれるように変更。
の2点です。後者はとある研究で必要になって機能追加しました。

大したアップデートではないですが、一応2つ溜まったので。

2024/02/07(水)g++-13と_Float64xとkv-0.4.56

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

今回は問題発生への対処です。ある人から、g++ version 13でkvがまったくコンパイルできない、clang++では問題ないと連絡が来ました。実際、ubuntu 23.10を入れて試してみると、test-rounding.ccとかtest-interval.ccみたいな基本的なプログラムすらコンパイルエラーになってしまいました。

調べてみると、_Float64x (IntelのFPUが持っている80bit拡張浮動小数点数) の関係する部分がエラーになっています。kvどころか、
#include <iostream>

int main()
{
        _Float64x a;

        std::cin >> a;
}
みたいな単純なプログラムがコンパイルエラーになってしまいます。また、
#include <iostream>
#include <limits>

int main()
{
        std::cout << std::numeric_limits<_Float64x>::epsilon() << std::endl;
        std::cout << std::numeric_limits<_Float64x>::max() << std::endl;
        std::cout << std::numeric_limits<_Float64x>::min() << std::endl;
        std::cout << std::numeric_limits<_Float64x>::infinity() << std::endl;

        std::cout << std::numeric_limits<__float80>::epsilon() << std::endl;
        std::cout << std::numeric_limits<__float80>::max() << std::endl;
        std::cout << std::numeric_limits<__float80>::min() << std::endl;
        std::cout << std::numeric_limits<__float80>::infinity() << std::endl;

        std::cout << std::numeric_limits<long double>::epsilon() << std::endl;
        std::cout << std::numeric_limits<long double>::max() << std::endl;
        std::cout << std::numeric_limits<long double>::min() << std::endl;
        std::cout << std::numeric_limits<long double>::infinity() << std::endl;
}
は、
0
0
0
0
1.0842e-19
1.18973e+4932
3.3621e-4932
inf
1.0842e-19
1.18973e+4932
3.3621e-4932
inf
のように、numeric_limits<_Float64x>はすべて0を返してきます。

調べてみると、「_Float64xをサポートしているのはgccのみで、もともとg++では_Float64xをまったくサポートしていなかった。しかしversion 12までは単なるlong doubleのtypedefだったので、たまたまg++でも動いてしまっていた。最新のg++ version 13では_Float64xを部分的にサポートするようになったがlibstdc++はサポートしていない状態になり、結果的にまったく使えなくなってしまった」ということのようです。やりとりはこの辺

gccのドキュメントの6.12 Additional Floating Typesでは、「As an extension, GNU C and GNU C++ support additional floating types, (中略) _float80 is available on the i386, x86_64, and IA-64 targets, and supports the 80-bit (XFmode) floating type. It is an alias for the type name _Float64x on these targets. 」とはっきり書いているので、g++でも_Float64xはサポートされているように見えるのですが。

IntelのCPUではlong double、__float80も_Float64xと同じ80bit拡張浮動小数点数です。こっちの2つはg++-13でも生き残っています。どちらかで全部書き換えることも考えたのですが、long doubleは異なるアーキテクチャだと異なる浮動小数点数フォーマットに対応していることが多く、トラブルを招く可能性が高い、__float80はgccのみのオリジナル拡張でclangで使えない、とそれぞれ問題があります。

もっともメジャーで、4月にはみんながubuntu 24.04を使い始めることを考えると、g++-13で動かないなんてのは使い物にならないも同然なので、緊急にkv-0.4.56をリリースしました。とりあえず最低限の変更で、g++-13では_Float64x関連の機能が全部働かないようにしました。g++-12以下ではすべて問題なく、g++-13では_Float64xを使う機能を除いてすべてコンパイルできるようになりました。他の変更は、ついでにoptimize.hppにちょっとした追加機能が入ったくらいです。

g++-14以降でちゃんと戻ることを祈るのみですが、こんなIntel CPUにしかない盲腸のような機能を使おうとするほうが悪い、という話でもあるのかな。

2022/09/16(金)kv-0.4.55

久しぶりに、kvを0.4.55にアップデートしました。

今年の1月に、M1 macの区間演算は遅い?という記事を書きました。このときに、ARM CPUでのinline assemblerによる丸めの向きの変更のコードを追加していたのですが、大して面白い変更では無かったので放置していました。気づいたら8ヶ月も経ってしまい、細かい修正も溜まってきたのでここでいったん公開することにしました。

というわけで、今回の変更はARM CPUで-DKV_FASTROUNDを付けるとinline assemblerによる丸めの向きの変更が使えるようになる、というものです。ARM 64bitだけでなく、一応ARM 32bitでも動くようにしたつもりですが、あまりテストされていません。

詳細は丸めモードの変え方とコンパイルオプションまとめの「6. ベンチマーク」に書きましたが、一般的にARM CPUでは丸めのモードを変更すると実行時間で大きなペナルティがあるようで、Intel CPUに比べてかなり遅いです。ARM CPUではハードウェアによる丸め変更を一切行わずにソフトウェアで方向付き丸めをエミュレートするのが一番速いという、残念な状況になっています。これを、精度保証に向かないCPUが普及しつつあって残念と見るか、エミュレーションのアルゴリズムを開発しておいて良かった、と見るか。

以下、とある区間演算を多用したプログラムを、端点がdouble精度の区間演算と、端点がdd精度の区間演算の場合について、オプションを変えながら実行したときの実行時間を上のページから抜き出しておきます。爆速のはずのM1 chipの計算速度が遅く、-DNOHWROUND (ソフトウェアエミュレートによる丸め変更) が一番まともになっているのが分かります。

core i5 11400

計算精度(端点の型) コンパイルオプション
(-O3 -DNDEBUGに加えて)
計算時間
-DKV_USE_TPFMAなし -DKV_USE_TPFMAあり
double なし 6.65 sec
-DKV_FASTROUND 3.74 sec
-DKV_NOHWROUND 7.88 sec 6.05 sec
-DKV_USE_AVX512 -mavx512f 2.74 sec
dd なし 37.20 sec 33.58 sec
-DKV_FASTROUND 25.08 sec 21.52 sec
-DKV_NOHWROUND 92.24 sec 71.81 sec
-DKV_USE_AVX512 -mavx512f 16.47 sec

M1 MacBook Air

計算精度(端点の型) コンパイルオプション
(-O3 -DNDEBUGに加えて)
計算時間
-DKV_USE_TPFMAなし -DKV_USE_TPFMAあり
double なし 24.45 sec
-DKV_FASTROUND 22.36 sec
-DKV_NOHWROUND 9.08 sec 6.03 sec
dd なし 108.8 sec 100.6 sec
-DKV_FASTROUND 102.4 sec 94.64 sec
-DKV_NOHWROUND 84.60 sec 53.76 sec

2022/01/13(木)M1 macの区間演算は遅い?

昨年3月に手に入れていたM1チップ搭載のMacBook Air、今頃開封して設定して遊んでいました。OSをMontereyにして、arm64なterminalでhomebrewをインストールしてarm64版のhomebrewを入れて、brew install gcc, brew install boostしたくらいです。とりあえずkvのプログラムは、
g++11 -O3 -I(kvのあるとこ) -I/opt/homebrew/include -L/opt/homebrew/lib hogehoge.c
とかすれば動くようになりました。

で、なぜか区間演算が異常に遅いことに気付きました。とりあえず区間演算のベンチマークとして、
#include <iostream>
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>

int main()
{
        int i;
        kv::interval<double> x;

        x = 1;

        for (i=0; i<1000000000; i++) {
                x = (x*x-3)/(x+x);
        }

        std::cout << x << std::endl;
}
みたいなのを使います。加減乗除を一回ずつ含んでいて、[1,1]と[-1,-1]を交互に繰り返します。本当は変なベンチマークテストみたいにいろんな値を取る写像を使いたかったのですが、区間演算だとすぐに発散して[-∞, ∞]になってしまい、ベンチマークとしてふさわしくないものになってしまいました。

これの実行時間を、ThinkPad X1 carbon (core i7 10510U, ubuntu 20.04, gcc 9.3)で計測すると、-O3で38.1秒程度でした。この状態では、丸めの向きの変更はfesetroundを使っています。kvライブラリには、Intelの64bit CPUのとき-DKV_FASTROUNDを付けるとfesetroundを使わずにSSE2の丸めの向きを制御するmxcsrレジスタに直接書き込むことによって高速化する機能があります。このときは、19.4秒と倍近くに高速化しました。更に、-DKV_NOHWROUNDを付けるとハードウェアによる丸めの向きの変更は一切行わず、twosumとtwoproductを用いて方向付き丸めのエミュレートを行わせることができます。このときは、35.7秒程度でした。

で、これをMacBook Air (M1, Monterey 12.1, gcc 11.2) で実行してみると、-O3でなんと142秒もかかってしまいました。fesetroundが重いのかと思い、aarch64なCPUに対しても-DKV_FASTROUNDが使えるようにkvを改造し(次期kv-0.4.55で入れる予定)実行してみると、114秒といくらか改善したものの、まだまだ遅いです。で、-DKV_NOHWROUNDだと、なんと23.8秒と劇的に改善しました。方向付き丸めのエミュレーションはかなり複雑な計算をしていて、こんなに速いはずはないのですが。まとめると、次のような感じです。
-O3-O3 -DKV_FASTROUND-O3 -DKV_NOHWROUUND
X1 carbon38.119.435.7
M1 MBA14211423.8
両PCの速度の比は、3つ目の35.7vs23.8が一番体感と一致します。丸めの変更でなにか大きなブレーキがかかっているのでは、と疑わざるを得ないような結果です。

そこで、次のような実験をしてみました。変なベンチマークテストのベンチマークに毎回丸めの向きの変更を挿入してみます。fesetroundの実装に影響されないように、アセンブリ埋め込みの最速と思われる実装にします。
#include <stdio.h>
#include <stdint.h>

int main()
{
	int i;
	double x;

#ifdef __aarch64__
	uint64_t reg;
	uint64_t regs[4];

	asm volatile ("mrs %0, fpcr" : "=r" (reg));
	regs[0] = (reg & ~(3ULL << 22)) | (0ULL << 22); // nearest
	regs[1] = (reg & ~(3ULL << 22)) | (2ULL << 22); // down
	regs[2] = (reg & ~(3ULL << 22)) | (1ULL << 22); // up
	regs[3] = (reg & ~(3ULL << 22)) | (3ULL << 22); // chop
#endif

#ifdef __x86_64__
	uint32_t reg;
	uint32_t regs[4];

	asm volatile ("stmxcsr %0" : "=m" (reg));
	regs[0] = (reg & ~(3UL << 13)) | (0UL << 13); // nearest
	regs[1] = (reg & ~(3UL << 13)) | (1UL << 13); // down
	regs[2] = (reg & ~(3UL << 13)) | (2UL << 13); // up
	regs[3] = (reg & ~(3UL << 13)) | (3UL << 13); // chop
#endif


	x = 0.5;
	for (i=0; i< 1000000000; i++) {
#ifdef __aarch64__
#ifdef ROUND_CHANGE1
		asm volatile ("msr fpcr, %0" : : "r" (regs[0]));
#endif
#ifdef ROUND_CHANGE2
		asm volatile ("msr fpcr, %0" : : "r" (regs[i%4]));
#endif
#endif

#ifdef __x86_64__
#ifdef ROUND_CHANGE1
		asm volatile ("ldmxcsr %0" : : "m" (regs[0]));
#endif
#ifdef ROUND_CHANGE2
		asm volatile ("ldmxcsr %0" : : "m" (regs[i%4]));
#endif
#endif
		x = 1 / (x * (x - 1)) + 4.6;
	}
	printf("%g\n", x);
}
埋め込んだ丸め変更命令に2種類あって、
  • -DROUND_CHANGE1 丸め変更命令を埋め込むが、丸めの向きはずっとnearest
  • -DROUND_CHANGE2 丸めの向きをnearest->down->up->chopと周期的に変更
としてみました。これの実行時間は、次のようになりました。
-O3-O3 -DROUND_CHANGE1-O3 -DROUND_CHANGE2
X1 carbon6.436.446.53
M1 MBA6.296.2915.04
X1 carbonだと、ほとんど速度は変わりません。M1だと、丸め変更の命令があっても丸めの向きがずっと一定ならば影響は認められませんが、丸めの向きが毎回変更される場合に顕著な速度低下が認められます。何らかのwaitが入ってしまっている感じです。ここまで遅いと、M1チップでの精度保証の方法を考え直さないといけないかも。線形計算のように丸めの向きの変更が頻繁でないケースは問題ないでしょうけど、非線形計算のように頻繁な丸めの向きの変更が必要な場合は、エミュレーションに頼る他ないかもしれません。

なお、これはaarch64なアーキテクチャ全般に見られる傾向かどうかが気になったので、3万円のchromebook (Lenovo IdeaPad Duet Chromebook, MediaTek Helio P60T) で試してみたら、
-O3-O3 -DROUND_CHANGE1-O3 -DROUND_CHANGE2
IdeaPad Duet15.820.720.9
とまた違った傾向を示しました。誰か、富岳で試してくれないかなあ。

ARMでのベンチマーク追加 (2022/6/30追記)

でのベンチマークを追加してみました。
-O3-O3 -DROUND_CHANGE1-O3 -DROUND_CHANGE2
X1 carbon6.436.446.53
M1 MBA6.296.2915.04
IdeaPad Duet15.820.720.9
Ampere Altra6.696.699.70
Gemini PDA19.723.028.8
Ampere AltraではM1 MBAほど遅くはなっていないのが分かります。

なお、とある方から富岳のA64FXでのベンチマークを聞かせてもらったのですが、少し異常な感じの結果(とても遅い)だったので何か測定ミスの可能性もあり、ここへの掲載はとりあえず控えます。

2022/01/06(木)kv-0.4.54

kvライブラリを0.4.54にアップデートしました。今回の変更は主に3点です。

Intel 80bit浮動小数点演算器の活用

1つ目は、Intelの80bit浮動小数点演算器の活用です。IntelのCPUの浮動小数点演算は32bit時代までは主にFPUで、64bit時代からはSSE2で行われ、現在ではハードウェアとしてのFPUは使われずに眠っている、盲腸のような存在です。FPUは、IEEE754よりも高精度な、全長80bit、指数部15bit、仮数部64bitの浮動小数点演算器を持っていて、それゆえIEEE754に完全に従わせるのが難しく、いろいろトラブルも起こしてきました。せっかく眠っているハードウェアを、逆にちょっと高精度な演算器として活用しようというのが、今回の話題です。

FPUの80bit演算は、長らくlong doubleという名前で使われることが多く、gcc/clangでは今でもその名前で使用することができますが、MSVCでは64bitのただのdoubleと同じ、他のアーキテクチャでは128bit floatだったりと、混乱を極めています。kvでは、最近規格化された、_Float64xという名前で80bit演算を行うことにしました。中田真秀先生もこの名前を推奨しています。

kvのintervalは、元々端点の型を自由にすげ替えられるように設計されていて、double, dd, mpfrを使うことができました。kv/rfloat64x.hppというファイルを作り、interval.hppの後にこれをincludeすることによって、interval<_Float64x>で区間演算を可能にしました。

次に、ちょうどdoubleを2つくっつけてdd型を作ったように、_Float64xを2つくっつけて、全長160bit、仮数部15bit、指数部128bit相当の、ddxという型を作りました。ddx.hppをincludeすると使うことができます。

更に、rddx.hppというファイルを作って、これをinterval.hppの後にincludeすることによって、ddx型を端点に持つような区間演算も行えるようにしました。

従来は、ddの指数部11bit、仮数部106bit相当で足りなければmpfrに頼るしか無かったのですが、ddxを使えば少しだけ粘ることができます。計算速度はddの1~2割増程度なので、まあまあ使えるかと思います。

i386 CPUの扱いの変更

今更サポートする意味は薄いですが、Intel CPUの32bit環境の扱いを変更しました。kvの機能のうち
  • double-double演算 (dd)
  • -DKV_NOHWROUNDを付けたときの、CPUによる丸めの向きの変更を一切使わない丸めエミュレーション
の2つは、twosumやtwoproductを使っていて、これらは演算がIEEE754に完全に従っていることを前提にしています。Intelの32bit環境はFPUの過剰精度のせいで完全にIEEE754に従っておらず、twosumやtwoproductが壊れてしまいます。これを防ぐため、従来はFPUの計算精度を53bitに制限することによって誤魔化していました。しかし、これでも完全にIEEE754に従うかは怪しく、また何よりもそれをやると今回の新機能である_Float64xを使った機能が一切使えなくなってしまいます。

そこで、float.hもしくはcfloatをincludeした後、FLT_EVAL_METHODマクロを調べ、これが0でなければddまたは-DKV_NOHWROUUNDを使ったプログラムはコンパイルできないようにしました。これが0であれば、数値の演算器上の長さとメモリ上の長さが一致します。Intelの32bit環境の場合そのままではコンパイルできず、-msse2 -mfpmath=sseとオプションを付けるなどしてSSE2を強制的に使わせるようにすればコンパイルできるようになります。FLT_EVAL_METHODが定義されない環境、あるいはSSE2を持たないIntel CPUは、古すぎるのでサポート外としました。

内部型の異なる区間同士の変換

内部型の異なる区間同士の変換に詳しく書きました。従来は、例えばmpfr型の区間xをdd型の区間yに代入するときは、
    impfrtoidd(x, y);
のようにしなければいけなかったのが、単に
    y = x;
と書けるようになりました。
OK キャンセル 確認 その他