JKフリップフロップに関する誤解

JKフリップフロップという便利な回路素子があります。 日本語で「JKフリップフロップ」と検索すると、多くの場合、それと異なるJKラッチの説明が出てきて混乱しがちです（JKラッチは便利ではないので使うことはあまりないと思います）。 これを整理するのがこの記事の目的です。

ラッチとフリップフロップ

ラッチやフリップフロップは、1ビットの情報が記憶できる素子で、順序回路を構成するために欠かせない基本素子です。 その名称の使い分けには三種類の流儀があります。 一つの流儀を一貫して使えばよいのですが、複数種類の流儀を混ぜると意味不明になります。 フリップフロップ - Wikipediaやラッチ回路 - Wikipediaはいろいろな流儀のものを寄せ集めて作ったのか、流儀が混在していて、それが原因で混乱が生じています。

	制御信号のないもの	レベルセンシティブな（＝Enable信号が1の時のみ透過になる）もの	エッジトリガな（＝Clock信号のエッジでのみ出力が変化する）もの
流儀1	セット・リセットタイプのラッチ	トランスペアレントタイプのラッチ	エッジトリガタイプのラッチ
流儀2	非同期式フリップフロップ	非同期式フリップフロップ	同期式フリップフロップ
流儀3	ラッチ	ゲート付きラッチ	フリップフロップ
例	SRラッチ、S̄R̄ラッチ、JKラッチ	ゲート付きSRラッチ、ゲート付きDラッチ、ゲート付きJKラッチ	Dフリップフロップ、JKフリップフロップ

本記事では順序回路 - WikipediaやFlip-flop (electronics) - Wikipediaで採用されている流儀3を使用していきます。

（セット・リセットタイプの）ラッチ

SRラッチ

最も基本的なラッチに、NORゲート二つを使ったSRラッチがあります。 SRラッチは、SとRという二つの入力と、QとQ̄という二つの出力を持ちます（Q̄はない場合もあるかもしれません）。 配線は、次の図のようになっています（引用元：File:RS Flip-flop (NOR).svg - Wikimedia Commons）。 上側のNORゲートはRとQ̄を入力としてその出力がQになります。 下側のNORゲートはSとQを入力としてその出力がQ̄になります。 このように、出力が入力に戻されるような、循環的な構造を持つことで情報を記憶します。

図1: SRラッチの回路図

SRラッチと呼ばれる理由は、片方の端子がSet（出力Qを1にする）、もう片方の端子がReset（出力Qを0にする）というわかりやすい意味を持つからです。 SとRが共に0の時、QとQ̄の値が現状維持となります。

動作例として、例えばQが1、Rが0で、Sが1である状況を考えます。 この時、下側のNORゲートの出力Q̄は、Qが1かつSが1なので0になります。 そして、上側のNORゲートの出力Qは、Rが0かつQ̄が0なので1になります。 このようにSRラッチ回路には循環がありますが、同じところに戻ってくるまでにNOT的な素子を二回通るようになっていて、循環が整合的に閉じているので発振しません（後述のメタステーブル状態になった場合を除く）。

この状態でSを0に変化させたとします。 この時、下側のNORゲートの出力Q̄は、Sが0であってもQが1なので0になります。 上側のNORゲートは、先ほどと入力が変わっていないので出力も同じで、出力Qは1になります。

下側のNORゲートの出力Q̄が0になる原因をたどってみると、Qが1であることが原因で、そしてそれはQ̄が0であることが原因です。 ここでQ̄が0になった原因を考えてみると、いつかの時点でSが1になったことが原因のはずです。

Sを1にすることはQを1にする効果がありますが、その後にSを0にしてもQは1のままです。 つまり、Qが1であるというのは、「以前にSが1になったことがある」という意味を持つ出力であり、この意味において情報を記憶する能力を持っています。

同様のことはRについても言えます。 RはSと逆に、Qを0にする効果を持ちます。 したがってより正確には、Qが1であるというのは「以前Rが1になった時よりも後に、Sが1になったことがある」を、Qが0であるというのは「以前Sが1になった時よりも後に、Rが1になったことがある」を、それぞれ意味します。

Sを1にすることは電池への充電、Rを1にすることは電池の放電、と考えるとわかりやすいかもしれません。 電池を充電しきった後、充電をやめても電池は充電された状態のままです。 電池を放電しきった後、放電をやめても電池は放電された状態のままです。 そして電池の今の状態は、充電と放電のどちらを最後に行ったかによります。

禁止入力(S,R)=(1,1)について

SRラッチにおいて、SとRをともに1にすることは禁止と言われます。 これは電池のアナロジーからも明らかで、充電しながら放電もするというのはおかしな話です。 NORゲート二つからなるSRラッチの場合には、ハードウェアが壊れるわけではありませんが、QもQ̄も0になってしまい、相補性（QとQ̄が常に異なる性質）が破れます。 この状態を積極的に活用することで回路が小さくなるみたいな曲芸がありますが、あまりやるべきではなさそうです。

他の用途としては早押し検出回路があります。 SとRがともに1の状態からともに0の状態に移行させようとしたとき、ふつうはSとRが変化するタイミングが完全に同一とはいかず、ほんのわずかな時間差があるはずです。 すると、(S,R)=(1,0)である瞬間か(S,R)=(0,1)である瞬間のどちらか一方があるはずです。 SRラッチはそれを見逃さず、前者であれば「以前Rが1になった時よりも後に、Sが1になったことがある」を示すQ=1を、後者であれば「以前Sが1になった時よりも後に、Rが1になったことがある」を示すQ=0を、それぞれ出力するようになります。

ただし、本当にほぼ同一タイミング、あるいはゲート遅延程度の時間差しかないような場合、相当長い時間出力が安定しなくなるメタステーブル状態になることがあります。 (S,R)=(1,1)が禁止入力と言われる背景には、潜在的にこのような問題を抱えることもあります。 なお、論理回路シミュレータだと循環を通して互いを追いかけるように状態が更新され続けるような挙動になったりしますが、実際の信号がそういう感じとは限りません。

S̄R̄ラッチ

NORゲートの代わりにNANDゲートを使って同じようなものを作ると、Qを1にしたいときには入力Aを0にする必要があり、Qを0にしたいときには入力Bを0にする必要があるようなラッチができます。 つまり、入力が負論理（何かをしたいときに信号を1にするのではなく、何かをしたいときに信号を0にする必要がある）になっています。 入力AはQを1にするというSet操作を負論理で行うためS̄と、入力BはQを0にするというReset操作を負論理で行うためR̄と、それぞれ呼ばれます。 この二つをとって、S̄R̄ラッチと呼ばれます。 配線は以下のようになっています（引用元：File:SR Flip-flop Diagram.svg - Wikimedia Commons）。

図2: S̄R̄ラッチの回路図

ちなみに、S̄R̄ラッチの禁止入力は(S̄,R̄)=(0,0)ですが、この時の出力はSRラッチと異なり(Q,Q̄)=(1,1)になります。 (S,R)=(1,1)が禁止入力であると言われる背景には、このように、中身をどう構成したかによって禁止入力に対する出力が変わってしまうことも挙げられます。

JKラッチ

SRラッチで禁止入力であった(S,R)=(1,1)に対して、出力を反転させるといった機能を割り当てたものがJKラッチです。 端子の名前はリセットする方をK、セットする方をJと呼んでいます。 端子の名前は、SRラッチの発展版がJKラッチであることと順番が対応していることから、Set→J、Reset→Kと覚えればよいです。 あるいは、0にする方＝killする方がKだと覚えればよいです。

JKラッチは単独では役に立ちません。 出力を反転させるという機能が直感的で回路を組むうえで便利だということ自体は正しい説明です。 しかしながら、JKラッチは(J,K)=(1,1)とした瞬間に出力が反転し続けます（発振します）。 つまり元の状態を「1回だけ」反転させる、という機能は実現できません。

図3: JKラッチの回路図

回路図はエレキ素人が何か考える（その９）：JKフリップフロップ #電子回路 - Qiitaを参考にして自分で描きました。

ゲート付きラッチ

単なるラッチは、入力信号が変化するとそれがすぐに出力に反映されてしまいます。 組み合わせ論理回路を適当に設計すると変なグリッチが立ったりしてしまうことがありますが、その際に意図せず記憶した情報が失われてしまうことにつながります。

ゲート付きラッチは、「入力の変化を反映させる」意味を持つEnable信号を明示的に持ちます。 この信号が0の時には、SやR（S̄・R̄・J・K）がどのように変わっても出力は変わりません。 この信号が1の時には、従来のラッチと同様に動作します。 このEnable信号をうまく使うと、以下のようにすることができて、意図しない情報消失を防げます。

• 計算期間。Enable信号を0にしておく。S入力やR入力につながる前段の組み合わせ論理回路でグリッチが立っても、Enableが0なので出力は変わらず、データの消失は起こらない。
• 更新期間。Enable信号を1にしておく。計算済みのS入力やR入力が取り込まれ、記憶状態が更新される。計算期間を十分長くとっておけば、更新期間中にグリッチが立つことはないはず。必要であれば(S,R)=(0,0)などとして記憶状態を維持することも可能。

ゲート付きラッチの作り方（ゲート付きSRラッチを例に）

ゲート付きラッチを作るのは簡単です。 S入力やR入力などにEnable信号とのANDゲートを追加するだけです。 例えば、ゲート付きSRラッチの回路図は以下のようになります（引用元：File:SR (Clocked) Flip-flop Diagram.svg - Wikimedia Commons）。

図4: ゲート付きSRラッチの回路図

（ゲート付き）Dラッチ

次の回路はゲート付きDラッチです（引用元：File:D-type Transparent Latch (NOR).svg - Wikimedia Commons）。 ゲート付きDラッチは、単にDラッチと呼ばれることが多いです（“ゲートなし”のDラッチはデータの保持ができなくて、記憶素子として役に立たないためです）。 この回路は、Enableが0の時はデータを保持、Enableが1の時は入力をそのまま出力、といった動作をします。 SRラッチが「1に変える」「0に変える」といった入力を持つのに対して、「入力をそのまま記憶する」という直感的な動作をする点がうれしいです。

図5: ゲート付きDラッチの回路図

出力を再利用するみたいな匠の技を使うと、以下のように作れるようです（引用元：File:D-Type Transparent Latch.svg - Wikimedia Commons）。

図6: ゲート付きDラッチの別の回路図

ゲート付きDラッチの回路シンボルは、Dフリップフロップと似たようなものが使われることもありますが、明示的に区別するために以下のようなものをよく見かけます。

図7: ゲート付きDラッチの回路図記号

回路図記号は私がhttps://infoscience.epfl.ch/record/178151/files/scm.pdfを参考にして自分で描きました。

この記号の出典はよくわかりませんでした。 standard cell-based memoryの論文でしか見かけないので、その分野に特有の記法なのかもしれません。 とはいえWikimedia commonsを見ると、手書きで「]」記号を使っている図（File:D latch flip-flop functional symbol.png - Wikimedia Commons・File:D master-slave flip-flop circuit.png - Wikimedia Commons）や、「)」記号を使っている図（File:Flip-flop synchronization types schematic.svg - Wikimedia Commons）を見つけることができました。 何かのツールで使われている記法が自然と広まっただけで標準化された記法ではないとかでしょうか。

ゲート付きJKラッチ

次の回路はゲート付きJKラッチです。

図8: ゲート付きJKラッチの回路図

ゲート付きJKラッチは次のようにも作れます。 後段のNORをNANDに変えたので、SRラッチからS̄R̄ラッチに変わりました。 それに伴い、上下が反転しました（JとKが入れ替わりました）。 SとRではなくS̄とR̄になったので、前段はANDではなくNANDになりました。

図9: ゲート付きJKラッチの別の回路図

Enable信号はクロックではない

これらの回路のEnable信号（ゲート信号）はフリップフロップのクロック信号とは性質が大きく異なります。 Enable信号はあくまで、入力（S・R・J・K・Dなど）の情報を後ろに通す通さないかを制御する信号です。 ゲート付きJKラッチは、(J,K)=(1,1)の時にEnable信号を1にしている間中ずっと出力が反転し続けるため、役に立つよう設計するためには工夫が必要です。 これに対してクロック信号を持つJKフリップフロップは、(J,K)=(1,1)の時、クロック信号の立ち上がり（または立ち下がり）の瞬間に一度だけ出力が反転するため、直感的に扱いやすいです。

そもそも、Enable信号が1の間は結局のところ、入力が素通しです。 ラッチの記憶状態が更新された時、その出力は先につながる組み合わせ論理回路に入力されますが、その計算結果がラッチの入力に戻ってくることがあります。 クロックに同期した回路のいいところはそういうことが起こらないことです。 そういったことが発生しうるということを明確にするためにも、ラッチ回路におけるゲート信号をクロックと呼ぶのはやめたほうがいいでしょう。

多くの日本語文献で混乱が見られるのは、JKラッチのことを流儀2でJKフリップフロップと呼びつつ、しかもEnable信号のことをクロック信号と呼んでいることに起因します。 これにより見た目上、流儀3の意味でのJKフリップフロップ（クロック信号がある）と区別がつかず、混乱が生じているようです。 流儀を明確にする・回路図を示す・Enableと呼び分ける、などで違いを明確にできるので、この記事ではそのようにしています。

Dラッチがデータを“ラッチする”様子

ゲート付きラッチ、特にDラッチにデータを覚えさせることを、（データを）“ラッチする”と言ったりします。 Enableが1から0になる（＝扉が閉まる）瞬間のデータを挟んで落とさないようにするようなイメージです。 次に新たにデータを覚えようとする場合、「Enableが1から0になる瞬間のデータを覚える」わけですから一旦Enableを1に戻す必要があり、この時には挟んでいたデータを落としてしまいます。 このイメージをアニメーションにしてみました。

動画1: Dラッチの動作のイメージ

Enableが1の場合は素通しで、Enableが0の時はDを無視して現状維持、という動作が確認できると思います。

パルスラッチ

ゲート付きラッチのEnable信号に、（0と1が一定周期で交代するようなクロック信号ではなく）一瞬だけ1になり他の時間はずっと0であるような信号を入れることで、上記の問題を解決することができます。 このような構成をパルスラッチと呼ぶことがあります。 ラッチはフリップフロップよりも小さく高速であるなどの利点がありますが、タイミング制約が厳しくなるなどの欠点もあり、素直にフリップフロップを使う設計が人気だそうです（参考：ウェスト＆ハリス CMOS VLSI回路設計 応用編 丸善出版 (2014) p.556）。

マスタースレーブ型構成のフリップフロップ

パルスラッチのようなことをしなくても、うまく回路を作ることで、クロック信号の立ち上がり（または立ち下がり）の瞬間のみで記憶状態が変わるようにすることができます。 そのような方法の一つに、マスタースレーブ型構成があります。 マスタースレーブ型構成では、ゲート付きラッチを直列に二つ並べます。 以下、ゲート付きラッチのことを単にラッチと呼びます。 前段（入力側）のラッチをマスターラッチ、後段（出力側）のラッチをスレーブラッチと呼ぶようです。 ここで、二つのラッチは互いに反転した信号をEnable信号として受け取るようにします。

このようにすると、入力が出力に素通しになることがなくなり、前述の問題が解決されます。 動作を順番に見ていきます。 マスターラッチがEnableの時、スレーブラッチはEnableではありません。 したがって、マスターラッチが素通しでも、スレーブラッチは現在の出力をそのまま続けます。 ここでクロックが変化すると、スレーブラッチがEnableで、マスターラッチがEnableでない状態に変わります。 この時スレーブラッチの状態が更新されます。 その後スレーブラッチは素通し状態になりますが、スレーブラッチの入力を作っているマスターラッチの出力は変化しないので、出力は固定されます。 その後クロックが再度変化すると、マスターラッチがEnableでスレーブラッチがEnableでない状態に戻ります。 この時、スレーブラッチの出力が固定されるので、やはり出力はそのままです。 このように二つのゲート付きラッチを、同時に開かない二重扉のごとく使うことで、入力の素通しを防ぎます。 このような構成は、マスターラッチがEnableでなくなった瞬間の入力を切り取って、それを一度だけ状態更新に使うような動作をします。

具体的な回路図をいくつか紹介します。

Dフリップフロップ

次の回路図は、Dラッチを二つ直列並べることで作ったDフリップフロップです（引用元：File:D-Type Flip-flop Diagram.svg - Wikimedia Commons）。 このDフリップフロップは、クロックが立ち上がった瞬間のD入力の値を、次にクロックが立ち上がるまでの間出力し続けます。

図10: マスタースレーブ型構成のDフリップフロップの回路図（ポジティブエッジトリガ）

以下のアニメーションを見ると、二重扉を使って入力の素通しを防いでいること、マスターラッチの扉が閉まるタイミングのみで出力Qが変わるようになっていることがわかると思います。

動画2: Dフリップフロップの動作のイメージ（ネガティブエッジトリガ）

Dフリップフロップの回路図記号は、以下のようなものが使われることが多いです。

図11: Dフリップフロップの回路図記号（ポジティブエッジトリガ）

回路図記号は私がhttps://infoscience.epfl.ch/record/178151/files/scm.pdfを参考にして自分で描きました。

JKフリップフロップ

次の回路図は、役に立つJKフリップフロップの回路図です。 ゲート付きJKラッチを前段に、ゲート付きS̄R̄ラッチを後段に、それぞれ配置して作っています。 ここで、ゲート付きJKラッチは「QやQ̄をフィードバックする」という本質部分を抜き出す形で採用している点に注意してください。 単にゲート付きJKラッチをゲート付きS̄R̄ラッチの前に配置するだけだとJKラッチ部分が発振する問題が解消されないのでダメです。

図12: マスタースレーブ型構成のJKフリップフロップの回路図（ポジティブエッジトリガ）

回路図はJKフリップフロップ - 電気主任技術者のナレッジノートを参考にして自分で描きました。

汎用ロジックICに入っているフリップフロップ

上記のマスタースレーブ構成はたくさんの素子が必要なため（でしょうか？）、汎用ロジックICの中にはより洗練された構成のフリップフロップが実装されています。

次の回路図は、汎用ロジックICに入っているDフリップフロップの回路を簡略化したものです。 S̄R̄ラッチの各入力の前にそれぞれS̄R̄ラッチがついていて、なんだか入力を整える働きをしているようです。 しかしながら回路構成がかなり非対称であり、しかもS̄R̄ラッチの禁止入力である(S̄,R̄)=(0,0)を有効利用することで後段のS̄R̄ラッチの保持入力(S̄,R̄)=(1,1)を作っているなど非常に技巧的な構成で、なぜこれでいいのかを直感的に説明できるような構成にはなっていません。 一体どうやったらこのような回路を思いつけるのでしょうか……。

図13: 洗練されたDフリップフロップの回路図（ポジティブエッジトリガ）

同様の技巧的な構成はJKフリップフロップにも応用できて、以下のような対称な形の回路になるようです。

図14: 洗練されたJKフリップフロップの回路図（ポジティブエッジトリガ）

回路図はhttps://ocw.u-tokyo.ac.jp/lecture_files/engin_04/7/notes/ja/07.pdfを参考にして自分で描きました。

汎用ロジックICに入っているJKフリップフロップの回路はこれとはちょっと違うようです。 四角で囲った部分はいずれもラッチになっていません。

図15: 洗練されたJKフリップフロップの回路図その2（ネガティブエッジトリガ）

回路図はTEXAS INSTRUMENTS社のデータシートを参考にして自分で描きました。

動作を確認してみます。

Clockが0の時、NANDゲートが常に1を出力しているため、外側のANDゲート経由で循環ができています（内側のANDゲートは常に0を出力していて循環に影響しません）。 Clockが0から1になる瞬間、外側のANDゲート経由の循環が途切れる前に内側のANDゲート経由の循環が完成します。 外側のAND経由の循環が途切れるまでにはNANDゲート一段分の遅延が存在するため、その短い時間を使ってバトンタッチが行われるのです。 これに加えて、残り2つのANDゲートのうちどちらか一方も循環に貢献します。 Qが0の場合、下のNANDは1を出力するため、下のANDゲート経由の循環が存在します。 Qが1の場合、上のNANDは1を出力するため、上のANDゲート経由の循環が存在します。

Clockが1から0になる瞬間が状態更新のタイミングです。 Clockが1から0になる瞬間は、内側のANDゲート経由の循環がまず先に途切れ、NANDゲート一段分の遅延の後、外側のANDゲート経由の循環が成立します。 この隙間が状態の更新の鍵です。 外側のANDゲート経由の循環が最初から成立している場合はバトンタッチできますが、そうでない場合にはバトンタッチに失敗します。

まず、(J,K)=(0,0)を考えます。 この場合、NANDゲートの出力はClock信号にかかわらず1であり、常に外側のANDゲート経由の循環が成立しています。 よって、Clockが1から0になるときも正しくバトンタッチが行われ、状態が保持されます。

次に、(J,K)=(0,1)を考えます。 まず、上のNANDゲートは1を出力しているため、上のANDゲート経由の循環は成立しています。 また、もともとQ=0であれば、下のANDゲート経由の循環はすでに成立しています。 よって、Q=0の場合は、状態は保持されます。 Q=1の場合、その情報を保持するための循環がNANDゲート一段分の遅延時間の間、完全に途切れます。 K=1により下のANDゲート経由の循環が有効になっていないためです。 これによりQ̄が1となり、上のANDゲート経由の循環を通じてQが0になります。

(J,K)=(1,0)の場合はこれと対称なので省略します。

最後に(J,K)=(1,1)を考えます。 もともとQ=1の時は、各素子の出力は(J,K)=(0,1)の時と同じであり、Q=1を伝える循環が一瞬途切れることにより状態がトグルされます。 もともとQ=0の時は、各素子の出力が(J,K)=(1,0)の時と同じであり、Q̄=1を伝える循環が一瞬途切れることにより状態がトグルされます。

ゲート一段分の遅延を有効利用してバトンタッチしたり、逆にそれにより生じるハザードを利用してバトンを落とさせて状態を更新したりする、賢い構成になっていることがわかりました。

CMOS VLSIで使われているフリップフロップ

LSI（CPUとか）で使われているDフリップフロップは以下のような、トランスミッションゲートを使用した回路になっているらしいです（transmission-gate flip-flop, TGFF）。

図16: TGFFの回路図

回路図はhttps://onlinelibrary.wiley.com/doi/10.1155/2016/8268917を参考にして自分で描きました。 実際には入出力を安定させるために追加のNOTゲートをDの後とQの前に入れるのが一般的なようです。

この図のうち、青い四角でくくった部分はトランスミッションゲート、黄色い四角で囲った部分はNOTゲートです。 緑の四角で囲った部分は見慣れない構成ですが、点線部分があればNOTゲートとトランスミッションゲートの組み合わせになっています。 点線がない場合は、クロックトトライステート/クロックトCMOS/C²MOSなどと呼ばれる素子になります（参考：ウェスト＆ハリス CMOS VLSI回路設計 応用編 丸善出版 (2014) p.541）。 この点線部分はあってもなくても論理的な動作には影響がありませんが、タイミング・面積・電力消費等に利点欠点があるため、点線がある構成とない構成の両方を見かけます。 この構成はマスタースレーブ構成ですが、通過可能・通過不能を直接実現するトランスミッションゲートを使うことでトランジスタ数を大幅に減らしています。

この構成は他の構成と大きく異なりますが、非常に直感的な設計になっています。 advanced-computer-organization/aco-shioya-02.pdf at master · shioyadan/advanced-computer-organization · GitHubの52ページに書いてあるような回路をそのままCMOS回路に落とし込んだものです。 ようするに、NOTゲート二つで作ったループ（1ビットの情報を記憶することができる）に値を取り込むかを単純にスイッチのON/OFFで制御しているということです。

まとめ

ラッチとフリップフロップの用語の混乱を整理し、代表的なラッチやフリップフロップの回路図を紹介しました。

回路図の参考文献

回路図を書く時に参考にしたと明示したもの以外で、配線が同じ回路図であることを確認したもの

• File:Edge triggered D flip flop.svg - Wikimedia Commons
  • 洗練されたD-FF
• https://www.ti.com/lit/ds/symlink/sn74s74.pdf
  • 洗練されたD-FF
• https://ocw.u-tokyo.ac.jp/lecture_files/engin_04/7/notes/ja/07.pdf
  • 洗練されたD-FF, マスタースレーブ型JK-FF, 洗練されたJK-FF
• 5. フリップフロップ — Digital Electronic Circuits 1.0 ドキュメント
  • マスタースレーブ型JK-FF
• 論理回路 / JKフリップフロップ | meyon's STUDY
  • マスタースレーブ型JK-FF
• エッジトリガ型フリップフロップ
  • 洗練されたD-FF, 洗練されたJK-FF

本文中で示したゲート付きDラッチの回路図記号が使われている文献

• https://ieeexplore.ieee.org/document/8050995🔒
• A Low Power 1024-Channels Spike Detector Using Latch-Based RAM for Real-Time Brain Silicon Interfaces🔓
• https://dl.acm.org/doi/pdf/10.1145/2890498🔓

去年の振り返りとか

2026年の最初の週が終わろうとしています。 今年もよろしくお願いします。

なんとなく去年はあまり何もしなかった一年だった気がしていたのですが、振り返って書き出してみるといろいろやっていたことがわかって意外でした。 今までは一年の振り返りとかはあまりやっていなかったのですが、書いておこうという気分になったので書いておきます。 本当はこういうのは年末にやるべきことのようにも思いますが、年末は忙しくて時間がなかったのと、年末ぎりぎりにいろいろ出来事もあったので、きっと今がまとめる時なのでしょう。

1月

GPU上で行列積を計算するコードと戦っていました。 あと、何度か東京に行った時に、中央線のグリーン車（お試し期間）に乗ってみました。

TODO: はねにトラッカーをつけてVRをやる

2月～3月

追っていた好きな作家さんが作っていたゲームが完成したらしいので、それを買って遊んでいました。 かなり長いこと、ゲームを買って遊ぶみたいなのはしていなかったので、なんだか新鮮でした。

コスプレ1回目

o3-mini-highに聞いて検索しづらい数学の疑問を解消したりしていました。 その他、5月ごろに量子計算の疑問を、10月ごろに生化学の疑問を、それぞれ解消していたようです。 このあたりの回答が正しいのかは私には判断できませんが、とりあえず分かった気にさせてくれるだけでうれしいです。

4月

能力が判明するみたいな話がありました。

5月

いろいろなところに行って、新しい知り合いをつくったり、旧交を温めたりしていました。

ChatGPTを使ってウェブサイトを作ってみたりもしていました。 「こんな感じにしたい」みたいな曖昧な要望で、それを実現する方法が出てくるのはさすがです。 私はJavascriptを書けませんが、雰囲気で読むことはできますし、コードを少し変えると挙動がどう変わるかを観察して何が起こっているかを把握するみたいなことはできます。 提示される複数種類の実装の中からなるべくコードが短そうなのを選ぶことで、把握しやすく、かつ定番の実装（と思われるもの）になるよう工夫しました。 私はデザインのセンスもなく、色を選ぶとかそういった簡単なことすらできないのですが、そのあたりも全部いい感じにやってくれて助かりました。 一から作る方法がわからなくても、良いか悪いかの判定と悪い点の指摘ができれば、少なくとも自分が満足するものは作れるという事実に感動します。

あと、ChatGPTに私の心理構造を分析してもらって遊んだりしていました。

6月

コスプレ2回目と3回目。 小道具の作り方に新たな技法が導入され、10年来*1の懸案が解決しました。 つくり方を褒められたりもしたのですが、なかなか素直に喜べない自分が恨めしいです。

Xで流れてきたゲームが面白そうだったので買って遊び、かなり良かったのでさらに同じジャンルのゲームを漁りました。 ゲームを買って遊ぶというのが習慣づいてきた気がします。 ただ、そのゲームのRTA動画を見てしまって、自由度の高いゲームなのに稼ぎ効率みたいなのが脳裏にちらつくようになったのは残念でした。 あとは、レベルが上がってもだんだん耐性が出てきて達成感が得られなくなるのが悲しい。

7月

Lean 4 を知ったので入門しました。 これは今年一年で一番大きな出来事だった気がします。 定理証明支援系は前からやってみたいとは思っていたのですがなかなかきっかけがなく、以下をきっかけに入門できたのは良いめぐりあわせでした。

https://t.co/rt2XwN1sns

定理証明支援系によるパーフェクトさんすう教室ありえなさすぎる pic.twitter.com/CxPxiF8dez

— SnO₂WMaN (@SnO2WMaN) 2025年7月9日

ただ、いきなり命題を証明するというのはかなり無謀で、証明を紙に書き下してからそれを形式化するという手順でやらないと厳しそうだと感じました。 体系だった数学の証明の訓練を受けていないので、その点に困難を感じました。 もともとできる場合にそれをブーストするものであって、もともとできない私が使うようなものではないのかもしれません。 あとは、私が浮動小数点数とかビット表現とかそのあたりの、低級であまりよい数学的構造を持たないものをいじっているのが原因かもしれません。 集合とか解析とかそのあたりの定理の証明には向いているように見えます（隣の芝生かもしれません）。

ひさしぶりにコスプレパフォーマンスをやりました。

8月

地元の（？）コスプレイベントであるコスサミに参加しました。

ChatGPTがバージョン5になって、前にできなかった数学の問題が解けたりで進化を感じました。 一番驚いたのは、単精度対数関数を精度よく実装する方法を聞いたときです。 とりあえず知っている高精度化のテクニックを列挙して、「これとこれは組み合わせることができない」「この組み合わせは試したがダメだった」という情報を与えたところ、私が考慮していなかった組み合わせで高精度になるような計算手順が出てきました（8/9）。 かなり高度な内容だと思っていたのですが、こういうことができるとすると、もうじきすべての分野で私を上回ることになりそうです。

夏コミは体調が悪かったのもあり、ものすごくひさしぶりにコスプレなしの参加でした。 友達がサークル参加しているのを見て、自分もサークル参加しようと思い立ち、申し込みました。 コスプレの時間を確保したいというややよこしまな目論見でもあります。 サークルカットの雰囲気はAIに提案してもらって、素材もAIに作ってもらいました（いわゆるAI絵ではないです）。

ゆりかごのそらが10周年でした。 もう10年も前のゲームなんですね。

あと、X上でのコミュニケーションに問題があったので改心しました……。

9月

TODO: mod_poppoさんに教えてもらった最適化を阻害するいろいろな技法を使ってみる

コスプレ8回目。 通常の反応を知ることができました。

10月～11月

メールの最後に、「～～も必要でしたらお申し付けください」みたいなことが書いてあって考えさせられました。 AIの出力の末尾にも、これと似た文面がついてくることが多いです。 しかし、その文面に対しては「できもしないことを書いてそれっぽくしている」「こちらがお願いすることがなさそうなタスクを提示している（ので実際にできなくても会話が破綻しづらいという読みがありそう）」と感じます。 「できもしない」というのはバージョン4を使ってた頃からの印象であって、バージョン5になったら実際にはできるのかもしれませんが。 一方でメールの方は、「こちらが依頼することがありそうなタスクを提示している（思いやりで提案している）」「現在の成果が不十分である可能性を考え謙遜している」と感じました。 これは果たしてAIへの差別的な態度なのか、それとも実際に思いやり（？）の程度に差があってこちらの受け取り方は正当なのか、よくわからなくなりました。

12月

冬コミで出す本に載せたかったので、空き時間を見つけてはひたすら、浮動小数点数加算器の正しさをLean 4で証明していました。 結局証明は全然完成せず、12/29に原稿を書いて12/30に表紙を作って印刷して東京に行って12/31にサークル参加する、みたいなひどい日程でやっていました。 よこしまな動機みたいな話もありましたが、自分の書いた本を目の前で手に取ってもらう、というのと、サークルスペースにコスプレした状態で居る、というのはやってみると良い体験でした。 コスプレ衣装は思い切って新しいものを仕立ててもらったのですが、細かい指定にもかかわらず良く再現されたものになっていて大変幸せな気持ちになりました。 2023年に環境が大きく変化して精神状態が良くなかったときに読んだ小説が心の支えになっていて、それがきっかけでそれ以来ずっとやりたいと思っていたコスプレだったのですが、今回ここで一つ区切りをつけることができました。 今回は開会2時間半くらいで準備した分がなくなってしまって、結果的に残り時間で会場を回れたのは良かったです（部数が読めなかったのと印刷が面倒だったのもありますが、よこしまな動機がなかったとも言えません）。 でもせっかくスペースを頂いているのだから、次回参加があれば、十分事前に告知して（今回は当日になってからの告知だった）、部数をもっと用意するようにしたいですね。 会場を回るためにはお手伝いさんを頼まないといけないかもしれません。 あとは、コスプレも準備が万端とは言えなかったので、再チャレンジしたいです。

その他、証明の息抜きで、Library checkerに上がっている、間違ったポラードのロー法の使い方をしているコードをいくつか撃墜しました。

あと、快適な執筆環境を求めて大きなディスプレイを買いましたが、スマートフォンが壊れました……。

おわり

特にまとまりはないですが振り返りました。

*1:記憶をたどってみたら、旧技法は少なくとも2013年3月以前から使っていたようです。

int8テンソルコアを用いたDGEMMエミュレーションライブラリは何を使えばよいか

最近のGPUは倍精度演算器が少なくDGEMM（倍精度行列積）を効率よく実行できませんが、代わりにたくさん載っているint8テンソルコアを使って実質的に倍精度行列積を計算する（＝エミュレーションする）方式が最近流行っているようです。 その先駆けは2023年6月に大友らがarXivに投稿したプレプリント（[2306.11975] DGEMM on Integer Matrix Multiplication Unit 🔓）で、その後内野らが改良手法を提案しています（Performance enhancement of the Ozaki Scheme on integer matrix multiplication unit - Yuki Uchino, Katsuhisa Ozaki, Toshiyuki Imamura, 2025 🔓）。 NVIDIAが2025年8月にリリースしたCUDA 13.0に含まれるcuBLASでは公式にint8テンソルコアを使用したエミュレーションが実装され、外部ライブラリに頼ることなく手軽に切り替えることができるようになりました。

しかしながらこれらの手法は、いずれも異なる精度と異なる速度を持っているため、どれを使うべきなのかよくわかりません。 内野らの論文には大友らの手法との比較が掲載されていますが、それより後に出たCUDA 13.0のcuBLASとの比較は当然ながら掲載されていません。 そこで、これらの手法の精度と速度を比較してみました。

環境

• WSL2 on Windows 11
• GeForce Game Readyドライバー581.80（2025-11-04）
• CUDA 13.0
• GeForce RTX 4090

記事を執筆しているうちにCUDA 13.1が出ましたが（2025-12-04）、もう一回実験するのも大変なのでCUDA 13.0でいきます。

また、以下では正方行列同士の行列積のみを取り扱います。 行列サイズNとは、一辺の長さのことを指します。 このとき、行列の要素数はN²、行列積の計算コストはΘ(N³)となります。

尾崎スキーム

尾崎スキームは、浮動小数点数の行列積結果を正確に求めるために固定小数点数の行列積を行う方法です。

尾崎スキームでは、まず入力の浮動小数点数行列を固定小数点数行列に変換します。 固定小数点数と言っても実際のところは、桁合わせされた複数の倍精度浮動小数点数の和として表されることが多いようです。 ここで、桁の区切りは53ビットとかではなく、例えば20ビットなどのあまり多くないビット数にするのがポイントです。 つまり、元の行列（M）を、2²⁰の位より上を持つ倍精度浮動小数点数の行列（A）＋残りのうち2⁰の位より上を持つ倍精度浮動小数点数の行列（B）＋残りのうち2⁻²⁰の位より上を持つ倍精度浮動小数点数の行列（C）＋残りのうち2⁻⁴⁰の位より上を持つ倍精度浮動小数点数の行列（D）＋残りのうち2⁻⁶⁰の位より上を持つ倍精度浮動小数点数の行列（E）＋……のように分割して表します。

桁の区切りを20ビットなどにしたのは、各行列同士の積を無誤差で計算できるようにするためです。 各行列（A・B・C・D・E・……）の各要素が20ビットしか仮数部を持っていないとき、その各積の仮数部は高々40ビットで表せます。 その場合、各行列のサイズが2¹³=4096以下であれば、各行列同士の積が無誤差になります。 各積は位が揃っているため、2¹³個足し合わせた結果が53ビットに収まるからです。 実際には論理が逆で、行列のサイズを見てから各行列が何ビット持つべきかを決めるようです。

さて、各行列同士の積は、通常のDGEMMライブラリが使えます。 求めたい行列積を M1M2 = (A1+B1+C1+D1+E1)(A2+B2+C2+D2+E2) と変換したとします。 このとき、右辺を展開すると A1A2 + A1B2 + A1C2 + A1D2 + A1E2 + B1A2 + …… + E1E2 となりますが、いずれの項も単独では行列積が無誤差で計算できます。 既存のチューニングされたDGEMMライブラリが使える点が尾崎スキームの重要な利点で、高精度行列積コードを一から作る労力から解放されます。

最後に、各項を足し合わせます。 ここはただ足し合わせるだけの処理で、特にメモリ律速なので、適当に書いたdouble-doubleの加算を使っておけばOKです。

尾崎スキームでは、精度と速度のトレードオフを柔軟に変更することができます。 例えば、E1E2のような小さい項が最終結果に与える影響は少なそうです。 その様な考えに基づいて、上から何ビットまでしか計算しない、のような最適化により、精度が低下する代わりに速度を改善することができます。 ただし、「上から何ビット」というのは、桁落ちにより「最終結果の上から何ビット」と大きく乖離する可能性があります。 なので、悪条件問題では、上から80ビット、とか計算しないと、最終結果を53ビット精度で求めることはできないかもしれません。 逆に言えば、問題の性質に応じて、精度に対して必要十分な計算量に設定できる点も尾崎スキームの利点の一つです。

GPUで使われる尾崎スキームは、無誤差計算部分でDGEMMの代わりにint8行列積を使う点が異なります。 ここまでで説明した方法は、普通にやってしまうと誤差が生じるDGEMMを、複数回の無誤差なDGEMMに分解して計算するものでした。 GPUで使われる方法は、普通にやってしまうと誤差が生じる――というか「遅い」――DGEMMを、複数回の無誤差なint8行列積に分解して計算するものです。 int8テンソルコア（int8行列積をint32に積算）を使う場合、上で説明したような技巧的な分割を行わなくとも、行列サイズが小さければ自動的に無誤差になります。 これは、もともと要素ごとの積の絶対値は2¹⁴以下なので、これを2¹⁷個くらい足してもint32に収まる、ということを利用しています。 最後に各項を足し合わせる時は、8nビットシフトして足し合わせるみたいなことが必要になります。

CUDA 13.0のcuBLASのDGEMMエミュレーションの観察

CUBLAS_EMULATE_DOUBLE_PRECISION

これを1に設定すると、DGEMMの呼び出しをint8テンソルコアを使った高速エミュレーションで置き換えることを許可できるようです。

CUBLAS_FIXEDPOINT_EMULATION_MANTISSA_BIT_COUNT

エミュレーション精度を指定することができます。 尾崎スキームでは普通、エミュレーション精度を変えたいときはスライス数を変化させますが、精度をビット数で指定するインタフェースのようです。 速度調査結果を見ると、

• 1～7を指定した時
• 8～15を指定した時
• 16～23を指定した時
• ……
• 8n～8n+7を指定した時
• ……

で速度が変化していないようでした。 このことからすると、1～7は1スライス、8～15は2スライス、16～23は3スライス、……使って計算しているのだと思います。 7ビットまでが1スライスというのはint8が符号を除いて7ビットであることと整合し、8～15ビットが2スライス・16～23ビットが3スライス・……というのは1スライスにつき8ビット分の表現力が追加されることと整合的です。 よく見たら、https://docs.nvidia.com/cuda/cublas/index.html#representation-and-mappingsにもそう書いてありました。

0を指定した場合、試した感じでは63を指定した場合（8スライス）と同じになっていました。 倍精度浮動小数点数の精度は53ビットで、これを正確に保持するためには55ビット分保持できる7スライスが必要ですが、絶対値が小さな値が出てくると固定小数点数に変換するときに誤差が生じてしまいます。 もう1スライス余分に計算しておけば、変換するときの誤差が最終結果に与える影響は十分小さいと言えるでしょう。 なので、「いい感じの」精度でやってくれるモードとして、8スライスでやっているのだと思います。

実験した範囲では9スライス以上にしても精度向上は見られなかったので、8スライスは必要十分なスライス数なのでしょう。 もしかすると、より悪条件な行列積において9スライス以上が有用なのかもしれません。 0をした場合の実行時間は8スライスの時より遅く9スライスの時より速いような感じだったので、行列の中身を見て条件数を推定してスライス数を切り替えたりする機能があるのかもしれません。

CUBLAS_EMULATION_STRATEGY

これはperformantまたはeagerに設定することができます。 eagerに設定すると、どんなDGEMMも指定した方式で実行するようです。 一方で、performantだと、

• N≦128の時、設定を無視して常にDGEMMを行う
• N>128の時、
  • CUBLAS_FIXEDPOINT_EMULATION_MANTISSA_BIT_COUNTが0～79の時はeagerの時と同じ
  • CUBLAS_FIXEDPOINT_EMULATION_MANTISSA_BIT_COUNTが80以上の場合は常にDGEMMを行う
    • CUDA_EMULATION_MANTISSA_CONTROL_DYNAMICのデフォルト値が79だからそうなる？（https://docs.nvidia.com/cuda/cublas/index.html#dynamic-mantissa-control）

のような動作をしているようでした。 けっこう雑な切り替えのようです。

私の環境では残念ながら、N=1024でさえDGEMMの方が速かったので、128で切り替えるのは微妙そうです。 逆にN=2048では、精度80ビット以上でDGEMMに切り替わると精度が落ちるわ速度も落ちるわで何もいいことがなさそうでした。

内野らの改良

Performance enhancement of the Ozaki Scheme on integer matrix multiplication unit - Yuki Uchino, Katsuhisa Ozaki, Toshiyuki Imamura, 2025を読み解いていきます。

errfree_sum

3.2節に書いてある方式です。 int8行列積の結果はint32行列ですが、ozIMMUではこれをdoubleに変換してから足していました。 一方でint32加算がオーバーフローしないなら、int32加算の方がdouble加算より高速なので、足してからdoubleに変換する方が高速なのでそうする、という方式のようです。 基本的には精度は変わらず、単に高速になる効果だけが得られるようです。

nearest_split

3.1節に書いてある方式です。 仮数部を分割する際、下位部分を切り捨てて分割するのではなく、下位部分を四捨五入丸めで分割することにより、スライス当たり1ビット表現力が高まるという方式のようです。 精度は上がるものの、四捨五入丸めを実現することには相応のコストがかかります。 そのため、論文のFigure14を見た感じでは、0.3スライス追加したみたいな感じの速度低下と精度向上になるようです。

最近接丸めを使うと128みたいなint8に収まらない数が出て来て困るのではないかなと思っていましたが、これは違うようでした。 分割結果は、-64～64しか使わないようになっていました。

hybrid

3.3節に書いてある方式です。 要するに上の二つは独立なので両方盛り込みましたという方式です。

実験

速度と精度の実験を行いました。 実験環境は「環境」節で書いた通りです。 入力行列は、各要素が独立に[-1.0,1.0]の一様分布に従うようなdouble行列としました。 より正確には、以下のような行列積計算を行いました。

    int N = atoi( argv[1] );
    std::mt19937_64 mt;
    std::uniform_real_distribution<double> dist( -1.0, 1.0 );

    double* a = (double*)malloc( N*N*sizeof(double) );
    double* b = (double*)malloc( N*N*sizeof(double) );
    double* c = (double*)malloc( N*N*sizeof(double) );

    for( int i = 0; i < N*N; ++i ) {
        a[i] = dist( mt );
        b[i] = dist( mt );
        c[i] = dist( mt );
    }

    double* device_a; cudaMalloc( (void**)&device_a, N*N*sizeof(double) );
    double* device_b; cudaMalloc( (void**)&device_b, N*N*sizeof(double) );
    double* device_c; cudaMalloc( (void**)&device_c, N*N*sizeof(double) );
    cudaMemcpy( device_a, a, N*N*sizeof(double), cudaMemcpyHostToDevice );
    cudaMemcpy( device_b, b, N*N*sizeof(double), cudaMemcpyHostToDevice );
    cudaMemcpy( device_c, c, N*N*sizeof(double), cudaMemcpyHostToDevice );

    cublasHandle_t handle; cublasCreate(&handle);
    double alpha = 1.0, beta = 1.0;

    cublasDgemm( handle, CUBLAS_OP_N, CUBLAS_OP_N, N, N, N, &alpha, device_b, N, device_a, N, &beta, device_c, N );

上のソースコードは、row-major形式の行列積を行うコードです。 cblasDgemmはcolumn-major形式で渡す必要があるので、行列積の順番を逆にすることで実質的にそれを実現しています。 row-majorとcolumn-majorの関係は転置に相当することを利用して、のように変換したということです。

実行時間は、cublasDgemmの前後にcudaEventRecord関数を挟んでcudaEventElapsedTimeで計測しました。 実行時間は、cublasDgemmを51回呼び出し、その中で最も高速だったものを採用しました。 これは、他のプロセスが動いていたりして理想的ではない環境で実行していることによる速度低下を取り除くためです。 また、初回実行はライブラリの遅延読み込みの時間がかかり、まともな実行時間計測結果にならないという問題を取り除くことも意図しています。

精度は、別途計算した“正解”データとの最大相対誤差（max_relative_error）を計測しました。 今回の実験の範囲では、並列性に起因して再現性がないなどの現象は確認されず、どの手法でも毎回同じ結果が得られました。

“正解”データは、CPU上でdouble-double行列積を行うことによって求めました（真の正解とはずれがあると思います）。 double-double行列積は、乗算にTwoProdを、加算にsloppy加算を用いました。 そのコードは以下の通りです。

struct HiLo { double hi; double lo; };
HiLo TwoProd( double a, double b ) {
    double hi = a * b;
    double lo = std::fma( a, b, -hi );
    return { hi, lo };
}
HiLo FastTwoSum( double a, double b ) {
    double hi = a + b;
    double lo = a - hi + b;
    return { hi, lo };
}
HiLo TwoSum( double a, double b ) {
    double hi = a + b;
    double tmp = hi - a;
    double lo = a - (hi - tmp) + (b - tmp);
    return { hi, lo };
}
HiLo DD_MulAdd( double a, double b, double dh, double dl ) {
    HiLo mul = TwoProd( a, b );
    HiLo add = TwoSum( dh, mul.hi );
    double lo = mul.lo + add.lo + dl;
    return FastTwoSum( add.hi, lo );
}

    for( size_t i = 0; i < N; ++i )
    for( size_t k = 0; k < N; ++k )
    for( size_t j = 0; j < N; ++j )
    {
        HiLo sum = DD_MulAdd( a[i*N+k], b[k*N+j], dh[i*N+j], dl[i*N+j] );
        dh[i*N+j] = sum.hi;
        dl[i*N+j] = sum.lo;
    }

測定対象のエミュレーション手法は以下の通りです。

• cuBLAS DGEMMそのまま
• cuBLAS DGEMMエミュレーション・performant・精度0～128
• cuBLAS DGEMMエミュレーション・eager・精度0～128
• ozIMMU スライス数3～18およびauto*1
• ozIMMU/errfree スライス数3～18およびauto
• ozIMMU/nearest スライス数3～18およびauto
• ozIMMU/nearest+errfree スライス数3～18およびauto

ozIMMUは、執筆時点で最新版の08eea9231729d54dbfd92955f2cbfc21ec236856（2025-12-02）を使いました。 三種の改良版は、d458fec19e29279556885a9da968d6fe0ac51058（2025-05-06）を使いましたが、これは動作しないため、ozIMMU_get_function_pointer() arguments are flipped · Issue #1 · RIKEN-RCCS/accelerator_for_ozIMMU · GitHubで示されているような修正を加えました。 なぜ本家にコントリビュートせず、差し替えソースコードを配布（しかも本家に追随できていない）みたいな運用になっているのでしょう？？

（追記：2025-12-07に三種の改良版がアップデートされたようですが、記事公開直前で気づかなかったため再実験はしていません）。

暖房を26度に設定し、Windowsからのディスプレイ出力がない状態で実験しました。

N = 1024

N = 1024の時の実験結果を図1に示します。 ただしozIMMU四種については、実験結果の図が見づらくなるのを防ぐため、精度が全く変わらず単に遅くなっただけのスライス数12～18はプロットしませんでした。 cuBLASのエミュレーションも同様に、精度が全く変わらず単に遅くなっただけのスライス数9以上はプロットしませんでした。 横軸はmax_relative_errorの二進対数に負号をつけたもので、右に行くほど高精度です。 おおよそ「何ビット精度か」というのを表していると考えてよいです。 例えば、横軸の値が27であれば、最大の相対誤差が2⁻²⁷であるということで、おおよそ27ビット精度ということになります。 縦軸は、本当にDGEMMを実行したと仮定した時に必要な2N³ FLOPを実行時間で割った、疑似的なFLOPS数で、上に行くほど高速です。 なお、RTX 4090の倍精度演算の理論性能は、1.29 TFLOPSです。

図1: N=1024の時の速度と精度の比較

cuBLASのエミュレーションは、DGEMMとほぼ同等の速度でした。 精度を下げると最大の相対誤差が増える一方でやや速くなるため、DGEMMに退化しているわけではなく、本当にint8テンソルコアを使っていると思います。

一方でozIMMU四種は、8スライス以上でDGEMMより速いのに精度が高い、という結果になりました。 9スライス以上では精度が変わっていませんが、これはdouble-doubleで計算した“正解”の方が間違っていることを示唆しているかもしれません。

N = 2048

N = 2048の時の実験結果を図2に示します。 図の見方は図1と同様です。

図2: N=2048の時の速度と精度の比較

cuBLASのエミュレーションは、DGEMMを超える性能となりました。 一方で、依然としてozIMMU四種の方がcuBLASのエミュレーションより高速です。 cuBLASのエミュレーションは、何かものすごく大きな固定オーバーヘッドがあるように見えます。

N = 4096

N = 4096の時の実験結果を図3に示します。 図の見方は図1と同様です。

図3: N=4096の時の速度と精度の比較

cuBLASのエミュレーションが、ozIMMU四種の性能とほぼ同じになりました。 固定オーバーヘッド（？）のせいでスライス数が小さい時は依然としてozIMMUよりも遅いようです。

N = 8192

N = 8192の時の実験結果を図4に示します。 図の見方は図1と同様です。

図4: N=8192の時の速度と精度の比較

スライス数4以上ではozIMMUとほぼ同等の精度-速度トレードオフが見られます。

N = 16384

N = 16384の時の実験結果を図5に示します。 図の見方は図1と同様です。

図5: N=16384の時の速度と精度の比較

やはり4スライス以上ではozIMMUとほぼ同等の精度-速度トレードオフが見られます。

なおcuBLASのエミュレーションでは、N=16384かつスライス数が2の時、Windowsが落ちる挙動が頻繁に発生しました。 おそらく電力を消費しすぎて電源の保護回路が働いたのだと思います。 火事になったりしても怖いので正確な条件は調べていませんが、そもそもこんな計算をコンシューマー向けGPUでやるべきではないのかもしれません。

スライス数あたりの精度

結果を見ると、cuBLASのエミュレーション6スライス版が、ozIMMU7スライス版と同等の精度を達成しています。 ozIMMUは6ビット毎、ozIMMU/nearestは7ビット毎に分割している一方、cuBLASのエミュレーションは8ビット毎に分割しているので、こういうことになっているのだと思います。

8ビット毎に分割するのは、おそらく下の方のビットから切り出して行っていると思います。 上の方のビットから切り出そうと思うと、特殊な丸めが必要になるからです。 ozIMMUは切り捨て丸めで6ビット毎の分割を、ozIMMU/nearestは最近接丸めで7ビットの分割を実現しています。 一方、8ビット毎に分割する場合、int8で表現できる数は-128～127と端（128）が含まれていないので、切り出すときの丸めに工夫が必要です。 下位桁を最大にすることで ULPを実現できる一方、下位桁を最小にすると ULPを実現できます。 よって、最近接丸め（0.5で切り上げか切り下げがが変わる）ではダメで、で切り上げか切り下げかが変わる丸めが必要とされます。 下の方のビットから切り出すときはこのような工夫は特に必要ないため、下の方のビットから切り出していると思います。

結局どれを使えばいいのか？

一様乱数で入力行列を使った今回の実験の範囲では、RTX 4090に関しては以下のようなことが言えると思います。←ChatGPTに入れておいた方がいいよって言われて追加した文

N=2048くらい以上で倍精度相当の行列積を行いたい場合、四種のozIMMUとcuBLASのエミュレーションは同等の速度と精度でした。 なのでこの場合には、外部ライブラリが不要で手軽なcuBLASのエミュレーションモードを使用すればよいでしょう。

一方でそれ以外の場合、つまりN=1024程度の小さな行列積を高速化したい場合や、倍精度だと過剰精度なので少ないスライス数で高速に計算したい場合は、四種のozIMMUはcuBLASのエミュレーションよりも速いです。 四種のozIMMUの中ではozIMMU/errfreeまたはozIMMU/nearest+errfreeが精度-速度トレードオフが優れているので、これを使うのがよさそうです。

また、四種のozIMMUはいずれもオープンソース（MITライセンス）であるという点も見逃せない点です。 特にozIMMU/errfreeまたはozIMMU/nearest+errfreeはcuBLASのエミュレーションと同等以上の性能を達成しているため、手法の研究や改良を行いたい場合はこれをもとにするのがよさそうです。

*1:変換誤差上限がデフォルトの0なのでDGEMMが呼び出されるだけだった

fatal: Fetched in submodule path 'dejagnu', but it did not contain c298959ef991b389b64a825f70094738c6a48780. Direct fetching of that commit failed.と出てRISC-Vのクロスコンパイラがビルドできない問題の対処

RISC-Vのクロスコンパイラが欲しい場合、もっとも楽なのはsudo apt install gcc-13-riscv64-linux-gnuなどとすることですが、これだとRV64GC向けのバイナリが生成されるため、困ることがあります。 -march=rv64gなどとすれば、ユーザーコードはRV64Gの範囲（C拡張なし）でコンパイルできるのですが、標準ライブラリはRV64GCになっているため、結局C拡張を含むバイナリができてしまいます。 そうなっては困る用途では、https://github.com/riscv/riscv-gnu-toolchainをクローンして自分でビルドしたクロスコンパイラを使うしかなさそうです。

ところが、このリポジトリはgit submoduleの仕組みでGitHubの外のリポジトリを参照していて、これが問題を引き起こします。 例えば、以前は動作していた以下の手順が、現在では動作しなくなっています。

git clone --recursive https://github.com/riscv/riscv-gnu-toolchain # すごく時間がかかる
cd riscv-gnu-toolchain
git checkout 2024.04.12 -f
git submodule update --init --recursive -f

上記手順は以下のようにエラーになります。 submoduleであるbinutilsはチェックアウトできますが、同じくsubmoduleであるdejagnuをチェックアウトできません。

Submodule path 'binutils': checked out 'c7f28aad0c99d1d2fec4e52ebfa3735d90ceb8e9'
error: Server does not allow request for unadvertised object c298959ef991b389b64a825f70094738c6a48780
fatal: Fetched in submodule path 'dejagnu', but it did not contain c298959ef991b389b64a825f70094738c6a48780. Direct fetching of that commit failed.

どうやらhttps://git.savannah.gnu.org/git/dejagnu.gitは、コミットハッシュで直指定してフェッチすることができない設定になっているようです。 GitHubはコミットハッシュで直指定して取ってくることができますが、普通のgitリポジトリはそれができないのがデフォルトです（uploadpack.allowReachableSHA1InWantとかuploadpack.allowAnySHA1InWantの規定値はfalseです。セキュリティ上の理由とか到達可能性判定の計算コストとかの理由によります。参考：Git - git-config Documentation）。 submodule先がそういう設定のリポジトリだとこういう問題が起きうるので、submodule先は同じorganizationにクローンして管理しておくのがよさそうです。

さて、解決方法ですが、少なくとも以下の二つが使えました。

ミラーを使う

Google Open Sourceのサーバー（dejagnu - Git at Google）は、少なくともc298959ef991b389b64a825f70094738c6a48780のコミットハッシュ直指定フェッチ要求を受け入れてくれます。 なので、

git config -f .gitmodules submodule.dejagnu.url https://gnu.googlesource.com/dejagnu
git submodule sync dejagnu

git submodule update --init --recursive -f

のようにやればdejagnuの必要なコミットをチェックアウトできて、ビルドに進めます。

深さ指定フェッチを行う

https://git.savannah.gnu.org/git/dejagnu.gitは、refのコミットしか（直には）取得できないとはいえ、その親コミットを取得することはできます。 そこで、

git -C dejagnu fetch origin master --depth=37

git submodule update --init --recursive -f

のようにやればc298959ef991b389b64a825f70094738c6a48780までたどり着けました（今後masterが伸びたら、37個では足りなくなるかもしれません）。 これにより、submoduleをすべてチェックアウトしてビルドに進めます。

まとめ

submoduleで記録されているコミットハッシュだけからではコミットを取得できない設定のリポジトリがあるので注意しましょう。

Golden Cove, Gracemont, Lion Cove, Skymontでのレイテンシとスループットの調査

先週の記事（Lion Coveにおける乗算命令のスループットの調査 - よーる）でLion Cove（Arrow LakeやLunar LakeのPコア）でのmulx命令のスループットを調査しました。 今後の課題として挙げていた、Gracemont（Alder Lake・Raptor Lake・Rapotor Lake RefreshのEコア）とSkymont（Arrow LakeやLunar LakeのEコア）における、mulx命令のスループットを調査します。 ついでに浮動小数点数命令や整数乗算命令のレイテンシとスループットを調査します。 また、Golden CoveとLion Coveでの値も載せました。

浮動小数点数演算命令のレイテンシ

Golden Cove（Pコア）

以前の記事（IntelのFast Adder (FADD)について - よーる）で調べた結果の流用です。

浮動小数点数加算命令のレイテンシは2 cycle、浮動小数点数乗算命令と浮動小数点数積和演算命令のレイテンシは4 cycleでした。 ただし、2 cycleが実現されるのは浮動小数点数加算の結果を浮動小数点数加算命令に渡すときだけで、その他の場合は3 cycleです。

Gracemont（Eコア）

実験環境：Intel Core i9 14900K (Raptor Lake-S Refresh)

Eコアの周波数は4.385 GHzでした。

浮動小数点数加算命令のレイテンシは3 cycle、浮動小数点数乗算命令のレイテンシは4 cycle、浮動小数点数積和演算命令のレイテンシは6 cycleでした。

Lion Cove（Pコア）

実験環境：Intel Core Ultra 9 285K (Arrow Lake-S)

Pコアの周波数は5.683 GHzでした。

浮動小数点数加算命令のレイテンシは2 cycle、浮動小数点数乗算命令のレイテンシは3 cycle、浮動小数点数積和演算命令のレイテンシは5 cycle、でした。 ただし、浮動小数点数加算から浮動小数点数乗算にフォワーディングする場合、1サイクル分のペナルティがあるようです。 これらは、以下の実験結果からの推論です。

• addの連鎖：2サイクル
• mulの連鎖：3サイクル
• fmaの連鎖：5サイクル
• addとmulを交互：6サイクル
• addとfmaを交互：7サイクル
• add-mul-fmaの順：11サイクル
• add-fma-mulの順：10サイクル

Lion Coveで浮動小数点数プログラムを動かすと妙に遅くて変だと思ったのですが、これはGolden Cove系列（Alder Lake, Raptor Lake, Rpotor Lake RefreshのPコア）でfmaが4 cycleだったのと比べて1 cycle遅くなったことが原因のようですね。

Skymont（Eコア）

実験環境：Intel Core Ultra 9 285K (Arrow Lake-S)

Eコアの周波数は4.586 GHzでした。

浮動小数点数加算命令のレイテンシは2 cycle、浮動小数点数乗算命令のレイテンシは3 cycle、浮動小数点数積和演算命令のレイテンシは4 cycleでした。 Golden CoveやLion Coveの場合、フォワーディング先次第で浮動小数点数加算命令のレイテンシが2 cycleにならない場合がありましたが、Skymontの「浮動小数点数加算命令のレイテンシは2 cycle」はフォワーディング先にかかわらず成り立つようです。

Skymontの方は、fmaのレイテンシが4 cycleです。 このため、積和演算の実時間で見たレイテンシは、Skymontが4 cycle÷4.6 GHz＝0.87 ns、Lion Coveが5 cycle÷5.7 GHz＝0.88 nsと、わずかではありますがEコアであるSkymontの方が短くなってしまっています。実際に測ってみるとこの逆転を確認することができます。

レイテンシまとめ表

	Golden Cove	Gracemont	Lion Cove	Skymont
add(→add)	2	3	2	2
add(→mul)	3	3	3	2
add(→fma)	3	3	2	2
mul	4	4	3	3
fma	4	6	5	4

浮動小数点数演算命令のスループット

Golden Cove（Pコア）

浮動小数点数加算命令のスループットは2（逆数スループットが0.5）、浮動小数点数乗算命令と浮動小数点数積和演算命令のスループットは2（逆数スループットが0.5）ですが、合計のスループットは3（逆数スループットが0.333）です。 これは、浮動小数点数加算をPort 1とPort 5で、浮動小数点数乗算と浮動小数点数積和演算をPort 0とPort 1で、それぞれ行うからです。 したがって、addと(mul/fma)の比が1:2～2:1であればスループット3、そうでなければ多い方の命令についてスループット2となります。

Gracemont（Eコア）

浮動小数点数加算命令と浮動小数点数乗算命令と浮動小数点数積和演算命令のスループットは2（逆数スループットが0.5）でした。 つまり、1 cycleに、addとmulとfmaを自由な組み合わせで2回することができます。

Lion Cove（Pコア）

浮動小数点数加算命令のスループットは2（逆数スループットが0.5）、浮動小数点数乗算命令と浮動小数点数積和演算命令のスループットは2（逆数スループットが0.5）、でした。 つまり、1 cycleに、addとmulをそれぞれ2回ずつ行ったりaddとfmaをそれぞれ2回ずつ行ったりすることはできますが、mulとfmaは合計で2回までしか行えません。

Skymont（Eコア）

浮動小数点数加算命令と浮動小数点数乗算命令と浮動小数点数積和演算命令のスループットは4（逆数スループットが0.25）でした。 つまり、1 cycleに、addとmulとfmaを自由な組み合わせで4回することができます。 おそらく、Javascriptのような、なんでもかんでも浮動小数点数で計算する言語で書かれたプログラムを高速に実行するためだと思いますが、それにしてもとんでもない構成です。

スループットまとめ表

	add	mul	fma
Golden Cove	addのスループット2	mul/fmaのスループット2
	ただし、合計でスループット3まで
Gracemont	add/mul/fmaのスループット2
Lion Cove	addのスループット2	mul/fmaのスループット2
Skymont	add/mul/fmaのスループット4

整数乗算命令

Golden Cove（Pコア）

先週の記事（Lion Coveにおける乗算命令のスループットの調査 - よーる）で調査しました。

整数乗算命令はどれも同じ演算器で実行され、その演算器のスループットは1のようでした。

Gracemont（Eコア）

実験環境：Intel Core i9 14900K (Raptor Lake-S Refresh)

mulx命令のスループットは1でした。

一方で、imul(64)命令のスループットは安定しません。 imul(64)を6個・8個・9個含むループで、逆数スループットがそれぞれ3.00, 4.00, 4.50のものを見つけることができたので、理論スループットは2（逆数スループットが0.5）となりそうですが、これを達成するのはかなり困難です。 NOPを詰めるとスループットが上がることがあるので、デコードグループの切れ目が影響していそうです。

混ぜてもスループットが改善する兆しはありませんでした。 混ぜたときの最良のスループットは以下のような感じでした。

• 9×imul(64), 1×mulxで逆数スループットが5.75
• 6×imul(64), 1×mulxで逆数スループットが4.16
• 6×imul(64), 4×mulxで逆数スループットが7.98

Lion Cove（Pコア）

先週の記事（Lion Coveにおける乗算命令のスループットの調査 - よーる）で調査しました。

64ビット乗算の上位部分を計算する演算器のスループットが1、下位部分を計算する演算器のスループットが3（逆数スループットが0.333）、となっているようです。

Skymont（Eコア）

実験環境：Intel Core Ultra 9 285K (Arrow Lake-S)

imul(64)命令のスループットは2（逆数スループットが0.5）でした。

一方で、mulx命令のスループットが安定しません。 mulx命令を2個含むループで、逆数スループットが1.5未満のものを見つけることができました。 mulx命令を増やすと1命令あたりの逆数スループット0.75は達成できなくなりますが、依然として1命令あたりの逆数スループットが0.8未満のものは見つけることができます。 NOPを詰めるとスループットが改善することがあるため、やはりデコードグループの切れ目に影響を受けているようです。

個々の命令のスループットが不明であるため、混ぜたときのスループットは測っていません。

また、レイテンシは、64ビット乗算の下位部が4 cycle、上位部が5 cycleでした。

まとめ

Golden Cove, Gracemont, Lion Cove, Skymontにおける、整数乗算命令のスループット、および浮動小数点数命令や整数乗算命令のレイテンシとスループットを測定しました。 Eコアは効率的を謳っており実際周波数が低いですが、スカラの浮動小数点数演算のスループットが高い上にレイテンシもPコア以上に短いなど、決して“遅い”コアではないことがわかりました。

Lion Coveにおける乗算命令のスループットの調査

x86における乗算命令としてはmul命令が古くからありますが、Haswell世代からmulx命令というのが追加されました。 これらの命令は、64bit符号なし整数と64bit符号なし整数の積を128bit符号なし整数として求める命令です。 一方で、64bit整数と64bit整数の積（の下64bitだけ）*1を求める命令には、imul命令があります*2。

ところで、最近のCPUには、1cycleに乗算命令を複数実行可能なものがあります。 いくつかのウェブサイトで公開されている解析結果を読むと、その命令だけを多数回実行した時のスループットは分かりますが、複数の異なる命令を混ぜて実行した時のスループットは分かりません。 今回はそれを調査してみました。

乗算命令の種類

今回計測する命令は以下の4種類の命令です。

• imul regA, regB：regAとregBの積をregAに格納（64bit乗算）
• mul regA：regAとraxの積の上位64bitをrdxに、下位64bitをraxに格納（符号なし128bit乗算）
• imul regA：regAとraxの積の上位64bitをrdxに、下位64bitをraxに格納（符号つき128bit乗算）
• mulx regA, regB, regC：rdxとregCの積の上位64bitをregAに、下位64bitをregBに格納（符号なし128bit乗算）

なお、mulx命令の2つのデスティネーション（regAとregB）は同じにすることができます。 そのようにした場合、上位64bitが格納されます。 下位64bitだけが欲しい場合はimul命令を使えばいいのでこのような仕様は合理的です。

測定するときに使う命令列は次の通りです。 破壊的な命令が多いので、命令間の依存ができないようにしてスループットを測ります。 実際、多倍長乗算などの場合は乗算命令間には依存がないので、このような条件でのベンチマークは適切だと思います。

• imul(64)とmulx(128)の混合を調べる時
  • imul(64)：mov reg1, rdx; imul reg1, rsiとすることでrdxとrsiの積をregに書き込む
  • mulx(128)：mulx reg1, reg2, rsiとすることでrdxとrsiの積をreg1とreg2に書き込む
  • mulx(128)（128bit積の上位）：mulx reg, reg, rsiとすることでrdxとrsiの積の上位64bitをregに書き込む
• imul(64)とmul(128)やimul(128)の混合を調べる時
  • imul(64)：mov reg1, r8; imul reg1, rsiとすることでr8とrsiの積をregに書き込む
  • mul(128)：mov rax, r8; mul rsiとすることでr8とrsiの積をrdxとraxに書き込む
  • imul(128)：mov rax, r8; imul rsiとすることでr8とrsiの積をrdxとraxに書き込む
    • 必要に応じて、mov命令で結果を特定のレジスタに移動させる場合も検証

Golden Cove (2021)

実験環境：Intel Core i9 12900K (Alder Lake-S)（のPコア）、Eコア無効状態

いずれの命令もスループットは等しく1で、混合しても変わりませんでした。 Golden CoveにはMul演算器が、MulHi演算器が1つあるとされています。 したがって、全ての命令がMul演算器を使うため、このような結果になったのでしょう。 mulx regA, regA, regCのようにしても、MulHi演算器を使うマイクロ命令だけが実行されるというわけではなくMul演算器を使うマイクロ命令も実行されるのか、混合してもスループットは上がりませんでした。

Zen 4 (2022)

実験環境：AMD Ryzen Threadripper 7980X

いずれの命令もスループットは等しく1で、混合しても変わりませんでした。 Zen 4の整数乗算器の個数はよくわかりませんでしたが、実験結果から察するに1つしかないのでしょう。

Lion Cove (2024)

実験環境：Intel Core Ultra 9 285K (Arrow Lake-S)（のPコア）

Lion Cove は、乗算器が3つあるらしいです。

以下、ループ内の乗算命令とmov命令の個数と逆数スループットを書いていきます。なお、mulx命令の結果の全体が必要か上位のみが必要かで逆数スループットはほぼ変わりませんでした。

• 5×mulx(128)……5.009
• 4×mulx(128)……4.003
• 1×imul(64)＋4×mulx(128)……4.005
• 4×imul(64)＋4×mulx(128)……4.009
• 8×imul(64)＋4×mulx(128)……4.117
• 9×imul(64)＋3×mulx(128)……4.005
• 6×imul(64)＋3×mulx(128)……3.052
• 6×imul(64)＋3×mulx(128)＋1×mov……3.104
• 6×imul(64)＋3×mul(128)＋3×mov……3.207
• 6×imul(64)＋3×mul(128)＋6×mov……3.341
• 6×imul(64)＋3×imul(128)＋6×mov……3.327

これらの実験結果からわかることは、

• mulx命令のスループットは1である
• mulx命令とimul命令の合計スループットは1/3である
• mulx命令とimul命令の比が2:1のとき、スループットが低めになる。おそらく、最適でないポート振り分けが起こっている
• mov命令は、1つあたり逆数スループットが0.05くらい増える

あたりでしょうか。

ここから推察される発行ポートの構造は、

• ポートX：128bit乗算も64bit乗算もできる
• ポートY：64bit乗算ができる
• ポートZ：64bit乗算ができる

ポート振り分けの戦略は、

• mulx(128)・mul(128)・imul(128)：すべてポートXに行く
• imul(64)：ポートX・ポートY・ポートZ、のどれか空いていそうなところに行く

という感じになっていると思います。

応用例としてモンゴメリ乗算を考えると、mulx(128)が2回とimul(64)が1回でできるので、乗算器が1つしかない時代はスループットが1/3だったのに対して、Lion Coveではスループットが1/2に向上したということになりそうです。 このような上位部が必要な計算は劇的には高速化しませんが、整数配列全体を何倍かするみたいな下位部だけで済む計算は3倍速くなりそうです。

TODO

Zen 5、Gracemont, Crestmont, Skymontは乗算器が複数搭載されているらしいので、これらも調査する必要がありそうです。

メモ

addsd xmm0, xmm1, xmm1を10個含むループを実行すると、予想されるスループット2よりも高くなる場合が20%くらいありました。 2.667くらいまで行くこともまれにありました。 かなり奇妙な挙動ですが、原因は分かりませんでした。

*1:下64bitだけであれば、符号つきでも符号なしでも結果（のビット表現）は変わりません。で考えれば明らかです。

*2:まぎらわしいことに、imul命令には複数の異なる形式があります。64bit符号つき整数と64bit符号つき整数の積を128bit符号つき整数として求める命令もimulです。

効率的な多倍長演算をC言語の範囲で書きたい（clang向け）

x86命令セットには、adc（キャリーつき加算）命令やsbb（ボローつき減算）命令があり、多倍長演算の際にはこれを使うと効率的です。 一方で、コンパイラにこれらの命令を出力してもらうのは、C言語の範囲での記述ではなかなか簡単ではありません。 なお、ここでは、「C言語の範囲」というのは、以下の要件を満たすことにします。

• __uint128_tや__int128_tは、加減算をadc命令やsbb命令にコンパイルするためには使わない
  • いわゆるgcc拡張
  • 二倍長乗算命令（x86のimul命令やRISC-Vのmulh命令など）を出力するのは、現状のコンパイラではどうやっても無理っぽいので、そのためにのみ使う
  • __int128_tや__uint128_tがない環境への移植コストの最小化が目的。二倍長乗算命令をエミュレートできることはすでに示した（64bit乗算の上位部分の計算 - よーる）ので、必要であればこれで差し替えればいいはず。一方で加減算は結構いろいろな落とし穴があるので非標準な拡張機能を使ったブラックボックスにしたくない
  • 192bitの計算がしたくなったりしたときは結局自分で書かないといけない、つまりキャリーやボローを取り扱う手法を勉強する必要がある
• _BitIntも使わない
  • BigIntかと見間違えた
  • C23で標準化された任意ビット数整数
  • 最大幅BITINT_MAXWIDTHが大きいかは処理系定義
  • 新しいので未対応な処理系も多い
  • C++は導入議論中みたい
• int64_tとuint64_tはあるとする
  • 実装されるかは処理系定義
  • x86_64向けの処理系でこれが定義されないということはないはず
  • これが定義されない処理系への移植では、int_least64_tやuint_least64_tをいちいちマスクしてエミュレートすればいいはず
  • _BitInt(64)は常に使えるので、それが最も移植性が高い……のかな？

また、x86向けのコンパイラは、かしこい最適化コンパイラであるclangを前提にします。

多語と一語の積（筆算の一行分）

以下のように書けば、意図通りキャリーつき加算になります。

void mul_256_64_320( uint64_t* A, uint64_t B, uint64_t* C) {
    __uint128_t p0 = (__uint128_t)A[0] * B;
    __uint128_t p1 = (__uint128_t)A[1] * B;
    __uint128_t p2 = (__uint128_t)A[2] * B;
    __uint128_t p3 = (__uint128_t)A[3] * B;
    uint64_t p0hi = p0 >> 64;
    uint64_t p0lo = p0;
    uint64_t p1hi = p1 >> 64;
    uint64_t p1lo = p1;
    uint64_t p2hi = p2 >> 64;
    uint64_t p2lo = p2;
    uint64_t p3hi = p3 >> 64;
    uint64_t p3lo = p3;

    uint64_t t4 =        p3lo;
    uint64_t t3 = p3hi + p2lo; p2hi += t3 < p3hi;
    uint64_t t2 = p2hi + p1lo; p1hi += t2 < p2hi;
    uint64_t t1 = p1hi + p0lo; p0hi += t1 < p1hi;
    uint64_t t0 = p0hi       ;

    C[0] = t0;
    C[1] = t1;
    C[2] = t2;
    C[3] = t3;
    C[4] = t4;
}

まず、z = x + yみたいな加算でキャリーが生じるかは、z < xで判定できます。 キャリーが生じるというのは要するに符号なしオーバーフローが発生するということです。 オーバーフローの発生は、「足したのに元の数よりも小さくなった」という事象と同値なので、これで判定できるわけです。 もちろん、加算は対称なので、z < yでも同じです。

さて、このコードで重要なのは、乗算結果の上位部にキャリーを足している（t2 = p2hi + p1lo; p1hi += t2 < p2hi）ことです。 この部分は、t1 < p2hi + (carry)はキャリー伝搬とみなされるのに対し、t2 < p1loはキャリー伝搬とみなされない点を考慮に入れたコードになっています。 「加算は対称なのでどちらでも同じ」と言っていましたが、これはキャリー生成の時だけのようで、キャリー伝搬の時は非対称のようです。 つまり、交換則は成り立っても結合則が成り立たないようです。 なお、下位部を最後以外で足す順番（c2 + p0lo + p1hiやp0lo + p1hi + c2）もやはりキャリー伝搬とみなされません。

一応、t1 = p1hi + c2 + p0lo; uint64_t c1 = t1 < p1hi + c2;のように書けば、変数の破壊的変更自体はなくせます。 ただ、同じ式が二回現れたり、比較の片方が加算結果になっていたりで、あまりきれいな形ではありません。 なお、変数の破壊的変更を許容する場合でも、乗算結果の下位にキャリーを足しこむ方式は、その時点でオーバーフローする可能性があってさらに記述が複雑化するのでうまくいきません。

多語と多語の和（多倍長加算）

乗算の場合は「乗算結果の上位は、1足してもオーバーフローしない」という性質をうまく利用することでうまくいったような雰囲気がありました。 したがって、多倍長加算はそれより難しいのかと思いましたが、意外と簡単に最適な命令列を出力させることができました。 コードを書いている時は気づきませんでしたが、加算の場合は「二つのuint64_tを足した結果」が「オーバーフローしたなら、さらに1足したときにオーバーフローしない数」になっている点をうまく活用するようになっているのかもしれません。

void add_256_256_256( uint64_t* A, uint64_t* B, uint64_t* C) {
    uint64_t a3 = A[3];
    uint64_t b3 = B[3];
    uint64_t a2 = A[2];
    uint64_t b2 = B[2];
    uint64_t a1 = A[1];
    uint64_t b1 = B[1];
    uint64_t a0 = A[0];
    uint64_t b0 = B[0];

    uint64_t t3 = a3 + b3;
    uint64_t s2 = a2 + b2;
    uint64_t t2 = s2 + (t3 < a3);
    uint64_t s1 = a1 + b1;
    uint64_t t1 = s1 + (t2 < s2 | s2 < a2);
    uint64_t s0 = a0 + b0;
    uint64_t t0 = s0 + (t1 < s1 | s1 < a1);

    C[0] = t0;
    C[1] = t1;
    C[2] = t2;
    C[3] = t3;
}

キャリーの部分は、OR（|）でも和（+）でも計算結果は同じになります。 これは、t2 < s2が成り立つときはs2 == 0xffffffffffffffffかつ(t3 < a3) == 1でt2 == 0になる場合のみであり、この時にs2 < a2は絶対成り立たないので、1OR1になることがないからです。 99999+99999+1=199999みたいなのを思い出せば、キャリー入力があったとしても、ある位から出力されるキャリーは高々1であり、上の位に2を足さないといけないパターンはあり得ない、という説明も可能です。

しかしながら、ORでないとclangはキャリー伝搬であると認識しません。 ORだとコードを局所的に眺めただけで上記推論ができますが、和で書かれていると大局的な推論の連鎖が必要であきらめてしまうということでしょうか。

変数の破壊的変更を許容する場合、以下のような短い書き方ができますが、まぁあまりやるべきではないタイプの書き方でしょう。

    uint64_t t3 = (b3 += a3);
    uint64_t t2 = (b2 += a2) + (b3 < a3);
    uint64_t t1 = (b1 += a1) + (b2 < a2 | t2 < b2);
    uint64_t t0 = (b0 += a0) + (b1 < a1 | t1 < b1);

多倍長加算の方式を利用すると、筆算一行分に相当する乗算は以下のように書けますが、これも特別見やすいというわけでもなさそうです。

    uint64_t t4 = p3lo;
    uint64_t t3 = p2lo + p3hi;
    uint64_t s2 = p1lo + p2hi;
    uint64_t t2 = s2 + (t3 < p2lo);
    uint64_t s1 = p0lo + p1hi;
    uint64_t t1 = s1 + (t2 < s2 | s2 < p1lo);
    uint64_t t0 = p0hi + (t1 < s1 | s1 < p0lo);

多語同士の乗算（多倍長乗算）

多倍長整数同士の乗算を効率よく記述する方法として、adcx命令とadox命令を使う方法が知られています。

zenn.dev

adcx命令は、キャリーつき加算を行いつつも他のフラグ、特にOF（オーバーフローフラグ）を書き換えない命令です。 一方、adox命令は、OFをキャリー入力として扱うキャリーつき加算を行い、生じたキャリーをOFに出力し、他のフラグは書き換えない命令です。 「他のフラグは書き換えない」みたいな命令はいかにもパーシャルレジスタの呪いの影響を受けそうですが、実はうまい設計になっているようです。 EFLAGSレジスタは、実質的にCFとそれ以外のフラグを分けて管理されているようです（最大公約数をもっと高速に求める（その2）【cmova命令は遅い】 - よーるで見たように、CFとそれ以外のフラグ（ZF）を両方読むcmova命令が遅くなるあたりから透けて見えます。EコアやZen系のコアではそのような性能低下は見られないので、「当時はそうだった」というのが適切かもしれません）。 そこで、adcx命令はCFだけ読み書きする、adox命令は「それ以外のフラグ」を全部読んでOFだけ書き換える、という実装が可能で、うまくやりたいことの実現と実装しやすさを両立させています。

この二つの命令は、Broadwell世代（後述のmulx命令が導入されたHaswellの次の世代。シュリンクするTick世代なのに命令が増えている……）で導入されたものですが、その価値は、キャリー用に使える“レジスタ”を二つ提供することにあります。 キャリー用の“レジスタ”が二つあれば、mulx命令（64bit×64bitを128bitで計算する、フラグを書き換えない命令）の出力として得られる二つの64bitレジスタをアキュムレータに即座に足しこんで、それらの64bitレジスタを解放することができます。 一方でキャリー用の“レジスタ”が一つしかなければ、最上位桁まで足しあげた後、また下の位に戻って足しあげる、という順番でやらないと、命令数が増えてしまいます。 そのような方法だと、乗算結果の片方を汎用64bitレジスタに保持しておく必要があります（メモリに追い出すならその時点で命令数が増えてしまいます）。 汎用64bitレジスタの数には限りがあるため、長い多倍長整数同士の乗算だとこの方法が取れなくて、キャリーをいったん汎用レジスタに出力する、みたいなことが必要になります。 なので、キャリー用の“レジスタ”が二つ提供されることが重要で、逆に三つ以上は特に必要ないのです。

さて、clangはadox命令を出力することはおそらくありません。 adcx命令も、adox命令と組み合わせていいかんじの多倍長乗算を実装するといった本来の目的に使用することはなさそうです（もしかするとフラグ保持の兼ね合いで出力することがなくはないのかもしれません）。

とりあえず、それほど長くない多倍長同士の（または、片方が長くない多倍長との）乗算であれば、上の二つを組み合わせることでおおよそ最適な命令列が生成されます。

void mul_192_192_384( uint64_t* A, uint64_t* B, uint64_t* C) {
    uint64_t pp[3][6] = {};
    for( int i = 0; i < 3; ++i ) {
        uint64_t hi = 0;
        for( int j = 2; j >= 0; --j ) {
            __uint128_t pij = (__uint128_t)B[j] * A[i];
            uint64_t lo = pij;
            pp[i][i+j+1] = hi + lo;
            hi = (pij >> 64) + (hi + lo < hi);
        }
        pp[i][i] = hi;
    }

    uint64_t sum[6] = {};
    for( int i = 2; i >= 0; --i ) {
        uint64_t carry = 0;
        for( int k = 5; k >= 0; --k ) {
            uint64_t t = pp[i][k];
            uint64_t u = t + sum[k];
            uint64_t v = u + carry;
            carry = (v < u | u < t);
            sum[k] = v;
        }
    }

    for( int k = 5; k >= 0; --k ) {
        C[k] = sum[k];
    }
}

これは次の命令列になったようです。

mul_192_192_384:
        push    rbp
        push    r15
        push    r14
        push    r13
        push    r12
        push    rbx
        mov     rax, rdx
        mov     r11, qword ptr [rsi + 16]
        mov     rdx, qword ptr [rdi]
        mov     rcx, qword ptr [rdi + 8]
        mulx    r10, r8, r11
        mov     rbx, qword ptr [rsi]
        mov     r14, qword ptr [rsi + 8]
        mulx    r15, r9, r14
        add     r9, r10
        mulx    rsi, r10, rbx
        mov     rdx, rcx
        mulx    r12, r13, r11
        adc     r10, r15
        adc     rsi, 0
        mulx    r15, rbp, r14
        add     rbp, r12
        mulx    rcx, r12, rbx
        adc     r12, r15
        adc     rcx, 0
        mov     rdx, qword ptr [rdi + 16]
        mulx    r15, rdi, r11
        mulx    r14, r11, r14
        add     r11, r15
        mulx    rbx, rdx, rbx
        adc     rdx, r14
        adc     rbx, 0
        add     r11, r13
        adc     rdx, rbp
        adc     rbx, r12
        adc     rcx, 0
        setb    bpl
        add     rdx, r8
        adc     rbx, r9
        adc     rcx, r10
        setb    r8b
        movzx   r8d, r8b
        add     bpl, 255
        adc     r8, rsi
        mov     qword ptr [rax], r8
        mov     qword ptr [rax + 8], rcx
        mov     qword ptr [rax + 16], rbx
        mov     qword ptr [rax + 24], rdx
        mov     qword ptr [rax + 32], r11
        mov     qword ptr [rax + 40], rdi
        pop     rbx
        pop     r12
        pop     r13
        pop     r14
        pop     r15
        pop     rbp
        ret

adc命令を二回やればいいだけなのに、キャリーをsetbで汎用レジスタに出してから、movzxしてから足したり、add 255でCFに戻したりしていて最適ではありませんが、頑張った方だと思います。 うまくいっていないのは最上位部分で、他の部分と形が異なる r64 + CF + CF みたいな形になっています。 おそらく、こういった形に対するパターンマッチがなくて万能な方法でやっているのでしょう。

なお、ソースコード上で変な書き方をしている部分とその理由は以下の通りです。

• 素直じゃないB[j] * A[i]の順番。これは、乗算の右オペランドの方がrdx（mulxの暗黙オペランド）に入りがちで、ループ内ではそれが固定の方がmov命令が減るからです。逆順だと、mulxのたびに毎回mov rdx, XXXみたいな命令が入ってしまいます。
• iのループが昇順だったり降順だったりする。これは、p22loに相当する部分（C[5]にストアされる）がなぜかスピルされてしまう問題を解決するためです。早くC[5]にストアすればいいのに、なぜか一回スタックにスピルしてから戻してストアする、みたいな命令列が出がちです。

これ以上長い乗算は、スピルなしでは難しいです。

ボローつきの符号付き引き算

以下の引き算を行いたいとします。

void sub_s128_s64( int64_t ah, uint64_t al, int64_t b, int64_t* ch, uint64_t* cl ) {
    __int128_t a = (__int128_t)ah << 64 | al;
    __int128_t c = a - b;
    *ch = c >> 64;
    *cl = c;
}

これをx86_64向けにコンパイルすると、次のような命令列になります。

sub_s128_s64:
        mov     rax, rdx
        sar     rax, 63
        sub     rsi, rdx
        sbb     rdi, rax
        mov     qword ptr [rcx], rdi
        mov     qword ptr [r8], rsi
        ret

これを__int128_tに頼らずに実装しようと思います。

void sub_s128_s64( int64_t ah, uint64_t al, int64_t b, int64_t* ch, uint64_t* cl ) {
    *cl = al - b;
    *ch = ah - (*cl > al) + (b < 0);
}

上のように実装しましたが、以下のように余分な命令が出てしまいます。 三項の足し引きの順番を変えてみても、うまくいく気配がありませんでした。 なお、ボローが出る条件である(*cl > al)（引いたはずなのに大きくなった）は(al < b)（小さい数から大きい数を引く）でも同じ結果になりました。 また、+ (b < 0)は- (b >> 63)や+ ((uint64_t)b >> 63)でも同じになりました（b >> 63は処理系定義動作で算術シフトになることを意図しています）。 この部分は、「bが負の時、bを（暗黙変換で）uint64_tにすると意図している値よりも大きくなるので、全体として引きすぎになる」のを補正しています。 - (b >> 63)で考えたほうがより直観的で、「bを符号拡張した時の上位部（bh相当）をahから引く」という操作に対応しています。

※ボローが出るというのは、x86においてはCFが1になる状況を意味します。記事の末尾に簡単な説明を付けました。

sub_s128_s64:
        mov     rax, rsi
        mov     rcx, rdx
        shr     rcx, 63
        add     rcx, rdi
        sub     rax, rdx
        sbb     rcx, 0
        mov     rdx, rcx
        ret

これは、キャリー/ボローに推論されないというよりは、(b < 0)を足す部分の問題のようです。

最終演算はボロー付きの命令（ボローが出ていたら1減る命令）にならないといけなくて、x86の範囲だとそのような命令はsbb命令しかありません。 したがって、最終演算は引き算にしたいです。 そこで、sar命令で-1/0を作ってそれを引くのがうまい方法になります。 実際、__int128_tを使った場合はそのような命令列にコンパイルされています。 しかしながら、+ (b < 0)や- (b >> 63)や+ ((uint64_t)b >> 63)はいずれも、clangが気を利かせたのか、shr命令で作った1/0の加算に変えられてしまうようです。 加算の方がなんとなくわかりやすいような気もしますし、コンパイラ制作者も混乱しないためにそのような正規化を行っているのかなと思いますが、残念ながら今回はそれが裏目に出てしまっています。

C言語の範囲で符号付き多倍長整数の演算を書くのは、clangの最適化能力をもってしてもなかなか難しいようです。 x86には符号なし整数と符号付き整数の積（の上位）を求める命令がなく（RISC-VにはMULHSU命令がある）、符号付き多倍長整数の乗算が難しいこともありますし、x86では素直に符号なしでやったほうがよさそうです。

二倍長乗算の最適化

キャリーつき加算を実現する方法が分かったので、__uint128_tを全く使わない二倍長乗算のフォールバック（64bit乗算の上位部分の計算 - よーる）を改善してみました。 実はキャリーつき加算はあまり関係なく、そもそも足してオーバーフローしない部分があるようです（コンパイラの出力を見て気づきました）。 かなりいいかんじに最適化されています。

uint64_t MULHUU( uint64_t x, uint64_t y ) {
    uint64_t xl = x & 0xffffffff;
    uint64_t xh = x >> 32;
    uint64_t yl = y & 0xffffffff;
    uint64_t yh = y >> 32;

    uint64_t z0 = xl * yl;
    uint64_t z1 = xl * yh;
    uint64_t z2 = xh * yl;
    uint64_t z3 = xh * yh;

    uint64_t z0h = z0 >> 32;
    uint64_t s = z0h + z1; // no overflow
    uint64_t t = s + z2;
    return z3 + (t >> 32) + ((uint64_t)(t < s) << 32);
}

ここで、z1は最大でもであるため、以下の数であるz0hを加えてもオーバーフローしません。 あるいは、96ビットあるxl * yの上位64ビット分であると考えても、オーバーフローしないことがわかります。

この関数は、以下のようにコンパイルされます。無駄な命令は一つもなさそうです。 movのうち2つは零拡張用に使われているので意味ある命令です。 もう1つのmovも消せそうにありません。 というのも、2×2で4回の乗算をやる必要がありますが、最初の1回だけは両オペランドともまだ後で使うので、破壊に備えてコピーが必要だからです（残りの3つの乗算は以降で使わないオペランドを破壊する形で実行されていることもわかります）。 32bit版の二倍長乗算命令を使えば>>32を一個消せる気がしますが、使用頻度の低そうな変な命令なので遅いのかもしれません（実際、uops.infoによれば、AlderLake-Pだとマイクロ命令数が増えるようです）。

MULHUU:
        mov     eax, edi
        shr     rdi, 32
        mov     ecx, esi
        shr     rsi, 32
        mov     rdx, rcx
        imul    rdx, rax
        imul    rax, rsi
        imul    rcx, rdi
        imul    rsi, rdi
        shr     rdx, 32
        add     rdx, rax
        xor     eax, eax
        add     rdx, rcx
        setb    al
        shr     rdx, 32
        add     rdx, rsi
        shl     rax, 32
        add     rax, rdx
        ret

xorでゼロ埋めしたあとsetbでキャリーを拾っているようです。xorでゼロ埋めした後adcでもいいですが、今回の目的ではこれで問題ないですね。 なんかたまにshld（ファネルシフト）命令にコンパイルされることがあるようですが、きまぐれでしょうか。

まとめ

なるべく移植性の高い方法で多倍長演算を記述しつつ、特定のコンパイラと特定の命令セット（clangでx86向け）を仮定した時に最適な命令列が出るような記述方法を模索しました。

おまけ：キャリーとボローの関係

キャリーは、次の桁への繰り上がり（0または1）です。 次の桁を計算する時、通常の加算はa + bですが、キャリーつき加算はa + b + carryとなります。 これにより、キャリーがあるときは1多く加算することができて、キャリーがない時は通常と同じ加算になって、いずれもOKです。

ボローは、次の桁からの繰り下がり（0または1）です。 次の桁を計算する時、通常の減算はa - bですが、ボローつき減算はa - b - borrowとなります。 これにより、ボローがあるときは1多く減算することができて、ボローがない時は通常と同じ減算になって、いずれもOKです。

これを二進加算器で実装することを考えます。 まず、キャリーつき加算の方は、二進加算器の最下位桁の余っているキャリー入力をキャリーフラグ読み出し口につなぐ＆二進加算器の最上位桁のキャリー出力をキャリーフラグ書き込み口につなぐ、とすれば実現できることはいいでしょう。 問題は、ボローにどう対応するかです。 まず一つの問題は、ボローフラグを作るかです。 これを作ってしまえばわかりやすく作ることができますが、キャリーフラグとボローフラグを同時に使う機会があるとは思えず、できれば共有したいです。 もう一つの問題は、どのように「追加で1引く」を実現するかです。 二進加算器には、「追加で1足す」を実現する入力信号（キャリー入力）はありますが、「追加で1引く」を実現できそうな入力信号は見当たりません。

この二つを同時に解決する、全くうまい方法があります。 二進加算器の性質をよく考えればある意味当たり前なのですが、命令セット定義と回路設計の境目みたいな話なのであまり考えたことがないかもしれません。

まず、引き算を二進加算器でどう実現していたかを思い出します。 a - bというのは、実はa + ~b + 1を実行していたのでした（ここで1はキャリー入力から入れています）。 したがって、a - bを実行する時、「追加で1を引く」は、このキャリー入力から1を入れているのをやめれば実現できることがわかります。 つまり、ボローが出ていれば二進加算器のキャリー入力に0を、ボローが出ていなければ（通常通り）1を、それぞれ入れたうえで、aと~bの加算を実行すればよいことになります。

次に、生成されるボローについて考えます。 ボローが出る時、つまりa - bの結果が符号なしオーバーフローしてaよりも大きくなる時、二進加算器のキャリー出力からは0が出ます。 逆にボローが出ない時、つまりa - bの結果がa以下になるとき、二進加算器のキャリー出力からは1が出ます。

この二つを合わせてみると、ボローが出るというのは「二進加算器のキャリー出力から0が出てきている＆次の桁では二進加算器のキャリー入力に0を入れたい」という状況であり、逆にボローが出ないというのは、「二進加算器のキャリー出力から1が出てきていて、次の桁では二進加算器のキャリー入力に1を入れたい」という状況です。 簡潔に言い換えると、ボローが出るというのは0を次回に受け渡したい、ボローが出ないというのは1を次回に受け渡したい、という状況です。

このことを考えると、ボローフラグを独立に作る必要はないことがわかります。 なぜなら、次回に受け渡したい値を保存するという、まさにそのものの役割を果たすキャリーフラグが既にあるからです。 この時、キャリーフラグが0＝ボローが出ている、キャリーフラグが1＝ボローが出ていない、という対応付けだと解釈することになります。 この方式では、追加回路は全く不要です。

ただ、このような極性の反転したフラグ方式は、デバッグ時に直観的に理解するのが困難という問題があります。 したがって、x86などではこのような方式はとられておらず、キャリーが出たときにCF=1、ボローが出た時もCF=1になります。 一方で、マイコンなどの回路を切り詰めたい用途向けの命令セットでは、極性反転方式が使われている印象です。 身近（？）なところでは、ARMのCフラグが、carry/no-borrowの時に1が立つという極性反転方式です。 とりあえず、今となっては、こんな細かいところでケチれるトランジスタ数などたかが知れており、二大派閥のどちらに属しているのかを覚える手間だけが余計に増えてややこしくなっているという感じがします……。