Intel 8087 チップ「FPTAN」命令のアルゴリズム逆解析

1980 年に発表された Intel 8087 は、IBM PC および他のシステムにおける浮動小数点演算の速度を大幅に向上させた歴史的なマイクロプロセッサです。本稿では、このチップが「接関（タンジェント）」計算を実行する

FPTAN

8087 は当時画期的な高速化を実現しており、接関計算にかかる時間を 13,000 マイクロ秒から 90 マイクロ秒 へと短縮しました。本記事では、チップ内部構造の調査結果と逆解析されたマイクロコードに基づき、この驚異的な性能を達成した仕組みを解説します。

チップ内部構造の解明

8087 の内部構造を可視化するためには、物理的にカバーを取り外し、顕微鏡で高精細な画像を取得しました。その結果、以下の主要なブロックが確認できました。

• マイクロコード ROM
  • チップの中心部に位置する大きな長方形領域です。
  • 内部に 1,648 個のマイクロ命令が格納されており、これらがチップ全体を制御しています。
• データパス（赤枠部）
  • チップの下半分（赤い四角）に位置する回路群です。
  • 80 ビットの浮動小数点値に対する計算を実行します。
• FPTAN に使用される機能ブロックの詳細
  • データパス内部に存在する主要なユニットが確認できます。

主要な機能ユニット一覧

アルゴリズム：CORDIC と多項式近似

三角関数を計算するための一般的なアプローチには、「CORDIC」と「多項式近似」の 2 つがあります。8087 は両者の手法を巧みに併用し、高い精度と高性能の両方を獲得しました。

1. CORDIC アルゴリズムとは

CORDIC (COordinate Rotation DIgital Computer) は、乗算や除法を行わずに、シフト、加算指令、およびテーブル参照（表引き）のみを用いて超越関数を高速に計算する巧妙なアルゴリズムです。

• 開発背景: 1956 年、マッハ 2 の爆撃機「B-58 ハスラー」向けのデジタルナビゲーションシステムとして開発されました。当時のトランジスタ技術では三角関数の生成が困難でしたが、簡素なハードウェアで高速計算を実現するために考案されました。
• 基本原理: 角度をベクトルに変換し、座標から三角関数値を導出します。鍵となるのは、角度を**特別な角度の列（$ \alpha_n = \arctan(2^{-n}) $）**に分解することです。これにより、複雑な乗算を行わずともビットシフトと加算だけで計算が可能です。
• 収束性: 反復処理（イテレーション）ごとに精度が向上し、急速に収束します。

2. 有理多項式近似の導入

CORDIC の精度は使用する項の数に依存します。64 ビットの完全な精度を得るには CORDIC 64 回分の計算が必要ですが、8087 はより高速な戦略を採用しました。

• ハイブリッド構成:
  • CORDIC: 入力角度の初期部分（約 16 ビット分）を処理。
  • 有理多項式 (Padé Approximant): CORDIC で残った微小な角度差に対して「$ \frac{3x}{3-x^2} $」というシンプルな式を使用。
• 精度: この近似の誤差は $ x^4 $ に比例し、入力が極めて小さいため最終的な 64 ビット精度要件を充足します。
• 高速化のポイント: 通常「分子 / 分母」として除算を行う必要がありますが、8087 はこれを独立した 2 つの値（X, Y）として返すため、実際の除算ステップを省略しています。これはコストゼロ（無料）の操作です。
  • 注釈: テイラー級数よりもこの比の方が優れているのは、tan 関数が $ \pi/2 $ で発散する特性に対し、多項式の比率であれば発散しないためです。

FPTAN 命令の実装フロー

8087 の

FPTAN

1. CORDIC の判定ビット決定（「疑似除算」）
2. 有理近似の計算
3. CORDIC の判定ビットに基づく回転適用（「疑似乗算」）

ステップごとの詳細解説

Step 1: 疑似除算（角度の分解）

入力角度を特別な角度の組み合わせへと分解します。目標角度より小さい特別角度を用いて減算し、その際に「1」を記録、しない場合は「0」と記録します。これにより、判定ビット列が生成されます。

• 実装手法: 長除法の原理に類似したシフト変換版を順次実行（「疑似除算」と呼ばれる）。
• 出力: 入力角度を特別角度によって減らし、最終的に微小な残差角度を残します。

Step 2: 有理近似計算（Padé Approximant）

Step 1 で得られた残差角度に対して、式 $ \frac{3x}{3-x^2} $ を適用して接関の値を計算します。

• 除算回避: マイクロコードは分子を初期ベクトルの Y 成分、分母を X 成分として扱うことで、高コストな除算ステップを実際に実行しません。
• 最高負荷演算: この段階で最もリソースを消費するのは角度の平方計算であり、完全な 64 ビット乗算が必要です。

Step 3: 疑似乗算（回転の適用）

Step 1 で決定された判定ビットに基づいて CORDIC の回転を適用します。これはある意味で「2 進数乗算」に似ており、判定ビットが 1 の場合に被乗数を加算するプロセスです。

• 順序: 逆順（最も小さな角度から先に）に適用し、丸め誤差を最小化します。
• ハードウェア効率: 個々の回転はシフト、加算、減算のみを使用するため非常に安価です。

マイクロコード実装の詳細分析

以下のリストには、

FPTAN

主要な処理パスの分岐

8087 は入力の指数（精度）に応じて 3 つのパスを自動選択します。

#1049

#1056

#1057

マイクロコードフローのハイライト

以下に FPTAN 命令の中核となるマイクロコードの一部を示します。

FPTAN:
#1039 st(0) -> tmpA          スタック上位の引数を tmpA レジスタに格納
#1040 stackPtr--             結果を格納するスタック位置に空きがあるか確認
#1041 stack overflow?         スタックオーバーフロー？
#1042 jmp #1022 if tmp empty/special/overflow/div exit if bad argument or stack overflow
                              もし引数不良、スタックアンダーフロー、スペシャル値、またはオーバーフローなら JMP #1022（除算例外終了）
#1043 jmp #1045 if tmpA:tag ZERO is argument 0?
                              もし tmpA のタグが ZERO で引数が 0 か？
#1044 except:precision       **精度例外**（0 を除く場合、他のすべての入力に対して発生）
#1045 expconst 0x3ff0        指数 >= -16 かテスト (CORDIC パスの閾値)
#1046 adder: tmpA:exp + 0 cin=1 argument exponent + 1
                              加算器：tmpA の指数に 0 を加算（キャリー入力で）、引数の指数 + 1 を得る
#1047 expConst -> Breg       定数 0x3ff0 を B レジスタへ
#1048 adder: sumreg:frac - Breg cin=1 subtract -15
                              加算器：和レジスタの有効数から B レジスタを減算（キャリー入力で、実質的に -15 を足す）
#1049 jmp #1061 if adder pos CORDIC if exponent >= -16
                              もし加算器の符号が正なら（指数 >= -16）、**CORDIC パス**へ JMP #1061

#--- CORDIC なパス (指数 < -16) ---
#1050 expconst 0x3fff        有理近似パスキープへの分岐指示
#1051 tmpA:exp -> Breg       オリジナルの指数をチェック
#1052 expConst -> tmpB:frac against exponent 0
                              定数を tmpB の有効数へ、指数 0 と比較
#1053 adder: tmpB:frac - Breg cin=1 subtract
                              加算器：tmpB の有効数から B レジスタを減算
#1054 sumreg:frac -> expConv 3fff - exp (i.e. -exp unbiased)
                              和レジスタの有効数を指数変換用に変換（バイアスなしの -exp）
#1055 sumreg:frac -> shiftcount
                              和レジスタの有効数をシフトカウンタへ
#1056 jmp #1082 if exponent[6:14] != 0 if exp > -64 goto rational approximation
                              もし指数の上位ビットが 0 でなければ（|exp| <= 64）、**有理近似**へ JMP #1082
#1057 high bit -> tmpB:frac otherwise return original argument
                              高位ビットを tmpB の有効数へ、さもなければ元の引数を返す

#--- 指数 >= -16 の CORDIC パス ---
#1049 (参照) ... CORDIC パス開始 ...

#1061 expconst 0x0010        CORDIC パス（回転適用準備）
#1062 sumreg:frac -> loopcounter loopcounter = exp + 16
                              和レジスタの有効数をループカウンタへ、ループ回数 = 指数 + 16
#1063 tmpA:frac -> tmpC      tmpC = オリジナルの角度（有効数部分）
#1064 trig const -> Breg/rnd 三角関数定数を B レジスタ・丸めモードへ
#1065 adder: tmpC - Breg cin=1, roundmode subtract first CORDIC angle
                              加算器：tmpC から最初の CORDIC 角度を減算（符号付き）
#1066 adder sign -> cordic, shl save bit in shift register
                              加算器の符号をシフトレジスタへ保存（決定ビットとして）

# ... CORDIC 擬似除算ループ (判定ビット記録) ...

#1079 jmp #1072 if not const latch zero bottom of CORDIC scan loop
                              もし定数ラッチがゼロでなければ、**CORDIC スキャンループの底部**で JMP #1072（ループ終了）
#1080 tmpC -> tmpA:frac      更新：tmpA, tmpB を tmpC から取得
#1081 expConst -> shiftcount 16 -> shiftcount (otherwise based on exponent)
                              シフトカウンタを 16 に設定（さもなければ指数に基づき）

#--- 有理近似計算 (Padé Approximation) ---
#1082 tmpA:frac -> tmpB:frac Padé approximation
                              **Padé近似計算開始**: tmpA の有効数を tmpB へ
#1083 tmpA:frac -> Breg      tmpA を ang（角度変数）として扱う
#1084 call SQUARE            angle^2 を計算（SQUARE サブルーチン呼出）
#1085 shift sumreg:frac R count byte bit
                              和レジスタの有効数を右シフト（平方値をスケール）
#1086 shift R -> sumreg:frac
                              シフト後を和レジスタへ
#1087 shift sumreg:frac R count byte bit
                              さらに右シフト
#1088 shift R -> Breg        B レジスタに angle^2 を保存
#1089 1.1 -> tmpB:frac       3 （指数を 1 にして）を表す定数（有効数 0x4000 の一部）
#1090 adder: tmpB:frac - Breg cin=0
                              加算器：tmpB から B レジスタ（3-ang^2 の計算）
#1091 sumreg:frac -> tmpB:frac tmpB (X) = 3-ang^2
                              結果を tmpB の有効数へ保存（分母 X）
#1092 shift tmpA:frac R 0 bytes, 1 bits
                              tmpA を右シフト
#1093 shift R -> Breg/rnd
                              シフト後の値を B レジスタ・丸emodeへ
#1094 adder: tmpA:frac + Breg cin=0, roundmode ang + ang/2
                              加算器：tmpA から B レジスタ（3×ang の計算）
#1095 sumreg:sign,frac,rnd -> tmpC tmpC (Y) = 3×ang with appropriate exponent
                              結果を tmpC へ保存（分子 Y）

# ... CORDIC 擬似乗算ループ (回転適用) ---

#1100 shift tmpC R 0 bytes, 1 bits CORDIC pseudo-multiplication path
                              **CORDIC 疑似乗算パス開始**: tmpC を右シフト
#1101 #15 -> loopcounter      ループカウンタを 15 に設定（CORDIC ステップ数）
#1102 shift R -> tmpC        シフト後の値を tmpC へ
#1103 jmp #1113 if not cordic[0] CORDIC: if saved decision bit...
                              もしシフトレジスタの上位ビット（決定ビット）が 1 でなければ、JMP #1113（この回転スキップ）

# ... 回転適用ロジック ...

#1115 jmp #1118 if cordic==0 done if shift register empty
                              シフトレジスタが全ゼロなら、**JMP #1118**（CORDIC ループ終了）

#--- 結果の正規化と出力 ---
#1118 expconst 0x3fff        CORDIC パス終了
#1119 tmpB:frac -> sumreg:frac Test top bit of tmpB (X)
                              tmpB の有効数の上位ビットをチェック
#1120 jmp #1124 if reg bit 63 If set, use exponent 0
                              もしビット 63 がセットなら、指数を 0 とする（正規化不要）
#1121 expconst 0x3ffe        さもなければ、指数を -1 と設定
#1122 shift sumreg:frac L 0 bytes, 1 bits shift left one bit to normalize
                              左シフトで正規化
#1123 shift L -> tmpB:frac   正規化後の値を tmpB へ
#1124 expConst -> tmpB:sign,exp store exponent in tmpB
                              符号・指数を tmpB へ保存（正规化済み）
#1125 tmpC -> tmpA:frac      tmpC (Y) を tmpA の有効数へ保存
#1132 tmpB -> st(0)          **tmpB (X:分母)** をスタック上位へ保存
#1134 tmpA -> st(0)          **tmpA (Y:分子)** をスタック次位に保存
#1136 RNI                    FPTAN 命令終了

マイクロコードの実装上の工夫

• 疑似除算の最適化:
  • コードは格納された 16 つの角度を使用して 16 回の CORDIC ステップを実行しますが、入力に対して大きな角度をスキップするなどの最適化が含まれています。
  • ループカウンタは入力角度の指数に基づき初期化されます（#1062）。
• ハードウェア効率:
  • シフトレジスタ: 内部有効数バスへのアクセスが必要なく、ビットの入出力には 2 つのラインのみで十分です。そのため、レイアウト上の未使用スペースを活用して配置できました。
  • コードの実行: コードは浮動小数点演算ではなく、64 ビット整数演算を使用して最適化されています。数値は帯有しない指数を持つ定点数として扱われ、ループ中で暗黙的に指数が変化（スケール）します。

性能と精度について

• 遅延命令:

  FPTAN

  は 8087 命令の中で比較的遅く、通常 450 クロックサイクルを必要とする場合があります（入力値により 30〜540 サイクル）。
• 時間分布の検証 (例: 入力 0.95):
  • CORDIC 疑似除算: 33%
  • 有理多項式計算: 15%（主に乗算部分）
  • CORDIC 疑似乗算: 47%
  • オーバーヘッド: 5%

結論と考察

Intel 8087 は、CORDIC という古典的なアルゴリズムと、 Padé 近似という現代的な数学的手法を巧みに組み合わせた独自のハイブリッド構成を採用しました。当時の技術水準において、これほど高精度かつ高速な浮動小数点演算を実現したのは画期的でした。

• 歴史的意義: Pentium など後継チップでは、ハードウェアの進化により多項式近似のみを使用する方向へ移行しましたが、8087 のこのハイブリッドアプローチは独自性と優位性を示しています。
• 現代との対比: 現在は MKL や SVML といったライブラリが SIMD 命令を活用した高性能実装を提供しており、x87 FPU (80 ビット浮動小数点) はほぼ廃れていますが、当時のエンジニアリングの知恵は依然として価値があります。

本記事では、物理的なチップ調査とマイクロコードの逆解析を通じて、Intel 8087 の

FPTAN