2026/04/11 18:03
マークの魔法倍増
RSS: https://news.ycombinator.com/rss
要約▶
Japanese 翻訳:
記事は、低消費電力のプロセッサ（完全な浮動小数点ユニット (FPU) を備えていないもの）において単精度浮動小数点乗算を高速化するための新しい RISC-V 拡張機能である Xh3sfx の導入について述べる。このアプローチは、遅いソフトウェアルーチンに依存するのではなく、特殊なハードウェアロジックを活用するものである。設計では、ARM アーキテクチャで以前に使用されていたビット操作のトリック（具体的には Mark Owen による最適化）に適応し、RISC-V 固有の命令と組み合わせて、乗算のレイテンシを 30 サイクル に達成している。これは、標準的な「教科書方式」の方法に必要な 33 サイクルよりも著しく速いものである。
現在の拡張機能は乗算に焦点を当てているが、14 サイクル以内に加算処理も行う。この設計は RP2350 チップ上で成功裏にタペアウトされており、Cortex-M0+ マイクロコントローラーのような微小デバイスにおける実用上の妥当性を示した。このアプローチは、完全なハードウェア FPU や純粋なソフトウェアエミュレーションと比較して、シリコンエリアを削減し、プロセッササイクル数を節約する。研究者たちは現在、同様の戦略を将来の倍精度計算へ拡張できるかどうかを探っている。
本文
この投稿の主題は、私が特に親しみ深かったり熱中したりしている分野です。つまり、組み込みプロセッサにおける単精度浮動小数点乗算の話です。
まず、なぜ最近このトピックにこれほど没頭してきたのかという背景から始め、私が独自に考案した実装手法を解説し、最後に 32 ビット組み込みコア向けの浮動小数点乗算に関するマーク・オーウェン（Mark Owen）氏の極めて大胆なトリックについて深掘りします。このトリックこそが、今回の投稿のそもそものインスピレーション源です。
⚠️ 注記：この記事には「浮動小数点」が含まれています。 カリフォルニア州当局によれば、浮動小数点は哺乳類二足動物（人間）に混乱と恐怖反応を引き起こすと知られています。標準的な推奨事項としては「Computer Scientists Who Should Know About Floating-Point Arithmetic Every」（『コンピューターサイエンティストが知っておくべき浮動小数点算術』）があります。また、実際の IEEE 754-2008 規格は、基数（radix）が 2 でないというファンタジー小説を無視すれば、意外に簡潔かつ読みやすいものです。より触覚的な体験をお探しであれば、IEEE 754 カルキュレーターでバイナリ値（binary16 から始めるのが良いでしょう）をポコポコ弄ってみてください。
ハードでもソフトでもない中間領域
最近、軟質浮動小数点ルーチンの高速化を目指したカスタム RISC-V 拡張命令セット「Xh3sfx」の開発に取り組んでいます。これは FPU を持つケースと持たないケースのちょうど中間地点にあたるもので、未だ十分に探求されていない空間だと感じています。「ファーム浮動小数点（firm floating point）」と呼ぶかもしれません。
C プログラムを、浮動小数点ハードウェアを持たないターゲット用にコンパイルする際、コンパイラは float 変数に対してリクエストされた操作を実行するために、
libgcccompiler-rtXh3sfx がカスタム拡張命令セットとはいえ、フォークしたコンパイラーを維持・配布するために名乗りを挙げるわけではありません。コンパイラの実行ルーチンを加速版のものに置換する方が簡単です。新しいルーチンは、浮動小数点フォーマットの難解で不細工な部分を処理するために数種の専用 ALU 操作を使い、実際の計算には通常の整数命令を混在させています。ランタイムライブラリは大部分が文書化され安定した API を有しており、プログラムへの追加支持は加速ライブラリンクやビルドソースへの配置のみで実現できます。これは組み込みファームウェアにとっては合理的なアプローチです。
Xh3sfx は数百ゲートのわずかなコストに対して、単精度の足し算を 14 サイクル、乗算を 16 サイクル（関数呼び出しオーバーヘッドを無視）を実現します。（他にもできることはありますが、これらは例にすぎません。）質的な面で言えば、これは浮動小数点を「ああ、神さま、なんでこんなに遅いのよ」という状況から、一般的なアプリケーションコードや軽量オーディオ DSP では「ただ動作する（Just Works™）」状態へと転換させます。元々は Mastodon でこれを投稿しましたが、命令セットの詳細についてはここ、ライブラリルーチンについてはこっちをご覧ください。
Xh3sfx を用いた乗算
Xh3sfx ライブラリのデフォルト単精度乗算実装は以下の手順を踏みます：
- 二つの浮動小数点入力から指数と有効桁（ significand ）をパッキング解除する。
- 正確な有効桁積を、32 ビット×32 ビット→64 ビットの乗算（
,mul
）を用いて計算する。mulh - 積を 32 ビットレジスタに圧縮し、丸め方向を保持した状態で格納する。
- 積を正規化する。
- 元の指数の和と組込みを行い、最終的な浮動小数点結果を得る。
具体的にはこのようになります：
__mulsf3: // 指数をパッキング解除 h3.bextmi a2, a0, 23, 8 h3.bextmi a3, a1, 23, 8 // 特殊ケース：入力片方が 0 / 劣化数/NaN/無限大である場合 h3.fcheck2e.s t0, a2, a3 bnez t0, __mulsf_special_exponent // 単精度有効桁を Q3.29 としてパッキング解除 h3.funpackq3.s a4, a0 h3.funpackq3.s a5, a1 // Q3.29 × Q3.29 → Q6.58 を計算 mul a1, a4, a5 mulh a4, a4, a5 // 下部小数節のスティッキービットを収集; a4 は実質的に Q6.26 となる snez a1, a1 or a4, a4, a1 // 新しい指数を計算； // -127 は冗長なバイアス用、+3 は Q6→Q3 のオフセット用 add a2, a2, a3 add a2, a2, a5 addi a2, a2, -127 + 3 // 丸め付き単精度パッキング h3.fpackrq3.s a0, a4, a2 ret
h3.*h3.bextmiこの実装は
mulmulhコンventional（学校教則）乗算
Hazard3 では、乗算用のハードウェア構成が 3 つ用意されています：
- Minimal（最小）：すべての除算・乗算命令はシークウェンシャルな乗算/除算回路で実行されます。通常、1 サイクルまたは 2 サイクルごとのビット処理に設定され、符号補正のために数サイクル追加がかかります。
- Intermediate（中間）：32×32→32 ビットの乗算は専用高速乗算器で実行され、残りの
はシークウェンシャル回路で実行されます。mulh/mulhu/mulhsu - Full（フル）：すべての乗算が専用高速の 32×32→64 ビット乗算器で実行されます。除算は引き続きシークウェンシャル回路で実行されます。
これらのオプションは、パフォーマンスと面積コストの両方において昇順に並べられています。RP2350 に Tape-out されたのはフル構成の 32×32→64 ビット乗算です。最小構成（シークウェンシャル乗算/除算のみ）の場合、前述の
mul; mulhmulこの構成に向けた最初の最適化試みでは、32×32→64 ビットの乗算を四個の 16×16→32 ビット乗算に分けることに成功しました。これらはそれぞれ
mul__mulsf3: // 指数のパッキング解除とチェック h3.bextmi a2, a0, 23, 8 h3.bextmi a3, a1, 23, 8 h3.fcheck2e.s t0, a2, a3 bnez t0, __mulsf_special_exponent // 乗算を容易化するために有効桁を符号なしとしてパッキング解除 h3.funpacku3.s a4, a0 h3.funpacku3.s a5, a1 // 符号と指数の計算；a2/a3 を解放する xor a0, a0, a1 add a1, a3, a2 addi a1, a1, -127 + 3 // U3.29 × U3.29 → U6.58 を計算 (32x32 -> 64) srli a2, a4, 16 srli a3, a5, 16 zext.h a4, a4 zext.h a5, a5 mul t0, a2, a3 // A_h * B_h mul t1, a2, a5 // A_h * B_l mul a2, a4, a3 // A_l * B_h mul a4, a4, a5 // A_l * B_l // 交差項の和を計算；hi[16] へのキャリー保存 add t1, a2, t1 sltu t2, t1, a2 // 交差項を加算して低単語へ；hi[0] へのキャリー保存 slli a2, t1, 16 add a4, a4, a2 sltu a3, a4, a2 // キャリーを正しい重み付けでパッキング pack a3, a3, t2 // 高単語：A_h * B_h + (交差項 >> 16) + ci0 + (ci1 << 16) srli t1, t1, 16 add a5, t0, t1 add a5, a5, a3 // 最終積完了：U6.58 が a5:a4 に格納。下部スティッキービットを収集。 snez a4, a4 or a5, a5, a4 // Q3 へ再正規化 h3.feadju3 a4, a5 h3.ssra a5, a5, a4 add a1, a1, a4 // パッキング前に符号を適用。有効なのは上部 U3.29 のみ。 h3.xorsign a5, a0, a5 h3.fpackrq3.s a0, a5, a1 ret
これは馴染みの深い二項積です：$(x_1 + x_0)(y_1 + y_0) = x_1y_1 + x_1y_0 + x_0y_1 + x_0y_0$。$A = a_1 \times 2^{16} + a_0$（ここで $a_1, a_0$ は 32 ビット A の 16 ビット半分）と置けば、積 $AB$ は四個の 16×16→32 ビット乗算の結果の和にいくつかの $2^n$ （即ちシフト）因子を伴って展開されます。
多くの時間は乗算そのものではなく、ビットを正しい場所へルーティングする処理で消費されています。特にキャリーの伝播は非常に複雑です。RISC-V で最も過小評価されている命令
packmul; mulh上記コードをさらに改善するのは困難です。もし試みようとすれば、文献は乗算回数を減らすことに没頭していることに直面するでしょうが、短いパイプラインにおいては
mulこの実装に満足していましたら、マーク・オーウェン氏が無作為にメールを持ってきてくれたこのトリック（Armv6-M）に出会いました。
lsls r0,#8 @ x のマンティサを左シフト lsls r1,#8 @ y のマンティサを左シフト lsrs r0,#9 // 9 ビット右シフト（符号なし仮定） lsrs r1,#9 adds r2,r0,r1 // 後続用 mov r12,r2 lsrs r2,r0,#7 @ x[22..7] Q16 lsrs r3,r1,#7 @ y[22..7] Q16 muls r2,r2,r3 @ 結果 [45..14] Q32:決して過大評価ではなく、最悪の誤差は // 2*(2^7-1)*(2^23-2^7)+(2^7-1)^2 = 2130690049 < 2^31 muls r0,r0,r1 @ 結果 [31..0] Q46 lsrs r2,#18 // 結果 [45..32] Q14 bcc 1f cmp r0,#0 bmi 1f adds r2,#1 // r2 の誤差を修正 1: lsls r3,r0,#9 @ 結果の下部からビットをずらす lsrs r0,#23 @ Q23 lsls r2,#9 adds r0,r2 // cut'n'shut（結合） add r0,r12 @ 欠落した 1 の挿入なしに対する補償 (x+y) * 1 ;@ 結果 - 1 が r3:r0 に格納 (Q23+32、つまり [0,3) 範囲)
マークは RP2040 ROM 浮動小数点ライブラリの著者であり、彼のメールは常に 80 文字で折り返されています。この関数は単精度乗算の正しく丸められた結果を返しますが、2 つの 32×32→32 ビットの部分積だけで実現されます。また、私のコンベンショナルな乗算に比べて、命令セットが格段に制限されているにもかかわらず、ビット操作や一般論的な議論（waffling）は遥かに少なく済んでいます。
核心となるトリックは、二回の乗算を使用して 23×23→46 ビットの積を計算することです。このトリックは 24 ビットの入力には適用できないため、彼は暗黙の 1 を除外し、単に 23 ビットの小数部分のみを乗算し、後に欠落した分の補償を行うのです。
まず第一段階：
muls r0,r0,r1 // 結果 [31..0] Q46
この下部の乗算が直接、結果の下位 32 ビットを与えてくれます。この部分は正確です。
第二段階：
lsrs r2,r0,#7 @ x[22..7] Q16 lsrs r3,r1,#7 @ y[22..7] Q16 muls r2,r2,r3 // 結果 [45..14] Q32:決して過大評価ではなく、最悪の誤差は // 2*(2^7-1)*(2^23-2^7)+(2^7-1)^2 = 2130690049 < 2^31
この上部の乗算は、我々の結果に近い数の上位ビット（45:14）を与えます。これらのビットに対する正確な結果は $(x * y) >> 14$（中間値の切断なし）となるべきですが、代わりに $(x >> 7) * (y >> 7)$ を計算しており、これは近似的で、各側面の 7 つの最低位ビット（LSB）からの寄与を欠いています。
もし高乗算から得られる上位 45:32 ビットを下乗算から得られる低位 31:0 ビットに接合すれば、完全な 46 ビットの結果が得られます。ただし、現時点ではこれは単なる近似的な結果（
approx[45:14]- 結果の上位ビット 45:31 は既に正しい。
- 結果の上位ビット 43:31 が間違っているが、$2^{31}$ を加えることで修正できる。
この修正は 30:14 ビットを修復しませんが、それは構いません：下乗算からはすでにそれらを得ているからです。修正論理は以下の通りです：
lsrs r2,#18 // 結果 [45..32] Q14 bcc 1f cmp r0,#0 bmi 1f adds r2,#1 // r2 の誤差を修正 1:
英語に訳せば、上記のコードは「
approx[31]exact[31]approx[32]approx[31:14]approx[31]exact[31]approx[31]残りのコードは、後の丸め処理のために便利な形にこれらのビットをシフトして結合し、欠落した暗黙の 1 の補償も部分的に行います。
マークにこれを尋ねた際、彼は以下のような説明を提供してくれました（これはエレガントですが、私の脳には少し形状が合いません）： (注：テキスト版では図は省略)
「したがって、計算した誤差の範囲を考えると、もし高乗算と下乗算で r[31] が異なれば、正しい修正は高乗算の r[31] をインクリメントすることです。あなたのコードでは、高乗算で r[31] がセットでかつ低乗算でクリアであるという場合にのみ r[32] にインクリメントを適用し、これは r[32] へのキャリーをもたらします。逆に、高乗算でクリアだが低乗算でセットになる場合は r[32] へのキャリーを生じないので、これを無視しているのですね。私の理解は正しいでしょうか？」
「はいです。私はこう考えています：下位 32 ビットは間違いなくすべて正確です。上位半分（高値部分）の積の上位 16 ビットをとる際の最悪の誤差は 2 だけです。したがって、もしこの二つの結果を 2 ビットだけオーバーラップさせることができれば、上位部分へのキャリーは最大一つに限定され、すべてを修正できます。」
私の乗算をより魔法的に
RISC-V ではかなり異なる結果になりますが、このトリックは実際に使用可能で有益であることが分かりました。
__mulsf3: h3.bextmi a2, a0, 23, 8 h3.bextmi a3, a1, 23, 8 h3.fcheck2e.s t0, a2, a3 bnez t0, __mulsf_special_exponent // ここでは、23 x 23 -> 46 ビットの積を二つの 32 ビット乗算で計算するための // マーク・オーウェン氏のトリックのバリエーションを使用します。 li t2, -1 << 23 andn a4, a0, t2 // x - 1, U0.23 andn a5, a1, t2 // y - 1, U0.23 bseti t2, a4, 23 add t2, t2, a5 // x + y - 1, U2.23 (後続で有用) xor a0, a0, a1 // 符号 add a1, a3, a2 // a2/a3 は現在解放済み addi a1, a1, -127 + 3 // 指数 mul a2, a4, a5 // exact[31:0] srli a4, a4, 7 srli a5, a5, 7 mul a4, a4, a5 // approx[45:14] // approx[45:0] の誤差は -2^31 < err <= 0 より少し良好。 // 二つの可能性のみ： // // * approx[45:31] が exact[45:31] に等しい // // * approx[31] と exact[31] が異なる；2^31 を加えると // 上位ビット 45:31 が修正される // // approx[32] へのキャリーは、approx[31] がセットでかつ exact [31] がクリアの場合に発生する。 // よって、approx[31] に !exact[31] を加えると、ビット 45:32（のみ）が修正される。 srli a4, a4, 17 bltz a2, 1f addi a4, a4, 1 1: srli a4, a4, 1 // 今や：exact[45:32] が a4 に、exact[31:0] が a2 に格納。U6.26 でスティッキー LSB と共にパッキング。 slli a4, a4, 12 li a5, 20 h3.ssrlsticky a2, a2, a5 or a4, a4, a2 // ここで得たのは：p = (x - 2^23)(y - 2^23) // = xy - (x + y) * 2^23 + 2^46 // 既に持っていた：s = x + y - 2^23 // // 望むもの： xy = p + (s + 2^23) * 2^23 - 2^46 // = p + s * 2^23 sh3add a4, t2, a4 // U3.29 へ正規化 h3.feadju3 a5, a4 h3.ssrlsticky a4, a4, a5 add a1, a1, a5 // パッキング前に符号を適用。 h3.xorsign a4, a0, a4 h3.fpackrq3.s a0, a4, a1 ret
いくつかの違いがあります。一つは、近似乗算の修正が少しずる賢く（cheeky）、ビット 31 までではなくビット 32 から修正される点です。また、暗黙の 1 の処理も完全に異なります。
これにより 3 サイクル節約となり、全体で 30 サイクルの
fmulこの技術を二重精度乗算における 52×52→104 ビットの内部積のための 32×32→64 ビット乗算に応用できる可能性を感じていますが、まだ証明していません。
⇥ wren.wtf に戻る