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

// FloatMath.hpp
// 多倍長浮動小数点数の数学関数ライブラリ
//
// このファイルでは、多倍長浮動小数点数（Float）のための
// 高精度な数学関数を提供します。各関数は、テイラー展開や
// 特殊なアルゴリズムを使用して実装されています。
//
// 特徴:
// - 解析的に必要な項数とガードビット数を計算
// - 高精度なテイラー級数計算
// - 効率的なアルゴリズム選択

#ifndef CALX_FLOAT_MATH_HPP
#define CALX_FLOAT_MATH_HPP

#include <math/core/mp/Float.hpp>

namespace calx {

/**
 * @brief 高精度加算のためのヘルパークラス
 * 
 * このクラスは、多くの項を加算する際の精度を維持するためにカハンの
 * 加算アルゴリズムを使用し、丸め誤差を補償します。
 */
class PrecisionSummer {
private:
    Float sum_;
    Float error_;
    int working_precision_;  // 作業精度
    int target_precision_;   // 目標精度

public:
    /**
     * @brief コンストラクタ
     * @param target_precision 目標精度（ビット単位）
     * @param guard_bits 追加のガードビット数（デフォルトは0）
     */
PrecisionSummer(int target_precision, int guard_bits = 0)
    : sum_(Float::zero(target_precision + guard_bits)),
      error_(Float::zero(target_precision + guard_bits)),
      working_precision_(target_precision + guard_bits),
      target_precision_(target_precision) {
      }

    /**
     * @brief 項を加算
     * @param term 加算する項
     */
    void add(const Float& term) {
        // カハン加算アルゴリズム
        Float y = term - error_;
        Float t = sum_ + y;
        error_ = (t - sum_) - y;
        sum_ = t;
    }

    /**
     * @brief 最終結果を取得
     * @return 目標精度に丸められた結果
     */
    Float result() const {
        Float final_result = sum_;
        final_result.setPrecision(target_precision_);
        return final_result;
    }
    
    /**
     * @brief 現在の合計を取得（作業精度のまま）
     * @return 現在の合計
     */
    const Float& current() const {
        return sum_;
    }
    
    /**
     * @brief 現在の合計値をリセット
     */
    void reset() {
        sum_ = Float::zero(working_precision_);
        error_ = Float::zero(working_precision_);
    }
};

/**
 * @brief 指数関数（e^x）を計算
 * 
 * テイラー展開を使用して e^x を計算します。
 * 必要な項数とガードビット数を解析的に計算します。
 * 
 * @param x 指数
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return e^x の値
 */
Float exp(const Float& x, int precision);

/**
 * @brief 自然対数（ln(x)）を計算
 * 
 * テイラー展開を使用して ln(x) を計算します。
 * 収束を高速化するために適切な変換を適用します。
 * 
 * @param x 入力値（正の数）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return ln(x) の値
 */
Float log(const Float& x, int precision);

/**
 * @brief 正の整数の自然対数（ln(n)）を計算
 *
 * 整数を小素数の積に分解し、キャッシュされた log(p) を
 * 利用して高速に計算します。分解できない場合は log(Float(n)) に
 * フォールバックします。
 *
 * @param n 正の整数
 * @param precision 精度（ビット単位）
 * @return ln(n) の値（n=0 の場合は -infinity）
 */
Float logUi(unsigned long long n, int precision);

/**
 * @brief 正弦関数（sin(x)）を計算
 * 
 * テイラー展開を使用して sin(x) を計算します。
 * 大きな引数の場合は周期性を利用して縮小します。
 * 
 * @param x 角度（ラジアン）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return sin(x) の値
 */
Float sin(const Float& x, int precision);

/**
 * @brief 余弦関数（cos(x)）を計算
 * 
 * テイラー展開を使用して cos(x) を計算します。
 * 大きな引数の場合は周期性を利用して縮小します。
 * 
 * @param x 角度（ラジアン）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return cos(x) の値
 */
Float cos(const Float& x, int precision);

/**
 * @brief 正接関数（tan(x)）を計算
 * 
 * sin(x)/cos(x) として計算されます。
 * 
 * @param x 角度（ラジアン）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return tan(x) の値
 */
Float tan(const Float& x, int precision);

/**
 * @brief 平方根（sqrt(x)）を計算
 * 
 * ニュートン法を使用して平方根を計算します。
 * 必要な反復回数を解析的に計算します。
 * 
 * @param x 入力値（非負）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return sqrt(x) の値
 */
Float sqrt(const Float& x, int precision);

/**
 * @brief 双曲線正弦関数（sinh(x)）を計算
 * 
 * @param x 入力値
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return sinh(x) の値
 */
Float sinh(const Float& x, int precision);

/**
 * @brief 双曲線余弦関数（cosh(x)）を計算
 * 
 * @param x 入力値
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return cosh(x) の値
 */
Float cosh(const Float& x, int precision);

/**
 * @brief 双曲線正接関数（tanh(x)）を計算
 * 
 * @param x 入力値
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return tanh(x) の値
 */
Float tanh(const Float& x, int precision);

/**
 * @brief 逆正弦関数（asin(x)）を計算
 * 
 * @param x 入力値（-1から1の範囲）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return asin(x) の値
 */
Float asin(const Float& x, int precision);

/**
 * @brief 逆余弦関数（acos(x)）を計算
 * 
 * @param x 入力値（-1から1の範囲）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return acos(x) の値
 */
Float acos(const Float& x, int precision);

/**
 * @brief 逆正接関数（atan(x)）を計算
 * 
 * @param x 入力値
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return atan(x) の値
 */
Float atan(const Float& x, int precision);

/**
 * @brief 2つの引数による逆正接関数（atan2(y, x)）を計算
 * 
 * @param y y座標
 * @param x x座標
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return atan2(y, x) の値
 */
Float atan2(const Float& y, const Float& x, int precision);

/**
 * @brief 累乗（x^y）を計算
 * 
 * @param x 底
 * @param y 指数
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return x^y の値
 */
Float pow(const Float& x, const Float& y, int precision);

/**
 * @brief 累乗（x^n）を計算（整数指数）
 * 
 * @param x 底
 * @param n 指数（整数）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return x^n の値
 */
Float pow(const Float& x, int n, int precision);

/**
 * @brief 逆双曲線正弦関数（asinh(x)）を計算
 *
 * asinh(x) = log(x + sqrt(x^2 + 1))
 *
 * @param x 入力値
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return asinh(x) の値
 */
Float asinh(const Float& x, int precision);

/**
 * @brief 逆双曲線余弦関数（acosh(x)）を計算
 *
 * acosh(x) = log(x + sqrt(x^2 - 1)), x >= 1
 *
 * @param x 入力値（1以上）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return acosh(x) の値
 */
Float acosh(const Float& x, int precision);

/**
 * @brief 逆双曲線正接関数（atanh(x)）を計算
 *
 * atanh(x) = 0.5 * log((1 + x) / (1 - x)), |x| < 1
 *
 * @param x 入力値（-1から1の範囲）
 * @param precision 精度（ビット単位、デフォルトはFloatのデフォルト精度）
 * @return atanh(x) の値
 */
Float atanh(const Float& x, int precision);

/**
 * @brief log(1 + x) を高精度に計算
 *
 * |x| が小さい場合に log(1+x) を直接計算するのではなく、
 * 桁落ちを回避した高精度な実装を使用します。
 *
 * @param x 入力値（-1より大きい値）
 * @param precision 精度（ビット単位）
 * @return log(1 + x) の値
 */
Float log1p(const Float& x, int precision);

/**
 * @brief 底2の対数を計算
 * @param x 入力値（正の数）
 * @param precision 精度（ビット単位）
 * @return log2(x) = log(x) / ln(2)
 */
Float log2(const Float& x, int precision);

/**
 * @brief 常用対数を計算
 * @param x 入力値（正の数）
 * @param precision 精度（ビット単位）
 * @return log10(x) = log(x) / ln(10)
 */
Float log10(const Float& x, int precision);

/**
 * @brief 2^x を計算
 * @param x 指数
 * @param precision 精度（ビット単位）
 * @return 2^x の値
 */
Float exp2(const Float& x, int precision);

/**
 * @brief 10^x を計算
 * @param x 指数
 * @param precision 精度（ビット単位）
 * @return 10^x の値
 */
Float exp10(const Float& x, int precision);

/**
 * @brief e^x - 1 を高精度に計算
 *
 * |x| が小さい場合の桁落ちを回避します。
 *
 * @param x 入力値
 * @param precision 精度（ビット単位）
 * @return e^x - 1 の値
 */
Float expm1(const Float& x, int precision);

/**
 * @brief 浮動小数点剰余を計算
 *
 * fmod(x, y) = x - trunc(x/y) * y
 *
 * @param x 被除数
 * @param y 除数
 * @return 剰余（x と同じ符号）
 */
Float fmod(const Float& x, const Float& y);
Float fmod(Float&& x, Float&& y);

/**
 * @brief IEEE 754 remainder を計算
 *
 * remainder(x, y) = x - round(x/y) * y
 *
 * @param x 被除数
 * @param y 除数
 * @return 剰余（|r| <= |y|/2）
 */
Float remainder(const Float& x, const Float& y);
Float remainder(Float&& x, Float&& y);

/**
 * @brief ユークリッドノルム（斜辺の長さ）を計算
 *
 * hypot(x, y) = sqrt(x^2 + y^2)（オーバーフロー/アンダーフロー回避）
 *
 * @param x 第1引数
 * @param y 第2引数
 * @param precision 精度（ビット単位）
 * @return sqrt(x^2 + y^2) の値
 */
Float hypot(const Float& x, const Float& y, int precision);

/**
 * @brief 立方根を計算
 * @param x 入力値
 * @param precision 精度（ビット単位）
 * @return x^(1/3) の値
 */
Float cbrt(const Float& x, int precision);

/**
 * @brief n乗根を計算
 * @param x 入力値
 * @param n 根の次数（正の整数）
 * @param precision 精度（ビット単位）
 * @return x^(1/n) の値
 */
Float nthRoot(const Float& x, int n, int precision);

/**
 * @brief 逆平方根を計算
 * @param x 入力値（正の数）
 * @param precision 精度（ビット単位）
 * @return 1/sqrt(x) の値
 */
Float recSqrt(const Float& x, int precision);

/**
 * @brief 融合乗加算: a * b + c
 * @param a 第1乗数
 * @param b 第2乗数
 * @param c 加数
 * @param precision 精度（ビット単位）
 * @return a * b + c の値（中間丸めなし）
 */
Float fma(const Float& a, const Float& b, const Float& c, int precision);

/**
 * @brief 融合乗減算: a * b - c
 * @param a 第1乗数
 * @param b 第2乗数
 * @param c 減数
 * @param precision 精度（ビット単位）
 * @return a * b - c の値（中間丸めなし）
 */
Float fms(const Float& a, const Float& b, const Float& c, int precision);

/**
 * @brief 二重積和: a * b + c * d (MPFR mpfr_fmma 相当)
 */
Float fmma(const Float& a, const Float& b, const Float& c, const Float& d, int precision);

/**
 * @brief 二重積差: a * b - c * d (MPFR mpfr_fmms 相当)
 */
Float fmms(const Float& a, const Float& b, const Float& c, const Float& d, int precision);

/**
 * @brief 専用二乗: x² (IntOps::square 利用で x*x より高速)
 */
Float sqr(const Float& x, int precision);
Float sqr(Float&& x, int precision);

/**
 * @brief sin と cos を同時に計算
 * @param x 角度（ラジアン）
 * @param[out] sin_result sin(x) の結果
 * @param[out] cos_result cos(x) の結果
 * @param precision 精度（ビット単位）
 */
void sinCos(const Float& x, Float& sin_result, Float& cos_result, int precision);

} // namespace calx

#endif // CALX_FLOAT_MATH_HPP