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

// Uncertain.hpp
// 誤差伝播つき不確かさ型 (P-28)
//
// Uncertain<T> は値と標準不確かさのペアを保持し、
// 算術演算・数学関数を通じて偏微分に基づく誤差伝播を行う。
//
// 使用例:
//   Uncertain<double> x(3.0, 0.1);
//   Uncertain<double> y(4.0, 0.2);
//   auto z = sqrt(x*x + y*y);  // z.value() ≈ 5.0, z.uncertainty() ≈ ...

#ifndef CALX_UNCERTAIN_HPP
#define CALX_UNCERTAIN_HPP

#include <cmath>
#include <ostream>
#include <type_traits>

namespace calx {

// =====================================================================
// Uncertain<T> — 値 ± 標準不確かさ
// =====================================================================

template<typename T>
class Uncertain {
    T value_;
    T uncertainty_;  // 常に非負

public:
    // -----------------------------------------------------------------
    // 構築
    // -----------------------------------------------------------------

    constexpr Uncertain() : value_(T(0)), uncertainty_(T(0)) {}

    constexpr explicit Uncertain(T val, T unc = T(0))
        : value_(val), uncertainty_(unc < T(0) ? -unc : unc) {}

    // -----------------------------------------------------------------
    // アクセサ
    // -----------------------------------------------------------------

    constexpr T value()       const { return value_; }
    constexpr T uncertainty() const { return uncertainty_; }

    // -----------------------------------------------------------------
    // 算術演算子 (Uncertain op Uncertain)
    // -----------------------------------------------------------------

    friend constexpr Uncertain operator+(const Uncertain& a, const Uncertain& b) {
        // δ(a+b) = sqrt(δa² + δb²)
        return Uncertain(a.value_ + b.value_,
                         std::sqrt(a.uncertainty_ * a.uncertainty_
                                 + b.uncertainty_ * b.uncertainty_));
    }

    friend constexpr Uncertain operator-(const Uncertain& a, const Uncertain& b) {
        // δ(a-b) = sqrt(δa² + δb²)
        return Uncertain(a.value_ - b.value_,
                         std::sqrt(a.uncertainty_ * a.uncertainty_
                                 + b.uncertainty_ * b.uncertainty_));
    }

    friend constexpr Uncertain operator*(const Uncertain& a, const Uncertain& b) {
        // δ(ab)/|ab| = sqrt((δa/a)² + (δb/b)²)
        T val = a.value_ * b.value_;
        T ra = (a.value_ != T(0)) ? (a.uncertainty_ / a.value_) : T(0);
        T rb = (b.value_ != T(0)) ? (b.uncertainty_ / b.value_) : T(0);
        T unc = std::abs(val) * std::sqrt(ra * ra + rb * rb);
        return Uncertain(val, unc);
    }

    friend constexpr Uncertain operator/(const Uncertain& a, const Uncertain& b) {
        // δ(a/b)/|a/b| = sqrt((δa/a)² + (δb/b)²)
        T val = a.value_ / b.value_;
        T ra = (a.value_ != T(0)) ? (a.uncertainty_ / a.value_) : T(0);
        T rb = (b.value_ != T(0)) ? (b.uncertainty_ / b.value_) : T(0);
        T unc = std::abs(val) * std::sqrt(ra * ra + rb * rb);
        return Uncertain(val, unc);
    }

    // -----------------------------------------------------------------
    // 算術演算子 (Uncertain op スカラー, スカラー op Uncertain)
    // -----------------------------------------------------------------

    friend constexpr Uncertain operator+(const Uncertain& a, const T& s) {
        return a + Uncertain(s);
    }
    friend constexpr Uncertain operator+(const T& s, const Uncertain& a) {
        return Uncertain(s) + a;
    }

    friend constexpr Uncertain operator-(const Uncertain& a, const T& s) {
        return a - Uncertain(s);
    }
    friend constexpr Uncertain operator-(const T& s, const Uncertain& a) {
        return Uncertain(s) - a;
    }

    friend constexpr Uncertain operator*(const Uncertain& a, const T& s) {
        return Uncertain(a.value_ * s, std::abs(s) * a.uncertainty_);
    }
    friend constexpr Uncertain operator*(const T& s, const Uncertain& a) {
        return Uncertain(s * a.value_, std::abs(s) * a.uncertainty_);
    }

    friend constexpr Uncertain operator/(const Uncertain& a, const T& s) {
        return Uncertain(a.value_ / s, a.uncertainty_ / std::abs(s));
    }
    friend constexpr Uncertain operator/(const T& s, const Uncertain& a) {
        return Uncertain(s) / a;
    }

    // -----------------------------------------------------------------
    // 単項演算子
    // -----------------------------------------------------------------

    constexpr Uncertain operator-() const {
        return Uncertain(-value_, uncertainty_);
    }

    constexpr Uncertain operator+() const {
        return *this;
    }

    // -----------------------------------------------------------------
    // 複合代入
    // -----------------------------------------------------------------

    Uncertain& operator+=(const Uncertain& rhs) { *this = *this + rhs; return *this; }
    Uncertain& operator-=(const Uncertain& rhs) { *this = *this - rhs; return *this; }
    Uncertain& operator*=(const Uncertain& rhs) { *this = *this * rhs; return *this; }
    Uncertain& operator/=(const Uncertain& rhs) { *this = *this / rhs; return *this; }

    // -----------------------------------------------------------------
    // 出力
    // -----------------------------------------------------------------

    friend std::ostream& operator<<(std::ostream& os, const Uncertain& u) {
        os << u.value_ << " \xC2\xB1 " << u.uncertainty_;  // UTF-8 ±
        return os;
    }
};

// =====================================================================
// 数学関数 — δf = |f'(a)| * δa
// =====================================================================

template<typename T>
Uncertain<T> sqrt(const Uncertain<T>& a) {
    T val = std::sqrt(a.value());
    // d/dx sqrt(x) = 1 / (2*sqrt(x))
    T unc = a.uncertainty() / (T(2) * val);
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> sin(const Uncertain<T>& a) {
    T val = std::sin(a.value());
    T unc = std::abs(std::cos(a.value())) * a.uncertainty();
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> cos(const Uncertain<T>& a) {
    T val = std::cos(a.value());
    T unc = std::abs(std::sin(a.value())) * a.uncertainty();
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> tan(const Uncertain<T>& a) {
    T val = std::tan(a.value());
    T c = std::cos(a.value());
    // d/dx tan(x) = 1 / cos²(x)
    T unc = a.uncertainty() / (c * c);
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> exp(const Uncertain<T>& a) {
    T val = std::exp(a.value());
    // d/dx exp(x) = exp(x)
    T unc = val * a.uncertainty();
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> log(const Uncertain<T>& a) {
    T val = std::log(a.value());
    // d/dx log(x) = 1/x
    T unc = a.uncertainty() / std::abs(a.value());
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> pow(const Uncertain<T>& base, const T& exponent) {
    T val = std::pow(base.value(), exponent);
    // d/dx x^n = n * x^(n-1)
    T unc = std::abs(exponent * std::pow(base.value(), exponent - T(1)))
          * base.uncertainty();
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> pow(const Uncertain<T>& base, const Uncertain<T>& exponent) {
    // f(a,b) = a^b, δf = sqrt((b*a^(b-1)*δa)² + (a^b*ln(a)*δb)²)
    T val = std::pow(base.value(), exponent.value());
    T da = exponent.value() * std::pow(base.value(), exponent.value() - T(1))
         * base.uncertainty();
    T db = val * std::log(base.value()) * exponent.uncertainty();
    T unc = std::sqrt(da * da + db * db);
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> abs(const Uncertain<T>& a) {
    return Uncertain<T>(std::abs(a.value()), a.uncertainty());
}

template<typename T>
Uncertain<T> asin(const Uncertain<T>& a) {
    T val = std::asin(a.value());
    // d/dx asin(x) = 1 / sqrt(1 - x²)
    T unc = a.uncertainty() / std::sqrt(T(1) - a.value() * a.value());
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> acos(const Uncertain<T>& a) {
    T val = std::acos(a.value());
    // d/dx acos(x) = -1 / sqrt(1 - x²)
    T unc = a.uncertainty() / std::sqrt(T(1) - a.value() * a.value());
    return Uncertain<T>(val, unc);
}

template<typename T>
Uncertain<T> atan(const Uncertain<T>& a) {
    T val = std::atan(a.value());
    // d/dx atan(x) = 1 / (1 + x²)
    T unc = a.uncertainty() / (T(1) + a.value() * a.value());
    return Uncertain<T>(val, unc);
}

} // namespace calx

#endif // CALX_UNCERTAIN_HPP
