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

// AffineArithmetic.hpp
//
// アフィン演算 (Affine Arithmetic, AA)
//
// 区間演算 (Interval Arithmetic) の過大評価を低減する手法。
// 各変数を x = x₀ + x₁ε₁ + x₂ε₂ + ... + xₙεₙ (εᵢ ∈ [-1,1]) の
// アフィン形式で表現し、変数間の依存性を追跡することで
// 区間幅の不必要な膨張を防ぐ。
//
// 例: 区間演算で x - x = [-2,2] - [-2,2] = [-4,4] となるが、
//     アフィン演算では x - x = 0 (依存性が相殺される)
//
// 主な API:
//   AffineForm(lo, hi)           — 区間からアフィン形式を構築
//   operator+, -, *, /           — 四則演算
//   toInterval()                 — 区間に変換
//   sqr(), sqrt(), exp(), log()  — 非線形関数 (Chebyshev 近似 + 誤差項)
//
// 参考文献:
//   - Stolfi & de Figueiredo, "Self-Validated Numerical Methods
//     and Applications" (1997)
//   - Comba & Stolfi, "Affine Arithmetic and its Applications
//     to Computer Graphics" (1993)

#ifndef CALX_CORE_AFFINE_ARITHMETIC_HPP
#define CALX_CORE_AFFINE_ARITHMETIC_HPP

#include <vector>
#include <cmath>
#include <algorithm>
#include <stdexcept>
#include <limits>
#include <numeric>
#include <cstddef>
#include <cstdint>

namespace calx {

/**
 * @brief アフィン形式 x = x₀ + Σ xᵢεᵢ
 *
 * x₀: 中心値 (central value)
 * xᵢ: ノイズ記号 εᵢ の係数 (partial deviation)
 * εᵢ ∈ [-1, 1]: ノイズ記号 (noise symbol)
 *
 * 区間に変換すると [x₀ - Σ|xᵢ|, x₀ + Σ|xᵢ|]
 */
class AffineForm {
public:
    /// 定数
    explicit AffineForm(double val = 0.0) : x0_(val) {}

    /// 区間 [lo, hi] からアフィン形式を構築 (新しいノイズ記号を割り当て)
    AffineForm(double lo, double hi) {
        if (lo > hi)
            throw std::invalid_argument("AffineForm: lo must be <= hi");
        x0_ = (lo + hi) * 0.5;
        double rad = (hi - lo) * 0.5;
        if (rad > 0.0) {
            std::size_t idx = nextSymbol();
            ensureSize(idx);
            noise_[idx] = rad;
        }
    }

    /// 中心値
    [[nodiscard]] double central() const { return x0_; }

    /// ノイズ係数
    [[nodiscard]] double noise(std::size_t i) const {
        return i < noise_.size() ? noise_[i] : 0.0;
    }

    /// ノイズ記号の数
    [[nodiscard]] std::size_t noiseCount() const { return noise_.size(); }

    /// 半径 (Σ|xᵢ|)
    [[nodiscard]] double radius() const {
        double r = 0.0;
        for (auto c : noise_) r += std::abs(c);
        return r;
    }

    /// 区間に変換
    [[nodiscard]] std::pair<double, double> toInterval() const {
        double r = radius();
        return {x0_ - r, x0_ + r};
    }

    // ════════════════════════════════════════════════════════════════
    // アフィン演算 (線形演算 — 誤差なしで正確)
    // ════════════════════════════════════════════════════════════════

    friend AffineForm operator+(const AffineForm& a, const AffineForm& b) {
        AffineForm r;
        r.x0_ = a.x0_ + b.x0_;
        std::size_t n = std::max(a.noise_.size(), b.noise_.size());
        r.noise_.resize(n, 0.0);
        for (std::size_t i = 0; i < a.noise_.size(); ++i) r.noise_[i] += a.noise_[i];
        for (std::size_t i = 0; i < b.noise_.size(); ++i) r.noise_[i] += b.noise_[i];
        return r;
    }

    friend AffineForm operator-(const AffineForm& a, const AffineForm& b) {
        AffineForm r;
        r.x0_ = a.x0_ - b.x0_;
        std::size_t n = std::max(a.noise_.size(), b.noise_.size());
        r.noise_.resize(n, 0.0);
        for (std::size_t i = 0; i < a.noise_.size(); ++i) r.noise_[i] += a.noise_[i];
        for (std::size_t i = 0; i < b.noise_.size(); ++i) r.noise_[i] -= b.noise_[i];
        return r;
    }

    friend AffineForm operator+(const AffineForm& a, double c) {
        AffineForm r = a;
        r.x0_ += c;
        return r;
    }

    friend AffineForm operator+(double c, const AffineForm& a) { return a + c; }

    friend AffineForm operator-(const AffineForm& a, double c) {
        AffineForm r = a;
        r.x0_ -= c;
        return r;
    }

    friend AffineForm operator-(double c, const AffineForm& a) {
        AffineForm r;
        r.x0_ = c - a.x0_;
        r.noise_.resize(a.noise_.size());
        for (std::size_t i = 0; i < a.noise_.size(); ++i) r.noise_[i] = -a.noise_[i];
        return r;
    }

    AffineForm operator-() const {
        AffineForm r;
        r.x0_ = -x0_;
        r.noise_.resize(noise_.size());
        for (std::size_t i = 0; i < noise_.size(); ++i) r.noise_[i] = -noise_[i];
        return r;
    }

    friend AffineForm operator*(const AffineForm& a, double c) {
        AffineForm r;
        r.x0_ = a.x0_ * c;
        r.noise_.resize(a.noise_.size());
        for (std::size_t i = 0; i < a.noise_.size(); ++i) r.noise_[i] = a.noise_[i] * c;
        return r;
    }

    friend AffineForm operator*(double c, const AffineForm& a) { return a * c; }

    // ════════════════════════════════════════════════════════════════
    // 準アフィン演算 (非線形 — 新しいノイズ記号で誤差を追跡)
    // ════════════════════════════════════════════════════════════════

    /// 乗算 (二次項は新しいノイズ記号に集約)
    friend AffineForm operator*(const AffineForm& a, const AffineForm& b) {
        // x*y = (x₀ + Σxᵢεᵢ)(y₀ + Σyⱼεⱼ)
        //     = x₀y₀ + Σ(x₀yᵢ + y₀xᵢ)εᵢ + Σᵢⱼ xᵢyⱼ εᵢεⱼ
        // 二次項の絶対値合計 ≤ rad(x) * rad(y) → 新しいノイズ記号
        AffineForm r;
        r.x0_ = a.x0_ * b.x0_;

        std::size_t n = std::max(a.noise_.size(), b.noise_.size());
        r.noise_.resize(n, 0.0);
        for (std::size_t i = 0; i < n; ++i) {
            double ai = i < a.noise_.size() ? a.noise_[i] : 0.0;
            double bi = i < b.noise_.size() ? b.noise_[i] : 0.0;
            r.noise_[i] = a.x0_ * bi + b.x0_ * ai;
        }

        // 二次項の上界
        double ra = a.radius();
        double rb = b.radius();
        double delta = ra * rb;
        if (delta > 0.0) {
            std::size_t idx = nextSymbol();
            r.ensureSize(idx);
            r.noise_[idx] = delta;
        }

        return r;
    }

    /// 除算 (b の区間が 0 を含まない場合のみ)
    friend AffineForm operator/(const AffineForm& a, const AffineForm& b) {
        // a / b = a * (1/b)
        return a * b.inverse();
    }

    /// 逆数 1/x (min-range Chebyshev 近似)
    [[nodiscard]] AffineForm inverse() const {
        auto [lo, hi] = toInterval();
        if (lo <= 0.0 && hi >= 0.0)
            throw std::domain_error("AffineForm::inverse: interval contains zero");
        return chebyshevApprox([](double x) { return 1.0 / x; },
                               [](double x) { return -1.0 / (x * x); },
                               lo, hi);
    }

    // ════════════════════════════════════════════════════════════════
    // 非線形関数 (Chebyshev 近似 + 新しいノイズ記号)
    // ════════════════════════════════════════════════════════════════

    /// x²
    [[nodiscard]] AffineForm sqr() const {
        // x² = x₀² + 2x₀·Σxᵢεᵢ + (Σxᵢεᵢ)²
        // 線形部: 2x₀xᵢ, 非線形部: ≤ rad²
        auto [lo, hi] = toInterval();
        return chebyshevApprox([](double x) { return x * x; },
                               [](double x) { return 2.0 * x; },
                               lo, hi);
    }

    /// √x (x > 0)
    [[nodiscard]] AffineForm sqrt() const {
        auto [lo, hi] = toInterval();
        if (lo < 0.0)
            throw std::domain_error("AffineForm::sqrt: interval contains negative values");
        lo = std::max(lo, 0.0);
        return chebyshevApprox([](double x) { return std::sqrt(x); },
                               [](double x) { return 0.5 / std::sqrt(x); },
                               std::max(lo, 1e-300), hi);
    }

    /// exp(x)
    [[nodiscard]] AffineForm exp() const {
        auto [lo, hi] = toInterval();
        return chebyshevApprox([](double x) { return std::exp(x); },
                               [](double x) { return std::exp(x); },
                               lo, hi);
    }

    /// log(x) (x > 0)
    [[nodiscard]] AffineForm log() const {
        auto [lo, hi] = toInterval();
        if (lo <= 0.0)
            throw std::domain_error("AffineForm::log: interval contains non-positive values");
        return chebyshevApprox([](double x) { return std::log(x); },
                               [](double x) { return 1.0 / x; },
                               lo, hi);
    }

    /// sin(x)
    [[nodiscard]] AffineForm sin() const {
        auto [lo, hi] = toInterval();
        return chebyshevApprox([](double x) { return std::sin(x); },
                               [](double x) { return std::cos(x); },
                               lo, hi);
    }

    /// cos(x)
    [[nodiscard]] AffineForm cos() const {
        auto [lo, hi] = toInterval();
        return chebyshevApprox([](double x) { return std::cos(x); },
                               [](double x) { return -std::sin(x); },
                               lo, hi);
    }

    /// ノイズ記号カウンタのリセット (テスト用)
    static void resetSymbolCounter() { symbolCounter() = 0; }

private:
    double x0_ = 0.0;                  // 中心値
    std::vector<double> noise_;         // ノイズ係数 (εᵢ の係数)

    void ensureSize(std::size_t idx) {
        if (idx >= noise_.size())
            noise_.resize(idx + 1, 0.0);
    }

    static std::size_t& symbolCounter() {
        static std::size_t counter = 0;
        return counter;
    }

    static std::size_t nextSymbol() {
        return symbolCounter()++;
    }

    /// Min-range Chebyshev 近似による非線形関数の適用
    ///
    /// f(x̂) ≈ α·x̂ + ζ + δ·εnew
    /// ここで α = (f(hi) - f(lo)) / (hi - lo) (secant slope)
    /// ζ, δ は残差の最小化から決定
    [[nodiscard]] AffineForm chebyshevApprox(
        double (*f)(double), double (*df)(double),
        double lo, double hi) const
    {
        if (hi - lo < 1e-15) {
            // 退化区間 — 点評価
            AffineForm r;
            r.x0_ = f(x0_);
            return r;
        }

        double fa = f(lo);
        double fb = f(hi);

        // secant slope
        double alpha = (fb - fa) / (hi - lo);

        // Chebyshev 残差の上下界を求める
        // r(x) = f(x) - α·x の最小値・最大値を [lo, hi] で求める
        // r'(x) = f'(x) - α = 0 の点が極値候補
        double rlo = fa - alpha * lo;
        double rhi = fb - alpha * hi;
        double rmin = std::min(rlo, rhi);
        double rmax = std::max(rlo, rhi);

        // 内部極値の探索 (f'(x) = alpha の点)
        // Newton 法で f'(x_c) = alpha を解く (最大 20 反復)
        double xc = (lo + hi) * 0.5;
        for (int iter = 0; iter < 20; ++iter) {
            double dfc = df(xc);
            // f''(x) の近似 (数値微分)
            double h = (hi - lo) * 1e-6;
            double dfp = df(std::min(xc + h, hi));
            double dfm = df(std::max(xc - h, lo));
            double d2f = (dfp - dfm) / (2.0 * h);
            if (std::abs(d2f) < 1e-30) break;
            double xn = xc - (dfc - alpha) / d2f;
            xn = std::max(lo, std::min(hi, xn));
            if (std::abs(xn - xc) < 1e-14 * (1.0 + std::abs(xc))) { xc = xn; break; }
            xc = xn;
        }
        if (xc > lo && xc < hi) {
            double rc = f(xc) - alpha * xc;
            rmin = std::min(rmin, rc);
            rmax = std::max(rmax, rc);
        }

        // ζ = (rmin + rmax) / 2, δ = (rmax - rmin) / 2
        double zeta = (rmin + rmax) * 0.5;
        double delta = (rmax - rmin) * 0.5;

        // 結果: α·x̂ + ζ + δ·εnew
        AffineForm r;
        r.x0_ = alpha * x0_ + zeta;
        r.noise_.resize(noise_.size());
        for (std::size_t i = 0; i < noise_.size(); ++i)
            r.noise_[i] = alpha * noise_[i];

        if (delta > 0.0) {
            std::size_t idx = nextSymbol();
            r.ensureSize(idx);
            r.noise_[idx] = delta;
        }

        return r;
    }
};

} // namespace calx

#endif // CALX_CORE_AFFINE_ARITHMETIC_HPP
