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

namespace calx {

    // Int型の前方宣言（循環参照を避けるため）
    class Int;

    // 数値状態管理のトレイト

    // 状態伝播ルールを定義するトレイト
    struct numeric_state_propagation_traits {
        // 二項演算の状態伝播ルール
        static NumericState propagate_binary_op(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // NaNは常に優先される
            if (NumericStateTraits::isNaN(lhs) || NumericStateTraits::isNaN(rhs)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // 無限大同士の演算で、符号が異なる場合はNaN
            if ((lhs == NumericState::PositiveInfinity && rhs == NumericState::NegativeInfinity) ||
                (lhs == NumericState::NegativeInfinity && rhs == NumericState::PositiveInfinity)) {
                resultError = NumericError::InfiniteIndeterminate;
                return NumericState::NaN;
            }

            // 無限大と他の値の演算では無限大が優先される
            if (lhs == NumericState::PositiveInfinity || lhs == NumericState::NegativeInfinity) {
                resultError = NumericError::None; // 無限大自体はエラーではない
                return lhs;
            }
            if (rhs == NumericState::PositiveInfinity || rhs == NumericState::NegativeInfinity) {
                resultError = NumericError::None; // 無限大自体はエラーではない
                return rhs;
            }

            // 複素無限大の処理
            if (lhs == NumericState::ComplexInfinity || rhs == NumericState::ComplexInfinity) {
                resultError = NumericError::None; // 無限大自体はエラーではない
                return NumericState::ComplexInfinity;
            }

            // オーバーフローは優先順位が高い
            if (lhs == NumericState::Overflow || rhs == NumericState::Overflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }

            // アンダーフロー
            if (lhs == NumericState::Underflow || rhs == NumericState::Underflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }

            // 発散状態
            if (NumericStateTraits::isNotConverged(lhs) || NumericStateTraits::isNotConverged(rhs)) {
                resultError = NumericError::DivergenceError;
                return NumericState::NotConverged;
            }

            // サブノーマル
            if (lhs == NumericState::Subnormal || rhs == NumericState::Subnormal) {
                resultError = NumericError::None;
                return NumericState::Subnormal;
            }

            // 通常値（ゼロを含む）
            resultError = NumericError::None;
            return NumericState::Normal;
        }

        // 加算の状態伝播ルール
        static NumericState propagate_add(NumericState lhs, NumericState rhs, NumericError& resultError) {
            return propagate_binary_op(lhs, rhs, resultError);
        }

        // 減算の状態伝播ルール
        static NumericState propagate_subtract(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // 減算では、rhs引数の符号を反転して加算のルールを適用
            if (rhs == NumericState::PositiveInfinity)
                rhs = NumericState::NegativeInfinity;
            else if (rhs == NumericState::NegativeInfinity)
                rhs = NumericState::PositiveInfinity;

            return propagate_add(lhs, rhs, resultError);
        }

        // 乗算の状態伝播ルール
        static NumericState propagate_multiply(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // NaNは常に優先される
            if (NumericStateTraits::isNaN(lhs) || NumericStateTraits::isNaN(rhs)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // 無限大 * 0 = NaN
            // ここでは、具体的な値ではなく、isZero()メソッドを呼び出す側で判定する必要がある
            // そのため、ここでのゼロチェックは不要となる（値の符号は状態から判断できない）

            // 無限大 * 無限大 = 符号に基づく無限大
            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;
            }

            // 無限大 * 通常値 = 符号に基づく無限大
            if (lhs == NumericState::PositiveInfinity) {
                if (rhs == NumericState::Normal) {
                    // ここで本来は値の符号を見るべきだが、状態だけではわからないため通常は正を仮定
                    resultError = NumericError::None;
                    return NumericState::PositiveInfinity;
                }
            }
            if (lhs == NumericState::NegativeInfinity) {
                if (rhs == NumericState::Normal) {
                    // 同上
                    resultError = NumericError::None;
                    return NumericState::NegativeInfinity;
                }
            }
            if (rhs == NumericState::PositiveInfinity || rhs == NumericState::NegativeInfinity) {
                // 同上
                if (lhs == NumericState::Normal) {
                    resultError = NumericError::None;
                    return rhs;
                }
            }

            // 複素無限大の処理
            if (lhs == NumericState::ComplexInfinity || rhs == NumericState::ComplexInfinity) {
                resultError = NumericError::None;
                return NumericState::ComplexInfinity;
            }

            // オーバーフロー
            if (lhs == NumericState::Overflow || rhs == NumericState::Overflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }

            // アンダーフロー（小さい値同士の乗算）
            if (lhs == NumericState::Underflow || rhs == NumericState::Underflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }

            // 発散状態
            if (NumericStateTraits::isDivergent(lhs) || NumericStateTraits::isDivergent(rhs)) {
                resultError = NumericError::DivergenceError;
                return NumericState::NotConverged;
            }

            // サブノーマル
            if (lhs == NumericState::Subnormal || rhs == NumericState::Subnormal) {
                resultError = NumericError::None;
                return NumericState::Subnormal;
            }

            // 通常値（ゼロを含む）
            resultError = NumericError::None;
            return NumericState::Normal;
        }

        // 除算の状態伝播ルール
        static NumericState propagate_divide(NumericState lhs, NumericState rhs, NumericError& resultError) {
            // NaNは常に優先される
            if (NumericStateTraits::isNaN(lhs) || NumericStateTraits::isNaN(rhs)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // 0除算チェックは、具体的な値でないとできないため、
            // ここでは実装せず、呼び出し側でisZero()メソッドを使って判定する

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

            // 通常値 / 無限大 = 0（これも値がわからないとサイズを判断できない）

            // 無限大 / 通常値 = 無限大（符号保持）
            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;
            }

            // 複素無限大の処理
            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; // 通常は0になるがここでは判断できない
            }

            // オーバーフロー
            if (lhs == NumericState::Overflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }
            if (rhs == NumericState::Underflow) { // 小さい値で割ると大きな値になる
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Overflow;
            }

            // アンダーフロー
            if (lhs == NumericState::Underflow) {
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }
            if (rhs == NumericState::Overflow) { // 大きな値で割ると小さな値になる
                resultError = NumericError::OutOfRangeInput;
                return NumericState::Underflow;
            }

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

            // 通常値（ゼロを含む）
            resultError = NumericError::None;
            return NumericState::Normal;
        }

        // 平方根演算の状態伝播ルール
        static NumericState propagate_sqrt(NumericState state, NumericError& resultError) {
            if (NumericStateTraits::isNaN(state)) {
                resultError = NumericError::NaNPropagation;
                return NumericState::NaN;
            }

            // 負の値の平方根 = 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;
            }

            // オーバーフローの平方根は通常値になる可能性がある
            if (state == NumericState::Overflow) {
                resultError = NumericError::None;
                return NumericState::Normal; // または適切なオーバーフロー状態
            }

            // アンダーフローの平方根はさらにアンダーフローになる
            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;
        }

        // 符号反転の状態伝播ルール
        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;
            }

            // その他の状態は変わらない
            resultError = NumericError::None;
            return state;
        }

        // 絶対値の状態伝播ルール
        static NumericState propagate_abs(NumericState state, NumericError& resultError) {
            if (state == NumericState::NegativeInfinity) {
                resultError = NumericError::None;
                return NumericState::PositiveInfinity;
            }

            // 絶対値は負の値を正にするだけなので、他の状態は変わらない
            resultError = NumericError::None;
            return state;
        }

        // 発散詳細の伝播
        static DivergenceDetail propagate_divergence(DivergenceDetail lhs, DivergenceDetail rhs) {
            // 振動が最も優先
            if (lhs == DivergenceDetail::Oscillating || rhs == DivergenceDetail::Oscillating)
                return DivergenceDetail::Oscillating;

            // 次に発散
            if (lhs == DivergenceDetail::Diverging || rhs == DivergenceDetail::Diverging)
                return DivergenceDetail::Diverging;

            // 収束速度が遅い場合
            if (lhs == DivergenceDetail::SlowConvergence || rhs == DivergenceDetail::SlowConvergence)
                return DivergenceDetail::SlowConvergence;

            // 打ち切り誤差
            if (lhs == DivergenceDetail::ConditionalConvergenceViolation || rhs == DivergenceDetail::ConditionalConvergenceViolation)
                return DivergenceDetail::ConditionalConvergenceViolation;

            // どちらも未定義ならNone
            return DivergenceDetail::None;
        }
    };

    // ベクトル/行列特性

    // ベクトル表現タグ
    struct vector_layout_tag {};
    struct dense_vector_tag : vector_layout_tag {};        // 密ベクトル
    struct sparse_vector_tag : vector_layout_tag {};       // 疎ベクトル

    // 行列表現タグ
    struct matrix_layout_tag {};
    struct row_major_tag : matrix_layout_tag {};           // 行優先
    struct column_major_tag : matrix_layout_tag {};        // 列優先
    struct dense_matrix_tag : matrix_layout_tag {};        // 密行列
    struct sparse_matrix_tag : matrix_layout_tag {};       // 疎行列
    struct banded_matrix_tag : matrix_layout_tag {};       // 帯行列
    struct diagonal_matrix_tag : matrix_layout_tag {};     // 対角行列
    struct symmetric_matrix_tag : matrix_layout_tag {};    // 対称行列
    struct hermitian_matrix_tag : matrix_layout_tag {};    // エルミート行列
    struct triangular_matrix_tag : matrix_layout_tag {};   // 三角行列

    // ベクトル特性の基本テンプレート
    template<typename VectorType>
    struct vector_traits {
        // デフォルトは不明
        static constexpr bool is_recognized = false;
    };

    // 行列特性の基本テンプレート
    template<typename MatrixType>
    struct matrix_traits {
        // デフォルトは不明
        static constexpr bool is_recognized = false;

        // 行列が対称行列かどうかを判定する関数
        // 修正：関数名を統一化（isSymmetric → is_symmetric）
        template<typename T>
        static bool is_symmetric(const MatrixType& m) {
            const auto rows = m.rows();
            const auto cols = m.cols();

            // 正方行列でないなら対称行列ではない
            if (rows != cols) {
                return false;
            }

            // すべての(i,j)と(j,i)の組について一致するか確認
            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;
        }
    };

    // 行列/ベクトル演算のためのタグディスパッチ関数

    // ベクトル加算用ディスパッチ
    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) {
            // デフォルト実装（要素ごとの加算）
            const auto size = a.size();
            for (std::size_t i = 0; i < size; ++i) {
                result[i] = a[i] + b[i];
            }
        }
    };

    // 行列加算用ディスパッチ
    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) {
            // デフォルト実装（要素ごとの加算）
            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);
                }
            }
        }
    };

    // 行列乗算用ディスパッチ
    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) {
            // デフォルト実装（三重ループによる乗算）
            const auto m = a.rows();
            const auto n = b.cols();
            const auto k = a.cols();

            // 結果を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();
                }
            }

            // 行列乗算: 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);
                    }
                }
            }
        }
    };

    // カスタム型のためのトレイト拡張メカニズム

    // 数値トレイトのカスタマイズポイント
    template<typename T, typename = void>
    struct custom_numeric_traits {};

    // カスタム数値トレイトを統合するためのヘルパー
    template<typename T>
    struct numeric_traits_base : custom_numeric_traits<T> {
        // カスタムトレイトがないか、サポートされていない場合のデフォルト
        static constexpr bool is_supported = false;
    };

    // 数値トレイトの特殊化の選択
    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>,     // 標準型の場合は組み込みトレイトを使用
            numeric_traits_base<T> // カスタム型の場合はカスタマイズポイントを使用
                    > ::type;
    };

    // トレイトのカスタマイズを簡単にするマクロ
#define CALX_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; \
};

// 状態管理トレイトのカスタマイズマクロ
#define CALX_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; \
};

// ヘルパー関数とメタ関数
//-----------------------------------------------------------------------------

// コンセプトチェック用のヘルパー
    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>;

    // 状態管理のチェック用ヘルパー
    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>;

    // ダックタイピングサポート
    // これらの関数は、特定の操作が型に対して定義されているかを検査します
    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>;
    };

    // 状態管理用ヘルパー
    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 calx