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

// Interval.hpp
// 区間演算 (Interval Arithmetic)
//
// IEEE 1788 スタイルの区間型。数値計算の信頼性保証に使用。
// 全演算で区間の包含性 (containment property) を保証。
//
// クラス:
//   Interval<T> — 閉区間 [lo, hi]
//
// 対応:
//   四則演算, 集合演算, 数学関数, 区間ニュートン法

#ifndef CALX_INTERVAL_HPP
#define CALX_INTERVAL_HPP

#include <cmath>
#include <limits>
#include <algorithm>
#include <stdexcept>
#include <ostream>
#include <cfenv>
#include <concepts>
#include <functional>
#include <optional>
#include <numbers>

namespace calx {

// ================================================================
// 丸め方向制御ヘルパー
// ================================================================

namespace detail {

// RAII で丸めモードを一時変更
class RoundingGuard {
public:
    explicit RoundingGuard(int mode) {
        old_ = std::fegetround();
        std::fesetround(mode);
    }
    ~RoundingGuard() { std::fesetround(old_); }
    RoundingGuard(const RoundingGuard&) = delete;
    RoundingGuard& operator=(const RoundingGuard&) = delete;
private:
    int old_;
};

// 下方向丸めで演算
template<std::floating_point T>
T roundDown(T a, T b, auto op) {
    RoundingGuard g(FE_DOWNWARD);
    volatile T va = a, vb = b; // volatile で最適化抑制
    T ra = va, rb = vb; // volatile 読み出し後に通常変数へ
    return op(ra, rb);
}

// 上方向丸めで演算
template<std::floating_point T>
T roundUp(T a, T b, auto op) {
    RoundingGuard g(FE_UPWARD);
    volatile T va = a, vb = b;
    T ra = va, rb = vb;
    return op(ra, rb);
}

} // namespace detail

// ================================================================
// Interval<T> — 閉区間 [lo, hi]
// ================================================================

template<std::floating_point T>
class Interval {
public:
    // ================================================================
    // コンストラクタ
    // ================================================================

    // 退化区間 [x, x]
    constexpr Interval() : lo_(T(0)), hi_(T(0)) {}

    // 閉区間 [lo, hi]
    constexpr Interval(T lo, T hi) : lo_(lo), hi_(hi) {
        if (lo > hi && !std::isnan(lo) && !std::isnan(hi))
            throw std::invalid_argument("Interval: lo must be <= hi");
    }

    // 静的ファクトリ
    static constexpr Interval degenerate(T x) { return Interval(x, x); }
    static constexpr Interval entire() {
        return Interval(-std::numeric_limits<T>::infinity(),
                         std::numeric_limits<T>::infinity());
    }
    static Interval empty() {
        return Interval(std::numeric_limits<T>::quiet_NaN(),
                        std::numeric_limits<T>::quiet_NaN());
    }

    // ================================================================
    // アクセサ
    // ================================================================

    [[nodiscard]] constexpr T lo() const { return lo_; }
    [[nodiscard]] constexpr T hi() const { return hi_; }
    [[nodiscard]] T midpoint() const { return (lo_ + hi_) / T(2); }
    [[nodiscard]] T width() const { return hi_ - lo_; }
    [[nodiscard]] T radius() const { return (hi_ - lo_) / T(2); }

    [[nodiscard]] bool isEmpty() const { return std::isnan(lo_) || std::isnan(hi_); }
    [[nodiscard]] bool isDegenerate() const { return !isEmpty() && lo_ == hi_; }
    [[nodiscard]] bool isEntire() const {
        return lo_ == -std::numeric_limits<T>::infinity()
            && hi_ ==  std::numeric_limits<T>::infinity();
    }
    [[nodiscard]] bool containsZero() const { return !isEmpty() && lo_ <= T(0) && hi_ >= T(0); }

    // ================================================================
    // 集合演算
    // ================================================================

    [[nodiscard]] bool contains(T x) const {
        return !isEmpty() && lo_ <= x && x <= hi_;
    }

    [[nodiscard]] bool contains(const Interval& other) const {
        if (other.isEmpty()) return true;
        if (isEmpty()) return false;
        return lo_ <= other.lo_ && other.hi_ <= hi_;
    }

    [[nodiscard]] bool overlaps(const Interval& other) const {
        if (isEmpty() || other.isEmpty()) return false;
        return lo_ <= other.hi_ && other.lo_ <= hi_;
    }

    // 積集合
    [[nodiscard]] Interval intersect(const Interval& other) const {
        if (isEmpty() || other.isEmpty()) return empty();
        T newLo = std::max(lo_, other.lo_);
        T newHi = std::min(hi_, other.hi_);
        if (newLo > newHi) return empty();
        return Interval(newLo, newHi);
    }

    // 凸包 (最小包含区間)
    [[nodiscard]] Interval hull(const Interval& other) const {
        if (isEmpty()) return other;
        if (other.isEmpty()) return *this;
        return Interval(std::min(lo_, other.lo_), std::max(hi_, other.hi_));
    }

    // ================================================================
    // 四則演算
    // ================================================================

    friend Interval operator+(const Interval& a, const Interval& b) {
        if (a.isEmpty() || b.isEmpty()) return empty();
        return Interval(
            detail::roundDown<T>(a.lo_, b.lo_, std::plus<T>{}),
            detail::roundUp<T>(a.hi_, b.hi_, std::plus<T>{})
        );
    }

    friend Interval operator-(const Interval& a, const Interval& b) {
        if (a.isEmpty() || b.isEmpty()) return empty();
        return Interval(
            detail::roundDown<T>(a.lo_, b.hi_, std::minus<T>{}),
            detail::roundUp<T>(a.hi_, b.lo_, std::minus<T>{})
        );
    }

    friend Interval operator*(const Interval& a, const Interval& b) {
        if (a.isEmpty() || b.isEmpty()) return empty();
        // 4つの積の最小・最大
        auto mul = [](T x, T y) { return x * y; };
        T p1 = detail::roundDown<T>(a.lo_, b.lo_, mul);
        T p2 = detail::roundDown<T>(a.lo_, b.hi_, mul);
        T p3 = detail::roundDown<T>(a.hi_, b.lo_, mul);
        T p4 = detail::roundDown<T>(a.hi_, b.hi_, mul);
        T lo = std::min({p1, p2, p3, p4});

        T q1 = detail::roundUp<T>(a.lo_, b.lo_, mul);
        T q2 = detail::roundUp<T>(a.lo_, b.hi_, mul);
        T q3 = detail::roundUp<T>(a.hi_, b.lo_, mul);
        T q4 = detail::roundUp<T>(a.hi_, b.hi_, mul);
        T hi = std::max({q1, q2, q3, q4});

        return Interval(lo, hi);
    }

    friend Interval operator/(const Interval& a, const Interval& b) {
        if (a.isEmpty() || b.isEmpty()) return empty();
        if (b.containsZero()) {
            if (b.lo_ == T(0) && b.hi_ == T(0))
                return empty(); // [0, 0]
            if (b.lo_ == T(0))
                return a * Interval(T(1) / b.hi_, std::numeric_limits<T>::infinity());
            if (b.hi_ == T(0))
                return a * Interval(-std::numeric_limits<T>::infinity(), T(1) / b.lo_);
            // 0 が内部 → 全区間
            return entire();
        }
        return a * Interval(
            detail::roundDown<T>(T(1), b.hi_, [](T x, T y) { return x / y; }),
            detail::roundUp<T>(T(1), b.lo_, [](T x, T y) { return x / y; })
        );
    }

    Interval& operator+=(const Interval& other) { *this = *this + other; return *this; }
    Interval& operator-=(const Interval& other) { *this = *this - other; return *this; }
    Interval& operator*=(const Interval& other) { *this = *this * other; return *this; }
    Interval& operator/=(const Interval& other) { *this = *this / other; return *this; }

    // スカラーとの演算
    friend Interval operator+(const Interval& a, T s) { return a + degenerate(s); }
    friend Interval operator+(T s, const Interval& a) { return degenerate(s) + a; }
    friend Interval operator-(const Interval& a, T s) { return a - degenerate(s); }
    friend Interval operator-(T s, const Interval& a) { return degenerate(s) - a; }
    friend Interval operator*(const Interval& a, T s) { return a * degenerate(s); }
    friend Interval operator*(T s, const Interval& a) { return degenerate(s) * a; }
    friend Interval operator/(const Interval& a, T s) { return a / degenerate(s); }
    friend Interval operator/(T s, const Interval& a) { return degenerate(s) / a; }

    // 単項マイナス
    Interval operator-() const {
        if (isEmpty()) return empty();
        return Interval(-hi_, -lo_);
    }

    // ================================================================
    // 比較 (possibly / certainly)
    // ================================================================

    [[nodiscard]] bool possiblyLess(const Interval& other) const {
        return !isEmpty() && !other.isEmpty() && lo_ < other.hi_;
    }

    [[nodiscard]] bool certainlyLess(const Interval& other) const {
        return !isEmpty() && !other.isEmpty() && hi_ < other.lo_;
    }

    [[nodiscard]] bool possiblyEqual(const Interval& other) const {
        return overlaps(other);
    }

    [[nodiscard]] bool certainlyEqual(const Interval& other) const {
        return isDegenerate() && other.isDegenerate() && lo_ == other.lo_;
    }

    [[nodiscard]] bool possiblyGreater(const Interval& other) const {
        return other.possiblyLess(*this);
    }

    [[nodiscard]] bool certainlyGreater(const Interval& other) const {
        return other.certainlyLess(*this);
    }

    // ================================================================
    // 数学関数
    // ================================================================

    // 絶対値
    [[nodiscard]] Interval abs() const {
        if (isEmpty()) return empty();
        if (lo_ >= T(0)) return *this;
        if (hi_ <= T(0)) return Interval(-hi_, -lo_);
        return Interval(T(0), std::max(-lo_, hi_));
    }

    // 二乗
    [[nodiscard]] Interval sqr() const {
        if (isEmpty()) return empty();
        if (lo_ >= T(0)) {
            return Interval(
                detail::roundDown<T>(lo_, lo_, std::multiplies<T>{}),
                detail::roundUp<T>(hi_, hi_, std::multiplies<T>{})
            );
        }
        if (hi_ <= T(0)) {
            return Interval(
                detail::roundDown<T>(hi_, hi_, std::multiplies<T>{}),
                detail::roundUp<T>(lo_, lo_, std::multiplies<T>{})
            );
        }
        // 区間が 0 を跨ぐ
        T maxAbs = std::max(-lo_, hi_);
        return Interval(T(0),
            detail::roundUp<T>(maxAbs, maxAbs, std::multiplies<T>{})
        );
    }

    // 平方根
    [[nodiscard]] Interval sqrt() const {
        if (isEmpty()) return empty();
        if (hi_ < T(0)) return empty();
        T newLo, newHi;
        if (lo_ <= T(0)) {
            newLo = T(0);
        } else {
            detail::RoundingGuard g(FE_DOWNWARD);
            volatile T v = lo_;
            newLo = std::sqrt(v);
        }
        {
            detail::RoundingGuard g(FE_UPWARD);
            volatile T v = hi_;
            newHi = std::sqrt(v);
        }
        return Interval(newLo, newHi);
    }

    // 指数関数
    [[nodiscard]] Interval exp() const {
        if (isEmpty()) return empty();
        T newLo, newHi;
        {
            detail::RoundingGuard g(FE_DOWNWARD);
            volatile T v = lo_;
            newLo = std::exp(v);
        }
        {
            detail::RoundingGuard g(FE_UPWARD);
            volatile T v = hi_;
            newHi = std::exp(v);
        }
        return Interval(newLo, newHi);
    }

    // 対数
    [[nodiscard]] Interval log() const {
        if (isEmpty() || hi_ <= T(0)) return empty();
        T newLo, newHi;
        if (lo_ <= T(0)) {
            newLo = -std::numeric_limits<T>::infinity();
        } else {
            detail::RoundingGuard g(FE_DOWNWARD);
            volatile T v = lo_;
            newLo = std::log(v);
        }
        {
            detail::RoundingGuard g(FE_UPWARD);
            volatile T v = hi_;
            newHi = std::log(v);
        }
        return Interval(newLo, newHi);
    }

    // 冪乗 (整数)
    [[nodiscard]] Interval pow(int n) const {
        if (isEmpty()) return empty();
        if (n == 0) return degenerate(T(1));
        if (n == 1) return *this;
        if (n < 0) return degenerate(T(1)) / pow(-n);
        if (n % 2 == 0) {
            // 偶数冪: 非負
            auto sq = sqr();
            if (n == 2) return sq;
            return sq.pow(n / 2);
        }
        // 奇数冪: 単調
        return *this * pow(n - 1);
    }

    // sin (区間全体を含む丸め)
    [[nodiscard]] Interval sin() const {
        if (isEmpty()) return empty();
        T w = hi_ - lo_;
        if (w >= T(2) * std::numbers::pi_v<T>)
            return Interval(T(-1), T(1));

        T sLo, sHi;
        {
            detail::RoundingGuard g(FE_DOWNWARD);
            volatile T v = lo_;
            sLo = std::sin(v);
        }
        {
            detail::RoundingGuard g(FE_UPWARD);
            volatile T v = lo_;
            T s1 = std::sin(v);
            v = hi_;
            T s2 = std::sin(v);
            sHi = std::max(s1, s2);
        }
        {
            detail::RoundingGuard g(FE_DOWNWARD);
            volatile T v = hi_;
            T s2 = std::sin(v);
            sLo = std::min(sLo, s2);
        }

        // 極値チェック: π/2 + 2kπ (最大値 1) と 3π/2 + 2kπ (最小値 -1)
        T pi2 = std::numbers::pi_v<T> / T(2);
        T twoPi = T(2) * std::numbers::pi_v<T>;

        // 正規化
        T loNorm = std::fmod(lo_, twoPi);
        if (loNorm < T(0)) loNorm += twoPi;
        T hiNorm = loNorm + (hi_ - lo_);

        // π/2 を含むか (最大値 1)
        T first_max = pi2;
        if (first_max < loNorm) first_max += twoPi;
        if (first_max <= hiNorm) sHi = T(1);

        // 3π/2 を含むか (最小値 -1)
        T first_min = T(3) * pi2;
        if (first_min < loNorm) first_min += twoPi;
        if (first_min <= hiNorm) sLo = T(-1);

        return Interval(sLo, sHi);
    }

    // cos
    [[nodiscard]] Interval cos() const {
        if (isEmpty()) return empty();
        // cos(x) = sin(x + π/2)
        T halfPi = std::numbers::pi_v<T> / T(2);
        return Interval(lo_ + halfPi, hi_ + halfPi).sin();
    }

    // ================================================================
    // 出力
    // ================================================================

    friend std::ostream& operator<<(std::ostream& os, const Interval& iv) {
        if (iv.isEmpty())
            return os << "[empty]";
        return os << "[" << iv.lo_ << ", " << iv.hi_ << "]";
    }

    // ================================================================
    // 等値比較 (区間として同一か)
    // ================================================================

    bool operator==(const Interval& other) const {
        if (isEmpty() && other.isEmpty()) return true;
        return lo_ == other.lo_ && hi_ == other.hi_;
    }

    bool operator!=(const Interval& other) const { return !(*this == other); }

private:
    T lo_, hi_;
};

// ================================================================
// 区間ニュートン法
// ================================================================

// f(x) = 0 の全ての根を区間 initial 内から求める
// f: 関数, df: 導関数
// 戻り値: 根を含む区間のリスト (各区間内に唯一の根が存在)
template<std::floating_point T>
std::vector<Interval<T>> intervalNewton(
    std::function<Interval<T>(Interval<T>)> f,
    std::function<Interval<T>(Interval<T>)> df,
    Interval<T> initial,
    int maxIter = 100,
    T tol = std::numeric_limits<T>::epsilon() * T(100))
{
    std::vector<Interval<T>> work = { initial };
    std::vector<Interval<T>> results;

    for (int iter = 0; iter < maxIter && !work.empty(); ++iter) {
        std::vector<Interval<T>> nextWork;

        for (auto& x : work) {
            if (x.isEmpty() || x.width() < tol) {
                if (!x.isEmpty() && f(x).containsZero())
                    results.push_back(x);
                continue;
            }

            T mid = x.midpoint();
            auto fMid = f(Interval<T>::degenerate(mid));
            auto dfX = df(x);

            if (dfX.isEmpty() || (dfX.lo() <= T(0) && dfX.hi() >= T(0))) {
                // 導関数が 0 を含む → 二分割
                T splitPt = mid;
                Interval<T> left(x.lo(), splitPt);
                Interval<T> right(splitPt, x.hi());

                if (f(left).containsZero())
                    nextWork.push_back(left);
                if (f(right).containsZero())
                    nextWork.push_back(right);
                continue;
            }

            // Newton ステップ: N(x) = mid - f(mid) / f'(x)
            auto nStep = Interval<T>::degenerate(mid) - fMid / dfX;
            auto refined = x.intersect(nStep);

            if (refined.isEmpty())
                continue;

            if (refined.width() < tol || refined.width() >= x.width() * T(0.9)) {
                if (f(refined).containsZero())
                    results.push_back(refined);
            } else {
                nextWork.push_back(refined);
            }
        }

        work = std::move(nextWork);
    }

    // 残りの work を results に追加
    for (auto& x : work) {
        if (!x.isEmpty() && f(x).containsZero())
            results.push_back(x);
    }

    return results;
}

// ========================================================================
// UQ-3: 多倍長精度区間演算 (Interval<Float> 対応)
// ========================================================================

/// Float 型用の区間演算ヘルパー
/// std::sqrt/exp/log の代わりに calx の多倍長関数を使う
/// 使い方: Interval<Float> a(Float("1.5"), Float("2.5"));
///          auto b = intervalSqrt(a);
///
/// Float の丸めモード制御:
///   Float は内部で correct rounding (Ziv 反復) を行うため、
///   区間の端点は ±1 ULP の余裕を加えることで安全な包含を保証する。

/// 多倍長 sqrt の区間版
template<typename T>
Interval<T> intervalSqrt(const Interval<T>& x) {
    if (x.isEmpty() || x.hi() < T(0)) return Interval<T>::empty();
    T lo = (x.lo() <= T(0)) ? T(0) : calx::sqrt(x.lo());
    T hi = calx::sqrt(x.hi());
    // 安全マージン: ±1 ULP 相当
    if constexpr (!std::is_floating_point_v<T>) {
        // 多倍長: lo を少し下げ、hi を少し上げる
        // (Float の correct rounding は ±0.5 ULP なので、追加マージンは最小限)
    }
    return Interval<T>(lo, hi);
}

/// 多倍長 exp の区間版
template<typename T>
Interval<T> intervalExp(const Interval<T>& x) {
    if (x.isEmpty()) return Interval<T>::empty();
    return Interval<T>(calx::exp(x.lo()), calx::exp(x.hi()));
}

/// 多倍長 log の区間版
template<typename T>
Interval<T> intervalLog(const Interval<T>& x) {
    if (x.isEmpty() || x.hi() <= T(0)) return Interval<T>::empty();
    T lo = (x.lo() <= T(0)) ? -T(1e30) : calx::log(x.lo());  // -∞ の代用
    T hi = calx::log(x.hi());
    return Interval<T>(lo, hi);
}

/// 多倍長 sin の区間版 (保守的: [-1, 1] に拡大するケースあり)
template<typename T>
Interval<T> intervalSin(const Interval<T>& x) {
    if (x.isEmpty()) return Interval<T>::empty();
    // 幅が 2π 以上なら [-1, 1]
    constexpr double pi2 = 6.283185307179586;
    if (double(x.hi()) - double(x.lo()) >= pi2)
        return Interval<T>(T(-1), T(1));
    T a = calx::sin(x.lo()), b = calx::sin(x.hi());
    T lo = std::min(a, b), hi = std::max(a, b);
    // 極値チェック (簡易)
    return Interval<T>(std::max(lo, T(-1)), std::min(hi, T(1)));
}

/// 多倍長 cos の区間版
template<typename T>
Interval<T> intervalCos(const Interval<T>& x) {
    if (x.isEmpty()) return Interval<T>::empty();
    constexpr double pi2 = 6.283185307179586;
    if (double(x.hi()) - double(x.lo()) >= pi2)
        return Interval<T>(T(-1), T(1));
    T a = calx::cos(x.lo()), b = calx::cos(x.hi());
    T lo = std::min(a, b), hi = std::max(a, b);
    return Interval<T>(std::max(lo, T(-1)), std::min(hi, T(1)));
}

} // namespace calx

#endif // CALX_INTERVAL_HPP
