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

// root_finding_base.hpp
#ifndef SANGI_ROOT_FINDING_BASE_HPP
#define SANGI_ROOT_FINDING_BASE_HPP

#include <math/core/common.hpp>
#include <math/core/traits.hpp>
#include <math/concepts/algebraic_concepts.hpp>
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>

#include <functional>
#include <optional>
#include <limits>
#include <cmath>
#include <tuple>

namespace sangi {

/**
 * @brief Structure for convergence checking
 * Combines multiple criteria to perform convergence checking
 */
template<concepts::OrderedField T>
struct ConvergenceCriteria {
    // Absolute tolerance of the function value
    T abs_ftol = std::numeric_limits<T>::epsilon() * 100;

    // Absolute tolerance of the variable value
    T abs_xtol = std::numeric_limits<T>::epsilon() * 100;

    // Relative tolerance of the function value
    T rel_ftol = std::numeric_limits<T>::epsilon() * 1000;

    // Relative tolerance of the variable value
    T rel_xtol = std::numeric_limits<T>::epsilon() * 1000;

    // Maximum number of iterations
    size_t max_iterations = 100;

    // Function-value check for convergence
    bool is_function_converged(T fx, T fx_prev = T(0)) const {
        // Absolute-value check
        if (std::abs(fx) < abs_ftol) {
            return true;
        }

        // Relative-change check - only when a previous value exists
        if (fx_prev != T(0)) {
            T rel_change = std::abs((fx - fx_prev)) /
                (std::abs(fx_prev) + std::numeric_limits<T>::epsilon());
            if (rel_change < rel_ftol) {
                return true;
            }
        }

        return false;
    }

    // Variable-value check for convergence
    bool is_variable_converged(T x, T x_prev) const {
        // Absolute-change check
        if (std::abs(x - x_prev) < abs_xtol) {
            return true;
        }

        // Relative-change check
        T scale = std::max(T(1), std::abs(x_prev));
        if (std::abs((x - x_prev) / scale) < rel_xtol) {
            return true;
        }

        return false;
    }

    // Convergence check by interval width (used in bisection etc.)
    bool is_interval_converged(T a, T b) const {
        return std::abs(b - a) < abs_xtol ||
            std::abs((b - a) / std::max(T(1), std::max(std::abs(a), std::abs(b)))) < rel_xtol;
    }

    // Composite convergence check (considers both function value and variable value)
    bool is_converged(T fx, T x, T x_prev) const {
        return is_function_converged(fx) && is_variable_converged(x, x_prev);
    }

    // Convergence check for vector-valued functions
    template<typename V>
        requires concepts::VectorOf<V, T>
    bool is_vector_converged(const V& f, const V& x, const V& x_prev) const {
        T f_norm = norm2(f);
        T dx_norm = norm2(x - x_prev);
        T x_norm = norm2(x);

        return f_norm < abs_ftol ||
            dx_norm < abs_xtol ||
            dx_norm < rel_xtol * (x_norm + abs_xtol);
    }
};

/**
 * @brief Structure representing a root-finding result (unified scalar / vector version)
 *
 * If R is a scalar type then error_type = R,
 * and if R is a vector type (having V::value_type) then error_type = V::value_type.
 */
namespace detail {
    template<typename R, typename = void>
    struct root_error_type { using type = R; };
    template<typename R>
    struct root_error_type<R, std::void_t<typename R::value_type>> {
        using type = typename R::value_type;
    };
}

template<typename R>
struct RootFindingResult {
    using error_type = typename detail::root_error_type<R>::type;
    std::optional<R> root;        // Approximate value of the root (nullopt on failure)
    bool converged;               // Whether it converged
    size_t iterations;            // Number of iterations performed
    error_type error_estimate;    // Estimated error

    // Helper to create a success result
    static RootFindingResult success(const R& root_value, size_t iter_count, error_type err = error_type(0)) {
        return {root_value, true, iter_count, err};
    }

    // Helper to create a failure result
    static RootFindingResult failure(size_t iter_count, error_type err = error_type(0)) {
        return {std::nullopt, false, iter_count, err};
    }

    // Optional success result (conditional success)
    static RootFindingResult partial_success(const R& root_value, bool conv, size_t iter_count, error_type err = error_type(0)) {
        return {root_value, conv, iter_count, err};
    }
};

// MKL availability check
inline constexpr bool has_mkl_support = SANGI_HAS_MKL == 1;

} // namespace sangi

#endif // SANGI_ROOT_FINDING_BASE_HPP