【Rust自作TUI電卓】ニュートン法×マジックナンバーで実装する汎用n乗根 #13

pub fn custom_cbrt(x: f64) -> f64 {
    if x == 0.0 {
        return 0.0;
    }

    let mut guess = x / 3.0; // 初期値は立方根の1/3
    let mut prev_guess;

    for _ in 0..20 {
        prev_guess = guess;
        guess = (2.0 * guess + x / (guess * guess)) / 3.0; // ニュートン法の更新式
        if (guess - prev_guess).abs() < f64::EPSILON {
            break;
        }
    }
    guess
}

pub fn custom_root_n(x: f64, n: u32) -> f64 {
    // エッジケースの処理
    if n == 0 { return f64::NAN; }
    if n == 1 { return x; }
    if x == 0.0 { return 0.0; }

    // 負の数に対する処理
    if x < 0.0 {
        if n % 2 == 0 {
            return f64::NAN; // 偶数乗根は実数解なし
        } else {
            // 奇数乗根なら絶対値で計算して最後に符号を戻す
            return -custom_root_n(-x, n);
        }
    }

    let n_f = n as f64;
    // 初期値の推定。ニュートン法準拠
    let mut guess = x / n_f; 
    let mut prev_guess;

    for _ in 0..30 { 
        prev_guess = guess;
        let x_pow_n_minus_1 = guess.powi(n as i32 - 1);
        guess = ((n_f - 1.0) * guess + x / x_pow_n_minus_1) / n_f;
        if (guess - prev_guess).abs() < f64::EPSILON {
            break;
        }
    }
    guess
}

// --- 3乗根専用 ---
pub fn custom_cbrt(x: f64) -> f64 {
    if x == 0.0 { return 0.0; }
    if x < 0.0 { return -custom_cbrt(-x); }

    let i = x.to_bits();
    let mut guess = f64::from_bits((i / 3) + 0x2a9f4795cfce2000);

    // ニュートン法
    for _ in 0..3 {
        let g2 = guess * guess;
        guess = (2.0 * guess + x / g2) / 3.0;
    }
    guess
}

// --- 汎用n乗根 ---
pub fn custom_root_n(x: f64, n: u32) -> f64 {
    // エッジケース
    if n == 0 { return f64::NAN; }
    if n == 1 { return x; }
    if x == 0.0 { return 0.0; }
    if x < 0.0 {
        if n % 2 == 0 { return f64::NAN; }
        else { return -custom_root_n(-x, n); }
    }

    // ステップ1: nに応じて初期値を推定
    let mut guess = match n {
        2 => sqrt(x)
        3 => custom_cbrt(x),
        4 => {
            let i = x.to_bits();
            f64::from_bits((i / 4) + 0x2ff6f08899d39000)
        },
        5 => {
            let i = x.to_bits();
            f64::from_bits((i / 5) + 0x3325da722c26d000)
        },
        6 => {
            let i = x.to_bits();
            f64::from_bits((i / 6) + 0x354760eed2f65000)
        },
        _ => {
            // 7乗根以上：定数を動的に計算
            let i = x.to_bits();
            let n_f = n as f64;
            let sigma = 0.0430; // 補正係数(大体の値なので最適化はできてない)
            let bias = 1023.0;
            let l = (1u64 << 52) as f64;
            
            let magic = (((n_f - 1.0) / n_f) * (bias - sigma) * l) as u64;
            f64::from_bits((i / n as u64) + magic)
        }
    };

    // すでに cbrt などで計算済みの場合はそのまま返す
    if n <= 3 { return guess; }

    // ニュートン法
    let n_f = n as f64;
    for _ in 0..5 {
        let prev_guess = guess;
        let x_pow_n_minus_1 = guess.powi(n as i32 - 1);
        if x_pow_n_minus_1 == 0.0 { break; }
        
        guess = ((n_f - 1.0) * guess + x / x_pow_n_minus_1) / n_f;

        if (guess - prev_guess).abs() < f64::EPSILON {
            break;
        }
    }
    guess
}