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

// IntSqrt.hpp
// Integer square root, n-th root, and perfect-power detection
//
// Algorithms:
// - sqrt: mpn-level Newton method (double precision for small values)
// - sqrtRem: compute square root and remainder simultaneously
// - nthRoot: n-th root computation via precision-doubling Newton method
// - nthRootRem: compute n-th root and remainder simultaneously
// - isSquare: GMP-style 3-layer filter (mod 256 + mod_34lsub1 + sqrtRem)
//
// References:
// - Brent & Zimmermann, "Modern Computer Arithmetic", section 1.5.2
// - GMP: mpz_sqrt, mpz_sqrtrem, mpz_root, mpz_rootrem, mpz_perfect_square_p

#ifndef SANGI_INT_SQRT_HPP
#define SANGI_INT_SQRT_HPP

#include "IntBase.hpp"
#include <utility>

namespace sangi {

class IntSqrt {
public:
    // Integer square root (floor(sqrt(value)))
    // Special states: NaN -> NaN, negative -> NaN (NegativeSqrt), +inf -> +inf, 0 -> 0
    // Algorithm: chosen automatically by size
    //   - MSB < 52: computed in double precision (exact)
    //   - otherwise: mpn-level Newton method
    static Int sqrt(const Int& value);

    // Return square root and remainder: value = sqrt^2 + remainder
    // remainder always satisfies 0 <= remainder < 2*sqrt + 1
    // GMP compatible: mpz_sqrtrem
    static Int sqrtRem(const Int& value, Int& remainder);

    // Perfect-square check.
    // If pSqrt != nullptr, the square root is stored in *pSqrt.
    // Acceleration: GMP-style 3-layer filter (mod 256, mod_34lsub1, sqrtRem)
    // GMP compatible: mpz_perfect_square_p
    static bool isSquare(const Int& value, Int* pSqrt = nullptr);

    // N-th root (floor(value^(1/n)))
    // Special states: NaN -> NaN, n<=0 -> NaN, negative with even n -> NaN, +inf -> +inf, 0 -> 0
    // Algorithm: precision-doubling Newton + remainder tracking
    //   x_{k+1} = ((n-1)*x_k + value/x_k^(n-1)) / n
    // GMP compatible: mpz_root
    static Int nthRoot(const Int& value, uint32_t n);

    // Return n-th root and remainder: value = root^n + remainder
    // remainder always satisfies 0 <= remainder
    // GMP compatible: mpz_rootrem
    static Int nthRootRem(const Int& value, uint32_t n, Int& remainder);

    // Perfect-power check: does there exist k such that value = b^k (k >= 2)?
    // Returns true if it exists, and stores b in pBase if pBase != nullptr.
    // Algorithm: try k = 2..bitLength, verify with nthRoot.
    static bool isPerfectPower(const Int& value, Int* pBase = nullptr, uint32_t* pExp = nullptr);

private:
    // Compute square root in double precision (for MSB < 52 bits)
    static Int sqrt_small(const Int& value);

    // Unified implementation of n-th root + optional remainder.
    // If pRemainder != nullptr, stores value - root^n into *pRemainder.
    // Since the verification block already computes R = value - x^n, this
    // eliminates the duplicate pow(root, n) call that nthRootRem would otherwise
    // perform (PLAN_NTHROOT_FOLLOWUP candidate A, follow-up to commit 4f86e68).
    static Int nthRoot_internal(const Int& value, uint32_t n, Int* pRemainder);
};

} // namespace sangi

#endif // SANGI_INT_SQRT_HPP