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

namespace sangi {

    // Forward declaration of the Int type (to avoid circular references)
    class Int;

    namespace concepts {
        // Forward declaration of basic concepts
        template<typename T>
        concept Field = requires(T a, T b) {
            { a + b } -> std::convertible_to<T>;
            { a - b } -> std::convertible_to<T>;
            { a* b } -> std::convertible_to<T>;
            { a / b } -> std::convertible_to<T>;
            { -a } -> std::convertible_to<T>;
            { T(0) } -> std::convertible_to<T>;
            { T(1) } -> std::convertible_to<T>;
        };

        template<typename T>
        concept OrderedField = Field<T> && requires(T a, T b) {
            { a < b } -> std::convertible_to<bool>;
            { a <= b } -> std::convertible_to<bool>;
            { a > b } -> std::convertible_to<bool>;
            { a >= b } -> std::convertible_to<bool>;
        };

        template<typename V, typename S>
        concept VectorSpace = requires(V v, V w, S s) {
            { v + w } -> std::convertible_to<V>;
            { v - w } -> std::convertible_to<V>;
            { s* v } -> std::convertible_to<V>;
            { v* s } -> std::convertible_to<V>;
            { -v } -> std::convertible_to<V>;
        };

        template<typename M, typename T>
        concept MatrixOf = requires(M m, std::size_t i, std::size_t j) {
            typename M::value_type;
            { m(i, j) } -> std::convertible_to<T>;
            { m.rows() } -> std::convertible_to<std::size_t>;
            { m.cols() } -> std::convertible_to<std::size_t>;
        };

        template<typename V, typename T>
        concept VectorOf = requires(V v, std::size_t i) {
            typename V::value_type;
            { v[i] } -> std::convertible_to<T>;
            { v.size() } -> std::convertible_to<std::size_t>;
        };

        // Support concept for numeric-state management
        template<typename T>
        concept HasDivergenceHandling = requires(T a) {
            { numeric_state_traits<T>::isNotConverged(a) } -> std::convertible_to<bool>;
            { numeric_state_traits<T>::getDivergenceDetail(a) } -> std::convertible_to<DivergenceDetail>;
        };

        template<typename T>
        concept HasOverflowDetection = requires(T a) {
            { numeric_state_traits<T>::isOverflow(a) } -> std::convertible_to<bool>;
        };

        // Types where numeric-state management can be detected and handled
        template<typename T>
        concept HasNumericState = requires(T a) {
            { a.isNormal() } -> std::convertible_to<bool>;
            { a.isNaN() } -> std::convertible_to<bool>;
            { a.isInfinite() } -> std::convertible_to<bool>;
            { a.getState() } -> std::convertible_to<typename T::state_type>;
        };
    }

    // Numeric-type traits

    // Identification of numeric-type categories
    struct numeric_category_tag {};
    struct real_tag : numeric_category_tag {};
    struct complex_tag : numeric_category_tag {};
    struct integer_tag : real_tag {};
    struct floating_point_tag : real_tag {};

    // Base template for numeric-type traits
    template<typename T>
    struct numeric_traits {
        // Undefined by default
        static constexpr bool is_supported = false;

        static constexpr T abs(T x) { return x < T(0) ? -x : x; }

        /// Pivot selection: is a a better pivot than b?
        /// Default: |a| > |b| (largest magnitude is stable)
        /// For exact types such as Rational, override to minimize height
        static bool pivotBetter(const T& a, const T& b) {
            return abs(a) > abs(b);
        }
    };

    // Specialization: floating-point types
    template<>
    struct numeric_traits<float> {
        using value_type = float;
        using category = floating_point_tag;

        static constexpr bool is_supported = true;
        static constexpr bool is_complex = false;
        static constexpr bool is_integer = false;
        static constexpr bool is_floating_point = true;

        static constexpr float zero() { return 0.0f; }
        static constexpr float one() { return 1.0f; }
        static constexpr float epsilon() { return std::numeric_limits<float>::epsilon(); }
        static constexpr float infinity() { return std::numeric_limits<float>::infinity(); }
        static constexpr float quiet_NaN() { return std::numeric_limits<float>::quiet_NaN(); }

        static constexpr float abs(float value) { return value < 0.0f ? -value : value; }
        static constexpr float conj(float value) { return value; }
        static constexpr float norm(float value) { return value * value; }
        static bool pivotBetter(float a, float b) { return abs(a) > abs(b); }

        // Numeric-state-related utilities - constexpr removed
        static bool isNaN(float value) { return std::isnan(value); }
        static bool isInfinite(float value) { return std::isinf(value); }
        static bool isNormal(float value) { return std::isnormal(value); }
        static bool isFinite(float value) { return std::isfinite(value); }

        static NumericState getState(float value) {
            if (std::isnan(value)) return NumericState::NaN;
            if (std::isinf(value)) {
                return value > 0 ? NumericState::PositiveInfinity : NumericState::NegativeInfinity;
            }
            if (std::isnormal(value) || value == 0.0f) return NumericState::Normal;
            return NumericState::Subnormal;
        }

        static NumericError getError(float value) {
            if (std::isnan(value)) return NumericError::ExplicitNaN;
            return NumericError::None;
        }
    };

    template<>
    struct numeric_traits<double> {
        using value_type = double;
        using category = floating_point_tag;

        static constexpr bool is_supported = true;
        static constexpr bool is_complex = false;
        static constexpr bool is_integer = false;
        static constexpr bool is_floating_point = true;

        static constexpr double zero() { return 0.0; }
        static constexpr double one() { return 1.0; }
        static constexpr double epsilon() { return std::numeric_limits<double>::epsilon(); }
        static constexpr double infinity() { return std::numeric_limits<double>::infinity(); }
        static constexpr double quiet_NaN() { return std::numeric_limits<double>::quiet_NaN(); }

        static constexpr double abs(double value) { return value < 0.0 ? -value : value; }
        static constexpr double conj(double value) { return value; }
        static constexpr double norm(double value) { return value * value; }
        static bool pivotBetter(double a, double b) { return abs(a) > abs(b); }

        // Numeric-state-related utilities - constexpr removed
        static bool isNaN(double value) { return std::isnan(value); }
        static bool isInfinite(double value) { return std::isinf(value); }
        static bool isNormal(double value) { return std::isnormal(value); }
        static bool isFinite(double value) { return std::isfinite(value); }

        static NumericState getState(double value) {
            if (std::isnan(value)) return NumericState::NaN;
            if (std::isinf(value)) {
                return value > 0 ? NumericState::PositiveInfinity : NumericState::NegativeInfinity;
            }
            if (std::isnormal(value) || value == 0.0) return NumericState::Normal;
            return NumericState::Subnormal;
        }

        static NumericError getError(double value) {
            if (std::isnan(value)) return NumericError::ExplicitNaN;
            return NumericError::None;
        }
    };

    // Extended-precision floating-point type (long double)
    template<>
    struct numeric_traits<long double> {
        using value_type = long double;
        using category = floating_point_tag;

        static constexpr bool is_supported = true;
        static constexpr bool is_complex = false;
        static constexpr bool is_integer = false;
        static constexpr bool is_floating_point = true;

        static constexpr long double zero() { return 0.0L; }
        static constexpr long double one() { return 1.0L; }
        static constexpr long double epsilon() { return std::numeric_limits<long double>::epsilon(); }
        static constexpr long double infinity() { return std::numeric_limits<long double>::infinity(); }
        static constexpr long double quiet_NaN() { return std::numeric_limits<long double>::quiet_NaN(); }

        static constexpr long double abs(long double value) { return value < 0.0L ? -value : value; }
        static constexpr long double conj(long double value) { return value; }
        static constexpr long double norm(long double value) { return value * value; }
        static bool pivotBetter(long double a, long double b) { return abs(a) > abs(b); }

        // Numeric-state-related utilities - constexpr removed
        static bool isNaN(long double value) { return std::isnan(value); }
        static bool isInfinite(long double value) { return std::isinf(value); }
        static bool isNormal(long double value) { return std::isnormal(value); }
        static bool isFinite(long double value) { return std::isfinite(value); }

        static NumericState getState(long double value) {
            if (std::isnan(value)) return NumericState::NaN;
            if (std::isinf(value)) {
                return value > 0 ? NumericState::PositiveInfinity : NumericState::NegativeInfinity;
            }
            if (std::isnormal(value) || value == 0.0L) return NumericState::Normal;
            return NumericState::Subnormal;
        }

        static NumericError getError(long double value) {
            if (std::isnan(value)) return NumericError::ExplicitNaN;
            return NumericError::None;
        }
    };

    // Integer-type specialization (for int)
    template<>
    struct numeric_traits<int> {
        using value_type = int;
        using category = integer_tag;

        static constexpr bool is_supported = true;
        static constexpr bool is_complex = false;
        static constexpr bool is_integer = true;
        static constexpr bool is_floating_point = false;

        static constexpr int zero() { return 0; }
        static constexpr int one() { return 1; }
        static constexpr int epsilon() { return 1; } // The smallest unit of an integer is 1

        // Integer types cannot represent "infinity"; return the maximum value
        static constexpr int infinity() { return std::numeric_limits<int>::max(); }

        // Integer types cannot represent NaN (return a special value as an error state)
        static constexpr int quiet_NaN() { return std::numeric_limits<int>::min(); }

        static constexpr int abs(int value) {
            return value < 0 ? (value == std::numeric_limits<int>::min() ?
                std::numeric_limits<int>::max() : -value) : value;
        }

        static constexpr int conj(int value) { return value; } // The conjugate of an integer is itself
        static constexpr int norm(int value) { return value * value; } // The norm of an integer is its square
        static bool pivotBetter(int a, int b) { return abs(a) > abs(b); }

        static NumericState getState(int) {
            return NumericState::Normal; // Integer types are always Normal
        }

        static NumericError getError(int) {
            return NumericError::None; // Default error for integer types is none
        }

        static int getSign(int value) {
            return value < 0 ? -1 : (value > 0 ? 1 : 0);
        }
    };

    // Specializations for other integer types
    template<>
    struct numeric_traits<long> : numeric_traits<int> {
        using value_type = long;

        static constexpr long zero() { return 0L; }
        static constexpr long one() { return 1L; }
        static constexpr long abs(long value) { return value < 0L ? -value : value; }

        static NumericState getState(long) {
            return NumericState::Normal;
        }

        static NumericError getError(long) {
            return NumericError::None;
        }

        static int getSign(long value) {
            return value < 0L ? -1 :
                value > 0L ? 1 : 0;
        }
    };

    template<>
    struct numeric_traits<long long> : numeric_traits<int> {
        using value_type = long long;

        static constexpr long long zero() { return 0LL; }
        static constexpr long long one() { return 1LL; }
        static constexpr long long abs(long long value) { return value < 0LL ? -value : value; }

        static NumericState getState(long long) {
            return NumericState::Normal;
        }

        static NumericError getError(long long) {
            return NumericError::None;
        }

        static int getSign(long long value) {
            return value < 0LL ? -1 :
                value > 0LL ? 1 : 0;
        }
    };

    // Forward definition of the Int class and the traits framework
    // The actual implementation lives in IntBase.hpp
    // A provisional implementation is provided here to avoid compile errors
    template<>
    struct numeric_traits<Int>;

    // Forward definition of the numeric-state traits specialization for the Int type
    template<>
    struct numeric_state_traits<Int>;

    // Specialization: complex-number types
    template<typename T>
    struct numeric_traits<std::complex<T>> {
        using value_type = std::complex<T>;
        using real_type = T;
        using category = complex_tag;

        static constexpr bool is_supported = numeric_traits<T>::is_supported;
        static constexpr bool is_complex = true;
        static constexpr bool is_integer = false;
        static constexpr bool is_floating_point = numeric_traits<T>::is_floating_point;

        static std::complex<T> zero() { return std::complex<T>(0, 0); }
        static std::complex<T> one() { return std::complex<T>(1, 0); }
        static T epsilon() { return numeric_traits<T>::epsilon(); }

        // constexpr removed (absolute value is the square root of the sum of squares of real and imaginary parts)
        static T abs(const std::complex<T>& value) {
            const T real = value.real();
            const T imag = value.imag();
            return std::sqrt(real * real + imag * imag);
        }

        static std::complex<T> conj(const std::complex<T>& value) { return std::conj(value); }
        static T norm(const std::complex<T>& value) { return std::norm(value); }
        static bool pivotBetter(const std::complex<T>& a, const std::complex<T>& b) {
            return abs(a) > abs(b);
        }

        // Numeric-state utilities for complex types - constexpr removed
        static bool isNaN(const std::complex<T>& value) {
            return numeric_traits<T>::isNaN(value.real()) || numeric_traits<T>::isNaN(value.imag());
        }

        static bool isInfinite(const std::complex<T>& value) {
            return numeric_traits<T>::isInfinite(value.real()) || numeric_traits<T>::isInfinite(value.imag());
        }

        static bool isNormal(const std::complex<T>& value) {
            return !isNaN(value) && !isInfinite(value) && (value.real() != 0 || value.imag() != 0);
        }

        static bool isFinite(const std::complex<T>& value) {
            return numeric_traits<T>::isFinite(value.real()) && numeric_traits<T>::isFinite(value.imag());
        }

        static NumericState getState(const std::complex<T>& value) {
            if (isNaN(value)) return NumericState::NaN;

            bool real_inf = numeric_traits<T>::isInfinite(value.real());
            bool imag_inf = numeric_traits<T>::isInfinite(value.imag());

            if (real_inf || imag_inf) {
                // Complex infinity has no sign; return a special state based on the sign of real/imaginary parts
                if (real_inf && value.real() > 0) return NumericState::PositiveInfinity;
                if (real_inf && value.real() < 0) return NumericState::NegativeInfinity;
                if (imag_inf) return NumericState::ComplexInfinity;
                return NumericState::ComplexInfinity;
            }

            return NumericState::Normal;
        }

        static NumericError getError(const std::complex<T>& value) {
            // For NaN, clarify the error cause
            if (isNaN(value)) {
                return NumericError::ExplicitNaN;
            }
            // Even for infinity, clarify the error cause
            if (isInfinite(value)) {
                return NumericError::OutOfRangeInput;
            }
            return NumericError::None;
        }

        static int getSign(const std::complex<T>& value) {
            // The sign of a complex number prefers the sign of the real part
            if (value.real() != 0) {
                return value.real() < 0 ? -1 :
                    value.real() > 0 ? 1 : 0;
            }

            // If the real part is zero, use the sign of the imaginary part
            if (value.imag() != 0) {
                return value.imag() < 0 ? -1 : 1;
            }

            // If both are zero, the sign is zero
            return 0;
        }
    };

} // namespace sangi