Payne-Hanek 法詳解：ビットスライドで読み解く高精度な剰余演算と Rust による最小実装

今回は、範囲還元や剰余計算に用いられる高精度な手法“Payne-Hanek法”について、具体例を挙げて実際に実装しながら、わかりやすく解説していこうと思います。

1. 理論の視覚化（ビットの「窓」と「スライド」）

まずは、Payne-Hanek法が「何を解決するためのものか」を整理しましょう。一言で言えば、「巨大な数に対しても、精度を1ビットも無駄にせず剰余（あまり）を求めるための技術」です。

1.1. 「桁落ち」問題

通常、ある数 \(x\) を周期 \(P\) で割った「あまり」は、以下の式で求められます。
\[x \pmod P = x – (n \cdot P)\]
（※ \(n\) は \(x/P\) の整数部）
しかし、ここに大きな罠があります。\(x\) が非常に巨大（例えば \(10^{20}\)）で、\(P\) が小さい（\(\pi \approx 3.14\)など）場合、\(x\) と \(n \cdot P\) はほとんど同じ値になってしまいます。

浮動小数点数で「ほぼ同じ大きさの数」同士を引き算すると、計算結果の有効数字がごっそり消えてしまう「桁落ち」が発生します。
これは、後述する範囲還元や、二次方程式の解の公式を計算させる場合にxの係数が特に大きい場合などで問題になります。

二次方程式解の公式を計算機で計算させるとき
\(ax^2 + bx + c = 0\) の解の公式 \(x = \frac{-b \pm \sqrt{b^2 – 4ac}}{2a}\) を考えます。
具体的な数値の例: \(a=1, b=10^8, c=1\)
計算の危うさ:
\(\sqrt{b^2 – 4ac}\) を計算すると、\(\sqrt{10^{16} – 4}\) となり、これはほぼ \(10^8\) に等しくなります。
桁落ち:
\[x_1 = \frac{-10^8 + \sqrt{10^{16} – 4}}{2}\]
ここでは「\(10^8\) に極めて近い数」から 「\(10^8\)」を引くことになるため、分子で深刻な桁落ちが発生し、解の精度が著しく低下します。

2. 状況に応じた「二つの手法」

Payne-Hanek法では、桁落ちを起こさずに剰余（あまり）を求めることができます。これが特に生かせるのが範囲還元（巨大な引数を小さな範囲に収める処理）です。これには、主に2つの手法を使い分けます。

• Cody-Waite法（軽量級）
  • 特徴: 高速。定数 \(\pi\) を少し多めの桁数で持っておき、通常の演算を行う。
  • 限界: \(x\) がある程度大きくなると精度が保てなくなり、無力化する。
• Payne-Hanek法（重いが正確）
  • 特徴: どれだけ \(x\) が巨大でも、数学的に完璧な精度で「あまり」を抽出できる。
  • 代償: 計算コストが高い。ビット列を「窓（ウィンドウ）」で覗き見るような特殊な処理を行う。

2. 巨大な数 x と高精度な定規

実際に巨大な数に対する余りの計算の実装をしながら、Payne-Hanek法の実装をしていこうと思います。
\(x\) が非常に大きい（例：\(2^{100}\)）場合、\(x \cdot (1/\pi)\) を浮動小数点数で計算すると、普通の精度では「整数部」だけで桁が埋まってしまい、肝心の「小数部（あまり）」が消えてしまいます。 Payne-Hanekの真髄は、「整数部を計算せずに、小数部が最初に出てくる場所（ビット位置）を計算で特定し、そこから数ワード分だけ抜き出して計算する」点にあります。

\[x \cdot \frac{1}{P} = \text{整数部（何周したか）} + \text{小数部（あまりの位置）}\]

2.1 仕組み解説その１ -小数部分を求める-

まず、私たちがやりたいことは \(x \cdot \frac{1}{P}\) の「小数部分」だけを完璧な精度で取り出すことです。

余り部分は、1/Pを掛けた結果の小数部分に出てくるので、小数部分だけ計算すれば剰余が求まることになります。

2.2 仕組み解説その２ -小数の桁を見つける-

上で説明したようにPayne-Hanek法では、1/Pを掛けて小数部分を見つけるという方法をとりますが、普通に計算しただけだとどのみち整数部分が桁数を占有し、小数部分の精度が出せません。
そこで賢く小数部分だけ求める手法として、浮動小数点数の指数部分を利用します。

簡易な例（\(x = 2^3\), \(1/P = 0.11011011…\) の場合）
わかりやすくするために、非常に小さな数でシミュレーションしてみます。

① \(x\) を準備する

\(x = 8\) （2進数で \(2^3\)）とします。このとき、浮動小数点の内部では「指数部」に 3 という情報が入っています。

② 定数テープから「窓」の位置を決める

\(1/P\) を \(0.1101101101…\) というビット列だとします。
通常、掛け算をするとこうなります：
\[2^3 \times 0.1101101101 = 110.1101101…\]
ここで、小数点より左側の 110（整数部）は「何周したか」を表すだけなので、余りの計算には一切不要です。

Payne-Hanek法はこう考えます。
「\(x\) の指数が \(3\) なんだから、\(1/P\)の最初の \(3\) ビットは整数になっちゃうな。じゃあ、最初から 4 ビット目以降だけを狙って掛け算すればいいじゃないか」

③ 抽出と掛け算

指数の分だけテープを「スライド」させて、必要な部分（小数部になるはずのビット列）だけを抜き出し、仮数部と掛け合わせます。ここが Payne-Hanek 法のコアの部分です。「指数が \(3\)」という情報を使って、定数テープから必要なビットだけをピンポイントで抜き出し計算ます。言葉だけでの説明だと理解しづらいと思うので以下(2.3項)に実際に計算してみた例を示します。

2.3 計算手順の例

簡易な例（\(x = 2^3\), \(1/P = 0.11011011…\) の場合）

1. 「窓」をセットする

\(x = 2^3\) なので、定数テープの 4ビット目（\(2^{-4}\) の位） から後ろを「窓」で切り取ります。

• 定数テープ (\(1/P\)): 0 . 1 1 0 [ 1 1 0 1 1 0 1 1 … ]
• 抜き出した値 (\(F_{win}\)): 0 . 1 1 0 1 1 0 1 1 ...
  （※窓の中身を便宜上 \(0\) から始まる小数として扱います）

浮動小数点数を \(M \times 2^E \)とすると指数部分はシフト演算

2. 固定小数点として掛け算する

次に、\(x\) の 仮数部（Mantissa） と、窓から抜き出した 定数ビット列 を掛け合わせます。
今回の例では \(x = 2^3\)（バイナリで \(1.000… \times 2^3\)）なので、仮数部はシンプルに 1 です。
\[\begin{aligned} (\text{仮数部}) &\times (\text{窓の値}) \\ 1 &\times 0.11011011… = 0.11011011… \end{aligned}\]

3. 結果の解釈

この掛け算の結果得られた 0.11011011... こそが、私たちが欲しかった 「あまり（小数部）」そのもの です！

Payne-Hanek法では、最初から「整数部になることが確定している3ビット（110）」を無視して、その次から掛け算を始めました。その結果、面倒な引き算をすることなく、ダイレクトに小数部のビット列だけが手元に残すことができます。(実際は、数値によっては整数部分が繰り下がるようなことがある為、少しの整数部分＋小数部分多めで窓をとる必要があります)

3. 3.14宇宙での最小実装

無理数で実装するとデバッグのために別途、もっと高精度な計算機を用意しないといけなくなってしまうため、円周率が3.14の別世界線の仮定の下、x[radian]を[2/3.14]の範囲に範囲還元する関数を実装してみます。 範囲還元が必要な理由：(記事は現在利用できません: ID 7311)

• 周期 \(P\): \(3.14/2 = 1.57\)
• 逆数 \(F = 1/P\): 約 \(0.636942675159235..\)
• 定数テーブル: \(F\) を \(2^{64}\) 倍した整数値を用いる。

3.14世界用の定数算出

Pythonなど高精度計算ができる環境であらかじめ計算します。

# pythonのコード例
from decimal import Decimal, getcontext

# 十分な精度（100ビット分以上）を確保します
getcontext().prec = 100
two_over_pi = Decimal(2) / (Decimal('3.14')) //Decimalで定義しないと誤差が大きくなるので注意

# 小数点以下をビット列として抽出
val = two_over_pi - int(two_over_pi)
bits = []
for _ in range(10):  # 64ビット × 10ワード分
    word = 0
    for i in range(64):
        val *= 2
        if val >= 1:
            word |= (1 << (63 - i))
            val -= 1
    bits.append(f"0x{word:016X}")

print(", ".join(bits))

16進数で10ワード分出力されます。

F = [0xA30EACD73C54CA30, 0xEACD73C54CA30EAC,
0xD73C54CA30EACD73, 0xC54CA30EACD73C54, 
0xCA30EACD73C54CA3, 0x0EC796DCE768B50E, 
0xBA7C6AECD378C700, 0x9A446D95C1F48D03, 
0x2ECF011D3759D081, 0x18630A7220343C2E]

最小実装コード (Rust)

// pi = 3.14 と仮定した世界の 1/pi (高精度ビット列の先頭)
const INV_PI_314: [u64; 10] = [
    0xA30EACD73C54CA30, 0xEACD73C54CA30EAC,
    0xD73C54CA30EACD73, 0xC54CA30EACD73C54, 
    0xCA30EACD73C54CA3, 0x0EC796DCE768B50E, 
    0xBA7C6AECD378C700, 0x9A446D95C1F48D03, 
    0x2ECF011D3759D081, 0x18630A7220343C2E
];


pub fn payne_hanek_sandbox(x: f64) -> (f64, u64) {
    let bits = x.to_bits();
    let exp = ((bits >> 52) & 0x7ff) as i32 - 1023;
    let mantissa = (bits & 0x000f_ffff_ffff_ffff) | 0x0010_0000_0000_0000;

    // 1. インデックスの特定
    // x = M * 2^exp, Table = F * 2^-64
    // 求める積は M * F * 2^(exp - 64) 
    // 1ワードの境界をまたぐため、exp/64 でベースとなるワードを特定する
    let idx = (exp.max(0) / 64) as usize;
    let shift_in_word = (exp.max(0) % 64) as u32;

    // 2. 必要なワードの取得
    let f0 = INV_PI_314[idx] as u128;
    let f1 = INV_PI_314[idx + 1] as u128;

    // 3. 掛け算 (128bit精度)
    // 仮数(53bit) * 定数(64bit) = 117bit
    let m = mantissa as u128;
    let p0 = m * f0; 
    let p1 = m * f1;

    // 4. ビットアライメントの核心
    // 数学的に、積の「小数点」の位置は 116 - shift_in_word ビット目になる
    let split_bit = 116 - shift_in_word;

    // 整数部 (何周したか / 象限)
    let integer_part = (p0 >> split_bit) as u64;

    // 小数部 (あまり) 
    // p0 の小数点以下と、p1 からの繰り上がりを合成して 64bit の小数ビットを作る
    let fractional_bits = if split_bit >= 64 {
        // 小数点が上位ワードの中にある場合
        (p0 >> (split_bit - 64)) as u64
    } else {
        // 小数点がワードをまたぐ場合
        (p0 << (64 - split_bit)) as u64 | (p1 >> split_bit) as u64
    };

    let reduced = (fractional_bits as f64) / (2.0f64.powi(64));

    (reduced, integer_part)
}

コメントはgeminiに書いてもらいました。

4. 検証（3.14宇宙の答え合わせ）

サンドボックス環境なので自明な値 3.14などで簡単にテストできます。以下が実際に検証してみた結果です。

テストケースと結果

計算結果を見ると、小さな誤差が含まれるようです。これは倍精度の計算機イプシロン約 \(2.22 \times 10^{-16}\)と同程度のものであり、浮動小数点数の仕組み上仕方ないものです。

また、1億は元の浮動小数点数に表現した時点で、仮数部の貴重な53ビット（約16桁分）のうち、1億の位から１までの表現に 27ビット分 を使い切っています。浮動小数点数で１億を入力した時点で小数部分は26ビットしかなく、これを10進数の精度に換算すると、約7.8桁になるため、余りの計算としては、10桁もあれば十分理想的なのでうまくできているといえます。

5. まとめ：本実装（本物の \(\pi\)）への復帰

3.14宇宙で「指数に応じてビットを抽出し、象限と小数部を切り分ける」ロジック(Payne-Hanek法)を実装してみました。今回10ワード分の1/Pを用意しましたが、これだけで理論上\(2^{640}\) 程度まで対応可能であるため倍精度の浮動小数点数ならすべて対応することができます。

私たちが住んでいる世界は円周率＝3.14ではないので、実際の計算では、あらかじめ用意するワードを3.14159265… = πで計算しておく必要がありますが、Payne-Hanek法の仕組みは今回の解説でご理解いただけたのではないかと思います。

ちゃんとした範囲還元関数を作る残り手順

1. テーブルの差し替え: 上の方で作成した \(\frac{2}{3.14}\)を\(\frac{2}{\pi}\)の正確なビット列テーブルに置き換えるだけで、本物の三角関数の範囲還元が可能になります。
2. 検証: 1億 [rad] のような巨大な値を入力し、標準ライブラリ（std::f64::consts::PI など）と比較して一致を確認すれば、それはもう「標準ライブラリ級」の実装と言えます。

TUI関数電卓シリーズの展望

このサイトで連載中の「TUI関数電卓を作る」シリーズでは、今回作成した基礎ロジックをもとに、sin や cos といった三角関数の範囲還元（Range Reduction）を実際に実装していきます。

三角関数編の公開は少し先になりますが、ここまでの「剰余計算の極意」を理解していれば、関数電卓の心臓部がどう動いているのかが理解しやすいと思います。

今回紹介したPayne-Hanek法を生かして三角関数を実装する：(記事は現在利用できません: ID 7334)