【Rustで作る電卓】Payne-Hanek法 Rustで挑む三角関数の実装(3) 高精度レンジリダクション #20

// 2/pi をビット列にしたもの(断片テーブル)を用意
const TWO_OVER_PI_TABLE: [u64; 10] = [
    0xA2F9836E4E441529, 0xFC271C4124204598, 0x23FB11315842E121,
    0xEE89B189FF1D91AD, 0x71245982399E4934, 0x61A2A83FB44A2941,
    0x6139824128084121, 0x124982399E493461, 0xA2A83FB44A294161,
    0x3982412808412161,
];

ub fn reduce_payne_hanek(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;

    let idx = (exp >> 6).max(0) as usize;
    let shift = (exp & 63) as u32;

    let m = mantissa as u128;
    // 3ワード分計算して精度を確保
    let p0 = m * (INV_PI_TABLE[idx] as u128);
    let p1 = m * (INV_PI_TABLE[idx + 1] as u128);

    // 小数点（1の位）がどこにあるかを計算
    // 理論上、m*T[idx]の結果において小数点位置は bit 116 から指数の余り分左にズレる
    // ここはマジックナンバーですが、下のコラムで説明します。
    let point_pos = 116 - shift;

    // 1. 象限（整数部の下位2ビット）の抽出
    let quadrant = (p0 >> point_pos) as u64 & 3;

    // 2. 小数部（剰余）の抽出
    // 小数点以下のビットを 128bit の固定小数として取得する
    let fraction_mask = (1u128 << point_pos) - 1;
    let frac_hi = p0 & fraction_mask;
    let frac_lo = p1;

    // 比率としての小数 [0, 1) を計算
    // frac_hi は 2^point_pos の重み、frac_lo はさらにその 2^-64 下
    let f_val = (frac_hi as f64) * 2.0f64.powi(-(point_pos as i32))
              + (frac_lo as f64) * 2.0f64.powi(-(point_pos as i32 + 64));

    // --- 精度向上のためのレンジ調整 ---
    // もし余りが 0.5 を超えていたら、次の象限に近いので反転させる
    // これにより reduced_x は常に [0, pi/4] に収まり、kernel の精度が最大化される
    let (final_frac, final_quadrant) = if f_val > 0.5 {
        (1.0 - f_val, quadrant + 1)
    } else {
        (f_val, quadrant)
    };

    let reduced_rad = final_frac * (std::f64::consts::PI / 2.0);

    (reduced_rad, final_quadrant)
}

/// Payne-Hanek法を統合した sin 関数 hanek_sin
pub fn sin(x: f64) -> f64 {
    let abs_x = x.abs();
    
    // 小さな数（1.0e6未満）は、これまでの安全な方法で計算する
    if abs_x < 1000000.0 {
        return  sin4small(x); // 前回実装したsin関数を名称変更しました。
    }

    let (reduced_x, quadrant) = reduce_payne_hanek(abs_x);
    let sign = if x < 0.0 { -1.0 } else { 1.0 };

    let res = match quadrant {
        0 => kernel_sin(reduced_x),
        1 => kernel_cos(reduced_x), // sin(pi/2 + x) = cos(x)
        2 => -kernel_sin(reduced_x), // sin(pi + x) = -sin(x)
        3 => -kernel_cos(reduced_x), // sin(3pi/2 + x) = -cos(x)
        _ => unreachable!(),
    };
    res * sign
}