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命令を使う方法が知られています。
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が立つという極性反転方式です。 とりあえず、今となっては、こんな細かいところでケチれるトランジスタ数などたかが知れており、二大派閥のどちらに属しているのかを覚える手間だけが余計に増えてややこしくなっているという感じがします……。
