PentiumのFDIVバグという有名な設計ミスがあります。 これは、PentiumプロセッサのFDIV命令（浮動小数点数除算命令）が正しくない値*1を返すという問題で、回路設計の誤りが原因です。 ブルン定数（双子素数の逆数の和の極限で、収束することが知られている）を数値的に計算しようとした際に発見されたそうです。
今日起きたら以下が流れてきたので、調べてみました。
ブルン定数 (双子素数の逆数和) を計算する過程で Pentium CPU の浮動小数点演算にバグが発見されてリコールになった話興味深い…!
— パジョカ (Pajoca)⁰😺Nya!😸 (@Pajoca_) 2025年3月22日
双子素数(824633702441, 824633702443)に差し掛かった時、1/824633702441 の数値が実際の値よりおよそ
0.000000000000000000004517509
だけ小さくなり計算がずれて発覚 pic.twitter.com/fFjzjSxRft
PentiumのFDIVアルゴリズム
Pentiumでは、浮動小数点数除算にSRT法を導入しました。 それまでの除算回路は、1cycleあたりに1bitしか商を求められませんでしたが、SRT法の導入により、1cycleあたりに2bitの商を求めることができるようになりました。 SRT法は、除数の上位ビットと余りの上位ビットを使って表引きし、商を求めるアルゴリズムです。 ここで、上位ビットしか見ていないので、下位ビットの具合によっては、商に誤差が含まれることがあります。 しかしながら、SRT法は商を冗長表現で求めているため、誤差に対して多少の猶予が存在します。 そのため、上位ビットを十分に（例えば6ビットくらい）見れば、その猶予内の商を求めることができます。
SRT法については、以前に記事にしました（SRT除算について - よーる）が、Pentiumで採用されているSRT法はより洗練されています。 記事で説明したような方法だと、毎サイクル余りを計算する必要があります。 この時、加算器の幅は53ビットとかで、キャリーの伝搬に結構な時間がかかります。 そこでPentiumでは、余りを冗長表現（carry-save表現、二つの二進数の和）で持っています。 冗長表現であれば、余りの計算は全加算器一段分の遅延で実行可能です。 表引きのための上位ビットは、冗長表現の上位7ビットを取り出して、それを幅7ビットの加算器に入力することで得ます。 このようにした場合、余りの上位ビットとは異なる値となることがあります（無視した下位ビットからの繰上りがある場合）。 しかし、この違いさえも、SRT法の猶予の範囲内であり、正しく表をつくれば、猶予内の商を求めることができます。
残念ながら、Pentiumの除算回路は、「正しく表をつくれば」の部分がうまくいっていません。 2という商を立てるべきところ参照されるエントリに、誤って0が記録されてしまっている部分が5か所あります。 AND-ORタイプのプログラマブルロジックアレイ(PLA)で表が実装されているのですが（Intel's $475 million error: the silicon behind the Pentium division bugという記事に回路の写真があります）、このエントリは参照されないと信じられていたらしく、謎の節約精神が働いたのかもしれません（上記記事では、修正版Pentiumプロセッサの回路の写真も載っています。修正版Pentiumプロセッサでは外側を全部2や-2にしていて、その方がかえってPLAを節約できるようです。境界ぎりぎりを攻めていないので当然ですね）。
この間違ったエントリが参照されるのは、かなり稀なことです。 もちろん、特定の入力に対しては必ず起こりますが、そのような入力は少ないという意味です。 これは、carry-save表現が特定の条件を満たしていないと参照されないエントリだからです。 除算開始時には(carry=0, sum=dividend)で初期化されますが、ここから最低でも8回のSRTアルゴリズムの反復を経ないと、その条件を満たすことはないそうです[1]。 さらに、除数が特定の条件を満たすことも必要です。 問題のある5か所のエントリはそれぞれ、除数の仮数部が0b1.0001xxxx...., 0b1.0100xxxx...., 0b1.0111xxxx...., 0b1.1010xxxx...., 0b1.1101xxxx....の時に使われうるエントリです。 ここで、xxxx....の部分の先頭に1が6連続以上並んでいないと、carry-save表現がその特定の条件を満たすことはありません[1]。 この部分は論文を流し読みしただけでわかっていないので、今度読んでみます*2。
具体的なアルゴリズム
表が間違っている、と言われても、carry-save表現のところをどのように扱っているのかを具体的に見ない限り、どう間違っているのかを理解することは困難です。 そこで、verilog HDLでアルゴリズムを書き出してみました。
module divider( input wire [55:0] x, input wire [55:0] y ); function[2:0] pentium_srt_table; input [6:0] rem; input [3:0] div; reg[6:0] cor = { 3'b001, div } % 3; reg[6:0] th1 = div < 2 ? 7'b0000011 : div < 8 ? 7'b0000100 : div < 14 ? 7'b0000101 : 7'b0000110; reg[6:0] th2 = div < 2 ? 7'b0001100 : div < 5 ? 7'b0001110 : div < 8 ? 7'b0010000 : div < 11 ? 7'b0010010 : div < 14 ? 7'b0010100 : 7'b0010110; reg[6:0] th3 = ( div < 2 ? 7'b0010110 : div < 5 ? 7'b0011010 : div < 8 ? 7'b0011110 : div < 11 ? 7'b0100010 : div < 14 ? 7'b0100111 : 7'b0101010 ) + { 6'b0, cor != 0 }; // `{ 6'b0, cor != 0 }` cause the FDIV bug; `cor` is correct reg[6:0] th4 = div < 2 ? 7'b1111011 : div < 8 ? 7'b1111010 : div < 14 ? 7'b1111001 : 7'b1111000; reg[6:0] th5 = div < 2 ? 7'b1110010 : div < 5 ? 7'b1110000 : div < 8 ? 7'b1101110 : div < 11 ? 7'b1101100 : div < 14 ? 7'b1101010 : 7'b1101000; reg[6:0] th6 = ( div < 2 ? 7'b1101001 : div < 5 ? 7'b1100101 : div < 8 ? 7'b1100001 : div < 11 ? 7'b1011101 : div < 14 ? 7'b1011001 : 7'b1010101 ) - cor; if( rem < th1 ) pentium_srt_table = 0; else if( rem < th2 ) pentium_srt_table = 1; else if( rem < th3 ) pentium_srt_table = 2; else if( rem > th4 ) pentium_srt_table = 0; else if( rem > th5 ) pentium_srt_table = -1; else if( rem > th6 ) pentium_srt_table = -2; else pentium_srt_table = 0; endfunction reg [5:0] iter; reg [55:0] rem1, rem2; // reminder in the carry-save representation always @* begin rem1 = 0; rem2 = x; for( iter = 0; iter < 28; ++iter ) begin reg[6:0] approx_rem = rem1[55:49] + rem2[55:49]; // avoid slow full-length addition reg[2:0] q = pentium_srt_table( approx_rem, y[51:48] ); // -q*d reg[55:0] mqd = $signed(q) == 1 ? ~y : $signed(q) == 2 ? ~( y << 1 ) : $signed(q) == -2 ? y << 1 : $signed(q) == -1 ? y : 0; // carry-save addition reg[55:0] c = (rem1 & rem2 | rem2 & mqd | mqd & rem1) << 1 | { 55'b0, q == 1 | q == 2 }; reg[55:0] s = rem1 ^ rem2 ^ mqd; rem1 = c << 2; rem2 = s << 2; if( y[51:48] == 1 && approx_rem == 7'b0010111 ) $display( "it = %d @ %x / %x", iter, x, y ); if( y[51:48] == 4 && approx_rem == 7'b0011011 ) $display( "it = %d @ %x / %x", iter, x, y ); if( y[51:48] == 7 && approx_rem == 7'b0011111 ) $display( "it = %d @ %x / %x", iter, x, y ); if( y[51:48] == 10 && approx_rem == 7'b0100011 ) $display( "it = %d @ %x / %x", iter, x, y ); if( y[51:48] == 13 && approx_rem == 7'b0100111 ) $display( "it = %d @ %x / %x", iter, x, y ); end end endmodule
間違った値が記録されているエントリを読み出したとき、$displayでその事実を出力します。
なお、テーブルの構成はintel pentium srt lookup table - Google スプレッドシートを参照しました。
ただし、SRT Division Visualizer - Lukas Kollmerと突き合わせたところ、th4は1ずれていそうだったので、修正しました。
また、二の補数の取り方などはSRT Division Visualizer - Lukas Kollmerを参考にしました。
間違いの発生する入力
このHDLを使って、問題が発生する除算にはどういったものがあるかを調べてみました。
まず、全探索によって、分母が整数であるものの中で最も小さいものを求めてみました。 その結果、1818617/2359287がそれに該当するものだとわかりました。 分母は0x23fff7であり、0b10001 111111 11111110111なので、要件を満たしていることもわかります。
他にも、このHDLを使っていろいろわかりました。
最も有名で、[1]の論文にも載っている4195835/3145727は、9回目のイテレーション（it=8）で問題のあるエントリを参照するものの中で、最も小さい整数の分母を持つもののようです。 8回目までのイテレーションはcarry-save表現が特定の条件を満たすことがなく必ず正しい商が立つため、間違った除算結果が返ってくる場合でも“そこそこ”正しい結果となることが知られています。 9回目のイテレーションで問題があるエントリを参照するのは最速であり、相対誤差が大きめになる例なのでしょう。
ブルン定数の計算に出てくるのは、824633702441（双子素数の片方）の逆数の計算、つまり1/824633702441らしいですが、逆数の計算の中で最も小さい整数の分母を持つものは、1/3221224323でした。 分母は0xbffffb83であり、0b10111 111111 111111111101110000011なので、やはり要件を満たしています。
なお、上記HDLがPentiumの動作と同じである保証はありませんが、すくなくともこれらの入力をSRT Division Visualizer - Lukas Kollmerで確認したところ、たしかに問題のあるエントリを使うようです。
参考文献
[1] Tim Coe and Ping Tak Peter Tang, It Takes Six Ones To Reach a Flaw, In 12th IEEE Symposium on Computer Arithmetic (ARITH-12), pp. 140–146, July 1995.
