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

// IntOps.cpp
// 多倍長整数の操作に関するユーティリティの実装

#include <math/core/mp/Int/IntOps.hpp>
#include <math/core/mp/Int/IntSpecialStates.hpp>
#include <math/core/mp/Int/IntDivision.hpp>
#include <math/core/mp/Int/IntMultiplication.hpp>
#include <math/core/mp/Int/IntSqrt.hpp>
#include <math/core/mp/Int/IntGCD.hpp>
#include <algorithm>
#include <stdexcept>

// C++20 <bit> により __builtin_clzll 互換マクロは不要になった
// std::countl_zero / std::countr_zero を使用

namespace calx {

    // 比較操作
    bool IntOps::compareAbsLess(const Int& lhs, const Int& rhs) {
        // 特殊状態の処理
        if (lhs.isSpecialState() || rhs.isSpecialState()) {
            return IntSpecialStates::compareLessThan(
                lhs.getSign() < 0 ? -lhs : lhs,
                rhs.getSign() < 0 ? -rhs : rhs
            );
        }

        // ワード数による比較
        if (lhs.m_words.size() != rhs.m_words.size()) {
            return lhs.m_words.size() < rhs.m_words.size();
        }

        // 同じワード数の場合は上位ワードから比較
        for (int i = static_cast<int>(lhs.m_words.size()) - 1; i >= 0; --i) {
            if (lhs.m_words[i] != rhs.m_words[i]) {
                return lhs.m_words[i] < rhs.m_words[i];
            }
        }

        // 完全に同じ
        return false;
    }

    bool IntOps::compareAbsGreater(const Int& lhs, const Int& rhs) {
        // 特殊状態の処理
        if (lhs.isSpecialState() || rhs.isSpecialState()) {
            return IntSpecialStates::compareLessThan(
                rhs.getSign() < 0 ? -rhs : rhs,
                lhs.getSign() < 0 ? -lhs : lhs
            );
        }

        // ワード数による比較
        if (lhs.m_words.size() != rhs.m_words.size()) {
            return lhs.m_words.size() > rhs.m_words.size();
        }

        // 同じワード数の場合は上位ワードから比較
        for (int i = static_cast<int>(lhs.m_words.size()) - 1; i >= 0; --i) {
            if (lhs.m_words[i] != rhs.m_words[i]) {
                return lhs.m_words[i] > rhs.m_words[i];
            }
        }

        // 完全に同じ
        return false;
    }

    bool IntOps::compareAbsEqual(const Int& lhs, const Int& rhs) {
        // 特殊状態の処理
        if (lhs.isSpecialState() || rhs.isSpecialState()) {
            return lhs.getState() == rhs.getState();
        }

        // ワード配列の直接比較
        return lhs.m_words == rhs.m_words;
    }

    // 加算操作
    void IntOps::addAbsolute(Int& result, const Int& other) {
        // 特殊状態の処理
        if (result.isSpecialState() || other.isSpecialState()) {
            result = IntSpecialStates::handleAddition(result, other);
            return;
        }

        // ゼロとの加算の特別処理
        if (other.getSign() == 0) {
            return;
        }
        if (result.getSign() == 0) {
            result.m_words = other.m_words;
            result.setSign(other.getSign());
            return;
        }

        // 加算処理
        size_t resultSize = result.m_words.size();
        size_t otherSize = other.m_words.size();
        size_t maxSize = std::max(resultSize, otherSize);

        // 結果のベクトルサイズを拡張
        result.m_words.resize(maxSize, 0);

        uint64_t carry = mpn::add(result.m_words.data(),
                                   result.m_words.data(), maxSize,
                                   other.m_words.data(), otherSize);
        if (carry) {
            result.m_words.push_back(carry);
        }
        // normalize 不要: 両入力は正規化済み → 最上位ワードは非ゼロ
    }

    // 3引数版 addAbsolute/subtractAbsolute は IntOps.hpp でインライン定義済み

    // インクリメント/デクリメント: mpn::add_1/sub_1 で Int(1) 構築を回避
    void IntOps::addDelta(Int& value, int delta) {
        if (value.m_sign == 0) {
            value.m_words.resize_uninitialized(1);
            value.m_words[0] = 1;
            value.m_sign = delta;
            value.m_state = NumericState::Normal;
            return;
        }
        // 正の delta と正の値 (または負の delta と負の値): 絶対値を +1
        if ((value.m_sign > 0) == (delta > 0)) {
            uint64_t carry = mpn::add_1(value.m_words.data(), value.m_words.size(), 1);
            if (carry) value.m_words.push_back(carry);
        } else {
            // 異符号: 絶対値を -1
            if (value.m_words.size() == 1 && value.m_words[0] == 1) {
                // |value| == 1 → 結果は 0
                value.m_words.clear();
                value.m_sign = 0;
                return;
            }
            mpn::sub_1(value.m_words.data(), value.m_words.size(), 1);
            if (value.m_words.back() == 0) value.m_words.pop_back();
        }
    }

    // 単一ワード加算: result += word (符号付き)
    void IntOps::addWord(Int& result, uint64_t word) {
        if (word == 0) return;
        if (result.m_sign == 0) {
            result.m_words.resize_uninitialized(1);
            result.m_words[0] = word;
            result.m_sign = 1;
            return;
        }
        if (result.m_sign > 0) {
            uint64_t carry = mpn::add_1(result.m_words.data(), result.m_words.size(), word);
            if (carry) result.m_words.push_back(carry);
        } else {
            // result < 0: -(|result| - word)
            if (result.m_words.size() == 1) {
                uint64_t v = result.m_words[0];
                if (v <= word) {
                    uint64_t diff = word - v;
                    if (diff == 0) { result.m_words.clear(); result.m_sign = 0; }
                    else { result.m_words[0] = diff; result.m_sign = 1; }
                    return;
                }
            }
            mpn::sub_1(result.m_words.data(), result.m_words.size(), word);
            if (result.m_words.back() == 0) result.m_words.pop_back();
            if (result.m_words.empty()) result.m_sign = 0;
        }
    }

    // 単一ワード減算: result -= word (符号付き)
    void IntOps::subWord(Int& result, uint64_t word) {
        if (word == 0) return;
        if (result.m_sign == 0) {
            result.m_words.resize_uninitialized(1);
            result.m_words[0] = word;
            result.m_sign = -1;
            return;
        }
        if (result.m_sign < 0) {
            // result < 0: -(|result| + word)
            uint64_t carry = mpn::add_1(result.m_words.data(), result.m_words.size(), word);
            if (carry) result.m_words.push_back(carry);
        } else {
            // result > 0: |result| - word
            if (result.m_words.size() == 1) {
                uint64_t v = result.m_words[0];
                if (v <= word) {
                    uint64_t diff = word - v;
                    if (diff == 0) { result.m_words.clear(); result.m_sign = 0; }
                    else { result.m_words[0] = diff; result.m_sign = -1; }
                    return;
                }
            }
            mpn::sub_1(result.m_words.data(), result.m_words.size(), word);
            if (result.m_words.back() == 0) result.m_words.pop_back();
            if (result.m_words.empty()) result.m_sign = 0;
        }
    }

    // 単一ワード除算: result /= word, 余りを返す
    uint64_t IntOps::divWord(Int& result, uint64_t word) {
        if (result.m_sign == 0) return 0;

        // 2のべき乗による除算 → 右シフト (divmod_1 より高速)
        if (word != 0 && (word & (word - 1)) == 0) {
            unsigned long shift;
            _BitScanForward64(&shift, word);
            // 余りは下位 shift ビット
            uint64_t rem = 0;
            if (shift > 0 && !result.m_words.empty()) {
                uint64_t mask = (uint64_t(1) << shift) - 1;
                rem = result.m_words[0] & mask;
            }
            rightShift(result, static_cast<int>(shift));
            return rem;
        }

        uint64_t rem = mpn::divmod_1(result.m_words.data(), result.m_words.data(),
                                      result.m_words.size(), word);
        while (!result.m_words.empty() && result.m_words.back() == 0)
            result.m_words.pop_back();
        if (result.m_words.empty()) result.m_sign = 0;
        return rem;
    }

    // 減算操作 (in-place: result -= |other|, 前提: |result| >= |other|)
    void IntOps::subtractAbsoluteInPlace(Int& result, const Int& other) {
        size_t resultSize = result.m_words.size();
        size_t otherSize = other.m_words.size();

        mpn::sub(result.m_words.data(),
                 result.m_words.data(), resultSize,
                 other.m_words.data(), otherSize);

        result.normalize();
    }

    // 乗算操作
    void IntOps::multiplyAbsolute(Int& result, const Int& other) {
        // 特殊状態の処理
        if (result.isSpecialState() || other.isSpecialState()) {
            result = IntSpecialStates::handleMultiplication(result, other);
            return;
        }

        // ゼロとの乗算の特別処理
        if (result.getSign() == 0 || other.getSign() == 0) {
            result.m_words.clear();
            result.setSign(0);
            return;
        }

        // 1との乗算の特別処理
        if (other.m_words.size() == 1 && other.m_words[0] == 1) {
            // 符号は呼び出し元で処理される
            return;
        }
        if (result.m_words.size() == 1 && result.m_words[0] == 1) {
            result.m_words = other.m_words;
            // 符号は呼び出し元で処理される
            return;
        }

        size_t an = result.m_words.size();
        size_t bn = other.m_words.size();

        // Small-size fast path: stack バッファで中間 Int 生成を回避
        // (ScratchScope, arena alloc, fromRawWords のオーバーヘッドを除去)
        if (an < mpn::KARATSUBA_THRESHOLD && bn < mpn::KARATSUBA_THRESHOLD) {
            size_t rn = an + bn;
            uint64_t rbuf[mpn::KARATSUBA_THRESHOLD * 2];
            mpn::multiply(rbuf, result.m_words.data(), an,
                          other.m_words.data(), bn, nullptr);
            // 乗算の積は最大 rn limbs、上位 0 は最上位 1 ワードのみ
            if (rbuf[rn - 1] == 0) --rn;
            result.m_words.assign(rbuf, rbuf + rn);
            return;
        }

        // 大サイズ: IntMultiplication 経由 (arena + Karatsuba/Toom/NTT)
        generalMultiply(result, other, result);
    }

    // 乗算操作 (3引数版: lhs, rhs を直接参照しコピー不要)
    // 前提: ゼロ・特殊状態は呼び出し元でチェック済み
    // 注意: 符号・状態は呼び出し元が設定する (オーバーヘッド削減)
    void IntOps::multiplyAbsolute(const Int& lhs, const Int& rhs, Int& result) {
        // ×1 の特別処理
        if (rhs.m_words.size() == 1 && rhs.m_words[0] == 1) {
            result.m_words = lhs.m_words;
            return;
        }
        if (lhs.m_words.size() == 1 && lhs.m_words[0] == 1) {
            result.m_words = rhs.m_words;
            return;
        }

        size_t an = lhs.m_words.size();
        size_t bn = rhs.m_words.size();

        // 1×1 特化: 128-bit 乗算を直接実行 (関数呼び出しチェーン全体を省略)
        if (an == 1 && bn == 1) {
            UInt128 prod = UInt128::multiply(lhs.m_words[0], rhs.m_words[0]);
            if (prod.high != 0) {
                result.m_words.resize_uninitialized(2);
                result.m_words[0] = prod.low;
                result.m_words[1] = prod.high;
            } else {
                result.m_words.resize_uninitialized(1);
                result.m_words[0] = prod.low;
            }
            return;
        }

        // N×1 特化: mul_1 直接呼び出し (multiply→mul_basecase のディスパッチを省略)
        if (an == 1 || bn == 1) {
            const uint64_t* long_data;
            size_t long_n;
            uint64_t short_val;
            if (bn == 1) {
                long_data = lhs.m_words.data(); long_n = an; short_val = rhs.m_words[0];
            } else {
                long_data = rhs.m_words.data(); long_n = bn; short_val = lhs.m_words[0];
            }
            result.m_words.resize_uninitialized(long_n + 1);
            result.m_words[long_n] = mpn::mul_1(result.m_words.data(), long_data, long_n, short_val);
            if (result.m_words.back() == 0) result.m_words.pop_back();
            return;
        }

        // Small-size: result.m_words に直接書き込み (中間コピー除去)
        if (an < mpn::KARATSUBA_THRESHOLD && bn < mpn::KARATSUBA_THRESHOLD) {
            size_t rn = an + bn;
            result.m_words.resize_uninitialized(rn);
            mpn::multiply(result.m_words.data(), lhs.m_words.data(), an,
                          rhs.m_words.data(), bn, nullptr);
            if (result.m_words.back() == 0) result.m_words.pop_back();
            return;
        }

        // Large-size: IntMultiplication 経由
        generalMultiply(lhs, rhs, result);
    }

    // 単一ワードとの乗算
    void IntOps::multiplyWord(Int& result, uint64_t word) {
        // 特殊状態の処理
        if (result.isSpecialState()) {
            return;
        }

        // ゼロとの乗算の特別処理
        if (result.getSign() == 0 || word == 0) {
            result.m_words.clear();
            result.setSign(0);
            return;
        }

        // 1との乗算の特別処理
        if (word == 1) {
            return;
        }

        // 2のべき乗との乗算 → 左シフト (mpn::mul_1 より高速)
        if ((word & (word - 1)) == 0) {
            // word は 2 のべき乗 (word >= 2 は上の word==0, word==1 チェックで保証)
            unsigned long shift;
            _BitScanForward64(&shift, word);
            leftShift(result, static_cast<int>(shift));
            return;
        }

        // mpn::mul_1 で単一ワード乗算 (MASM 4x アンロール対応)
        uint64_t carry = mpn::mul_1(result.m_words.data(),
                                     result.m_words.data(),
                                     result.m_words.size(), word);
        if (carry > 0) {
            result.m_words.push_back(carry);
        }

        result.normalize();
    }

    // 3引数版 mul: result のバッファを再利用 (GMP mpz_mul 相当)
    void IntOps::mulUnchecked(const Int& lhs, const Int& rhs, Int& result) {
        if (lhs.getSign() == 0 || rhs.getSign() == 0) [[unlikely]] {
            result.m_words.clear();
            result.m_sign = 0;
            result.m_state = NumericState::Normal;
            return;
        }

        int resultSign = lhs.getSign() * rhs.getSign();

        if (&lhs == &rhs) {
            square(lhs, result);
        } else {
            multiplyAbsolute(lhs, rhs, result);
        }
        result.m_sign = resultSign;
        result.m_state = NumericState::Normal;
    }

    void IntOps::mul(const Int& lhs, const Int& rhs, Int& result) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] {
            result = lhs * rhs;
            return;
        }
        mulUnchecked(lhs, rhs, result);
    }

    // 平方計算
    void IntOps::square(const Int& value, Int& result) {
        // 特殊状態の処理
        if (value.isSpecialState()) {
            result = IntSpecialStates::handleMultiplication(value, value);
            return;
        }

        // ゼロの平方の特別処理
        if (value.getSign() == 0) {
            result.m_words.clear();
            result.m_sign = 0;
            result.m_state = NumericState::Normal;
            return;
        }

        size_t n = value.m_words.size();

        // 小サイズ: ASM mul_basecase のほうが ASM sqr_basecase より速い
        // (sqr の 3-pass 構造 + zeroing のオーバーヘッドが n(n-1)/2 の削減を相殺)
        // crossover: n=8 MUL 22%有利, n=9 TIE, n=11 SQR 13%有利 (bench_sqr_fallback.cpp)
        constexpr size_t SQR_MUL_FALLBACK = 8;
        if (n <= SQR_MUL_FALLBACK) {
            multiplyAbsolute(value, value, result);
            result.m_sign = 1;
            result.m_state = NumericState::Normal;
            return;
        }

        // 大サイズ: mpn::square による専用自乗 (対称性利用で高速化)
        result.m_words.resize_uninitialized(2 * n);

        size_t scratch_sz = mpn::square_scratch_size(n);
        if (scratch_sz > 0) {
            ScratchScope scope;
            uint64_t* scratch = getThreadArena().alloc_limbs(scratch_sz);
            mpn::square(result.m_words.data(), value.m_words.data(), n, scratch);
        } else {
            mpn::square(result.m_words.data(), value.m_words.data(), n, nullptr);
        }

        // 自乗の結果: 最上位 word が 0 なら除去 (normalize ループ不要)
        if (result.m_words.back() == 0) result.m_words.pop_back();
        result.m_sign = 1;
        result.m_state = NumericState::Normal;
    }

    // 3引数版 divUnchecked: isSpecialState チェック省略版
    void IntOps::divUnchecked(const Int& dividend, const Int& divisor, Int& result) {
        // ゼロ除算 (特殊状態ではないが安全のため保持)
        if (divisor.m_sign == 0) [[unlikely]] {
            result = dividend / divisor;
            return;
        }

        // 被除数がゼロ
        if (dividend.m_sign == 0) {
            result.m_words.clear();
            result.m_sign = 0;
            result.m_state = NumericState::Normal;
            return;
        }

        size_t an = dividend.m_words.size();
        size_t bn = divisor.m_words.size();

        // |dividend| < |divisor| → 商 = 0
        int abscmp = mpn::cmp(dividend.m_words.data(), an, divisor.m_words.data(), bn);
        if (abscmp < 0) {
            result.m_words.clear();
            result.m_sign = 0;
            result.m_state = NumericState::Normal;
            return;
        }
        if (abscmp == 0) {
            result.m_words.resize(1);
            result.m_words[0] = 1;
            result.m_sign = dividend.m_sign * divisor.m_sign;
            result.m_state = NumericState::Normal;
            return;
        }

        int resultSign = dividend.m_sign * divisor.m_sign;

        // 1-limb 除数: divmod_1 fast path
        if (bn == 1) {
            size_t qn = an;
            result.m_words.resize_uninitialized(qn);
            mpn::divmod_1(result.m_words.data(), dividend.m_words.data(), an, divisor.m_words[0]);
            while (qn > 0 && result.m_words[qn - 1] == 0) qn--;
            result.m_words.resize(qn);
            result.m_sign = (qn == 0) ? 0 : resultSign;
            result.m_state = NumericState::Normal;
            return;
        }

        // マルチワード: mpn::divide を直接呼ぶ
        {
            size_t scratch_sz = mpn::divide_scratch_size(an, bn);
            size_t qn = an - bn + 1;

            // スタックバッファで小サイズの heap 確保を回避
            constexpr size_t STACK_LIMIT = 128;
            uint64_t q_stack[STACK_LIMIT], r_stack[STACK_LIMIT], scratch_stack[STACK_LIMIT * 8];
            uint64_t* qp;
            uint64_t* rp;
            uint64_t* scratch;

            if (qn <= STACK_LIMIT && bn <= STACK_LIMIT && scratch_sz <= STACK_LIMIT * 8) {
                qp = q_stack;
                rp = r_stack;
                scratch = scratch_stack;
            } else {
                qp = new uint64_t[qn + 1];
                rp = new uint64_t[bn];
                scratch = new uint64_t[scratch_sz];
            }

            mpn::divide(qp, rp, dividend.m_words.data(), an,
                        divisor.m_words.data(), bn, scratch);

            while (qn > 0 && qp[qn - 1] == 0) qn--;
            result.m_words.assign(qp, qp + qn);
            result.m_sign = (qn == 0) ? 0 : resultSign;
            result.m_state = NumericState::Normal;

            if (qp != q_stack) { delete[] qp; delete[] rp; delete[] scratch; }
        }
    }

    // 3引数版 div: チェック付きラッパー
    void IntOps::div(const Int& dividend, const Int& divisor, Int& result) {
        if (dividend.isSpecialState() || divisor.isSpecialState()) [[unlikely]] {
            result = dividend / divisor;
            return;
        }
        divUnchecked(dividend, divisor, result);
    }

    // ================================================================
    // addmul / submul: rop ±= a * b
    // ================================================================

    // rop += a * b
    void IntOps::addmul(Int& rop, const Int& a, const Int& b) {
        // 特殊状態
        if (a.isSpecialState() || b.isSpecialState() || rop.isSpecialState()) [[unlikely]] {
            rop = rop + a * b;
            return;
        }
        // a または b がゼロなら何もしない
        if (a.getSign() == 0 || b.getSign() == 0) return;

        // rop がゼロなら単純に rop = a * b
        if (rop.getSign() == 0) {
            mul(a, b, rop);
            return;
        }

        // a * b を一時バッファに計算
        Int product;
        mul(a, b, product);

        // rop += product (3引数 add でバッファ再利用)
        add(rop, product, rop);
    }

    // rop -= a * b
    void IntOps::submul(Int& rop, const Int& a, const Int& b) {
        if (a.isSpecialState() || b.isSpecialState() || rop.isSpecialState()) [[unlikely]] {
            rop = rop - a * b;
            return;
        }
        if (a.getSign() == 0 || b.getSign() == 0) return;

        if (rop.getSign() == 0) {
            mul(a, b, rop);
            rop.m_sign = -rop.m_sign;
            return;
        }

        Int product;
        mul(a, b, product);
        sub(rop, product, rop);
    }

    // rop += a * word (単一ワード版 — mpn_addmul_1 を直接活用)
    void IntOps::addmul(Int& rop, const Int& a, uint64_t word) {
        if (a.isSpecialState() || rop.isSpecialState()) [[unlikely]] {
            rop = rop + a * Int(word);
            return;
        }
        if (a.getSign() == 0 || word == 0) return;

        if (rop.getSign() == 0) {
            rop = a;
            multiplyWord(rop, word);
            return;
        }

        // |a| * word を一時バッファに計算
        size_t an = a.m_words.size();
        size_t pn = an + 1;  // 最大ワード数
        constexpr size_t BUF_LIMIT = 128;
        uint64_t stack_buf[BUF_LIMIT];
        uint64_t* pbuf = (pn <= BUF_LIMIT) ? stack_buf : new uint64_t[pn];

        uint64_t carry = mpn::mul_1(pbuf, a.m_words.data(), an, word);
        if (carry) {
            pbuf[an] = carry;
        } else {
            pn = an;
        }

        // product の符号: sign(a) (word は常に正)
        int prodSign = a.getSign();

        // rop ±= product を mpn レベルで実行
        size_t rn = rop.m_words.size();
        // product を Int に変換して add/sub に委譲
        // (mpn レベルでのエイリアシング問題を回避)
        Int product;
        product.m_words.resize_uninitialized(pn);
        std::copy_n(pbuf, pn, product.m_words.data());
        product.m_sign = prodSign;
        product.m_state = NumericState::Normal;

        if (pbuf != stack_buf) delete[] pbuf;

        add(rop, product, rop);
    }

    // rop -= a * word (単一ワード版)
    void IntOps::submul(Int& rop, const Int& a, uint64_t word) {
        if (a.isSpecialState() || rop.isSpecialState()) [[unlikely]] {
            rop = rop - a * Int(word);
            return;
        }
        if (a.getSign() == 0 || word == 0) return;

        // submul を addmul で実装: rop -= a*w == rop += (-a)*w
        // ただし a は const → 符号反転して addmul を呼ぶ代わりに
        // prodSign を反転して直接処理

        if (rop.getSign() == 0) {
            rop = a;
            multiplyWord(rop, word);
            rop.m_sign = -rop.m_sign;
            return;
        }

        // |a| * word を一時バッファに計算
        size_t an = a.m_words.size();
        size_t pn = an + 1;
        constexpr size_t BUF_LIMIT = 128;
        uint64_t stack_buf[BUF_LIMIT];
        uint64_t* pbuf = (pn <= BUF_LIMIT) ? stack_buf : new uint64_t[pn];

        uint64_t carry = mpn::mul_1(pbuf, a.m_words.data(), an, word);
        if (carry) {
            pbuf[an] = carry;
        } else {
            pn = an;
        }

        // product の実効符号 (submul なので反転): -sign(a)
        int prodSign = -a.getSign();

        Int product;
        product.m_words.resize_uninitialized(pn);
        std::copy_n(pbuf, pn, product.m_words.data());
        product.m_sign = prodSign;
        product.m_state = NumericState::Normal;

        if (pbuf != stack_buf) delete[] pbuf;

        add(rop, product, rop);
    }

    // 除算操作
    void IntOps::divideAbsolute(Int& result, const Int& divisor) {
        // 特殊状態の処理
        if (result.isSpecialState() || divisor.isSpecialState()) {
            result = IntSpecialStates::handleDivision(result, divisor);
            return;
        }

        // ゼロによる除算チェック
        if (divisor.getSign() == 0) {
            result.setState(NumericState::NaN, NumericError::DivideByZero);
            return;
        }

        // ゼロの除算
        if (result.getSign() == 0) {
            return; // 0 / x = 0 (x ≠ 0)
        }

        // 同じ値による除算
        if (compareAbsEqual(result, divisor)) {
            result.m_words.resize(1);
            result.m_words[0] = 1;
            // 符号は呼び出し元で処理される
            return;
        }

        // |result| < |divisor| の場合
        if (compareAbsLess(result, divisor)) {
            result.m_words.clear(); // 0
            result.setSign(0);
            return;
        }

        // グローバル設定に基づいてアルゴリズムを選択
        switch (g_division_algorithm) {
        case DivisionAlgorithm::Knuth: {
            // Knuth Algorithm D を使用
            // divKnuth は符号付き Int を扱うので、符号を保存しておく
            int original_sign = result.getSign();
            result.setSign(1);  // 絶対値として扱う
            Int abs_divisor = divisor;
            abs_divisor.setSign(1);

            Int remainder;
            Int quotient = IntDivision::divKnuth(result, abs_divisor, remainder);
            result = std::move(quotient);

            // 符号は呼び出し元で処理される（絶対値除算なので正のまま）
            if (!result.isZero()) {
                result.setSign(1);
            }
            return;
        }

        case DivisionAlgorithm::BitByBit:
            // ビット単位除算へ直接ジャンプ
            goto bit_by_bit_division;

        case DivisionAlgorithm::SingleWordOnly:
            // 単一ワード除数のみ最適化、それ以外はビット単位
            if (divisor.m_words.size() == 1) {
                goto single_word_division;
            } else {
                goto bit_by_bit_division;
            }

        case DivisionAlgorithm::Auto:
        default:
            // 自動選択（ベンチマーク結果に基づく最適化）

            // 1. 2のべき乗検出 → PowerOfTwo (28-168x speedup)
            {
                bool is_power_of_two = true;
                size_t power = 0;
                bool found_nonzero = false;

                // 除数が2のべき乗かチェック
                for (size_t i = 0; i < divisor.m_words.size(); ++i) {
                    if (divisor.m_words[i] != 0) {
                        // 非ゼロワードが複数あれば2のべき乗ではない
                        if (found_nonzero) {
                            is_power_of_two = false;
                            break;
                        }
                        found_nonzero = true;

                        // ビットが1つだけ立っているかチェック
                        uint64_t w = divisor.m_words[i];
                        if ((w & (w - 1)) != 0) {
                            is_power_of_two = false;
                            break;
                        }

                        // ビット位置を計算
                        for (int b = 0; b < 64; ++b) {
                            if (w & (1ULL << b)) {
                                power = i * 64 + b;
                                break;
                            }
                        }
                    }
                }

                if (is_power_of_two && power > 0) {
                    // PowerOfTwo を使用
                    Int remainder_tmp;
                    Int quotient_tmp = IntDivision::divPowerOfTwo(result, power, remainder_tmp);
                    result = std::move(quotient_tmp);
                    return;
                }
            }

            // 2. 単一ワード除数 → SingleWordOnly (1.67x speedup)
            if (divisor.m_words.size() == 1) {
                goto single_word_division;
            }

            // 3. マルチワード除数 → mpn::divide (BZ / schoolbook)
            if (divisor.m_words.size() >= 2) {
                const uint64_t* a = result.m_words.data();
                size_t an = result.m_words.size();
                const uint64_t* b = divisor.m_words.data();
                size_t bn = divisor.m_words.size();

                // --- 2-word 除数 × 2-word 被除数 特化 (128÷128 bit) ---
                // lshift/div_basecase/assign のオーバーヘッドを回避し
                // Knuth 3÷2 除算をインラインで実行 (商は最大 1 word)
                if (bn == 2 && an == 2) {
                    unsigned shift = std::countl_zero(b[1]);
                    uint64_t nb1, nb0, u2, u1, u0;
                    if (shift > 0) {
                        nb1 = (b[1] << shift) | (b[0] >> (64 - shift));
                        nb0 = b[0] << shift;
                        u2 = a[1] >> (64 - shift);
                        u1 = (a[1] << shift) | (a[0] >> (64 - shift));
                        u0 = a[0] << shift;
                    } else {
                        nb1 = b[1]; nb0 = b[0];
                        u2 = 0; u1 = a[1]; u0 = a[0];
                    }
                    // u2 < nb1 は常に成立 (u2 < 2^shift, nb1 >= 2^63)
                    auto [q, re] = UInt128::divmod_fast(u2, u1, nb1);
                    // Knuth adjustment (最大 2 回)
                    auto prod = UInt128::multiply(q, nb0);
                    if (prod.high > re || (prod.high == re && prod.low > u0)) {
                        q--;
                        re += nb1;
                        if (re >= nb1) { // overflow なし → 再チェック
                            prod = UInt128::multiply(q, nb0);
                            if (prod.high > re || (prod.high == re && prod.low > u0)) {
                                q--;
                            }
                        }
                    }
                    if (q == 0) {
                        result.m_words.clear();
                        result.setSign(0);
                    } else {
                        result.m_words.resize(1);
                        result.m_words[0] = q;
                    }
                    return;
                }

                // --- 小サイズ fast-path: スタックバッファで直接 div_basecase ---
                // ScratchScope/Arena 確保のオーバーヘッドを回避
                constexpr size_t SMALL_DIV_THRESHOLD = 16;  // 1024 bits
                if (an <= SMALL_DIV_THRESHOLD && bn < mpn::BZ_THRESHOLD) {
                    uint64_t na[SMALL_DIV_THRESHOLD + 2];
                    uint64_t nb[SMALL_DIV_THRESHOLD + 1];
                    uint64_t qq[SMALL_DIV_THRESHOLD + 1];

                    unsigned shift = std::countl_zero(b[bn - 1]);
                    if (shift > 0) {
                        mpn::lshift(nb, b, bn, shift);
                        na[an] = mpn::lshift(na, a, an, shift);
                    } else {
                        std::memcpy(nb, b, bn * sizeof(uint64_t));
                        std::memcpy(na, a, an * sizeof(uint64_t));
                        na[an] = 0;
                    }
                    size_t nan = an + (na[an] ? 1 : 0);
                    na[nan] = 0;  // sentinel

                    mpn::div_basecase(qq, na, nan, nb, bn);

                    size_t qn = mpn::normalized_size(qq, nan - bn + 1);
                    if (qn == 0) {
                        result.m_words.clear();
                        result.setSign(0);
                    } else {
                        result.m_words.assign(qq, qq + qn);
                        result.setSign(1);
                    }
                    return;
                }

                // --- 通常パス: Arena + mpn::divide (BZ / schoolbook) ---
                ScratchScope scope;
                size_t qn_max = an - bn + 1;
                uint64_t* q_buf = getThreadArena().alloc_limbs(qn_max + 1);
                size_t scratch_size = mpn::divide_scratch_size(an, bn);
                uint64_t* scratch = getThreadArena().alloc_limbs(scratch_size);
                uint64_t* r_buf = getThreadArena().alloc_limbs(bn);

                std::memset(q_buf, 0, (qn_max + 1) * sizeof(uint64_t));
                size_t qn = mpn::divide(q_buf, r_buf, a, an, b, bn, scratch);

                if (qn == 0) {
                    result = Int::Zero();
                } else {
                    result = Int::fromRawWords(std::span<const uint64_t>(q_buf, qn), 1);
                }
                return;
            }

            // 4. デフォルト → BitByBit (single-word のフォールバック)
            break;
        }

        // 単一ワードの除数の場合の最適化 (in-place, heap alloc 不要)
single_word_division:
        if (divisor.m_words.size() == 1) {
            uint64_t d = divisor.m_words[0];
            uint64_t rem = 0;

            // 上位ワードから処理 — 読んだ位置に直接商を書く
            for (int i = static_cast<int>(result.m_words.size()) - 1; i >= 0; --i) {
                auto [q, r] = UInt128::divmod_fast(rem, result.m_words[i], d);
                result.m_words[i] = q;
                rem = r;
            }

            // normalize: 先頭が 0 になるのは最大 1 word
            if (!result.m_words.empty() && result.m_words.back() == 0) {
                result.m_words.pop_back();
            }
            if (result.m_words.empty()) {
                result.setSign(0);
            }
            return;
        }

        // ビット単位ロングディビジョン（DivisionAlgorithm::BitByBit 指定時のフォールバック）
        // 注: Auto モードでは Knuth Algorithm D を使用（Phase 0 でバグ修正済み）
bit_by_bit_division:
        Int remainder = Int::Zero();
        Int quotient = Int::Zero();

        // 桁数の差
        int bitDiff = static_cast<int>((result.m_words.size() - divisor.m_words.size()) * 64);

        // 最上位ワードの最上位ビット位置を計算
        uint64_t resultMsw = result.m_words.back();
        uint64_t divisorMsw = divisor.m_words.back();

        int resultMsBit = 63;
        if (resultMsw > 0) {
            // resultMsw の最上位ビット位置を求める
            for (int i = 63; i >= 0; --i) {
                if (resultMsw & (1ULL << i)) {
                    resultMsBit = i;
                    break;
                }
            }
        }

        int divisorMsBit = 63;
        if (divisorMsw > 0) {
            // divisorMsw の最上位ビット位置を求める
            for (int i = 63; i >= 0; --i) {
                if (divisorMsw & (1ULL << i)) {
                    divisorMsBit = i;
                    break;
                }
            }
        }

        bitDiff += resultMsBit - divisorMsBit;

        // 除数をシフトして位置合わせ
        Int shiftedDivisor = divisor;
        leftShift(shiftedDivisor, bitDiff);

        // 位置合わせした除数が被除数より大きくなった場合は調整
        if (compareAbsGreater(shiftedDivisor, result)) {
            rightShift(shiftedDivisor, 1);
            bitDiff--;
        }

        // 除算処理
        remainder = result;
        quotient.m_words.resize((bitDiff / 64) + 1, 0);

        while (bitDiff >= 0) {
            if (compareAbsGreater(remainder, shiftedDivisor) || compareAbsEqual(remainder, shiftedDivisor)) {
                subtractAbsolute(remainder, shiftedDivisor, remainder);
                // 商のビットを立てる
                quotient.m_words[bitDiff / 64] |= (1ULL << (bitDiff % 64));
            }

            rightShift(shiftedDivisor, 1);
            bitDiff--;
        }

        result = std::move(quotient);
        result.normalize();
        // quotient は Int::Zero() から構築されるため sign=0 のまま。
        // 正規化後に非ゼロなら sign を設定する（絶対値除算なので正）。
        if (!result.m_words.empty()) {
            result.setSign(1);
        }
    }

    // 剰余操作
    void IntOps::moduloAbsolute(Int& result, const Int& divisor) {
        // 特殊状態の処理
        if (result.isSpecialState() || divisor.isSpecialState()) {
            result = IntSpecialStates::handleModulo(result, divisor);
            return;
        }

        // ゼロによる剰余チェック
        if (divisor.getSign() == 0) {
            result.setState(NumericState::NaN, NumericError::DivideByZero);
            return;
        }

        // ゼロの剰余
        if (result.getSign() == 0) {
            return; // 0 % x = 0 (x ≠ 0)
        }

        // 同じ値による剰余
        if (compareAbsEqual(result, divisor)) {
            result.m_words.clear();
            result.setSign(0);
            return;
        }

        // |result| < |divisor| の場合
        if (compareAbsLess(result, divisor)) {
            return; // 絶対値がすでに小さいので変更なし
        }

        // 単一ワードの除数の場合の最適化
        if (divisor.m_words.size() == 1) {
            uint64_t d = divisor.m_words[0];
            uint64_t remainder = 0;

            // 上位ワードから処理
            for (int i = static_cast<int>(result.m_words.size()) - 1; i >= 0; --i) {
                std::pair<uint64_t, uint64_t> qr = UInt128::divmod_fast(remainder, result.m_words[i], d);
                remainder = qr.second;
            }

            if (remainder == 0) {
                result.m_words.clear();
                result.setSign(0);
            }
            else {
                result.m_words.resize(1);
                result.m_words[0] = remainder;
                // 符号は呼び出し元で処理される
            }

            return;
        }

        // 通常の剰余処理
        divmod(result, divisor, result);
    }

    // 除算と剰余の同時計算
    Int IntOps::divmod(const Int& dividend, const Int& divisor, Int& remainder) {
        // 特殊状態の処理
        if (dividend.isSpecialState() || divisor.isSpecialState()) {
            remainder = IntSpecialStates::handleModulo(dividend, divisor);
            return IntSpecialStates::handleDivision(dividend, divisor);
        }

        // ゼロによる除算チェック
        if (divisor.getSign() == 0) {
            remainder.setState(NumericState::NaN, NumericError::DivideByZero);
            Int quotient = Int::NaN();
            quotient.setState(NumericState::NaN, NumericError::DivideByZero);
            return quotient;
        }

        // ゼロの除算
        if (dividend.getSign() == 0) {
            remainder = Int::Zero();
            return Int::Zero();
        }

        // 符号の決定
        int quotientSign = (dividend.getSign() == divisor.getSign()) ? 1 : -1;

        // 絶対値を使用
        Int absDividend = dividend;
        if (absDividend.getSign() < 0) {
            absDividend.setSign(1);
        }

        Int absDivisor = divisor;
        if (absDivisor.getSign() < 0) {
            absDivisor.setSign(1);
        }

        // |dividend| < |divisor| の場合
        if (compareAbsLess(absDividend, absDivisor)) {
            remainder = dividend; // 絶対値がすでに小さいので剰余は被除数自体
            return Int::Zero();
        }

        // 単一ワードの除数の場合の最適化
        if (absDivisor.m_words.size() == 1) {
            uint64_t d = absDivisor.m_words[0];
            uint64_t r = 0;

            // 上位ワードから処理
            std::vector<uint64_t> quotient(absDividend.m_words.size(), 0);
            for (int i = static_cast<int>(absDividend.m_words.size()) - 1; i >= 0; --i) {
                UInt128 current(r, absDividend.m_words[i]);
                std::pair<uint64_t, uint64_t> qr = UInt128::divmod_fast(r, absDividend.m_words[i], d);
                quotient[i] = qr.first;
                r = qr.second;
            }

            // 剰余を設定
            if (r == 0) {
                remainder = Int::Zero();
            }
            else {
                remainder.m_words.resize(1);
                remainder.m_words[0] = r;
                remainder.setSign(dividend.getSign()); // 剰余の符号は被除数と同じ
                remainder.setState(NumericState::Normal);
            }

            // 商を設定
            Int result;
            result.m_words.assign(quotient.data(), quotient.data() + quotient.size());
            result.setSign(quotientSign);
            result.setState(NumericState::Normal);
            result.normalize();

            return result;
        }

        // 通常の除算処理: mpn::divide (BZ / schoolbook)
        if (absDivisor.m_words.size() >= 2) {
            const uint64_t* a = absDividend.m_words.data();
            size_t an = absDividend.m_words.size();
            const uint64_t* b = absDivisor.m_words.data();
            size_t bn = absDivisor.m_words.size();

            ScratchScope scope;
            size_t qn_max = an - bn + 1;
            uint64_t* q_buf = getThreadArena().alloc_limbs(qn_max + 1);
            uint64_t* r_buf = getThreadArena().alloc_limbs(bn);
            size_t scratch_size = mpn::divide_scratch_size(an, bn);
            uint64_t* scratch = getThreadArena().alloc_limbs(scratch_size);

            std::memset(q_buf, 0, (qn_max + 1) * sizeof(uint64_t));
            size_t qn = mpn::divide(q_buf, r_buf, a, an, b, bn, scratch);

            // 商を Int に変換
            Int quotient;
            if (qn == 0) {
                quotient = Int::Zero();
            } else {
                quotient = Int::fromRawWords(std::span<const uint64_t>(q_buf, qn), quotientSign);
            }

            // 余りを Int に変換
            size_t rn = mpn::normalized_size(r_buf, bn);
            if (rn == 0) {
                remainder = Int::Zero();
            } else {
                remainder = Int::fromRawWords(std::span<const uint64_t>(r_buf, rn), dividend.getSign());
            }

            return quotient;
        }

        // フォールバック: BitByBit（single-word 除数が上で処理されるため、ここには来ないはず）
        Int quotient = Int::Zero();
        remainder = absDividend;

        // 桁数の差
        int bitDiff = static_cast<int>((absDividend.m_words.size() - absDivisor.m_words.size()) * 64);

        // 最上位ワードの最上位ビット位置を計算
        uint64_t dividendMsw = absDividend.m_words.back();
        uint64_t divisorMsw = absDivisor.m_words.back();

        int dividendMsBit = 63;
        if (dividendMsw > 0) {
            // dividendMsw の最上位ビット位置を求める
            for (int i = 63; i >= 0; --i) {
                if (dividendMsw & (1ULL << i)) {
                    dividendMsBit = i;
                    break;
                }
            }
        }

        int divisorMsBit = 63;
        if (divisorMsw > 0) {
            // divisorMsw の最上位ビット位置を求める
            for (int i = 63; i >= 0; --i) {
                if (divisorMsw & (1ULL << i)) {
                    divisorMsBit = i;
                    break;
                }
            }
        }

        bitDiff += dividendMsBit - divisorMsBit;

        // 除数をシフトして位置合わせ
        Int shiftedDivisor = absDivisor;
        leftShift(shiftedDivisor, bitDiff);

        // 位置合わせした除数が被除数より大きくなった場合は調整
        if (compareAbsGreater(shiftedDivisor, absDividend)) {
            rightShift(shiftedDivisor, 1);
            bitDiff--;
        }

        // 除算処理
        quotient.m_words.resize((bitDiff / 64) + 1, 0);

        while (bitDiff >= 0) {
            if (compareAbsGreater(remainder, shiftedDivisor) || compareAbsEqual(remainder, shiftedDivisor)) {
                subtractAbsolute(remainder, shiftedDivisor, remainder);
                // 商のビットを立てる
                quotient.m_words[bitDiff / 64] |= (1ULL << (bitDiff % 64));
            }

            rightShift(shiftedDivisor, 1);
            bitDiff--;
        }

        // 符号を調整
        quotient.setSign(quotientSign);
        if (!remainder.isZero()) {
            remainder.setSign(dividend.getSign()); // 剰余の符号は被除数と同じ
        }

        quotient.normalize();
        remainder.normalize();

        return quotient;
    }

    // =========================================================================
    // ビット操作（2の補数セマンティクス）
    //
    // 符号付き整数のビット演算は、2の補数表現として解釈する。
    // 内部表現は符号-絶対値方式なので、以下の恒等式で変換する:
    //
    //   正の数 +m: ビット列 = m（上位は 0 で無限延長）
    //   負の数 -m: ビット列 = ~(m-1)（上位は 1 で無限延長）
    //
    // 各演算の符号ごとの公式:
    //   AND:  (+a)&(+b) = +(a&b)
    //         (+a)&(-b) = +(a & ~(|b|-1))
    //         (-a)&(+b) = +(~(|a|-1) & b)
    //         (-a)&(-b) = -(((|a|-1)|(|b|-1)) + 1)
    //
    //   OR:   (+a)|(+b) = +(a|b)
    //         (+a)|(-b) = -(((|b|-1) & ~a) + 1)
    //         (-a)|(+b) = -(((|a|-1) & ~b) + 1)
    //         (-a)|(-b) = -(((|a|-1)&(|b|-1)) + 1)
    //
    //   XOR:  (+a)^(+b) = +(a^b)
    //         (+a)^(-b) = -((a^(|b|-1)) + 1)
    //         (-a)^(+b) = -(((|a|-1)^b) + 1)
    //         (-a)^(-b) = +((|a|-1)^(|b|-1))
    //
    //   NOT:  ~n = -(n+1)
    // =========================================================================

    // ヘルパー: ワード配列から 1 を引く（|x|-1 の計算用）
    // 前提: 配列は非ゼロ
    template<typename Container>
    static void words_sub1(Container& w) {
        mpn::sub_1(w.data(), w.size(), 1ULL);
    }

    // ヘルパー: ワード配列に 1 を足す（結果 + 1 の計算用）
    template<typename Container>
    static void words_add1(Container& w) {
        uint64_t carry = mpn::add_1(w.data(), w.size(), 1ULL);
        if (carry) w.push_back(1);
    }

    // ヘルパー: normalize — 上位のゼロワードを除去
    template<typename Container>
    static void words_normalize(Container& w) {
        while (!w.empty() && w.back() == 0) {
            w.pop_back();
        }
    }

    void IntOps::bitwiseAnd(Int& result, const Int& other) {
        // 特殊状態の処理
        if (result.isSpecialState() || other.isSpecialState()) {
            result.setState(NumericState::NaN, NumericError::InvalidBitOperation);
            return;
        }

        int sa = result.getSign();
        int sb = other.getSign();

        // ゼロとの AND は常にゼロ方向に吸収
        if (sa == 0) { return; }
        if (sb == 0) { result = Int::Zero(); return; }

        if (sa > 0 && sb > 0) {
            // (+a) & (+b) = +(a & b)
            size_t minSize = std::min(result.m_words.size(), other.m_words.size());
            result.m_words.resize(minSize);
            for (size_t i = 0; i < minSize; ++i) {
                result.m_words[i] &= other.m_words[i];
            }
            result.normalize();
        }
        else if (sa > 0 && sb < 0) {
            // (+a) & (-b) = +(a & ~(|b|-1))
            // -b の2の補数ビット: ~(|b|-1)。上位は 1 で延長。
            SboWords bm = other.m_words;
            words_sub1(bm); // |b| - 1
            size_t aSize = result.m_words.size();
            size_t bSize = bm.size();
            // 結果は a のサイズ以下（上位の -b ビットは 1 なので a がそのまま残る）
            for (size_t i = 0; i < aSize; ++i) {
                uint64_t b_bits = (i < bSize) ? ~bm[i] : UINT64_MAX;
                result.m_words[i] &= b_bits;
            }
            // 結果は正
            result.normalize();
        }
        else if (sa < 0 && sb > 0) {
            // (-a) & (+b) = +(~(|a|-1) & b)
            SboWords am = result.m_words;
            words_sub1(am); // |a| - 1
            size_t aSize = am.size();
            size_t bSize = other.m_words.size();
            // 結果は b のサイズ以下
            result.m_words.resize(bSize);
            for (size_t i = 0; i < bSize; ++i) {
                uint64_t a_bits = (i < aSize) ? ~am[i] : UINT64_MAX;
                result.m_words[i] = a_bits & other.m_words[i];
            }
            result.setSign(1);
            result.normalize();
        }
        else {
            // (-a) & (-b) = -(((|a|-1) | (|b|-1)) + 1)
            SboWords am = result.m_words;
            SboWords bm = other.m_words;
            words_sub1(am);
            words_sub1(bm);
            size_t maxSize = std::max(am.size(), bm.size());
            am.resize(maxSize, 0);
            bm.resize(maxSize, 0);
            result.m_words.resize(maxSize);
            for (size_t i = 0; i < maxSize; ++i) {
                result.m_words[i] = am[i] | bm[i];
            }
            words_add1(result.m_words);
            result.setSign(-1);
            result.normalize();
        }
    }

    void IntOps::bitwiseOr(Int& result, const Int& other) {
        // 特殊状態の処理
        if (result.isSpecialState() || other.isSpecialState()) {
            result.setState(NumericState::NaN, NumericError::InvalidBitOperation);
            return;
        }

        int sa = result.getSign();
        int sb = other.getSign();

        // ゼロとの OR
        if (sa == 0) { result = other; return; }
        if (sb == 0) { return; }

        if (sa > 0 && sb > 0) {
            // (+a) | (+b) = +(a | b)
            size_t maxSize = std::max(result.m_words.size(), other.m_words.size());
            result.m_words.resize(maxSize, 0);
            for (size_t i = 0; i < other.m_words.size(); ++i) {
                result.m_words[i] |= other.m_words[i];
            }
            result.normalize();
        }
        else if (sa > 0 && sb < 0) {
            // (+a) | (-b) = -(((|b|-1) & ~a) + 1)
            SboWords bm = other.m_words;
            words_sub1(bm);
            size_t aSize = result.m_words.size();
            size_t bSize = bm.size();
            size_t maxSize = std::max(aSize, bSize);
            result.m_words.resize(maxSize, 0);
            bm.resize(maxSize, 0);
            for (size_t i = 0; i < maxSize; ++i) {
                uint64_t a_val = result.m_words[i];
                result.m_words[i] = bm[i] & ~a_val;
            }
            words_normalize(result.m_words);
            words_add1(result.m_words);
            result.setSign(-1);
            result.normalize();
        }
        else if (sa < 0 && sb > 0) {
            // (-a) | (+b) = -(((|a|-1) & ~b) + 1)
            SboWords am = result.m_words;
            words_sub1(am);
            size_t aSize = am.size();
            size_t bSize = other.m_words.size();
            size_t maxSize = std::max(aSize, bSize);
            am.resize(maxSize, 0);
            result.m_words.resize(maxSize);
            for (size_t i = 0; i < maxSize; ++i) {
                uint64_t b_val = (i < bSize) ? other.m_words[i] : 0;
                result.m_words[i] = am[i] & ~b_val;
            }
            words_normalize(result.m_words);
            words_add1(result.m_words);
            result.setSign(-1);
            result.normalize();
        }
        else {
            // (-a) | (-b) = -(((|a|-1) & (|b|-1)) + 1)
            SboWords am = result.m_words;
            SboWords bm = other.m_words;
            words_sub1(am);
            words_sub1(bm);
            size_t minSize = std::min(am.size(), bm.size());
            result.m_words.resize(minSize);
            for (size_t i = 0; i < minSize; ++i) {
                result.m_words[i] = am[i] & bm[i];
            }
            words_normalize(result.m_words);
            words_add1(result.m_words);
            result.setSign(-1);
            result.normalize();
        }
    }

    void IntOps::bitwiseXor(Int& result, const Int& other) {
        // 特殊状態の処理
        if (result.isSpecialState() || other.isSpecialState()) {
            result.setState(NumericState::NaN, NumericError::InvalidBitOperation);
            return;
        }

        int sa = result.getSign();
        int sb = other.getSign();

        // ゼロとの XOR
        if (sa == 0) { result = other; return; }
        if (sb == 0) { return; }

        if (sa > 0 && sb > 0) {
            // (+a) ^ (+b) = +(a ^ b)
            size_t maxSize = std::max(result.m_words.size(), other.m_words.size());
            result.m_words.resize(maxSize, 0);
            for (size_t i = 0; i < other.m_words.size(); ++i) {
                result.m_words[i] ^= other.m_words[i];
            }
            result.normalize();
        }
        else if (sa > 0 && sb < 0) {
            // (+a) ^ (-b) = -((a ^ (|b|-1)) + 1)
            SboWords bm = other.m_words;
            words_sub1(bm);
            size_t aSize = result.m_words.size();
            size_t bSize = bm.size();
            size_t maxSize = std::max(aSize, bSize);
            result.m_words.resize(maxSize, 0);
            bm.resize(maxSize, 0);
            for (size_t i = 0; i < maxSize; ++i) {
                result.m_words[i] ^= bm[i];
            }
            words_normalize(result.m_words);
            words_add1(result.m_words);
            result.setSign(-1);
            result.normalize();
        }
        else if (sa < 0 && sb > 0) {
            // (-a) ^ (+b) = -(((|a|-1) ^ b) + 1)
            SboWords am = result.m_words;
            words_sub1(am);
            size_t aSize = am.size();
            size_t bSize = other.m_words.size();
            size_t maxSize = std::max(aSize, bSize);
            am.resize(maxSize, 0);
            result.m_words.resize(maxSize);
            for (size_t i = 0; i < maxSize; ++i) {
                uint64_t b_val = (i < bSize) ? other.m_words[i] : 0;
                result.m_words[i] = am[i] ^ b_val;
            }
            words_normalize(result.m_words);
            words_add1(result.m_words);
            result.setSign(-1);
            result.normalize();
        }
        else {
            // (-a) ^ (-b) = +((|a|-1) ^ (|b|-1))
            SboWords am = result.m_words;
            SboWords bm = other.m_words;
            words_sub1(am);
            words_sub1(bm);
            size_t maxSize = std::max(am.size(), bm.size());
            am.resize(maxSize, 0);
            bm.resize(maxSize, 0);
            result.m_words.resize(maxSize);
            for (size_t i = 0; i < maxSize; ++i) {
                result.m_words[i] = am[i] ^ bm[i];
            }
            result.setSign(1);
            result.normalize();
        }
    }

    void IntOps::bitwiseNot(Int& result) {
        // 特殊状態の処理
        if (result.isSpecialState()) {
            result.setState(NumericState::NaN, NumericError::InvalidBitOperation);
            return;
        }

        // ~n = -(n + 1) (2の補数の恒等式)
        if (result.getSign() == 0) {
            // ~0 = -1
            result = Int(-1);
        }
        else if (result.getSign() > 0) {
            // ~(+m) = -(m + 1)
            words_add1(result.m_words);
            result.setSign(-1);
        }
        else {
            // ~(-m) = m - 1
            words_sub1(result.m_words);
            result.setSign(1);
            result.normalize(); // m-1 が 0 になる場合あり（m=1 のとき）
        }
    }

    void IntOps::leftShift(Int& value, int shift) {
        // 特殊状態の処理
        if (value.isSpecialState()) {
            return;
        }

        // 0のシフトや0シフトは何もしない
        if (value.getSign() == 0 || shift == 0) {
            return;
        }

        // ワード単位のシフト（64ビット単位）
        int wordShift = shift / 64;
        int bitShift = shift % 64;

        // 結果のサイズ計算
        size_t newSize = value.m_words.size() + wordShift + (bitShift > 0 ? 1 : 0);
        std::vector<uint64_t> result(newSize, 0);

        // ビット単位のシフト
        if (bitShift > 0) {
            uint64_t carry = 0;
            for (size_t i = 0; i < value.m_words.size(); ++i) {
                uint64_t word = value.m_words[i];
                result[i + wordShift] = (word << bitShift) | carry;
                carry = word >> (64 - bitShift);
            }
            result[value.m_words.size() + wordShift] = carry;
        }
        else {
            // ビットシフトなし（ワード単位のシフトのみ）
            for (size_t i = 0; i < value.m_words.size(); ++i) {
                result[i + wordShift] = value.m_words[i];
            }
        }

        value.m_words.assign(result.data(), result.data() + result.size());
        value.normalize();
    }

    void IntOps::rightShift(Int& value, int shift) {
        // 特殊状態の処理
        if (value.isSpecialState()) {
            return;
        }

        // 0のシフトや0シフトは何もしない
        if (value.getSign() == 0 || shift == 0) {
            return;
        }

        // 2の補数セマンティクス: 負の値の右シフトは floor 除算
        // -7 >> 1 = -4 (floor(-7/2) = -4), not -3 (truncation)
        // 負の値でシフトアウトされるビットが非ゼロなら、絶対値に 1 を加算してからシフト
        bool needRoundUp = false;
        if (value.getSign() < 0) {
            // シフトアウトされるビットに 1 が含まれるか検査
            for (int i = 0; i < std::min(static_cast<int>(value.m_words.size()), shift / 64); ++i) {
                if (value.m_words[i] != 0) {
                    needRoundUp = true;
                    break;
                }
            }
            if (!needRoundUp && (shift % 64) > 0) {
                int wordIdx = shift / 64;
                if (wordIdx < static_cast<int>(value.m_words.size())) {
                    uint64_t mask = (1ULL << (shift % 64)) - 1;
                    if (value.m_words[wordIdx] & mask) {
                        needRoundUp = true;
                    }
                }
            }
        }

        // ワード単位のシフト（64ビット単位）
        int wordShift = shift / 64;
        int bitShift = shift % 64;

        // 完全にシフトアウトされる場合
        if (wordShift >= static_cast<int>(value.m_words.size())) {
            if (needRoundUp) {
                // 負の値で全ビットシフトアウト → floor = -1
                value.m_words.resize(1);
                value.m_words[0] = 1;
                value.setSign(-1);
            } else {
                value.m_words.clear();
                value.setSign(0);
            }
            return;
        }

        if (bitShift == 0) {
            // ワード単位のみ: erase で in-place 削除 (alloc なし)
            value.m_words.erase(value.m_words.begin(),
                                value.m_words.begin() + wordShift);
        } else {
            // ビット+ワードシフト: 2段階で処理
            // Step 1: ワード単位の erase (memmove)
            if (wordShift > 0) {
                value.m_words.erase(value.m_words.begin(),
                                    value.m_words.begin() + wordShift);
            }
            // Step 2: in-place ビットシフト (wordShift 済みなので overlap なし)
            size_t n = value.m_words.size();
            uint64_t carry = 0;
            for (int i = static_cast<int>(n) - 1; i >= 0; --i) {
                uint64_t w = value.m_words[i];
                value.m_words[i] = (w >> bitShift) | carry;
                carry = w << (64 - bitShift);
            }
            // MSB ゼロワード除去
            while (!value.m_words.empty() && value.m_words.back() == 0) {
                value.m_words.pop_back();
            }
        }

        // 負の値の floor 補正: 絶対値に 1 を加算
        if (needRoundUp) {
            uint64_t c = mpn::add_1(value.m_words.data(), value.m_words.size(), 1ULL);
            if (c) {
                value.m_words.push_back(1);
            }
        }

        value.normalize();
    }

    // 3引数版 leftShift: result = value << shift (バッファ再利用)
    void IntOps::leftShift(const Int& value, int shift, Int& result) {
        if (value.isSpecialState() || value.getSign() == 0 || shift == 0) {
            if (&result != &value) result = value;
            return;
        }

        size_t n = value.m_words.size();
        int wordShift = shift / 64;
        int bitShift = shift % 64;
        size_t newSize = n + wordShift + (bitShift > 0 ? 1 : 0);

        result.m_words.resize_uninitialized(newSize);
        uint64_t* rp = result.m_words.data();
        const uint64_t* ap = value.m_words.data();

        // 下位ワードをゼロ埋め
        std::memset(rp, 0, wordShift * sizeof(uint64_t));

        if (bitShift == 0) {
            std::memcpy(rp + wordShift, ap, n * sizeof(uint64_t));
        } else {
            // mpn::lshift で直接シフト (ASM 対応)
            uint64_t overflow = mpn::lshift(rp + wordShift, ap, n, bitShift);
            rp[wordShift + n] = overflow;
        }

        // 正規化 (上位のゼロ除去)
        while (newSize > 0 && rp[newSize - 1] == 0) newSize--;
        result.m_words.resize(newSize);

        result.m_sign = value.m_sign;
        result.m_state = NumericState::Normal;
    }

    // 3引数版 rightShift: result = value >> shift (バッファ再利用)
    void IntOps::rightShift(const Int& value, int shift, Int& result) {
        if (value.isSpecialState() || value.getSign() == 0 || shift == 0) {
            if (&result != &value) result = value;
            return;
        }

        int wordShift = shift / 64;
        int bitShift = shift % 64;
        size_t n = value.m_words.size();

        // 完全にシフトアウト
        if (wordShift >= static_cast<int>(n)) {
            if (value.getSign() < 0) {
                // 負の値の floor 除算: -1
                result.m_words.resize(1);
                result.m_words[0] = 1;
                result.m_sign = -1;
            } else {
                result.m_words.clear();
                result.m_sign = 0;
            }
            result.m_state = NumericState::Normal;
            return;
        }

        // 負の値の floor 補正チェック
        bool needRoundUp = false;
        if (value.getSign() < 0) {
            for (int i = 0; i < std::min(static_cast<int>(n), wordShift); ++i) {
                if (value.m_words[i] != 0) { needRoundUp = true; break; }
            }
            if (!needRoundUp && bitShift > 0 && wordShift < static_cast<int>(n)) {
                uint64_t mask = (1ULL << bitShift) - 1;
                if (value.m_words[wordShift] & mask) needRoundUp = true;
            }
        }

        size_t srcN = n - wordShift;
        const uint64_t* ap = value.m_words.data() + wordShift;

        result.m_words.resize_uninitialized(srcN);
        uint64_t* rp = result.m_words.data();

        if (bitShift == 0) {
            std::memcpy(rp, ap, srcN * sizeof(uint64_t));
        } else {
            mpn::rshift(rp, ap, srcN, bitShift);
        }

        // 正規化
        while (srcN > 0 && rp[srcN - 1] == 0) srcN--;
        result.m_words.resize(srcN);

        if (srcN == 0) {
            result.m_sign = needRoundUp ? -1 : 0;
            if (needRoundUp) {
                result.m_words.resize(1);
                result.m_words[0] = 1;
            }
        } else {
            result.m_sign = value.m_sign;
            if (needRoundUp) {
                uint64_t c = mpn::add_1(result.m_words.data(), result.m_words.size(), 1ULL);
                if (c) result.m_words.push_back(1);
            }
        }
        result.m_state = NumericState::Normal;
    }

    // 2の累乗係数
    Int IntOps::pow2(uint64_t exponent) {
        // ワード単位のシフト（64ビット単位）
        int wordShift = (int)(exponent / 64);
        int bitShift = exponent % 64;

        // 結果の初期化
        Int result;
        result.m_words.resize(wordShift + 1, 0);
        result.m_words[wordShift] = 1ULL << bitShift;
        result.setSign(1);

        return result;
    }

    // GCD (最大公約数)
    Int IntOps::gcd(const Int& a, const Int& b) {
        // Lehmer GCD (IntGCD) に委譲
        return IntGCD::gcd(a, b);
    }

    // 因子除去の実装
    uint64_t IntOps::removeFactor(const Int& value, const Int& factor, Int& result) {
        // 特殊状態のチェック
        if (value.isNaN() || factor.isNaN()) {
            result = Int::NaN();
            return 0;
        }

        if (value.isInfinite() || factor.isInfinite()) {
            result = Int::NaN();
            return 0;
        }

        // factor が 0, ±1 の場合は不正
        if (factor.isZero() || abs(factor).isOne()) {
            result = Int::NaN();
            return 0;
        }

        // value が 0 の場合
        if (value.isZero()) {
            result = Int::Zero();
            return 0;
        }

        // GMP 互換: 因子は常に絶対値で扱う
        Int abs_factor = abs(factor);

        // 因子を除去していく
        result = value;
        uint64_t count = 0;

        // result が abs_factor で割り切れる限り除去を続ける
        while (true) {
            Int quotient = result / abs_factor;
            Int remainder = result % abs_factor;

            if (!remainder.isZero()) {
                // 割り切れなくなったら終了
                break;
            }

            result = quotient;
            count++;

            // カウントオーバーフロー防止
            if (count == UINT64_MAX) {
                break;
            }
        }

        return count;
    }

    // 平方根 (IntSqrt に委譲)
    Int IntOps::sqrt(const Int& value) {
        return IntSqrt::sqrt(value);
    }

    // カラツバ法による乗算
    void IntOps::karatsubaMultiply(const Int& a, const Int& b, Int& result) {
        // IntMultiplication クラスの Karatsuba 実装を使用
        result = IntMultiplication::karatsubaMultiply(a, b);
    }

    // Toom-Cook-3 による乗算
    void IntOps::toomCookMultiply(const Int& a, const Int& b, Int& result) {
        // IntMultiplication クラスの Toom-Cook-3 実装を使用
        result = IntMultiplication::toomCook3Multiply(a, b);
    }

    // Toom-Cook-4 による乗算
    void IntOps::toomCook4Multiply(const Int& a, const Int& b, Int& result) {
        // IntMultiplication クラスの Toom-Cook-4 実装を使用
        result = IntMultiplication::toomCook4Multiply(a, b);
    }

    // 汎用乗算 (アンバランス対応)
    // result.m_words に直接書き込み、fromRawWords の中間コピーを除去
    void IntOps::generalMultiply(const Int& a, const Int& b, Int& result) {
        size_t an = a.m_words.size();
        size_t bn = b.m_words.size();

        ScratchScope scope;
        size_t scratch_size = mpn::multiply_scratch_size(an, bn);
        uint64_t* scratch = (scratch_size > 0)
            ? getThreadArena().alloc_limbs(scratch_size) : nullptr;

        size_t rn = an + bn;
        result.m_words.resize_uninitialized(rn);
        mpn::multiply(result.m_words.data(), a.m_words.data(), an,
                      b.m_words.data(), bn, scratch);
        // normalize: 上位のゼロワードを除去
        while (!result.m_words.empty() && result.m_words.back() == 0)
            result.m_words.pop_back();
        result.setSign(1);
        result.setState(NumericState::Normal);
    }

    // NTT キャッシュ付き絶対値乗算: rhs の forward NTT をキャッシュして再利用
    void IntOps::mulAbsCached(const Int& lhs, const Int& rhs, Int& result,
                               prime_ntt::NttCache& cache) {
        size_t an = lhs.m_words.size();
        size_t bn = rhs.m_words.size();

        // NTT 閾値未満はキャッシュ不要 → 通常乗算にフォールバック
        if (an == 0 || bn == 0) {
            result.m_words.clear();
            result.setSign(0);
            result.setState(NumericState::Normal);
            return;
        }
        if (std::min(an, bn) < mpn::PRIME_NTT_THRESHOLD) {
            multiplyAbsolute(lhs, rhs, result);
            result.setSign(1);
            result.setState(NumericState::Normal);
            return;
        }

        size_t rn = an + bn;
        result.m_words.resize_uninitialized(rn);
        prime_ntt::mul_prime_ntt_cached(result.m_words.data(),
                                         lhs.m_words.data(), an,
                                         rhs.m_words.data(), bn, cache);
        while (!result.m_words.empty() && result.m_words.back() == 0)
            result.m_words.pop_back();
        result.setSign(1);
        result.setState(NumericState::Normal);
    }

    // YC-2: Fused multiply-add: result = a*b + c*d (符号処理込み)
    void IntOps::mulAdd(const Int& a, const Int& b,
                         const Int& c, const Int& d, Int& result) {
        int sa = a.getSign(), sb = b.getSign();
        int sc = c.getSign(), sd = d.getSign();
        int sign_ab = sa * sb;  // a*b の符号
        int sign_cd = sc * sd;  // c*d の符号

        // いずれかの積がゼロ → 単純な乗算
        if (sign_ab == 0 && sign_cd == 0) {
            result = Int(0); return;
        }
        if (sign_ab == 0) { mul(c, d, result); return; }
        if (sign_cd == 0) { mul(a, b, result); return; }

        size_t an = a.m_words.size();
        size_t bn = b.m_words.size();
        size_t cn = c.m_words.size();
        size_t dn = d.m_words.size();

        // 符号が異なる場合: NTT 融合不可 → 通常の乗算 + 加算にフォールバック
        // (NTT 内で mod p の引き算は可能だが、CRT 復元時に符号付き大数を扱えない)
        if (sign_ab != sign_cd) {
            result = a * b;
            result += c * d;
            return;
        }

        // 以降、sign_ab == sign_cd (両方正 or 両方負)
        // |a*b| + |c*d| を NTT 融合で計算し、共通符号を適用

        // 両方の小さい側が NTT 閾値以上のときのみ融合
        size_t min_ab = std::min(an, bn);
        size_t min_cd = std::min(cn, dn);
        if (min_ab < mpn::PRIME_NTT_THRESHOLD || min_cd < mpn::PRIME_NTT_THRESHOLD) {
            // フォールバック
            result = a * b;
            result += c * d;
            return;
        }

        // NTT 融合パス: |a|*|b| + |c|*|d|
        size_t rn = std::max(an + bn, cn + dn) + 1;  // +1 for carry from addition
        result.m_words.resize_uninitialized(rn);
        prime_ntt::mul_add_prime_ntt(result.m_words.data(), rn,
                                      a.m_words.data(), an,
                                      b.m_words.data(), bn,
                                      c.m_words.data(), cn,
                                      d.m_words.data(), dn);
        while (!result.m_words.empty() && result.m_words.back() == 0)
            result.m_words.pop_back();
        result.setSign(result.m_words.empty() ? 0 : sign_ab);
        result.setState(NumericState::Normal);
    }

    // FFTを使用した乗算
    void IntOps::fftMultiply(const Int& a, const Int& b, Int& result) {
        // IntMultiplication クラスの FFT 実装を使用
        result = IntMultiplication::fftMultiply(a, b);
    }

    // ニュートン法による除算
    void IntOps::newtonDivision(const Int& dividend, const Int& divisor, Int& result) {
        // ニュートン法による高速除算を使用（超大型数向け）
        Int remainder;
        result = IntDivision::divNewton(dividend, divisor, remainder);
    }

    // 床除算の実装
    Int IntOps::floorDiv(const Int& dividend, const Int& divisor) {
        // 特殊状態のチェック
        if (dividend.isNaN() || divisor.isNaN()) {
            return Int::NaN();
        }
        if (divisor.isZero()) {
            return Int::NaN();  // ゼロ除算
        }
        if (dividend.isInfinite() || divisor.isInfinite()) {
            // ∞ / 有限数 = ∞, 有限数 / ∞ = 0, ∞ / ∞ = NaN
            if (dividend.isInfinite() && !divisor.isInfinite()) {
                return (dividend.getSign() * divisor.getSign() > 0) ?
                    Int::PositiveInfinity() : Int::NegativeInfinity();
            } else if (!dividend.isInfinite() && divisor.isInfinite()) {
                return Int::Zero();
            } else {
                return Int::NaN();
            }
        }

        // 通常の除算
        Int quotient = dividend / divisor;
        Int remainder = dividend % divisor;

        // 余りが 0 なら通常の除算と同じ
        if (remainder.isZero()) {
            return quotient;
        }

        // 符号が異なる場合、商を 1 減らす（負の無限大方向に丸める）
        if ((dividend.getSign() > 0 && divisor.getSign() < 0) ||
            (dividend.getSign() < 0 && divisor.getSign() > 0)) {
            quotient = quotient - Int::One();
        }

        return quotient;
    }

    // 床除算と剰余を同時に返す
    Int IntOps::floorDivMod(const Int& dividend, const Int& divisor, Int& remainder) {
        if (dividend.isNaN() || divisor.isNaN()) {
            remainder = Int::NaN();
            return Int::NaN();
        }
        if (divisor.isZero()) {
            remainder = Int::NaN();
            return Int::NaN();
        }

        Int quotient = divmod(dividend, divisor, remainder);

        // 余りが非ゼロで符号が異なる場合、調整
        if (!remainder.isZero() &&
            ((dividend.getSign() > 0 && divisor.getSign() < 0) ||
             (dividend.getSign() < 0 && divisor.getSign() > 0))) {
            quotient = quotient - Int::One();
            remainder = remainder + divisor;
        }

        return quotient;
    }

    // 床剰余の実装
    Int IntOps::floorMod(const Int& dividend, const Int& divisor) {
        Int remainder;
        floorDivMod(dividend, divisor, remainder);
        return remainder;
    }

    // 天井除算の実装
    Int IntOps::ceilDiv(const Int& dividend, const Int& divisor) {
        // 特殊状態のチェック
        if (dividend.isNaN() || divisor.isNaN()) {
            return Int::NaN();
        }
        if (divisor.isZero()) {
            return Int::NaN();  // ゼロ除算
        }
        if (dividend.isInfinite() || divisor.isInfinite()) {
            // ∞ / 有限数 = ∞, 有限数 / ∞ = 0, ∞ / ∞ = NaN
            if (dividend.isInfinite() && !divisor.isInfinite()) {
                return (dividend.getSign() * divisor.getSign() > 0) ?
                    Int::PositiveInfinity() : Int::NegativeInfinity();
            } else if (!dividend.isInfinite() && divisor.isInfinite()) {
                return Int::Zero();
            } else {
                return Int::NaN();
            }
        }

        // 通常の除算
        Int quotient = dividend / divisor;
        Int remainder = dividend % divisor;

        // 余りが 0 なら通常の除算と同じ
        if (remainder.isZero()) {
            return quotient;
        }

        // 符号が同じ場合、商を 1 増やす（正の無限大方向に丸める）
        if ((dividend.getSign() > 0 && divisor.getSign() > 0) ||
            (dividend.getSign() < 0 && divisor.getSign() < 0)) {
            quotient = quotient + Int::One();
        }

        return quotient;
    }

    // 正確な除算の実装
    Int IntOps::divExact(const Int& dividend, const Int& divisor) {
        // 特殊状態のチェック
        if (dividend.isNaN() || divisor.isNaN()) {
            return Int::NaN();
        }
        if (divisor.isZero()) {
            return Int::NaN();  // ゼロ除算
        }
        if (dividend.isInfinite() || divisor.isInfinite()) {
            if (dividend.isInfinite() && !divisor.isInfinite()) {
                return (dividend.getSign() * divisor.getSign() > 0) ?
                    Int::PositiveInfinity() : Int::NegativeInfinity();
            } else if (!dividend.isInfinite() && divisor.isInfinite()) {
                return Int::Zero();
            } else {
                return Int::NaN();
            }
        }

        // 割り切れることが保証されているので、単純に除算
        // （高速化のため余りチェックは省略）
        return dividend / divisor;
    }

} // namespace calx