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

// IntDivision.hpp
// Division algorithms for multi-precision integers
//
// This file provides division algorithms for the multi-precision integer
// library. The implementation is based on Knuth's division algorithm
// (Algorithm D).
//
// Main features:
// - Knuth's division algorithm
// - Efficient quotient and remainder computation
// - Special-case optimizations
// - Integration with state management

#ifndef SANGI_INT_DIVISION_HPP
#define SANGI_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 sangi {

    /**
     * @brief Selection of division algorithm
     *
     * An enum used to switch between different division algorithms in tests
     * and benchmarks. The algorithm chosen via the global variable
     * g_division_algorithm is used at runtime.
     */
    enum class DivisionAlgorithm {
        Auto,              // Automatic selection (chooses the optimal algorithm based on size)
        BitByBit,          // Bit-by-bit long division (O(n^2 * w), reliable but slow)
        SingleWordOnly,    // Only optimize the single-word divisor case (otherwise BitByBit)
        Knuth,             // Knuth Algorithm D (O(n^2), fast and accurate)
        Newton,            // Newton's method based division (reserved for future extension)
        BurnikelZiegler,   // Burnikel-Ziegler method (reserved for future extension)
        PowerOfTwo,        // Power-of-two specialized (right shift)
        PowerOfTen         // Power-of-ten specialized (internally uses Knuth)
    };

    /**
     * @brief Global division algorithm setting
     *
     * The default is Auto (automatic selection).
     * By changing this variable in tests or benchmarks, the division
     * algorithm in use can be switched.
     *
     * Example:
     *   sangi::g_division_algorithm = DivisionAlgorithm::Knuth;
     *   Int result = a / b;  // Uses Knuth Algorithm D
     *
     *   sangi::g_division_algorithm = DivisionAlgorithm::BitByBit;
     *   Int result2 = a / b; // Uses bit-by-bit division
     */
    inline DivisionAlgorithm g_division_algorithm = DivisionAlgorithm::Auto;

    /**
     * @brief Class providing division algorithms for multi-precision integers
     */
    class IntDivision {
    public:
        // ================================================================
        // Exact division (fast division assuming divisibility)
        // Used in Toom-Cook interpolation, etc. Complexity O(n).
        // ================================================================

        /**
         * @brief Exact division by 3 (assumes divisibility by 3)
         *
         * Uses a Montgomery inverse: 3^(-1) mod 2^64 = 0xAAAAAAAAAAAAAAAB
         * Complexity: O(n) -- an alternative to general O(n^2) division
         *
         * @param x Dividend (must be divisible by 3)
         * @return x / 3
         */
        static Int divexactBy3(const Int& x) {
            if (x.isZero()) return Int(0);
            if (x.isSpecialState()) return x;

            // Save the sign and operate on the absolute value
            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 (with borrow detection)
                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
                // Note: since 3*yi == xi_minus_carry (mod 2^64), prod.low == xi_minus_carry
                UInt128 prod = UInt128::multiply(yi, 3);
                carry = prod.high + borrow;
            }

            // Construct the result
            Int result;
            result.m_words.assign(out.data(), out.data() + out.size());
            result.setSign(result_sign);
            result.normalize();
            return result;
        }

        /**
         * @brief Exact division by 24 (assumes divisibility by 24)
         *
         * Decomposition: x/24 = (x >> 3) / 3 -- shift + divexactBy3
         * Complexity: O(n)
         *
         * @param x Dividend (must be divisible by 24)
         * @return x / 24
         */
        static Int divexactBy24(const Int& x) {
            if (x.isZero()) return Int(0);
            if (x.isSpecialState()) return x;

            // Save the sign and operate on the absolute value
            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;

            // Convert the shifted result into an 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's division algorithm
         * Reference: D. E. Knuth, "The Art of Computer Programming", Volume 2, 4.3.1
         *
         * @param dividend Dividend
         * @param divisor Divisor
         * @param remainder Output location for the remainder
         * @return Quotient
         */
        static Int divKnuth(const Int& dividend, const Int& divisor, Int& remainder) {
            // Handle special states
            if (divisor.isZero()) {
                // Division by zero handling
                if (dividend.isZero()) {
                    // 0/0 is NaN
                    remainder = Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                    return Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                }
                else {
                    // x/0 (x != 0) is also NaN (no integer quotient is defined)
                    remainder = Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                    return Int::fromState(NumericState::NaN, NumericError::DivideByZero);
                }
            }

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

            // Handle NaN, infinities, and similar special states
            if (dividend.getState() != NumericState::Normal || divisor.getState() != NumericState::Normal) {
                return IntSpecialStates::handleDivision(dividend, divisor);
            }

            // Special handling for divisors with absolute value 1
            if (abs(divisor).isOne()) {
                remainder = Int(0);
                Int q = abs(dividend);
                if (dividend.getSign() != divisor.getSign() && !q.isZero()) {
                    q.negate();
                }
                return q;
            }

            // Sign handling
            bool quotient_negative = dividend.getSign() != divisor.getSign();
            bool remainder_negative = dividend.getSign() < 0;
            Int abs_dividend = abs(dividend);
            Int abs_divisor = abs(divisor);

            // If the dividend is smaller than the divisor
            if (abs_dividend < abs_divisor) {
                remainder = dividend; // The remainder is the dividend itself
                return Int(0);
            }

            // For small values, compute directly
            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;

                // Build the result
                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;
            }

            // From here on, Knuth's division algorithm
            // Prepare input data (little-endian layout)
            std::vector<uint64_t> U = abs_dividend.words(); // Dividend
            std::vector<uint64_t> V = abs_divisor.words();  // Divisor

            // Normalization - adjust top digit so divisor[n-1] >= 2^32
            int normalization_shift = 0;
            if (!V.empty() && V.back() != 0) {
                // Normalization shift amount: shift until the MSB is set
                // Knuth Algorithm D requires the MSB of the top divisor word to be set
                normalization_shift = static_cast<int>(std::countl_zero(V.back()));
            }

            if (normalization_shift > 0) {
                // Normalize both dividend and divisor (left shift)
                Int normalized_dividend = abs_dividend << normalization_shift;
                Int normalized_divisor = abs_divisor << normalization_shift;

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

            // Compute digit counts
            size_t m = U.size();
            size_t n = V.size();

            // Allocate storage for the quotient
            // Maximum quotient length is dividend length - divisor length + 1
            std::vector<uint64_t> Q;
            if (m >= n) {
                Q.resize(m - n + 1, 0);
            }
            else {
                Q.resize(1, 0);
            }

            // Working array (copy of U)
            std::vector<uint64_t> U_work = U;
            // Append 0 to the top to prevent overflow
            U_work.push_back(0);

            // Implementation of Knuth's Algorithm D (multi-precision integer division)
            for (int j = static_cast<int>(m - n); j >= 0; j--) {
                // Bounds checks
                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");
                }

                // Compute the trial quotient estimate
                uint64_t dividend_high = U_work[j + n];
                uint64_t dividend_mid = U_work[j + n - 1];
                uint64_t divisor_high = V[n - 1];

                // Trial quotient computation
                uint64_t q_hat;
                uint64_t r_hat;
                bool skip_correction = false;

                // Improved single-word division
                if (dividend_high == 0) {
                    // Easy case - ordinary single-word division
                    q_hat = dividend_mid / divisor_high;
                    r_hat = dividend_mid % divisor_high;
                }
                else if (dividend_high >= divisor_high) {
                    // Trial quotient saturates at the maximum value (B-1)
                    q_hat = UINT64_MAX;

                    // Correct r_hat computation:
                    // r_hat = (dividend_high * B + dividend_mid) - (B-1) * divisor_high
                    //       = (dividend_high - divisor_high) * B + dividend_mid + divisor_high
                    //
                    // When dividend_high > divisor_high: r_hat >= B -> does not fit in 64 bits
                    // -> trial quotient correction is unnecessary
                    //    (r_hat >= B implies q_hat*v2 < r_hat*B always holds)
                    if (dividend_high > divisor_high) {
                        // r_hat >= B -> skip the correction loop
                        r_hat = 0;
                        skip_correction = true;
                    } else {
                        // dividend_high == divisor_high
                        // r_hat = dividend_mid + divisor_high
                        r_hat = dividend_mid + divisor_high;

                        // If r_hat overflowed (r_hat >= B), correction is also unnecessary
                        if (r_hat < dividend_mid) {
                            skip_correction = true;
                        }
                    }
                }
                else {
                    // Two-word / one-word division
                    // divmod_fast correctly handles the overflow flag
                    // (UInt128::divide has a 64-bit overflow bug in its intermediate variable r)
                    auto div_result = UInt128::divmod_fast(dividend_high, dividend_mid, divisor_high);
                    q_hat = div_result.first;
                    r_hat = div_result.second;
                }

                // Trial quotient correction (when a second word is present)
                // skip_correction: no correction needed when r_hat >= B (does not fit in 64 bits)
                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;

                    // Correction loop
                    bool need_adjustment = true;
                    int adjustment_count = 0;
                    while (need_adjustment) {
                        // Check 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;

                            // If r_hat overflowed, stop the loop
                            if (r_hat < old_r_hat) {
                                need_adjustment = false;
                            }

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

                // Subtract q_hat * divisor from the dividend (single-pass, Knuth Algorithm D / GMP style)
                // Fixes an overflow bug in the old "product.low + borrow" code
                uint64_t mul_carry = 0;

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

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

                    // product.low + mul_carry (carry from the previous iteration)
                    uint64_t sum_low = product.low + mul_carry;
                    uint64_t carry_from_add = (sum_low < product.low) ? 1 : 0;

                    // Subtract sum_low from U_work[j+i]
                    uint64_t prev = U_work[j + i];
                    U_work[j + i] = prev - sum_low;
                    uint64_t borrow_from_sub = (prev < sum_low) ? 1 : 0;

                    // Next carry = product.high + carry_from_add + borrow_from_sub
                    // Note: product.high <= 2^64-2 (maximum of a 64-bit x 64-bit multiplication)
                    // carry_from_add + borrow_from_sub <= 2, but
                    // when both are 1 we have product.high < 2^64-2, so the
                    // total is always at most 2^64-1 (no overflow).
                    mul_carry = product.high + carry_from_add + borrow_from_sub;
                }

                // Subtract the remaining carry from the top word
                {
                    uint64_t prev = U_work[j + n];
                    U_work[j + n] = prev - mul_carry;
                    // If prev < mul_carry, the result went negative (add-back needed)
                    mul_carry = (prev < mul_carry) ? 1 : 0;
                }

                // Adjustment when the subtraction result is negative (q_hat was too large)
                if (mul_carry != 0) {
                    // Decrement q_hat by 1 and add the divisor back
                    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];

                        // Two-stage addition to track the carry accurately
                        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;
                    }

                    // Add the carry into the top digit
                    U_work[j + n] += carry;
                }

                // Set the quotient digit
                Q[j] = q_hat;
            }

            // Set the remainder (the first n digits of U_work)
            std::vector<uint64_t> R;
            for (size_t i = 0; i < n && i < U_work.size(); i++) {
                R.push_back(U_work[i]);
            }

            // Restore the remainder from before the normalization shift
            Int normalizedRemainder = Int::fromRawWords(R, +1);
            normalizedRemainder.normalize();

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

            // Build the quotient
            Int quotient = Int::fromRawWords(Q, +1);
            quotient.normalize();

            // Apply signs
            if (remainder_negative && !remainder.isZero()) {
                remainder.negate();
            }

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

            return quotient;
        }

        /**
         * @brief Division implementation (quotient only)
         * @param dividend Dividend
         * @param divisor Divisor
         * @return Quotient
         */
        static Int divide(const Int& dividend, const Int& divisor) {
            Int remainder;
            return divKnuth(dividend, divisor, remainder);
        }

        /**
         * @brief Modulo implementation
         * @param dividend Dividend
         * @param divisor Divisor
         * @return Remainder
         */
        static Int modulo(const Int& dividend, const Int& divisor) {
            Int remainder;
            divKnuth(dividend, divisor, remainder);
            return remainder;
        }

        /**
         * @brief Fast division by a power of two
         *
         * Division by a power of two can be performed quickly as a right shift.
         *
         * @param dividend Dividend
         * @param power Power-of-two exponent (2^power)
         * @param remainder Output location for the remainder
         * @return Quotient
         */
        static Int divPowerOfTwo(const Int& dividend, size_t power, Int& remainder) {
            if (power == 0) {
                // Division by 2^0 = 1
                remainder = Int(0);
                return dividend;
            }

            // Save the sign
            int sign = dividend.getSign();

            // Handle special states
            if (dividend.getState() != NumericState::Normal) {
                if (dividend.isNaN() || NumericStateTraits::isInfinite(dividend.getState())) {
                    remainder = Int::fromState(NumericState::NaN, NumericError::NaNPropagation);
                    return dividend; // NaN or infinity passes through
                }
            }

            // Take the absolute value
            Int abs_dividend = abs(dividend);

            // Compute the remainder (extract the low `power` bits)
            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()) {
                // Copy whole words
                for (size_t i = 0; i < std::min(words.size(), word_shift); ++i) {
                    remainder_words.push_back(words[i]);
                }
                // Mask the leftover bits
                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);

            // Compute the quotient (right shift)
            Int quotient = abs_dividend >> static_cast<int>(power);

            // Adjust the sign
            if (sign < 0 && !quotient.isZero()) {
                    quotient.negate();
            }

            return quotient;
        }

        /**
         * @brief Fast division by a power of ten
         *
         * Used when shifting digits in the decimal representation
         *
         * @param dividend Dividend
         * @param power Power-of-ten exponent (10^power)
         * @param remainder Output location for the remainder
         * @return Quotient
         */
        static Int divPowerOfTen(const Int& dividend, size_t power, Int& remainder) {
            // Compute 10^power
            Int divisor(10);

            // Use the unsigned-int overload of pow
            Int ten_power = sangi::pow(divisor, static_cast<unsigned int>(power));

            // Use ordinary division
            return divKnuth(dividend, ten_power, remainder);
        }

        /**
         * @brief Newton-Raphson approximation of the reciprocal
         *
         * Uses double precision to obtain an initial approximation, then
         * refines it through Newton-Raphson iteration. Even for huge
         * integers, the exponent is adjusted so the value fits in the
         * double range before computation.
         *
         * Algorithm:
         * 1. Normalize divisor into the double range (adjust by bit shift)
         * 2. Compute an initial approximation of 1/divisor in double precision
         * 3. Newton-Raphson iteration: X_{n+1} = X_n * (2 - D * X_n / 2^k)
         * 4. Return the result in fixed-point: floor(2^scale / divisor)
         *
         * @param divisor Divisor (positive integer)
         * @param scale Scale factor (result is 2^scale / divisor)
         * @return floor(2^scale / divisor) - fixed-point reciprocal
         */
        static Int newtonRaphsonReciprocal(const Int& divisor, size_t scale) {
            // Handle special cases
            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();

            // Small divisors are computed directly
            if (divisor_bits <= 64) {
                Int numerator = Int(1) << static_cast<int>(scale);
                return numerator / divisor;
            }

            // Step 1: normalize the divisor into the double range
            // double's mantissa is 53 bits
            const size_t DOUBLE_MANTISSA_BITS = 53;

            // Extract the top ~53 bits of the divisor
            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;

            // Convert to 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);
                // Approximate high * 2^64 + low as a 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;
            }

            // Step 2: compute the initial reciprocal in double precision
            double recip_double = 1.0 / d_approx;

            // Convert recip_double into a fixed-point Int
            // recip_double is the reciprocal of divisor_normalized
            // Since divisor_normalized = divisor >> shift_down,
            // 1/divisor_normalized = (1/divisor) * 2^shift_down
            // i.e. recip_double ~= (1/divisor) * 2^shift_down

            // working_scale is the same as scale
            // The reciprocal X needs ~2^(scale-divisor_bits) bits of precision
            // for 2^working_scale / D. If working_scale is too small, the
            // integer representation of X has too few bits and the quotient
            // loses precision (the correction loop spins forever).
            const size_t working_scale = scale;

            // Initial approximation: 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);
            }

            // Step 3: improve precision quadratically via Newton-Raphson iteration
            // Iteration: X_{n+1} = X_n * (2*2^k - D * X_n) / 2^k
            // Each iteration doubles the number of accurate bits (quadratic convergence)
            //
            // Initial precision: 53-bit double mantissa
            // Required precision: working_scale bits (until diff < D)

            // Compute the number of iterations: 53 -> 106 -> 212 -> ... -> working_scale
            size_t initial_precision = 50;  // Conservative estimate from double precision
            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) {
                // Compute D * X_n (scale 2^k)
                Int DX = divisor * X;

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

                // Correction: 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);
            }

            // Since working_scale == scale, no scale adjustment is needed
            return X;
        }

        /**
         * @brief Newton-based division (fast algorithm for very large numbers)
         *
         * Computes 1/divisor via Newton-Raphson iteration and multiplies it
         * by the dividend to perform the division.
         * Complexity is O(M(n)) where M(n) is multiplication time
         * (O(n^1.585) for Karatsuba, O(n log n) for FFT).
         *
         * Reference: "Modern Computer Arithmetic" Algorithm 1.6
         *
         * @param dividend Dividend
         * @param divisor Divisor
         * @param remainder Output location for the remainder
         * @return Quotient
         */
        static Int divNewton(const Int& dividend, const Int& divisor, Int& remainder) {
            // Handle special cases
            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);
            }

            // Save signs
            int dividend_sign = dividend.getSign();
            int divisor_sign = divisor.getSign();
            int quotient_sign = dividend_sign * divisor_sign;

            // Operate on absolute values
            Int abs_dividend = abs(dividend);
            Int abs_divisor = abs(divisor);

            // If the divisor is larger than the dividend
            if (abs_dividend < abs_divisor) {
                remainder = dividend;  // Preserve the sign
                return Int(0);
            }

            // Newton division algorithm
            // Compute Q = floor(A / D)
            // 1. Compute R = 1/D (fixed-point: 2^k / D)
            // 2. Q_approx = A * R / 2^k
            // 3. Compute the remainder and correct the quotient

            size_t dividend_bits = abs_dividend.bitLength();

            // Choose the scale: dividend_bits + headroom
            // For the reciprocal R = 2^scale / D to carry quotient_bits of
            // precision, we need scale >= dividend_bits + margin.
            // divisor_bits + 64 is insufficient precision, causing the
            // correction loop to spin forever.
            size_t scale = dividend_bits + 64;

            // Step 1: compute 1/divisor
            Int reciprocal = newtonRaphsonReciprocal(abs_divisor, scale);

            // Step 2: compute the quotient approximation
            // Q ~= A * (2^scale / D) / 2^scale = A * reciprocal / 2^scale
            Int quotient = (abs_dividend * reciprocal) >> static_cast<int>(scale);

            // Step 3: compute the remainder
            Int rem = abs_dividend - (quotient * abs_divisor);

            // If the remainder is negative, decrement the quotient
            // (usually only 1-2 adjustments are required)
            while (rem.getSign() < 0) {
                quotient -= 1;
                rem = rem + abs_divisor;
            }

            // If the remainder is >= divisor, increment the quotient
            // (usually only 1-2 adjustments are required)
            while (rem >= abs_divisor) {
                quotient += 1;
                rem = rem - abs_divisor;
            }

            // Apply signs
            if (quotient_sign < 0 && !quotient.isZero()) {
                quotient.negate();
            }
            if (dividend_sign < 0 && !rem.isZero()) {
                rem.negate();
            }

            remainder = rem;
            return quotient;
        }

    private:
        // Internal implementation helper functions
    };

} // namespace sangi

#endif // SANGI_INT_DIVISION_HPP