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

namespace sangi {

    // Forward declaration (to avoid a circular include)
    class Int;

    // Traits for numeric-state management

    // Traits defining the state-propagation rules
    struct numeric_state_propagation_traits {
        // State-propagation rule for binary operations
        static NumericState propagate_binary_op(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // NaN always takes priority
            if (NumericStateTraits::isNaN(lhs) || NumericStateTraits::isNaN(rhs)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // Infinity vs infinity with opposite signs yields NaN
            if ((lhs == NumericState::PositiveInfinity && rhs == NumericState::NegativeInfinity) ||
                (lhs == NumericState::NegativeInfinity && rhs == NumericState::PositiveInfinity)) {
                resultError = NumericError::InfiniteIndeterminate;
                return NumericState::NaN;
            }

            // For operations between an infinity and another value, infinity takes priority
            if (lhs == NumericState::PositiveInfinity || lhs == NumericState::NegativeInfinity) {
                resultError = NumericError::None; // Infinity itself is not an error
                return lhs;
            }
            if (rhs == NumericState::PositiveInfinity || rhs == NumericState::NegativeInfinity) {
                resultError = NumericError::None; // Infinity itself is not an error
                return rhs;
            }

            // Handle complex infinity
            if (lhs == NumericState::ComplexInfinity || rhs == NumericState::ComplexInfinity) {
                resultError = NumericError::None; // Infinity itself is not an error
                return NumericState::ComplexInfinity;
            }

            // Overflow has high priority
            if (lhs == NumericState::Overflow || rhs == NumericState::Overflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }

            // Underflow
            if (lhs == NumericState::Underflow || rhs == NumericState::Underflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }

            // Divergent state
            if (NumericStateTraits::isNotConverged(lhs) || NumericStateTraits::isNotConverged(rhs)) {
                resultError = NumericError::DivergenceError;
                return NumericState::NotConverged;
            }

            // Subnormal
            if (lhs == NumericState::Subnormal || rhs == NumericState::Subnormal) {
                resultError = NumericError::None;
                return NumericState::Subnormal;
            }

            // Normal value (including zero)
            resultError = NumericError::None;
            return NumericState::Normal;
        }

        // State-propagation rule for addition
        static NumericState propagate_add(NumericState lhs, NumericState rhs, NumericError& resultError) {
            return propagate_binary_op(lhs, rhs, resultError);
        }

        // State-propagation rule for subtraction
        static NumericState propagate_subtract(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // For subtraction, flip the sign of the rhs argument and apply the addition rule
            if (rhs == NumericState::PositiveInfinity)
                rhs = NumericState::NegativeInfinity;
            else if (rhs == NumericState::NegativeInfinity)
                rhs = NumericState::PositiveInfinity;

            return propagate_add(lhs, rhs, resultError);
        }

        // State-propagation rule for multiplication
        static NumericState propagate_multiply(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // NaN always takes priority
            if (NumericStateTraits::isNaN(lhs) || NumericStateTraits::isNaN(rhs)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // Infinity * 0 = NaN
            // Since we cannot test the actual value here, the caller must check isZero()
            // Therefore no zero check is needed here (the value's sign cannot be inferred from state alone)

            // Infinity * infinity = signed infinity
            if (NumericStateTraits::isInfinite(lhs) && NumericStateTraits::isInfinite(rhs)) {
                bool positive = (lhs == NumericState::PositiveInfinity && rhs == NumericState::PositiveInfinity) ||
                    (lhs == NumericState::NegativeInfinity && rhs == NumericState::NegativeInfinity);
                resultError = NumericError::None;
                return positive ? NumericState::PositiveInfinity : NumericState::NegativeInfinity;
            }

            // Infinity * normal = signed infinity
            if (lhs == NumericState::PositiveInfinity) {
                if (rhs == NumericState::Normal) {
                    // We would ideally inspect the value's sign, but it cannot be inferred from state alone; assume positive
                    resultError = NumericError::None;
                    return NumericState::PositiveInfinity;
                }
            }
            if (lhs == NumericState::NegativeInfinity) {
                if (rhs == NumericState::Normal) {
                    // Same as above
                    resultError = NumericError::None;
                    return NumericState::NegativeInfinity;
                }
            }
            if (rhs == NumericState::PositiveInfinity || rhs == NumericState::NegativeInfinity) {
                // Same as above
                if (lhs == NumericState::Normal) {
                    resultError = NumericError::None;
                    return rhs;
                }
            }

            // Handle complex infinity
            if (lhs == NumericState::ComplexInfinity || rhs == NumericState::ComplexInfinity) {
                resultError = NumericError::None;
                return NumericState::ComplexInfinity;
            }

            // Overflow
            if (lhs == NumericState::Overflow || rhs == NumericState::Overflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }

            // Underflow (multiplication of small values)
            if (lhs == NumericState::Underflow || rhs == NumericState::Underflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }

            // Divergent state
            if (NumericStateTraits::isDivergent(lhs) || NumericStateTraits::isDivergent(rhs)) {
                resultError = NumericError::DivergenceError;
                return NumericState::NotConverged;
            }

            // Subnormal
            if (lhs == NumericState::Subnormal || rhs == NumericState::Subnormal) {
                resultError = NumericError::None;
                return NumericState::Subnormal;
            }

            // Normal value (including zero)
            resultError = NumericError::None;
            return NumericState::Normal;
        }

        // State-propagation rule for division
        static NumericState propagate_divide(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // NaN always takes priority
            if (NumericStateTraits::isNaN(lhs) || NumericStateTraits::isNaN(rhs)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // The zero-divisor check cannot be done from state alone, so it is not
            // implemented here; the caller is expected to use isZero() instead

            // Infinity / infinity = NaN
            if (NumericStateTraits::isInfinite(lhs) && NumericStateTraits::isInfinite(rhs)) {
                resultError = NumericError::InfiniteIndeterminate;
                return NumericState::NaN;
            }

            // Normal / infinity = 0 (its magnitude is also indeterminate without the value)

            // Infinity / normal = infinity (sign preserved)
            if (lhs == NumericState::PositiveInfinity && rhs == NumericState::Normal) {
                resultError = NumericError::None;
                return NumericState::PositiveInfinity;
            }
            if (lhs == NumericState::NegativeInfinity && rhs == NumericState::Normal) {
                resultError = NumericError::None;
                return NumericState::NegativeInfinity;
            }

            // Handle complex infinity
            if (lhs == NumericState::ComplexInfinity) {
                if (rhs == NumericState::ComplexInfinity) {
                    resultError = NumericError::InfiniteIndeterminate;
                    return NumericState::NaN;
                }
                resultError = NumericError::None;
                return NumericState::ComplexInfinity;
            }
            if (rhs == NumericState::ComplexInfinity) {
                resultError = NumericError::None;
                return NumericState::Normal; // Typically 0, but indeterminate here
            }

            // Overflow
            if (lhs == NumericState::Overflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }
            if (rhs == NumericState::Underflow) { // Dividing by a small value yields a large one
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }

            // Underflow
            if (lhs == NumericState::Underflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }
            if (rhs == NumericState::Overflow) { // Dividing by a large value yields a small one
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }

            // Divergent state
            if (NumericStateTraits::isDivergent(lhs) || NumericStateTraits::isDivergent(rhs)) {
                resultError = NumericError::DivergenceError;
                return NumericState::Divergent;
            }

            // Normal value (including zero)
            resultError = NumericError::None;
            return NumericState::Normal;
        }

        // State-propagation rule for square root
        static NumericState propagate_sqrt(NumericState state, NumericError& resultError) {
            if (NumericStateTraits::isNaN(state)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // Square root of a negative value = NaN
            if (state == NumericState::NegativeInfinity) {
                resultError = NumericError::NegativeSqrt;
                return NumericState::NaN;
            }

            if (state == NumericState::PositiveInfinity) {
                resultError = NumericError::None;
                return NumericState::PositiveInfinity;
            }

            if (state == NumericState::ComplexInfinity) {
                resultError = NumericError::None;
                return NumericState::ComplexInfinity;
            }

            // Square root of an overflow value may return to a normal value
            if (state == NumericState::Overflow) {
                resultError = NumericError::None;
                return NumericState::Normal; // Or an appropriate overflow state
            }

            // Square root of an underflow value remains underflowed
            if (state == NumericState::Underflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }

            if (NumericStateTraits::isDivergent(state)) {
                resultError = NumericError::DivergenceError;
                return NumericState::Divergent;
            }

            if (state == NumericState::Subnormal) {
                resultError = NumericError::None;
                return NumericState::Subnormal;
            }

            resultError = NumericError::None;
            return NumericState::Normal;
        }

        // State-propagation rule for negation
        static NumericState propagate_negate(NumericState state, NumericError& resultError) {
            if (state == NumericState::PositiveInfinity) {
                resultError = NumericError::None;
                return NumericState::NegativeInfinity;
            }
            if (state == NumericState::NegativeInfinity) {
                resultError = NumericError::None;
                return NumericState::PositiveInfinity;
            }

            // All other states are unchanged
            resultError = NumericError::None;
            return state;
        }

        // State-propagation rule for absolute value
        static NumericState propagate_abs(NumericState state, NumericError& resultError) {
            if (state == NumericState::NegativeInfinity) {
                resultError = NumericError::None;
                return NumericState::PositiveInfinity;
            }

            // abs() only flips negative to positive; other states are unchanged
            resultError = NumericError::None;
            return state;
        }

        // Propagation of divergence detail
        static DivergenceDetail propagate_divergence(DivergenceDetail lhs, DivergenceDetail rhs) {
            // Oscillation has top priority
            if (lhs == DivergenceDetail::Oscillating || rhs == DivergenceDetail::Oscillating)
                return DivergenceDetail::Oscillating;

            // Then divergence
            if (lhs == DivergenceDetail::Diverging || rhs == DivergenceDetail::Diverging)
                return DivergenceDetail::Diverging;

            // Slow convergence
            if (lhs == DivergenceDetail::SlowConvergence || rhs == DivergenceDetail::SlowConvergence)
                return DivergenceDetail::SlowConvergence;

            // Truncation error
            if (lhs == DivergenceDetail::ConditionalConvergenceViolation || rhs == DivergenceDetail::ConditionalConvergenceViolation)
                return DivergenceDetail::ConditionalConvergenceViolation;

            // None if both are undefined
            return DivergenceDetail::None;
        }
    };

    // Vector/matrix traits

    // Vector representation tags
    struct vector_layout_tag {};
    struct dense_vector_tag : vector_layout_tag {};        // dense vector
    struct sparse_vector_tag : vector_layout_tag {};       // sparse vector

    // Matrix representation tags
    struct matrix_layout_tag {};
    struct row_major_tag : matrix_layout_tag {};           // row-major
    struct column_major_tag : matrix_layout_tag {};        // column-major
    struct dense_matrix_tag : matrix_layout_tag {};        // dense matrix
    struct sparse_matrix_tag : matrix_layout_tag {};       // sparse matrix
    struct banded_matrix_tag : matrix_layout_tag {};       // banded matrix
    struct diagonal_matrix_tag : matrix_layout_tag {};     // diagonal matrix
    struct symmetric_matrix_tag : matrix_layout_tag {};    // symmetric matrix
    struct hermitian_matrix_tag : matrix_layout_tag {};    // Hermitian matrix
    struct triangular_matrix_tag : matrix_layout_tag {};   // triangular matrix

    // Base template for vector traits
    template<typename VectorType>
    struct vector_traits {
        // Unknown by default
        static constexpr bool is_recognized = false;
    };

    // Base template for matrix traits
    template<typename MatrixType>
    struct matrix_traits {
        // Unknown by default
        static constexpr bool is_recognized = false;

        // Function that tests whether a matrix is symmetric
        // Fix: unify the function name (isSymmetric → is_symmetric)
        template<typename T>
        static bool is_symmetric(const MatrixType& m) {
            const auto rows = m.rows();
            const auto cols = m.cols();

            // Not square => not symmetric
            if (rows != cols) {
                return false;
            }

            // Check equality of every (i,j) and (j,i) pair
            for (std::size_t i = 0; i < rows; ++i) {
                for (std::size_t j = i + 1; j < cols; ++j) {
                    if (m(i, j) != m(j, i)) {
                        return false;
                    }
                }
            }

            return true;
        }
    };

    // Tag-dispatch functions for matrix/vector operations

    // Dispatch for vector addition
    template<typename VectorType, typename Tag = typename vector_traits<VectorType>::layout_tag>
    struct vector_addition_impl {
        static void apply(const VectorType& a, const VectorType& b, VectorType& result) {
            // Default implementation (element-wise addition)
            const auto size = a.size();
            for (std::size_t i = 0; i < size; ++i) {
                result[i] = a[i] + b[i];
            }
        }
    };

    // Dispatch for matrix addition
    template<typename MatrixType, typename Tag = typename matrix_traits<MatrixType>::layout_tag>
    struct matrix_addition_impl {
        static void apply(const MatrixType& a, const MatrixType& b, MatrixType& result) {
            // Default implementation (element-wise addition)
            const auto rows = a.rows();
            const auto cols = a.cols();
            for (std::size_t i = 0; i < rows; ++i) {
                for (std::size_t j = 0; j < cols; ++j) {
                    result(i, j) = a(i, j) + b(i, j);
                }
            }
        }
    };

    // Dispatch for matrix multiplication
    template<typename MatrixTypeA, typename MatrixTypeB, typename MatrixTypeC,
        typename TagA = typename matrix_traits<MatrixTypeA>::layout_tag,
        typename TagB = typename matrix_traits<MatrixTypeB>::layout_tag>
    struct matrix_multiplication_impl {
        static void apply(const MatrixTypeA& a, const MatrixTypeB& b, MatrixTypeC& result) {
            // Default implementation (triple-loop multiplication)
            const auto m = a.rows();
            const auto n = b.cols();
            const auto k = a.cols();

            // Initialize the result with 0
            for (std::size_t i = 0; i < m; ++i) {
                for (std::size_t j = 0; j < n; ++j) {
                    result(i, j) = numeric_traits<typename MatrixTypeC::value_type>::zero();
                }
            }

            // Matrix multiplication: C = A * B
            for (std::size_t i = 0; i < m; ++i) {
                for (std::size_t j = 0; j < n; ++j) {
                    for (std::size_t l = 0; l < k; ++l) {
                        result(i, j) += a(i, l) * b(l, j);
                    }
                }
            }
        }
    };

    // Mechanism for extending traits to custom types

    // Customization point for numeric traits
    template<typename T, typename = void>
    struct custom_numeric_traits {};

    // Helper that integrates custom numeric traits
    template<typename T>
    struct numeric_traits_base : custom_numeric_traits<T> {
        // Default when no custom traits exist or are not supported
        static constexpr bool is_supported = false;
    };

    // Specialization selector for numeric traits
    template<typename T>
    struct numeric_traits_selector {
        using type = typename std::conditional<
            std::is_same_v<T, float> ||
            std::is_same_v<T, double> ||
            std::is_same_v<T, long double> ||
            std::is_same_v<T, int> ||
            std::is_same_v<T, long> ||
            std::is_same_v<T, long long> ||
            std::is_same_v<T, std::complex<float>> ||
            std::is_same_v<T, std::complex<double>>,
            numeric_traits<T>,     // For standard types, use the built-in traits
            numeric_traits_base<T> // For custom types, use the customization point
                    > ::type;
    };

    // Macro that simplifies customizing the traits
#define SANGI_DECLARE_NUMERIC_TRAITS(Type, IsSupported, IsComplex, IsInteger, IsFloatingPoint) \
template<> \
struct custom_numeric_traits<Type> { \
   using value_type = Type; \
   static constexpr bool is_supported = IsSupported; \
   static constexpr bool is_complex = IsComplex; \
   static constexpr bool is_integer = IsInteger; \
   static constexpr bool is_floating_point = IsFloatingPoint; \
};

// Customization macro for state-management traits
#define SANGI_DECLARE_STATE_TRAITS(Type) \
template<> \
struct custom_numeric_traits<Type> { \
   using value_type = Type; \
   using state_type = NumericState; \
   using error_type = NumericError; \
   using divergence_detail_type = DivergenceDetail; \
   static constexpr bool has_state_management = true; \
};

// Helper functions and metafunctions
//-----------------------------------------------------------------------------

// Helpers for concept checking
    template<typename T>
    constexpr bool is_field_v = concepts::Field<T>;

    template<typename T>
    constexpr bool is_ordered_field_v = concepts::OrderedField<T>;

    template<typename V, typename S>
    constexpr bool is_vector_space_v = concepts::VectorSpace<V, S>;

    template<typename M, typename T>
    constexpr bool is_matrix_of_v = concepts::MatrixOf<M, T>;

    template<typename V, typename T>
    constexpr bool is_vector_of_v = concepts::VectorOf<V, T>;

    // Helpers for state-management checks
    template<typename T>
    constexpr bool has_numeric_state_management_v = concepts::HasNumericState<T>;

    template<typename T>
    constexpr bool has_divergence_handling_v = concepts::HasDivergenceHandling<T>;

    template<typename T>
    constexpr bool has_overflow_detection_v = concepts::HasOverflowDetection<T>;

    // Duck-typing support
    // These constants test whether a given operation is defined for the type
    template<typename T>
    constexpr bool has_addition_v = requires(T a, T b) { { a + b } -> std::convertible_to<T>; };

    template<typename T>
    constexpr bool has_multiplication_v = requires(T a, T b) { { a* b } -> std::convertible_to<T>; };

    template<typename T>
    constexpr bool has_division_v = requires(T a, T b) { { a / b } -> std::convertible_to<T>; };

    template<typename T>
    constexpr bool has_negation_v = requires(T a) { { -a } -> std::convertible_to<T>; };

    template<typename T>
    constexpr bool has_equality_v = requires(T a, T b) { { a == b } -> std::convertible_to<bool>; };

    template<typename T>
    constexpr bool has_comparison_v = requires(T a, T b) { { a < b } -> std::convertible_to<bool>; };

    template<typename T>
    constexpr bool has_abs_v = requires(T a) {
        { std::abs(static_cast<double>(a)) } -> std::convertible_to<double>;
    };

    template<typename T>
    constexpr bool has_sqrt_v = requires(T a) {
        { std::sqrt(static_cast<double>(a)) } -> std::convertible_to<double>;
    };

    // Helpers for state management
    template<typename T>
    constexpr bool has_is_normal_v = requires(T a) { { a.isNormal() } -> std::convertible_to<bool>; };

    template<typename T>
    constexpr bool has_is_nan_v = requires(T a) { { a.isNaN() } -> std::convertible_to<bool>; };

    template<typename T>
    constexpr bool has_is_infinite_v = requires(T a) { { a.isInfinite() } -> std::convertible_to<bool>; };

    template<typename T>
    constexpr bool has_is_divergent_v = requires(T a) { { a.isDivergent() } -> std::convertible_to<bool>; };

} // namespace sangi