2024/08/27(火)kv-0.4.57
今回は、
- 以前にるふぁさんにご指摘頂いた、defint_autostepで使用する級数の次数を1にしたときに正しく計算できていなかったバグの修正。
- ode_maffineに自動微分型を指定し、なおかつcallback関数を指定したとき、 callback関数が初期値に関する微分の情報をも渡してくれるように変更。
大したアップデートではないですが、一応2つ溜まったので。
2024/02/07(水)g++-13と_Float64xと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
今年の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の区間演算は遅い?
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 carbon | 38.1 | 19.4 | 35.7 |
M1 MBA | 142 | 114 | 23.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 carbon | 6.43 | 6.44 | 6.53 |
M1 MBA | 6.29 | 6.29 | 15.04 |
なお、これはaarch64なアーキテクチャ全般に見られる傾向かどうかが気になったので、3万円のchromebook (Lenovo IdeaPad Duet Chromebook, MediaTek Helio P60T) で試してみたら、
-O3 | -O3 -DROUND_CHANGE1 | -O3 -DROUND_CHANGE2 | |
---|---|---|---|
IdeaPad Duet | 15.8 | 20.7 | 20.9 |
ARMでのベンチマーク追加 (2022/6/30追記)
- OCI (Oracle Cloud Infrastructure)のAlways Free ARMで提供されている、Ampere Altraプロセッサ。Ubuntu 20.04、gcc 9.4.0。
- Planexのスマートフォン、Gemini PDA (CPUはMediatek MT6797X Helio X27)にTermuxで入れたclang 14.0.1
-O3 | -O3 -DROUND_CHANGE1 | -O3 -DROUND_CHANGE2 | |
---|---|---|---|
X1 carbon | 6.43 | 6.44 | 6.53 |
M1 MBA | 6.29 | 6.29 | 15.04 |
IdeaPad Duet | 15.8 | 20.7 | 20.9 |
Ampere Altra | 6.69 | 6.69 | 9.70 |
Gemini PDA | 19.7 | 23.0 | 28.8 |
なお、とある方から富岳のA64FXでのベンチマークを聞かせてもらったのですが、少し異常な感じの結果(とても遅い)だったので何か測定ミスの可能性もあり、ここへの掲載はとりあえず控えます。
2022/01/06(木)kv-0.4.54
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による丸めの向きの変更を一切使わない丸めエミュレーション
そこで、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;と書けるようになりました。