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

// BallArithmetic.hpp
//
// Ball 演算 (Ball Arithmetic / Arb 相当)
//
// 実数を [mid ± rad] の球で表現し、誤差を自動追跡する。
// 区間演算より精密: midpoint + radius で丸め誤差を累積。
//
// 主な API:
//   RealBall(mid, rad)           — 実数球
//   ComplexBall(re, im)          — 複素数球
//   operator+, -, *, /           — 四則演算 (誤差伝播)
//   sqrt, exp, log, sin, cos     — 初等関数
//   contains(ball, x)            — 包含判定
//   overlaps(a, b)               — 重なり判定
//   isFinite, isExact            — 性質判定

#ifndef CALX_CORE_BALL_ARITHMETIC_HPP
#define CALX_CORE_BALL_ARITHMETIC_HPP

#include <cstdint>
#include <cmath>
#include <stdexcept>
#include <algorithm>
#include <limits>
#include <numbers>

namespace calx {

/**
 * @brief 実数球 [mid ± rad]
 *
 * mid: 中心値, rad: 半径 (≥ 0)
 * 真の値 x は mid - rad ≤ x ≤ mid + rad に含まれる。
 */
class RealBall {
public:
    constexpr RealBall() : mid_(0.0), rad_(0.0) {}
    constexpr explicit RealBall(double mid) : mid_(mid), rad_(0.0) {}
    constexpr RealBall(double mid, double rad)
        : mid_(mid), rad_(rad >= 0 ? rad : -rad) {}

    constexpr double mid() const { return mid_; }
    constexpr double rad() const { return rad_; }

    /// 下限
    constexpr double lower() const { return mid_ - rad_; }
    /// 上限
    constexpr double upper() const { return mid_ + rad_; }

    /// 厳密値かどうか (rad == 0)
    constexpr bool isExact() const { return rad_ == 0.0; }

    /// 有限かどうか
    bool isFinite() const {
        return std::isfinite(mid_) && std::isfinite(rad_);
    }

    /// ゼロを含むか
    constexpr bool containsZero() const {
        return mid_ - rad_ <= 0.0 && 0.0 <= mid_ + rad_;
    }

    /// 正と確定できるか
    constexpr bool isPositive() const { return mid_ - rad_ > 0.0; }
    /// 負と確定できるか
    constexpr bool isNegative() const { return mid_ + rad_ < 0.0; }

    // 四則演算
    friend RealBall operator+(const RealBall& a, const RealBall& b) {
        double mid = a.mid_ + b.mid_;
        double rad = a.rad_ + b.rad_ + eps(mid);
        return {mid, rad};
    }

    friend RealBall operator-(const RealBall& a, const RealBall& b) {
        double mid = a.mid_ - b.mid_;
        double rad = a.rad_ + b.rad_ + eps(mid);
        return {mid, rad};
    }

    RealBall operator-() const { return {-mid_, rad_}; }

    friend RealBall operator*(const RealBall& a, const RealBall& b) {
        double mid = a.mid_ * b.mid_;
        // |a·b - mid_a·mid_b| ≤ |mid_a|·rad_b + |mid_b|·rad_a + rad_a·rad_b
        double rad = std::abs(a.mid_) * b.rad_ +
                     std::abs(b.mid_) * a.rad_ +
                     a.rad_ * b.rad_ + eps(mid);
        return {mid, rad};
    }

    friend RealBall operator/(const RealBall& a, const RealBall& b) {
        if (b.containsZero())
            throw std::domain_error("RealBall: division by interval containing zero");
        double mid = a.mid_ / b.mid_;
        // |a/b - mid_a/mid_b| の上界
        double bLower = std::abs(b.mid_) - b.rad_;
        double rad = (std::abs(a.mid_) * b.rad_ +
                      std::abs(b.mid_) * a.rad_ +
                      a.rad_ * b.rad_) / (bLower * std::abs(b.mid_)) + eps(mid);
        return {mid, rad};
    }

    // スカラー演算
    friend RealBall operator+(const RealBall& a, double b) {
        return a + RealBall(b);
    }
    friend RealBall operator+(double a, const RealBall& b) {
        return RealBall(a) + b;
    }
    friend RealBall operator*(const RealBall& a, double b) {
        return {a.mid_ * b, a.rad_ * std::abs(b) + eps(a.mid_ * b)};
    }
    friend RealBall operator*(double a, const RealBall& b) { return b * a; }

private:
    double mid_;
    double rad_;

    /// 浮動小数点丸め誤差の上界
    static double eps(double x) {
        return std::abs(x) * std::numeric_limits<double>::epsilon();
    }
};

/// 包含判定: ball が値 x を含むか
[[nodiscard]]
constexpr bool contains(const RealBall& ball, double x) {
    return ball.lower() <= x && x <= ball.upper();
}

/// 重なり判定
[[nodiscard]]
constexpr bool overlaps(const RealBall& a, const RealBall& b) {
    return a.lower() <= b.upper() && b.lower() <= a.upper();
}

// =====================================================================
// 初等関数
// =====================================================================

/// 平方根
[[nodiscard]]
inline RealBall sqrt(const RealBall& x) {
    if (x.isNegative())
        throw std::domain_error("RealBall::sqrt: negative argument");
    double mid = std::sqrt(std::max(0.0, x.mid()));
    // |√(m+r) - √m| ≤ r / (2√m) for m > 0
    double rad;
    if (mid > 0) {
        rad = x.rad() / (2.0 * mid) + RealBall(0).rad();
    } else {
        rad = std::sqrt(x.rad());
    }
    rad += std::abs(mid) * std::numeric_limits<double>::epsilon();
    return {mid, rad};
}

/// 指数関数
[[nodiscard]]
inline RealBall exp(const RealBall& x) {
    double mid = std::exp(x.mid());
    // |exp(m+r) - exp(m)| ≤ exp(m) · (exp(r) - 1) ≤ exp(m) · r · exp(r)
    double expRad = std::exp(x.rad());
    double rad = mid * (expRad - 1.0) +
                 std::abs(mid) * std::numeric_limits<double>::epsilon();
    return {mid, rad};
}

/// 自然対数
[[nodiscard]]
inline RealBall log(const RealBall& x) {
    if (!x.isPositive())
        throw std::domain_error("RealBall::log: argument must be positive");
    double mid = std::log(x.mid());
    // |log(m+r) - log(m)| ≤ r / (m - r) for m > r
    double denom = x.mid() - x.rad();
    double rad = x.rad() / denom +
                 std::abs(mid) * std::numeric_limits<double>::epsilon();
    return {mid, rad};
}

/// 正弦
[[nodiscard]]
inline RealBall sin(const RealBall& x) {
    double mid = std::sin(x.mid());
    // |sin(m+r) - sin(m)| ≤ r (Lipschitz 定数 = 1)
    double rad = x.rad() + std::abs(mid) * std::numeric_limits<double>::epsilon();
    return {mid, rad};
}

/// 余弦
[[nodiscard]]
inline RealBall cos(const RealBall& x) {
    double mid = std::cos(x.mid());
    double rad = x.rad() + std::abs(mid) * std::numeric_limits<double>::epsilon();
    return {mid, rad};
}

/// 正接
[[nodiscard]]
inline RealBall tan(const RealBall& x) {
    return sin(x) / cos(x);
}

/// 絶対値
[[nodiscard]]
inline RealBall abs(const RealBall& x) {
    if (x.isPositive()) return x;
    if (x.isNegative()) return -x;
    // ゼロを含む場合
    double mid = (x.upper() + std::max(0.0, x.lower())) / 2.0;
    double rad = x.upper() / 2.0;
    return {std::abs(x.mid()), x.rad()};
}

/// べき乗 (整数指数)
[[nodiscard]]
inline RealBall pow(const RealBall& x, int n) {
    if (n == 0) return RealBall(1.0);
    if (n == 1) return x;
    if (n < 0) return RealBall(1.0) / pow(x, -n);
    // 二乗反復
    if (n % 2 == 0) {
        auto half = pow(x, n / 2);
        return half * half;
    }
    return x * pow(x, n - 1);
}

/// π の ball 表現
[[nodiscard]]
inline RealBall ballPi() {
    return {std::numbers::pi, std::numeric_limits<double>::epsilon() * 4};
}

/// e の ball 表現
[[nodiscard]]
inline RealBall ballE() {
    return {std::numbers::e, std::numeric_limits<double>::epsilon() * 3};
}

// =====================================================================
// 複素数球
// =====================================================================

/**
 * @brief 複素数球 (re ± rad_re) + (im ± rad_im)i
 */
class ComplexBall {
public:
    ComplexBall() = default;
    ComplexBall(const RealBall& re, const RealBall& im = RealBall())
        : re_(re), im_(im) {}

    const RealBall& real() const { return re_; }
    const RealBall& imag() const { return im_; }

    bool isFinite() const { return re_.isFinite() && im_.isFinite(); }

    friend ComplexBall operator+(const ComplexBall& a, const ComplexBall& b) {
        return {a.re_ + b.re_, a.im_ + b.im_};
    }

    friend ComplexBall operator-(const ComplexBall& a, const ComplexBall& b) {
        return {a.re_ - b.re_, a.im_ - b.im_};
    }

    friend ComplexBall operator*(const ComplexBall& a, const ComplexBall& b) {
        return {a.re_ * b.re_ - a.im_ * b.im_,
                a.re_ * b.im_ + a.im_ * b.re_};
    }

    ComplexBall conj() const { return {re_, -im_}; }

    /// |z|² の ball
    RealBall normSq() const { return re_ * re_ + im_ * im_; }

private:
    RealBall re_, im_;
};

} // namespace calx

#endif // CALX_CORE_BALL_ARITHMETIC_HPP
