[未掲載分] 「除算が遅い」の補足 (6)

アプリを作る際、レジスタの本数が足りず演算結果を一旦メモリ上に保存したり、読み込み直すことが速度向上の妨げとなる。「除算が遅い」の補足 (3)で綴ったように、32ビット版のアプリを組む際に顕著である。
「除算が遅い」の補足 (4)や「除算が遅い」の補足 (5)では、お手軽に新しく追加された機能を活かすアプリを作る方法を触れた。
しかし、お手軽な方では、期待したほど高速化されないケースも生じる。理想的な命令とは異なる命令を用いた実行コードが生成されるケースも考えられる。

XMM レジスタの幅は128ビット分。C/C ++言語で言えば、double 型の数値を２つ分、float 型の数値を４つ分同時に扱える。ゆえに、その分高速化されるハズ。ところが、同時に複数の加算や乗算が行われないような実行コードに翻訳されてしまうこともある。この辺はソースコードの組み方だけでなく、コンパイラの賢さにも左右される。もちろん、コンパイラのバージョンが上がる毎に改善されている。

前話で述べたように、高速化をコンパイラに任せる副作用として、要件を満たさない環境で実行した際のトラブルも希に発生する。
そこで、今回は
/arch:SSE2 オプションを適用しない。
高速化したい箇所のみ新機能を活かすようなソースコードに書き換え
について・・・

---

お知らせ
活動休止にともない、この記事を事前に予約投稿してあります。
トップ記事の固定を目的としています

---

実数 ( 浮動小数点数 ) と整数で状況は異なります。今回は、お手軽ではない方の浮動小数を中心としたお話。

C/C++ 言語で浮動小数を扱いたい場合、一般的には float や double といった変数の型を用いる。倍精度浮動小数点数型のdouble 型 を軸に・・・ double 型の数値は64ビット幅うんぬん・・・と文字では判り難いヒトもいることでしょう。簡単な図で示します。

まずは FPU と 倍精度浮動小数点数。



FPU ( 浮動小数点演算処理装置 ) は 80ビット幅の数値まで扱えます。
西暦1993年頃から、CPU にFPU が標準搭載された。それまで、FPU はオプション扱いで、搭載されていない PC も多かった。
今後、FPU は廃止へと向かっています。古くからの書籍やWeb サイトでは FPU を前提とした記事が載っているのでは・・・

次に、(SSE 登場時より追加された) XMM レジスタでのイメージを図で表すと下のよう。



レジスタの幅が128ビットあり、64ビット幅の数値であれば２つ分を扱える。命令ひとつで複数の演算が可能となっている。FPU が除算を１回行うのに比べ、XMM レジスタ内の２組分の除算を行えるのがメリット。クロック数、つまり、待ち時間はほぼ一緒。
※ 一部のCPUを除く。

簡単に言えば、トラックの荷台が広くなり、今までの倍の荷物を積めるようになったようなもの・・・

さて、XMM レジスタが下図のような状況もありうる。



左半分のグレー「未使用」と記した部分は、「使われない」「利用されない」という意味。

冒頭に述べた理想的な命令とは異なる命令に翻訳された場合はこのような状況。XMM レジスタのうち、下位の半分しか「利用されない」。ひとつの命令で同時に２つの演算ではなく、ひとつの演算しか行わない。
原因としては、連続的に配置されていない、ハードウェア的に備わっていない機能を補うため ( 三角関数に関するライブラリを呼び出す ) などが挙げられる。

簡単に言えば、トラックの荷台がまだ半分空いているのに出発するようなもの・・・

---

高速さの話と反してしまうが、不慣れなヒト向けに少し補足。
演算結果の一貫性、つまり、精度に関して頭に入れておいたほうが良い。
大雑把な言い方をすれば、正しい値からのズレをどの程度妥協できるか・・・

FPU は 80ビット幅の数値まで扱えます。
Visual Studio で作成したアプリにおいて、double 型の数値は64ビット幅としてメモリ上に置かれたり、ファイルに書き込まれる。これらは読み出された後、FPU の中で80ビットの精度で演算が行われる。演算結果を FPU からメモリに書き出す際には64ビット幅の値に丸められる。
一方、XMM レジスタでは64ビット幅の倍精度のまま演算が行われる。もろん、ビット数が高い (桁数が多い) ほうが、結果は誤差が少なくなる。誤差という表記は混乱するかもしれません。

コンピュータは正しい結果を示すと思い込んでいるヒトが多い。
たいていの PC において浮動小数点数の計算は、正確な答えに一番近い数を探し出しているに過ぎません。
誤解が生じないように加えると、PC のはじき出した数字が信用できないということではありません。正確さが重視されるアプリ等は、人間からみてより正確な答えに近づくよう補正する仕組みを備えています。

人間の目線で考えて「1.0」が正しい答えであるとしても、途中で誤差が生じて
「0.99999999～～～」や「1.0000 ～～～たくさんのゼロ ～～～ 03」等の答えが返ってくることもあります。これはバグや欠陥ではありません。

例えば「Ａ ÷ Ｂ × ○△□ × Ｃ」といった式において、変数 Ａ、Ｂ、Ｃ の内容は実行時まで不明だとします。実行時にＡ が 1.0 、Ｂ が 3.0 、Ｃ が 3.0 の数が代入されたならば、
「1.0 ÷ 3.0 × ○△□ × 3.0」となります。

人間の場合、分数の概念があるので「1.0 ÷ 3.0」は３分の１、さらにそれを３倍して通分して～～～と考えることができます。最終的には Ａ、Ｂ、Ｃ 以外の部分を計算した結果と答えが等しいと判断できます。
一方、PC の内部では分数ではなく、「1.0 ÷ 3.0」の部分は「0.33333333～～～」といった割り切れない数として扱われます。よって、末尾の「× 3.0」で３倍して「0.99999999～～～」等の数値が返ってきます。
先ほど「正しい答え」ではなく「正確な答えに一番近い数」と記したのはこのようなケースがあるから。

しばしば、高い精度のほうが好まれる。それは割り切れない数などから生じる誤差をなるべく少なくしたいからである。突き詰めれば、より正確な値を扱う手法もあるのだが高速ではない。
とはいえ、画像や音を処理に関して言えば、よほど敏いヒトでないかぎり感じないレベル。
高速化を目指すなら、ある程度の近い値で妥協するといった割り切りも必要である・・・

---

例として、
Ａ ＝ Ｂ × Ｃ ÷ Ｄ を10000回繰り返すとしましょう。流れとしては、ロード・演算・ストアのひたすら繰り返し。

過去記事除算が遅いなどで
・浮動小数点の除算ならば、「逆数の掛け算」に置き換えると良いと綴りました。

Ａ ＝ Ｂ × Ｃ ÷ Ｄ は
Ａ ＝ Ｂ × ( Ｃ ÷ Ｄ ) と同じ。
基のＢ は毎回メモリから読み出すため異なる値、Ｃ、Ｄ の値が決まった値、つまり固定としよう。( Ｃ ÷ Ｄ ) の部分を予め求めておく。それを Ｅ ＝ Ｃ ÷ Ｄ と置き換えることにより、演算の部分 は Ａ ＝ Ｂ × Ｅ となる。単純に見て、演算のステップがひとつ減る。
ちょっと判り難いかも。何らかの数を当てはめてみれば判り易くなるだろうか。
Ａ ＝ Ｂ × Ｃ ÷ Ｄ の Ｃ は ４、Ｄ は ５としよう。
Ａ ＝ Ｂ × ４ ÷ ５となる。
Ａ ＝ Ｂ × ０．８
・・・といった具合に、除数 ( 割るほうの数 ) が変わらないならば、遅くなりがちな除算を避けることも可能。

今回は意図的に除算を避けない形で話を進めます。さらに、前提条件として、基のＢは 10000個の値が連続的に配置されているものとします。

高速さを求めるならば、
読み出す値がどのように格納されているかも重要なカギ。



連続的に配置されているならば、データを読み出しがスムーズ。
そうでない場合、先に述べた XMM レジスタの下位半分しか使わない命令に翻訳される確率が高くなります。

簡単に言えば、何かを作ろうとして原料・材料があちこちバラバラに散っている場合よりも、ひとまとめに揃えてある方が速やかに捗るというのと同じ。
もし、基となる数値がバラバラに散っているのであれば、その辺の見直しすだけでも効果があります。
演算ループに入る前に暫定的なメモリ空間を確保し、基の数値を連続的して配置。つまり、一列に整列させてから演算ループに入ると良いでしょう。ループから脱したら暫定的なメモリ空間の解放をお忘れなく。
メモリ空間を確保や解放については、過去記事車輪の再発明 (6)や車輪の再発明 (8)で触れていますので省きます。
また、それらの過去記事で触れたように、メモリ空間のアドレスが16の倍数に揃っていれば若干速度が向上します・・・

予め記しておきます。汎用レジスタが対象の場合、メモリからデータを「ロード」、結果をメモリに格納する際には「ストア」の語を用います。FPU に対する操作は 「ロード」、「ストア」ではなく「プッシュ ( 積む )」「ポップ 」を用います。Push の反対はPullではないかとの声も聞こえてきそうですが・・・

Ａ ＝ Ｂ × Ｃ ÷ Ｄ を10000回繰り返す。
まずは、旧来からのFPU ( 浮動小数点演算処理装置 ) を経由して演算を繰り返すとして
(A-1) 64ビット幅の浮動小数Ｂをスタックへプッシュ
(A-2) FPU のスタックの値とＣを掛け
(A-3) スタックにある値をとＤで割り
(A-4) FPU スタックの値をメモリへポップ
(A-5) 浮動小数Ｂが格納されているポインタを64ビット幅の浮動小数ひとつ分後方へ・・・
といった感じでしょうか。 (A-1)から(A-4) を10000 回繰り返します。

次に、XMM レジスタを使った、つまり、新しい機能を活かした流れは
(B-1) 64ビット幅の浮動小数Ｃを２組 XMM レジスタ(0) へロード
(B-2) 64ビット幅の浮動小数Ｄを２組 XMM レジスタ(1) へロード
(B-3-1) 64ビット幅の浮動小数Ｂをふたつ、Ｂ[0] と Ｂ[1] を XMM レジスタ(2) へロード
(B-3-2) XMM レジスタ(2) と XMM レジスタ(0)を掛け
(B-3-3) XMM レジスタ(2) と XMM レジスタ(1)で割り
(B-3-4) XMM レジスタ(2) の値をメモリへストア
(B-3-5) 浮動小数Ｂが格納されているポインタを64ビット幅の浮動小数ふたつ分後方へ・・・
といった感じで (B-3-1)から(B-3-5) を5000 回繰り返す。

10000回繰り返すのでは！？！？！？
と誤解するヒトがいるかもしれませんので補足。
一見、(A-2) と (B-3-2)、(A-3) と (B-3-3) が同じように見えます。
前者は一度に一組の乗算や除算、後者は一度に二組の乗算や除算を行います。同等のクロック数 ( 待ち時間 ) で２組分の除算が行える点が後者のメリット、アドバンテージ。

流れ (A)、(B) に即したコードを比べてみましょう。まず、流れ (A)

for(i = 0 ; i < 10000 ; i += 1){ double b = *p;// (A-2) b *= c;// (A-3) b /= d;// (A-4) p += 1;// (A-5)

・・・のようになる。
簡単に説明すると、変数 C、D は for(分の前に値を代入、変数 i は ループ回数、変数 p は浮動小数Ｂを読み出す位置、つまりポインタ。

次に、流れ (B) に従い SSE2 命令を使うように変更を加えるならば、

__m128d xmm0 = _mm_set1_pd(4.0);// (B-1) __m128d xmm1 = _mm_set1_pd(5.0);// (B-2) for(i = 0 ; i < 10000 ; i += 2){ __m128d xmm2 = _mm_loadu_pd( p );// (B-3-1) xmm2 = _mm_mul_pd(xmm2 , xmm0);// (B-3-2) xmm2 = _mm_div_pd(xmm2 , xmm1);// (B-3-3) p += 2;//(B-3-5)

・・・のような感じ。

__m128d データ型の部分は、２組のdouble 型をひとつに扱うための変更。流れ(A) と間違わないように変数名を xmm ～と付けてあります。
ほか、パッと見、_mm_ ○△□ _pdがたくさん!?!?!?
この辺の読み方は過去記事車輪の再発明 (8)で触れています。
頭に「_mm_」が付くのが共通として、
○△□ の部分が掛け算なら MUL (Multiply の略) 、 割り算ならDIV ( Divide の略 )、
末尾「_pd」は同時に２組の64ビット幅の数値演算を行うことを意味する。ちなみに、末尾が「_sd」ならば１組の64ビット幅しか扱わない、つまり、XMM レジスタのうち下位の半分を使うという命令。



数値をロードする部分は _mm_loadu_pd と _mm_load_pd の選択肢があります。
loadの後ろに「u」が付く・付かないの違いがあります。
前者は16バイトの倍数以外のアドレスからデータを読み出せますが、後者は16バイトの倍数以外のアドレスからデータを読み出そうとすると例外が発生し、アプリの不正終了へとつながります。
データを読み出すアドレスが16バイトの倍数に整っているのが不確実なので前者を用いました。データのアドレスが16バイトの倍数であることが確実ならば、後者に切り替えることで若干速度が向上します。

たしかに、_mm_ ○△□ _pdを覚えるのは面倒。Visual Studio にはそれらを軟らげる便利なクラスが用意されている。dvec.h や fvec.h 等のヘッダーファイルをインクルードすることで便利なクラスを使うことができる。
なかでも __m128d 型 に対応するのはF64vec2 クラス。Float , 64ビット幅、２組分の略と考えれば覚え易いでしょう。
これを用いれば

F64vec2 xmm0(4.0 , 4.0);// (B-1) F64vec2 xmm1(5.0 , 5.0); // (B-2)

F64vec2 xmm2(*(double *) p);// (B-3-1) xmm2 *= xmm0;// (B-3-2) xmm2 /= xmm1;// (B-3-3)

・・・のような感じ。
_mm_ ○△□ _pdなどのコンパイラ組み込み関数を用いる場合に比べ、流れ(A) から変更箇所が少なく済む・・・


※ これは2012年7月 5日に掲載した記事の画像です

まだまだ、速度向上の余地がありそう。「除算が遅い」の補足 (4)で取り上げた
並列に実行できる、できない
などを練り直すことで、さらなる速度向上を狙えるかも！！！

32ビット環境では 8本、64ビット環境では 16本 の XMM レジスタが備わっています。XMM レジスタをさらに２組使うとして

(C-1-1) 64ビット幅の浮動小数Ｃを２組 XMM レジスタ(0) へロード
(C-1-2) 64ビット幅の浮動小数Ｃを２組 XMM レジスタ(1) へロード
(C-2-1) 64ビット幅の浮動小数Ｄを２組 XMM レジスタ(2) へロード
(C-2-2) 64ビット幅の浮動小数Ｄを２組 XMM レジスタ(3) へロード
(C-3-1-1) 64ビット幅の浮動小数Ｂをふたつ、Ｂ[0] と Ｂ[1] を XMM レジスタ(4) へロード
(C-3-1-2) 64ビット幅の浮動小数Ｂをふたつ、Ｂ[2] と Ｂ[3] を XMM レジスタ(5) へロード
(C-3-2-1) XMM レジスタ(4) と XMM レジスタ(0)を掛け
(C-3-2-1) XMM レジスタ(5) と XMM レジスタ(1)を掛け
(C-3-3-1) XMM レジスタ(4) と XMM レジスタ(2)で割り
(C-3-2-1) XMM レジスタ(5) と XMM レジスタ(3)を掛け
(C-3-4-1) XMM レジスタ(4) の値をメモリへストア
(C-3-4-1) XMM レジスタ(5) の値をメモリへストア
(C-3-5) 浮動小数Ｂが格納されているポインタを64ビット幅の浮動小数 ４つ分後方へ・・・
といった感じにすれば (C-3-1-1)から(C-3-5) は 2500回となります。演算を繰り返す回数も 1/4。
ひとつ前の実行結果を待つ必要が無い部分は並列に実行されるハズ。

理想としては、当初の FPU を経由した演算に比べ４倍に向上してほしいところです。
「なぜ４倍なの!?!?!?」と疑問に感じるヒトもいることでしょう。トラックの荷台で言えば、荷台の広さが２倍、利用できるトラックの台数も２倍に増えたのに等しいから。

実際は、並列処理を手助けする回路に左右されます。昨今の CPU では 実行ユニットを複数備えています。
実行ユニットの数が多いほど並列処理も速やか。体感から言えば、これを書いている時点の CPU では ２～３命令が並列に実行されているようです。
単純に４倍の速度向上には至らないまでも、当初に比べ 2.5 ～ 3 倍は速くなるのでは・・・

細かいことを言えば、XMM レジスタは並列処理が得意な反面、EAX ～ EDX 等の汎用レジスタに比べロードやストアに比べ3 ～ 5 倍のクロック数 ( 待ち時間 )。
※ なお、これを下書きしたのは2012年頃。ご覧いただいている時点では CPU の進化とともに、ロード、ストアは遅さは解消されているかもしれません。

かつてはこの辺を人間の手で調整( ロードや演算の手順を修正 ) することでさらなるパフォーマンス改善も期待できた。同時に、ソースコードが曲芸的に偏りがちで、読み難くなる傾向にあった。
昨今は CPU の進化にともない、メモリアクセスの効率化、先に実行できる命令か否かなどを能動的に判断するような仕組み等も強化されている。
後日の保守作業を想定するならば、曲芸的なコードを書くことは避けたいものである・・・

長くなりましたので今回はこの辺で・・・