// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// IntDivision.hpp
// 多倍長整数の除算アルゴリズム
//
// このファイルでは、多倍長整数ライブラリの除算アルゴリズムを提供します。
// Knuthの除算アルゴリズム（Algorithm D）に基づいた実装です。
//
// 主な機能:
// - Knuthの除算アルゴリズム
// - 効率的な商と余りの計算
// - 特殊ケースの最適化
// - 状態管理の統合

#ifndef CALX_INT_DIVISION_HPP
#define CALX_INT_DIVISION_HPP

#include <math/core/basic_types.hpp>
#include <math/core/mp/Int.hpp>
#include <math/core/numeric_state.hpp>
#include <utility>
#include <algorithm>
#include <vector>
#include <cstdint>
#include <bit>

namespace calx {

    /**
     * @brief 除算アルゴリズムの選択
     *
     * テストやベンチマークで異なる除算アルゴリズムを切り替えるための列挙型。
     * グローバル変数 g_division_algorithm で選択されたアルゴリズムが使用される。
     */
    enum class DivisionAlgorithm {
        Auto,              // 自動選択（サイズに応じて最適なアルゴリズムを選択）
        BitByBit,          // ビット単位ロングディビジョン（O(n²·w), 確実だが遅い）
        SingleWordOnly,    // 単一ワード除数のみ最適化（それ以外はBitByBit）
        Knuth,             // Knuth Algorithm D（O(n²), 高速で正確）
        Newton,            // ニュートン法による除算（将来の拡張用）
        BurnikelZiegler,   // Burnikel-Ziegler法（将来の拡張用）
        PowerOfTwo,        // 2のべき乗専用（右シフト）
        PowerOfTen         // 10のべき乗専用（内部でKnuthを使用）
    };

    /**
     * @brief グローバル除算アルゴリズム設定
     *
     * デフォルトは Auto（自動選択）。
     * テストやベンチマークでこの変数を変更することで、
     * 使用する除算アルゴリズムを切り替えることができる。
     *
     * 例:
     *   calx::g_division_algorithm = DivisionAlgorithm::Knuth;
     *   Int result = a / b;  // Knuth Algorithm D を使用
     *
     *   calx::g_division_algorithm = DivisionAlgorithm::BitByBit;
     *   Int result2 = a / b; // ビット単位除算を使用
     */
    inline DivisionAlgorithm g_division_algorithm = DivisionAlgorithm::Auto;

    /**
     * @brief 多倍長整数の除算アルゴリズムをサポートするクラス
     */
    class IntDivision {
    public:
        // ================================================================
        // Exact division（割り切れることが前提の高速除算）
        // Toom-Cook 補間などで使用。計算量 O(n)。
        // ================================================================

        /**
         * @brief 3 での exact division（3 で割り切れることが前提）
         *
         * Montgomery inverse を使用: 3^(-1) mod 2^64 = 0xAAAAAAAAAAAAAAAB
         * 計算量: O(n) — 汎用除算 O(n²) の代替
         *
         * @param x 被除数（3 で割り切れること）
         * @return x / 3
         */
        static Int divexactBy3(const Int& x) {
            if (x.isZero()) return Int(0);
            if (x.isSpecialState()) return x;

            // 符号を保存して絶対値で計算
            int result_sign = x.getSign();
            const auto& words = x.words();
            size_t n = words.size();

            // Montgomery inverse: 3 * INV3 ≡ 1 (mod 2^64)
            static constexpr uint64_t INV3 = 0xAAAAAAAAAAAAAAABULL;

            std::vector<uint64_t> out(n);
            uint64_t carry = 0;

            for (size_t i = 0; i < n; i++) {
                // x[i] - carry（borrow 検出）
                uint64_t xi = words[i];
                uint64_t xi_minus_carry = xi - carry;
                uint64_t borrow = (xi < carry) ? 1 : 0;

                // y[i] = (x[i] - carry) * INV3 mod 2^64
                uint64_t yi = xi_minus_carry * INV3;
                out[i] = yi;

                // carry = (yi * 3).high + borrow
                // 注: 3*yi ≡ xi_minus_carry (mod 2^64) なので prod.low == xi_minus_carry
                UInt128 prod = UInt128::multiply(yi, 3);
                carry = prod.high + borrow;
            }

            // 結果を構築
            Int result;
            result.m_words.assign(out.data(), out.data() + out.size());
            result.setSign(result_sign);
            result.normalize();
            return result;
        }

        /**
         * @brief 24 での exact division（24 で割り切れることが前提）
         *
         * 分解: x/24 = (x >> 3) / 3 — shift + divexactBy3
         * 計算量: O(n)
         *
         * @param x 被除数（24 で割り切れること）
         * @return x / 24
         */
        static Int divexactBy24(const Int& x) {
            if (x.isZero()) return Int(0);
            if (x.isSpecialState()) return x;

            // 符号を保存して絶対値で計算
            int original_sign = x.getSign();
            const auto& words = x.words();
            size_t n = words.size();

            // Phase 1: right shift by 3 bits (exact division by 8)
            std::vector<uint64_t> shifted(n);
            for (size_t i = 0; i < n - 1; i++) {
                shifted[i] = (words[i] >> 3) | (words[i + 1] << 61);
            }
            shifted[n - 1] = words[n - 1] >> 3;

            // シフト結果を Int に変換
            Int shifted_int;
            shifted_int.m_words.assign(shifted.data(), shifted.data() + shifted.size());
            shifted_int.setSign(original_sign);
            shifted_int.normalize();

            // Phase 2: divexact by 3
            return divexactBy3(shifted_int);
        }

        // ================================================================
        // Knuth Algorithm D
        // ================================================================

        /**
         * @brief Knuthの除算アルゴリズム
         * 参考: D. E. Knuth, "The Art of Computer Programming", Volume 2, 4.3.1
         *
         * @param dividend 被除数
         * @param divisor 除数
         * @param remainder 余りの格納先
         * @return 商
         */
        static Int divKnuth(const Int& dividend, const Int& divisor, Int& remainder) {
            // 特殊状態の処理
            if (divisor.isZero()) {
                // ゼロ除算の処理
                if (dividend.isZero()) {
                    // 0/0 はNaN
                    remainder = Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                    return Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                }
                else {
                    // x/0 (x≠0) もNaN（整数除算では商が定義できない）
                    remainder = Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                    return Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                }
            }

            if (dividend.isZero()) {
                remainder = Int(0);
                return Int(0);
            }

            // NaNや無限大などの処理
            if (dividend.getState() != NumericState::Normal || divisor.getState() != NumericState::Normal) {
                return IntSpecialStates::handleDivision(dividend, divisor);
            }

            // 絶対値が1の除数の特別処理
            if (abs(divisor).isOne()) {
                remainder = Int(0);
                Int q = abs(dividend);
                if (dividend.getSign() != divisor.getSign() && !q.isZero()) {
                    q.negate();
                }
                return q;
            }

            // 符号の処理
            bool quotient_negative = dividend.getSign() != divisor.getSign();
            bool remainder_negative = dividend.getSign() < 0;
            Int abs_dividend = abs(dividend);
            Int abs_divisor = abs(divisor);

            // 被除数が除数より小さい場合
            if (abs_dividend < abs_divisor) {
                remainder = dividend; // 余りは被除数そのもの
                return Int(0);
            }

            // 小さな値の場合、直接計算を行う
            if (abs_dividend.size() == 1 && abs_divisor.size() == 1) {
                uint64_t a = abs_dividend.word(0);
                uint64_t b = abs_divisor.word(0);
                uint64_t q = a / b;
                uint64_t r = a % b;

                // 結果の構築
                remainder = Int(static_cast<int64_t>(r));
                if (remainder_negative && r != 0) {
                    remainder.negate();
                }

                Int quotient(static_cast<int64_t>(q));
                if (quotient_negative && q != 0) {
                    quotient.negate();
                }

                return quotient;
            }

            // ここからKnuthの除算アルゴリズム
            // 入力データの準備（リトルエンディアン形式）
            std::vector<uint64_t> U = abs_dividend.words(); // 被除数
            std::vector<uint64_t> V = abs_divisor.words();  // 除数

            // 正規化（normalize）- 最上位桁を調整してdivisor[n-1] >= 2^32にする
            int normalization_shift = 0;
            if (!V.empty() && V.back() != 0) {
                // 正規化シフト量: 最上位ビットが立つまでシフト
                // Knuth Algorithm D では divisor の最上位ワードの MSB が立つ必要がある
                normalization_shift = static_cast<int>(std::countl_zero(V.back()));
            }

            if (normalization_shift > 0) {
                // 被除数と除数を正規化（左シフト）
                Int normalized_dividend = abs_dividend << normalization_shift;
                Int normalized_divisor = abs_divisor << normalization_shift;

                U = normalized_dividend.words();
                V = normalized_divisor.words();
            }

            // 桁数の計算
            size_t m = U.size();
            size_t n = V.size();

            // 商の格納領域確保
            // 商の長さは被除数の長さ - 除数の長さ + 1（最大）
            std::vector<uint64_t> Q;
            if (m >= n) {
                Q.resize(m - n + 1, 0);
            }
            else {
                Q.resize(1, 0);
            }

            // 作業用配列（Uをコピー）
            std::vector<uint64_t> U_work = U;
            // オーバーフローを防ぐため最上位に0を追加
            U_work.push_back(0);

            // Knuthのアルゴリズム D（多倍長整数除算）の実装
            for (int j = static_cast<int>(m - n); j >= 0; j--) {
                // 境界チェック
                size_t idx_high = j + n;
                size_t idx_mid = j + n - 1;
                if (idx_high >= U_work.size() || idx_mid >= U_work.size()) {
                    throw std::out_of_range("U_work index out of bounds");
                }

                // 商の推定値を計算
                uint64_t dividend_high = U_work[j + n];
                uint64_t dividend_mid = U_work[j + n - 1];
                uint64_t divisor_high = V[n - 1];

                // 商の推定計算
                uint64_t q_hat;
                uint64_t r_hat;
                bool skip_correction = false;

                // 1ワード除算の改善された計算
                if (dividend_high == 0) {
                    // 簡単なケース - 通常の単純な除算
                    q_hat = dividend_mid / divisor_high;
                    r_hat = dividend_mid % divisor_high;
                }
                else if (dividend_high >= divisor_high) {
                    // 商は最大値 (B-1)
                    q_hat = UINT64_MAX;

                    // 正しい r_hat の計算:
                    // r_hat = (dividend_high * B + dividend_mid) - (B-1) * divisor_high
                    //       = (dividend_high - divisor_high) * B + dividend_mid + divisor_high
                    //
                    // dividend_high > divisor_high の場合: r_hat >= B → 64ビットに収まらない
                    // → trial quotient correction は不要（r_hat >= B なら q_hat*v2 < r_hat*B は常に成立）
                    if (dividend_high > divisor_high) {
                        // r_hat >= B → correction ループをスキップ
                        r_hat = 0;
                        skip_correction = true;
                    } else {
                        // dividend_high == divisor_high
                        // r_hat = dividend_mid + divisor_high
                        r_hat = dividend_mid + divisor_high;

                        // r_hat がオーバーフローした場合（r_hat >= B）も correction 不要
                        if (r_hat < dividend_mid) {
                            skip_correction = true;
                        }
                    }
                }
                else {
                    // 2ワード÷1ワードの除算
                    // divmod_fast は overflow フラグで正しく処理する
                    // (UInt128::divide は中間変数 r の 64-bit overflow バグあり)
                    auto div_result = UInt128::divmod_fast(dividend_high, dividend_mid, divisor_high);
                    q_hat = div_result.first;
                    r_hat = div_result.second;
                }

                // 商の修正（2ワード目が存在する場合）
                // skip_correction: r_hat >= B（64ビットに収まらない）の場合は修正不要
                if (n >= 2 && !skip_correction) {
                    uint64_t divisor_next = V[n - 2];
                    uint64_t dividend_low = j + n - 2 < U_work.size() ? U_work[j + n - 2] : 0;

                    // 修正ループ
                    bool need_adjustment = true;
                    int adjustment_count = 0;
                    while (need_adjustment) {
                        // q_hat * divisor_next > r_hat * 2^64 + dividend_low を確認
                        UInt128 qd_product = UInt128::multiply(q_hat, divisor_next);
                        UInt128 r_dividend(r_hat, dividend_low);

                        if (qd_product.high < r_dividend.high ||
                            (qd_product.high == r_dividend.high && qd_product.low <= r_dividend.low)) {
                            need_adjustment = false;
                        }
                        else {
                            q_hat--;

                            uint64_t old_r_hat = r_hat;
                            r_hat += divisor_high;

                            // r_hatがオーバーフローしたら終了
                            if (r_hat < old_r_hat) {
                                need_adjustment = false;
                            }

                            adjustment_count++;
                            if (adjustment_count > 2) {
                                need_adjustment = false;
                            }
                        }
                    }
                }

                // q_hat * divisor を被除数から引く (single-pass, Knuth Algorithm D / GMP style)
                // 旧コードの product.low + borrow オーバーフローバグを修正
                uint64_t mul_carry = 0;

                for (size_t i = 0; i < n; i++) {
                    // 境界チェック
                    size_t idx = j + i;
                    if (idx >= U_work.size()) {
                        throw std::out_of_range("U_work index out of bounds in mul-sub");
                    }

                    // q_hat * V[i] を計算
                    UInt128 product = UInt128::multiply(q_hat, V[i]);

                    // product.low + mul_carry（前回からの繰り上がり）
                    uint64_t sum_low = product.low + mul_carry;
                    uint64_t carry_from_add = (sum_low < product.low) ? 1 : 0;

                    // U_work[j+i] から sum_low を引く
                    uint64_t prev = U_work[j + i];
                    U_work[j + i] = prev - sum_low;
                    uint64_t borrow_from_sub = (prev < sum_low) ? 1 : 0;

                    // 次の繰り上がり = product.high + carry_from_add + borrow_from_sub
                    // 注: product.high <= 2^64-2（64bit×64bit乗算の最大）
                    // carry_from_add + borrow_from_sub <= 2 だが、
                    // 両方 1 になるケースでは product.high < 2^64-2 なので
                    // 合計は常に 2^64-1 以下（オーバーフローしない）
                    mul_carry = product.high + carry_from_add + borrow_from_sub;
                }

                // 最上位ワードから残りの carry を引く
                {
                    uint64_t prev = U_work[j + n];
                    U_work[j + n] = prev - mul_carry;
                    // prev < mul_carry なら結果が負（add-back が必要）
                    mul_carry = (prev < mul_carry) ? 1 : 0;
                }

                // 引き算結果が負になった場合の調整（q_hat が大きすぎた）
                if (mul_carry != 0) {
                    // q_hatを1減らし、divisorを足し戻す
                    q_hat--;

                    uint64_t carry = 0;

                    for (size_t i = 0; i < n; i++) {
                        uint64_t u = U_work[j + i];
                        uint64_t v = V[i];

                        // 2段階加算でキャリーを正確に追跡
                        uint64_t sum = u + v;
                        uint64_t c1 = (sum < u) ? 1ULL : 0ULL;
                        uint64_t sum2 = sum + carry;
                        uint64_t c2 = (sum2 < sum) ? 1ULL : 0ULL;
                        U_work[j + i] = sum2;
                        carry = c1 + c2;
                    }

                    // 最上位桁のキャリーを加算
                    U_work[j + n] += carry;
                }

                // 商を設定
                Q[j] = q_hat;
            }

            // 余りを設定（U_workの先頭n桁）
            std::vector<uint64_t> R;
            for (size_t i = 0; i < n && i < U_work.size(); i++) {
                R.push_back(U_work[i]);
            }

            // 正規化シフト前の余りを復元
            Int normalizedRemainder = Int::fromRawWords(R, +1);
            normalizedRemainder.normalize();

            if (normalization_shift > 0) {
                remainder = normalizedRemainder >> normalization_shift;
            }
            else {
                remainder = normalizedRemainder;
            }

            // 商を構築
            Int quotient = Int::fromRawWords(Q, +1);
            quotient.normalize();

            // 符号の設定
            if (remainder_negative && !remainder.isZero()) {
                remainder.negate();
            }

            if (quotient_negative && !quotient.isZero()) {
                quotient.negate();
            }

            return quotient;
        }

        /**
         * @brief 除算の実装（商のみ）
         * @param dividend 被除数
         * @param divisor 除数
         * @return 商
         */
        static Int divide(const Int& dividend, const Int& divisor) {
            Int remainder;
            return divKnuth(dividend, divisor, remainder);
        }

        /**
         * @brief 剰余の実装
         * @param dividend 被除数
         * @param divisor 除数
         * @return 余り
         */
        static Int modulo(const Int& dividend, const Int& divisor) {
            Int remainder;
            divKnuth(dividend, divisor, remainder);
            return remainder;
        }

        /**
         * @brief 2のべき乗で除算する高速アルゴリズム
         *
         * 2のべき乗による除算は右シフトで高速に処理できる
         *
         * @param dividend 被除数
         * @param power 2の累乗値（2^power）
         * @param remainder 余りの格納先
         * @return 商
         */
        static Int divPowerOfTwo(const Int& dividend, size_t power, Int& remainder) {
            if (power == 0) {
                // 2^0 = 1 で除算する場合
                remainder = Int(0);
                return dividend;
            }

            // 符号の保存
            int sign = dividend.getSign();

            // 特殊状態の処理
            if (dividend.getState() != NumericState::Normal) {
                if (dividend.isNaN() || NumericStateTraits::isInfinite(dividend.getState())) {
                    remainder = Int::fromState(NumericState::NaN, NumericError::NaNPropagation);
                    return dividend; // NaNや無限大はそのまま
                }
            }

            // 絶対値を取得
            Int abs_dividend = abs(dividend);

            // 余りを計算（下位 power ビットを抽出）
            size_t word_shift = power / 64;
            size_t bit_shift = power % 64;

            std::vector<uint64_t> remainder_words;
            std::vector<uint64_t> words = abs_dividend.words();

            if (!words.empty()) {
                // 完全なワードをコピー
                for (size_t i = 0; i < std::min(words.size(), word_shift); ++i) {
                    remainder_words.push_back(words[i]);
                }
                // 端数ビットのマスク処理
                if (bit_shift > 0 && word_shift < words.size()) {
                    uint64_t mask = (1ULL << bit_shift) - 1;
                    remainder_words.push_back(words[word_shift] & mask);
                }
            }

            remainder = Int::fromRawWords(remainder_words, sign);

            // 商を計算（右シフト）
            Int quotient = abs_dividend >> static_cast<int>(power);

            // 符号の調整
            if (sign < 0 && !quotient.isZero()) {
                    quotient.negate();
            }

            return quotient;
        }

        /**
         * @brief 10のべき乗で除算する高速アルゴリズム
         *
         * 10進数表現での桁移動の際に使用
         *
         * @param dividend 被除数
         * @param power 10の累乗値（10^power）
         * @param remainder 余りの格納先
         * @return 商
         */
        static Int divPowerOfTen(const Int& dividend, size_t power, Int& remainder) {
            // 10^power の計算
            Int divisor(10);

            // unsigned int版のpow関数を使用
            Int ten_power = calx::pow(divisor, static_cast<unsigned int>(power));

            // 通常の除算を使用
            return divKnuth(dividend, ten_power, remainder);
        }

        /**
         * @brief ニュートン・ラフソン法による逆数の近似計算
         *
         * doubleの精度を使って初期値を求め、Newton-Raphson反復で精度を向上させる。
         * 巨大な整数でも、指数部を調整してdouble範囲に収めてから計算。
         *
         * アルゴリズム:
         * 1. divisor を double 範囲に正規化（ビットシフトで調整）
         * 2. double精度で 1/divisor の初期近似値を計算
         * 3. Newton-Raphson反復: X_{n+1} = X_n * (2 - D * X_n / 2^k)
         * 4. 固定小数点表現で結果を返す: floor(2^scale / divisor)
         *
         * @param divisor 除数（正の整数）
         * @param scale スケールファクター（結果は 2^scale / divisor）
         * @return floor(2^scale / divisor) - 固定小数点での逆数
         */
        static Int newtonRaphsonReciprocal(const Int& divisor, size_t scale) {
            // 特殊ケースの処理
            if (divisor.isZero()) {
                return Int::fromState(NumericState::NaN, NumericError::DivideByZero);
            }

            if (divisor.isNaN() || NumericStateTraits::isInfinite(divisor.getState())) {
                return Int::fromState(NumericState::NaN, NumericError::NaNPropagation);
            }

            size_t divisor_bits = divisor.bitLength();

            // 小さい除数は直接計算
            if (divisor_bits <= 64) {
                Int numerator = Int(1) << static_cast<int>(scale);
                return numerator / divisor;
            }

            // ステップ1: divisor を double 範囲に正規化
            // double の仮数部は53ビット
            const size_t DOUBLE_MANTISSA_BITS = 53;

            // divisor の最上位53ビット付近を取り出す
            int shift_down = static_cast<int>(divisor_bits) - static_cast<int>(DOUBLE_MANTISSA_BITS);
            if (shift_down < 0) shift_down = 0;

            Int divisor_normalized = divisor >> shift_down;

            // double に変換
            double d_approx;
            if (divisor_normalized.size() >= 2) {
                uint64_t high = divisor_normalized.word(divisor_normalized.size() - 1);
                uint64_t low = divisor_normalized.word(divisor_normalized.size() - 2);
                // high * 2^64 + low を double で近似
                d_approx = static_cast<double>(high) * 18446744073709551616.0 + static_cast<double>(low);
            } else if (divisor_normalized.size() == 1) {
                d_approx = static_cast<double>(divisor_normalized.word(0));
            } else {
                d_approx = 1.0;
            }

            // ステップ2: double精度で初期逆数を計算
            double recip_double = 1.0 / d_approx;

            // recip_double を固定小数点 Int に変換
            // recip_double は divisor_normalized の逆数
            // divisor_normalized = divisor >> shift_down なので
            // 1/divisor_normalized = (1/divisor) * 2^shift_down
            // つまり recip_double ≈ (1/divisor) * 2^shift_down

            // working_scaleはscaleと同じにする
            // 逆数Xは2^working_scale/D ≈ 2^(scale-divisor_bits)ビットの精度が必要
            // working_scaleが小さすぎるとXの整数表現のビット数が不足し、
            // 商の精度が足りなくなる（修正ループが無限に回る）
            const size_t working_scale = scale;

            // 初期近似: X_0 = floor(2^working_scale / divisor)
            // recip_double ≈ (1/divisor) * 2^shift_down = 2^shift_down / divisor
            // We want X_0 = floor(2^working_scale / divisor)
            // X_0 = recip_double * 2^(working_scale - shift_down)

            // Extract mantissa and exponent from recip_double using frexp
            int exponent;
            double mantissa = std::frexp(recip_double, &exponent);
            // Now recip_double = mantissa * 2^exponent, where 0.5 <= |mantissa| < 1.0

            // Convert mantissa to fixed-point (multiply by 2^53 to get all 53 bits)
            uint64_t recip_int = static_cast<uint64_t>(mantissa * (1ULL << 53));
            Int X = Int(recip_int);

            // Total shift = (working_scale - shift_down) + exponent - 53
            int total_shift = static_cast<int>(working_scale) - shift_down + exponent - 53;
            if (total_shift > 0) {
                X = X << total_shift;
            } else if (total_shift < 0) {
                X = X >> (-total_shift);
            }

            // ステップ3: Newton-Raphson反復で精度を倍々に向上
            // 反復式: X_{n+1} = X_n * (2*2^k - D * X_n) / 2^k
            // 各反復で精度ビット数が倍増する（二次収束）
            //
            // 初期精度: doubleの53ビット仮数部
            // 必要精度: working_scale ビット（diff < D となるまで）

            // 反復回数を計算: 53 → 106 → 212 → ... → working_scale
            size_t initial_precision = 50;  // double精度から保守的に見積もり
            size_t iterations_needed = 0;
            {
                size_t current = initial_precision;
                while (current < working_scale) {
                    current *= 2;
                    iterations_needed++;
                }
            }

            Int R = Int(1) << static_cast<int>(working_scale);

            for (size_t iter = 0; iter < iterations_needed; ++iter) {
                // D * X_n を計算（スケール2^k）
                Int DX = divisor * X;

                // 2 * 2^k
                Int two_R = R << 1;

                // 補正: E = 2*2^k - D*X_n
                Int E = two_R - DX;

                // X_{n+1} = X_n * E / 2^k
                X = (X * E) >> static_cast<int>(working_scale);
            }

            // working_scale == scale なのでスケール調整は不要
            return X;
        }

        /**
         * @brief ニュートン法による除算（超大型数向け高速アルゴリズム）
         *
         * ニュートン・ラフソン法で 1/divisor を計算し、それを dividend に掛けることで除算を行う。
         * 計算量は O(M(n)) where M(n) is multiplication time
         * (Karatsuba なら O(n^1.585), FFT なら O(n log n))
         *
         * 参考: "Modern Computer Arithmetic" Algorithm 1.6
         *
         * @param dividend 被除数
         * @param divisor 除数
         * @param remainder 余りの格納先
         * @return 商
         */
        static Int divNewton(const Int& dividend, const Int& divisor, Int& remainder) {
            // 特殊ケースの処理
            if (divisor.isZero()) {
                remainder = Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                return Int::fromState(NumericState::NaN, NumericError::DivideByZero);
            }

            if (dividend.isZero()) {
                remainder = Int(0);
                return Int(0);
            }

            if (dividend.isNaN() || divisor.isNaN() ||
                NumericStateTraits::isInfinite(dividend.getState()) ||
                NumericStateTraits::isInfinite(divisor.getState())) {
                remainder = Int::fromState(NumericState::NaN, NumericError::NaNPropagation);
                return Int::fromState(NumericState::NaN, NumericError::NaNPropagation);
            }

            // 符号を保存
            int dividend_sign = dividend.getSign();
            int divisor_sign = divisor.getSign();
            int quotient_sign = dividend_sign * divisor_sign;

            // 絶対値で計算
            Int abs_dividend = abs(dividend);
            Int abs_divisor = abs(divisor);

            // 除数が被除数より大きい場合
            if (abs_dividend < abs_divisor) {
                remainder = dividend;  // 符号を保持
                return Int(0);
            }

            // Newton除算アルゴリズム
            // Q = floor(A / D) を求める
            // 1. R = 1/D を計算（固定小数点: 2^k / D）
            // 2. Q_approx = A * R / 2^k
            // 3. 余りを計算して商を補正

            size_t dividend_bits = abs_dividend.bitLength();
            size_t divisor_bits = abs_divisor.bitLength();

            // スケールを決定: dividend_bits + 余裕
            // 逆数 R = 2^scale / D が quotient_bits ≈ dividend_bits - divisor_bits の
            // 精度を持つためには、scale >= dividend_bits + margin が必要
            // scale = divisor_bits + 64 では精度不足で補正ループが無限に回る
            size_t scale = dividend_bits + 64;

            // ステップ1: 1/divisor を計算
            Int reciprocal = newtonRaphsonReciprocal(abs_divisor, scale);

            // ステップ2: 商の近似を計算
            // Q ≈ A * (2^scale / D) / 2^scale = A * reciprocal / 2^scale
            Int quotient = (abs_dividend * reciprocal) >> static_cast<int>(scale);

            // ステップ3: 余りを計算
            Int rem = abs_dividend - (quotient * abs_divisor);

            // 余りが負の場合、商を減らす（通常は1-2回の調整で済む）
            while (rem.getSign() < 0) {
                quotient -= 1;
                rem = rem + abs_divisor;
            }

            // 余りが除数以上の場合、商を増やす（通常は1-2回の調整で済む）
            while (rem >= abs_divisor) {
                quotient += 1;
                rem = rem - abs_divisor;
            }

            // 符号を設定
            if (quotient_sign < 0 && !quotient.isZero()) {
                quotient.negate();
            }
            if (dividend_sign < 0 && !rem.isZero()) {
                rem.negate();
            }

            remainder = rem;
            return quotient;
        }

    private:
        // 内部実装ヘルパー関数
    };

} // namespace calx

#endif // CALX_INT_DIVISION_HPP