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

// FloatMath.hpp
// Mathematical function library for multi-precision floating-point numbers
//
// This file provides high-precision mathematical functions
// for multi-precision floating-point numbers (Float). Each function is implemented
// using Taylor expansions or specialized algorithms.
//
// Features:
// - Analytically computes the required number of terms and guard bits
// - High-precision Taylor series computation
// - Efficient algorithm selection

#ifndef SANGI_FLOAT_MATH_HPP
#define SANGI_FLOAT_MATH_HPP

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

namespace sangi {

/**
 * @brief Helper class for high-precision addition
 * 
 * This class uses Kahan's summation algorithm
 * to maintain precision when summing many terms by compensating for rounding error.
 */
class PrecisionSummer {
private:
    Float sum_;
    Float error_;
    int working_precision_;  // Working precision
    int target_precision_;   // Target precision

public:
    /**
     * @brief Constructor
     * @param target_precision Target precision (in bits)
     * @param guard_bits Additional guard bits (default 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 Add a term
     * @param term Term to add
     */
    void add(const Float& term) {
        // Kahan summation algorithm
        Float y = term - error_;
        Float t = sum_ + y;
        error_ = (t - sum_) - y;
        sum_ = t;
    }

    /**
     * @brief Get the final result
     * @return Result rounded to the target precision
     */
    Float result() const {
        Float final_result = sum_;
        final_result.setPrecision(target_precision_);
        return final_result;
    }
    
    /**
     * @brief Get the current sum (at working precision)
     * @return Current sum
     */
    const Float& current() const {
        return sum_;
    }
    
    /**
     * @brief Reset the current sum
     */
    void reset() {
        sum_ = Float::zero(working_precision_);
        error_ = Float::zero(working_precision_);
    }
};

/**
 * @brief Compute the exponential e^x
 * 
 * Compute e^x using a Taylor expansion.
 * Analytically compute the required number of terms and guard bits.
 * 
 * @param x Exponent
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of e^x
 */
Float exp(const Float& x, int precision);

/**
 * @brief Compute the natural logarithm ln(x)
 * 
 * Compute ln(x) using a Taylor expansion.
 * Apply appropriate transformations to accelerate convergence.
 * 
 * @param x Input value (positive)
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of ln(x)
 */
Float log(const Float& x, int precision);

/**
 * @brief Compute the natural logarithm of a positive integer (ln(n))
 *
 * Factor the integer into small primes and use cached log(p) for
 * fast computation. If factorization is unavailable, fall back to log(Float(n)).
 * (fallback)
 *
 * @param n Positive integer
 * @param precision Precision (in bits)
 * @return Value of ln(n) (returns -infinity when n=0)
 */
Float logUi(unsigned long long n, int precision);

/**
 * @brief Compute the sine sin(x)
 * 
 * Compute sin(x) using a Taylor expansion.
 * For large arguments, reduce using periodicity.
 * 
 * @param x Angle (in radians)
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of sin(x)
 */
Float sin(const Float& x, int precision);

/**
 * @brief Compute the cosine cos(x)
 * 
 * Compute cos(x) using a Taylor expansion.
 * For large arguments, reduce using periodicity.
 * 
 * @param x Angle (in radians)
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of cos(x)
 */
Float cos(const Float& x, int precision);

/**
 * @brief Compute the tangent tan(x)
 * 
 * Computed as sin(x)/cos(x).
 * 
 * @param x Angle (in radians)
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of tan(x)
 */
Float tan(const Float& x, int precision);

/**
 * @brief Compute the square root sqrt(x)
 * 
 * Compute the square root using Newton's method.
 * Analytically determine the number of iterations required.
 * 
 * @param x Input value (non-negative)
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of sqrt(x)
 */
Float sqrt(const Float& x, int precision);

/**
 * @brief Compute the hyperbolic sine sinh(x)
 * 
 * @param x Input value
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of sinh(x)
 */
Float sinh(const Float& x, int precision);

/**
 * @brief Compute the hyperbolic cosine cosh(x)
 * 
 * @param x Input value
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of cosh(x)
 */
Float cosh(const Float& x, int precision);

/**
 * @brief Compute the hyperbolic tangent tanh(x)
 * 
 * @param x Input value
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of tanh(x)
 */
Float tanh(const Float& x, int precision);

/**
 * @brief Compute the inverse sine asin(x)
 * 
 * @param x Input value (in [-1, 1])
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of asin(x)
 */
Float asin(const Float& x, int precision);

/**
 * @brief Compute the inverse cosine acos(x)
 * 
 * @param x Input value (in [-1, 1])
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of acos(x)
 */
Float acos(const Float& x, int precision);

/**
 * @brief Compute the inverse tangent atan(x)
 * 
 * @param x Input value
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of atan(x)
 */
Float atan(const Float& x, int precision);

/**
 * @brief Compute the two-argument inverse tangent atan2(y, x)
 * 
 * @param y y coordinate
 * @param x x coordinate
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of atan2(y, x)
 */
Float atan2(const Float& y, const Float& x, int precision);

/**
 * @brief Compute the power x^y
 * 
 * @param x Base
 * @param y Exponent
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of x^y
 */
Float pow(const Float& x, const Float& y, int precision);

/**
 * @brief Compute the power x^n (integer exponent)
 * 
 * @param x Base
 * @param n Exponent (integer)
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of x^n
 */
Float pow(const Float& x, int n, int precision);

/**
 * @brief Compute the inverse hyperbolic sine asinh(x)
 *
 * asinh(x) = log(x + sqrt(x^2 + 1))
 *
 * @param x Input value
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of asinh(x)
 */
Float asinh(const Float& x, int precision);

/**
 * @brief Compute the inverse hyperbolic cosine acosh(x)
 *
 * acosh(x) = log(x + sqrt(x^2 - 1)), x >= 1
 *
 * @param x Input value (>= 1)
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of acosh(x)
 */
Float acosh(const Float& x, int precision);

/**
 * @brief Compute the inverse hyperbolic tangent atanh(x)
 *
 * atanh(x) = 0.5 * log((1 + x) / (1 - x)), |x| < 1
 *
 * @param x Input value (in [-1, 1])
 * @param precision Precision (in bits; defaults to the Float default precision)
 * @return Value of atanh(x)
 */
Float atanh(const Float& x, int precision);

/**
 * @brief Compute log(1 + x) with high precision
 *
 * Instead of computing log(1+x) directly when |x| is small,
 * use a high-precision implementation that avoids cancellation.
 *
 * @param x Input value (> -1)
 * @param precision Precision (in bits)
 * @return Value of log(1 + x)
 */
Float log1p(const Float& x, int precision);

/**
 * @brief Compute the base-2 logarithm
 * @param x Input value (positive)
 * @param precision Precision (in bits)
 * @return log2(x) = log(x) / ln(2)
 */
Float log2(const Float& x, int precision);

/**
 * @brief Compute the common (base-10) logarithm
 * @param x Input value (positive)
 * @param precision Precision (in bits)
 * @return log10(x) = log(x) / ln(10)
 */
Float log10(const Float& x, int precision);

/**
 * @brief Compute 2^x
 * @param x Exponent
 * @param precision Precision (in bits)
 * @return Value of 2^x
 */
Float exp2(const Float& x, int precision);

/**
 * @brief Compute 10^x
 * @param x Exponent
 * @param precision Precision (in bits)
 * @return Value of 10^x
 */
Float exp10(const Float& x, int precision);

/**
 * @brief Compute e^x - 1 with high precision
 *
 * Avoids cancellation when |x| is small.
 *
 * @param x Input value
 * @param precision Precision (in bits)
 * @return Value of e^x - 1
 */
Float expm1(const Float& x, int precision);

/**
 * @brief Stable computation of (1+x)^n (equivalent to MPFR mpfr_compound_si)
 *
 * Common in financial (compound interest) and numerical analysis.
 * For small x, maintain precision via log1p / expm1.
 *
 * @param x Small quantity added to the base (x > -1)
 * @param n Integer exponent
 * @param precision Computation precision
 * @return (1+x)^n
 */
Float compound(const Float& x, long n, int precision);

/**
 * @brief Sine integral Si(x) = integral from 0 to x of sin(t)/t dt
 * @note Multi-precision Float version (Taylor expansion, for moderate |x|)
 */
Float sinIntegral(const Float& x, int precision);

/**
 * @brief Cosine integral Ci(x) = gamma + ln(x) + integral from 0 to x of (cos(t)-1)/t dt (x>0)
 */
Float cosIntegral(const Float& x, int precision);

/**
 * @brief Fresnel sine integral S(x) = integral from 0 to x of sin(pi*t^2/2) dt
 */
Float fresnelS(const Float& x, int precision);

/**
 * @brief Fresnel cosine integral C(x) = integral from 0 to x of cos(pi*t^2/2) dt
 */
Float fresnelC(const Float& x, int precision);

/**
 * @brief Lambert W function (principal branch W_0): W(x)*exp(W(x)) = x, x >= -1/e
 * @note Halley iteration; initial guess from the double version
 */
Float lambertW(const Float& x, int precision);

/**
 * @brief Compute the floating-point remainder
 *
 * fmod(x, y) = x - trunc(x/y) * y
 *
 * @param x Dividend
 * @param y Divisor
 * @return Remainder (same sign as x)
 */
Float fmod(const Float& x, const Float& y);
Float fmod(Float&& x, Float&& y);

/**
 * @brief Compute the IEEE 754 remainder
 *
 * remainder(x, y) = x - round(x/y) * y
 *
 * @param x Dividend
 * @param y Divisor
 * @return Remainder (|r| <= |y|/2)
 */
Float remainder(const Float& x, const Float& y);
Float remainder(Float&& x, Float&& y);

/**
 * @brief Compute the Euclidean norm (hypotenuse length)
 *
 * hypot(x, y) = sqrt(x^2 + y^2) (avoids overflow/underflow)
 *
 * @param x First argument
 * @param y Second argument
 * @param precision Precision (in bits)
 * @return Value of sqrt(x^2 + y^2)
 */
Float hypot(const Float& x, const Float& y, int precision);

/**
 * @brief Compute the cube root
 * @param x Input value
 * @param precision Precision (in bits)
 * @return Value of x^(1/3)
 */
Float cbrt(const Float& x, int precision);

/**
 * @brief Compute the n-th root
 * @param x Input value
 * @param n Degree of the root (positive integer)
 * @param precision Precision (in bits)
 * @return Value of x^(1/n)
 */
Float nthRoot(const Float& x, int n, int precision);

/**
 * @brief Compute the inverse square root
 * @param x Input value (positive)
 * @param precision Precision (in bits)
 * @return Value of 1/sqrt(x)
 */
Float recSqrt(const Float& x, int precision);

/**
 * @brief Fused multiply-add: a * b + c
 * @param a First multiplicand
 * @param b Second multiplicand
 * @param c Addend
 * @param precision Precision (in bits)
 * @return Value of a * b + c (no intermediate rounding)
 */
Float fma(const Float& a, const Float& b, const Float& c, int precision);

/**
 * @brief Fused multiply-subtract: a * b - c
 * @param a First multiplicand
 * @param b Second multiplicand
 * @param c Subtrahend
 * @param precision Precision (in bits)
 * @return Value of a * b - c (no intermediate rounding)
 */
Float fms(const Float& a, const Float& b, const Float& c, int precision);

/**
 * @brief Sum of two products: a * b + c * d (equivalent to MPFR mpfr_fmma)
 */
Float fmma(const Float& a, const Float& b, const Float& c, const Float& d, int precision);

/**
 * @brief Difference of two products: a * b - c * d (equivalent to MPFR mpfr_fmms)
 */
Float fmms(const Float& a, const Float& b, const Float& c, const Float& d, int precision);

/**
 * @brief Dedicated square: x^2 (faster than x*x via IntOps::square)
 */
Float sqr(const Float& x, int precision);
Float sqr(Float&& x, int precision);

/**
 * @brief Compute sin and cos together
 * @param x Angle (in radians)
 * @param[out] sin_result Result of sin(x)
 * @param[out] cos_result Result of cos(x)
 * @param precision Precision (in bits)
 */
void sinCos(const Float& x, Float& sin_result, Float& cos_result, int precision);

} // namespace sangi

#endif // SANGI_FLOAT_MATH_HPP