/*********************************************************************************************************
 *  ------------------------------------------------------------------------------------------------------
 *  file description
 *  ------------------------------------------------------------------------------------------------------
 *         \file  flmath.c
 *         \unit  flmath
 *        \brief  This is a simple large float number math calculate module for C language
 *       \author  Lamdonn
 *      \version  v1.0.0
 *      \license  GPL-2.0
 *    \copyright  Copyright (C) 2023 Lamdonn.
 ********************************************************************************************************/
#include "flmath.h"

// Number of terms in the Taylor series for trigonometric functions
static const int trig_terms = 15 * sizeof(floatl) / 8; 
// Number of terms in the Taylor series for the exponential function
static const int exp_terms = 15 * sizeof(floatl) / 8; 
// Number of terms in the Taylor series for the natural logarithm function
static const int ln_terms = 15 * sizeof(floatl) / 8; 

/**
 * \brief Calculates the factorial of an integer.
 * \param n: The integer for which the factorial is to be calculated.
 * \return A 'floatl' number representing the factorial of 'n'.
 */
static floatl floatl_int_fact(int n) 
{
    floatl result = FLOATL_CONST_1;
    for (int i = 1; i <= n; i++) 
    {
        result = floatl_mul(result, floatl(i));
    }
    return result;
}

/**
 * \brief Calculates the power of a 'floatl' number with an integer exponent.
 * \param base: The base 'floatl' number.
 * \param exponent: The integer exponent.
 * \return A 'floatl' number representing 'base' raised to the power of 'exponent'.
 */
static floatl floatl_int_power(floatl base, int exponent) 
{
    floatl result = FLOATL_CONST_1;
    for (int i = 0; i < exponent; i++) 
    {
        result = floatl_mul(result, base);
    }
    return result;
}

/**
 * \brief Calculates the sine of a 'floatl' number using the Taylor series expansion.
 * \param x: The 'floatl' number for which the sine is to be calculated.
 * \return A 'floatl' number representing the sine of 'x'.
 */
floatl floatl_sin(floatl x) 
{
    floatl result = FLOATL_CONST_0;
    // Normalize x to the range [-π, π]
    while (floatl_gt(x, FLOATL_PI)) 
    {
        x = floatl_sub(x, FLOATL_PI);
        x = floatl_sub(x, FLOATL_PI);
    }
    while (floatl_lt(x, floatl_neg(FLOATL_PI))) 
    {
        x = floatl_add(x, FLOATL_PI);
        x = floatl_add(x, FLOATL_PI);
    }
    for (int n = 0; n < trig_terms; n++) 
    {
        floatl term = floatl_int_power(floatl_neg(FLOATL_CONST_1), n);
        term = floatl_mul(term, floatl_int_power(x, 2 * n + 1));
        term = floatl_div(term, floatl_int_fact(2 * n + 1));
        result = floatl_add(result, term);
    }
    return result;
}

/**
 * \brief Calculates the cosine of a 'floatl' number using the sine function.
 * \param x: The 'floatl' number for which the cosine is to be calculated.
 * \return A 'floatl' number representing the cosine of 'x'.
 */
floatl floatl_cos(floatl x) 
{
    x = floatl_add(x, floatl_mul(floatl(0.5), FLOATL_PI));
    return floatl_sin(x);
}

/**
 * \brief Calculates the secant of a 'floatl' number.
 * \param x: The 'floatl' number for which the secant is to be calculated.
 * \return A 'floatl' number representing the secant of 'x'.
 */
floatl floatl_sec(floatl x)
{
    return floatl_div(FLOATL_CONST_1, floatl_cos(x));
}

/**
 * \brief Calculates the cosecant of a 'floatl' number.
 * \param x: The 'floatl' number for which the cosecant is to be calculated.
 * \return A 'floatl' number representing the cosecant of 'x'.
 */
floatl floatl_csc(floatl x)
{
    return floatl_div(FLOATL_CONST_1, floatl_sin(x));
}

/**
 * \brief Calculates the tangent of a 'floatl' number.
 * \param x: The 'floatl' number for which the tangent is to be calculated.
 * \return A 'floatl' number representing the tangent of 'x'. Returns infinity or negative infinity
 *         if the cosine of 'x' is zero.
 */
floatl floatl_tan(floatl x)
{
    floatl s = floatl_sin(x);
    floatl c = floatl_cos(x);
    if (floatl_eq(c, FLOATL_CONST_0)) 
    {
        return (floatl_gt(s, FLOATL_CONST_0)) ? FLOATL_INF : floatl_neg(FLOATL_INF);
    }
    return floatl_div(s, c);
}

/**
 * \brief Calculates the cotangent of a 'floatl' number.
 * \param x: The 'floatl' number for which the cotangent is to be calculated.
 * \return A 'floatl' number representing the cotangent of 'x'.
 */
floatl floatl_cot(floatl x)
{
    x = floatl_sub(floatl_mul(floatl(0.5), FLOATL_PI), x);
    return floatl_tan(x);
}

/**
 * \brief Calculates the exponential function of a 'floatl' number using the Taylor series expansion.
 * \param x: The 'floatl' number for which the exponential is to be calculated.
 * \return A 'floatl' number representing e raised to the power of 'x'.
 */
floatl floatl_exp(floatl x) 
{
    floatl result = FLOATL_CONST_0;
    for (int n = 0; n < exp_terms; n++) 
    {
        floatl term = floatl_int_power(x, n);
        term = floatl_div(term, floatl_int_fact(n));
        result = floatl_add(result, term);
    }
    return result;
}

/**
 * \brief Calculates the natural logarithm of a 'floatl' number using the Taylor series expansion.
 * \param x: The 'floatl' number for which the natural logarithm is to be calculated.
 * \return A 'floatl' number representing the natural logarithm of 'x'. Returns -1 if 'x' is less than or equal to 0.
 */
floatl floatl_ln(floatl x)
{
    if (floatl_le(x, FLOATL_CONST_0)) return floatl_neg(FLOATL_CONST_1);
    floatl result = FLOATL_CONST_0;
    // Transform x to the range (0, 2]
    floatl y = floatl_sub(x, FLOATL_CONST_1);
    y = floatl_div(y, floatl_add(x, FLOATL_CONST_1));
    for (int n = 0; n < ln_terms; n++) 
    {
        floatl term = floatl_int_power(y, 2 * n + 1);
        term = floatl_div(term, floatl(2 * n + 1));
        result = floatl_add(result, term);
    }
    result = floatl_mul(result, floatl(2.0));
    return result;
}

/**
 * \brief Calculates the logarithm of a 'floatl' number with a given base.
 * \param x: The base 'floatl' number.
 * \param y: The 'floatl' number for which the logarithm is to be calculated.
 * \return A 'floatl' number representing the logarithm of 'y' with base 'x'. Returns negative infinity
 *         if 'x' is zero or 'y' is less than or equal to 0.
 */
floatl floatl_log(floatl x, floatl y)
{
    if ((floatl_eq(x, FLOATL_CONST_0)) || (floatl_le(y, FLOATL_CONST_0))) return floatl_neg(FLOATL_INF);
    return floatl_div(floatl_ln(y), floatl_ln(x));
}

/**
 * \brief Calculates the logarithm of a 'floatl' number with base 2.
 * \param x: The 'floatl' number for which the logarithm is to be calculated.
 * \return A 'floatl' number representing the logarithm of 'x' with base 2. Returns negative infinity
 *         if 'x' is zero.
 */
floatl floatl_log2(floatl x)
{
    if (floatl_eq(x, FLOATL_CONST_0)) return floatl_neg(FLOATL_INF);
    return floatl_div(floatl_ln(x), floatl_ln(floatl(2.0)));
}

/**
 * \brief Calculates the logarithm of a 'floatl' number with base 10.
 * \param x: The 'floatl' number for which the logarithm is to be calculated.
 * \return A 'floatl' number representing the logarithm of 'x' with base 10. Returns negative infinity
 *         if 'x' is zero.
 */
floatl floatl_log10(floatl x)
{
    if (floatl_eq(x, FLOATL_CONST_0)) return floatl_neg(FLOATL_INF);
    return floatl_div(floatl_ln(x), floatl_ln(floatl(10.0)));
}

/**
 * \brief Calculates the sigmoid function of a 'floatl' number.
 * \param x: The 'floatl' number for which the sigmoid is to be calculated.
 * \return A 'floatl' number representing the sigmoid of 'x'.
 */
floatl floatl_sigmoid(floatl x)
{
    return floatl_div(FLOATL_CONST_1, floatl_add(FLOATL_CONST_1, floatl_exp(floatl_neg(x))));
}

/**
 * \brief Calculates the derivative of the sigmoid function of a 'floatl' number.
 * \param x: The 'floatl' number for which the derivative of the sigmoid is to be calculated.
 * \return A 'floatl' number representing the derivative of the sigmoid of 'x'.
 */
floatl floatl_dsigmoid(floatl x)
{
    floatl s = floatl_sigmoid(x);
    return floatl_mul(s, floatl_sub(FLOATL_CONST_1, s));
}

/**
 * \brief Calculates the hyperbolic sine of a 'floatl' number.
 * \param x: The 'floatl' number for which the hyperbolic sine is to be calculated.
 * \return A 'floatl' number representing the hyperbolic sine of 'x'.
 */
floatl floatl_sinh(floatl x)
{
    floatl a = floatl_exp(x);
    return floatl_div(floatl_sub(a, floatl_div(FLOATL_CONST_1, a)), floatl(2.0));
}

/**
 * \brief Calculates the hyperbolic cosine of a 'floatl' number.
 * \param x: The 'floatl' number for which the hyperbolic cosine is to be calculated.
 * \return A 'floatl' number representing the hyperbolic cosine of 'x'.
 */
floatl floatl_cosh(floatl x)
{
    floatl a = floatl_exp(x);
    return floatl_div(floatl_add(a, floatl_div(FLOATL_CONST_1, a)), floatl(2.0));
}

/**
 * \brief Calculates the hyperbolic tangent of a 'floatl' number.
 * \param x: The 'floatl' number for which the hyperbolic tangent is to be calculated.
 * \return A 'floatl' number representing the hyperbolic tangent of 'x'.
 */
floatl floatl_tanh(floatl x)
{
    floatl a = floatl_exp(x); 
    floatl b = floatl_div(FLOATL_CONST_1, a);
    return floatl_div(floatl_sub(a, b), floatl_add(a, b));
}

/**
 * \brief Calculates the hyperbolic cotangent of a 'floatl' number.
 * \param x: The 'floatl' number for which the hyperbolic cotangent is to be calculated.
 * \return A 'floatl' number representing the hyperbolic cotangent of 'x'.
 */
floatl floatl_coth(floatl x)
{
    floatl a = floatl_exp(x); 
    floatl b = floatl_div(FLOATL_CONST_1, a);
    return floatl_div(floatl_add(a, b), floatl_sub(a, b));
}

/**
 * \brief Calculates the hyperbolic secant of a 'floatl' number.
 * \param x: The 'floatl' number for which the hyperbolic secant is to be calculated.
 * \return A 'floatl' number representing the hyperbolic secant of 'x'.
 */
floatl floatl_sech(floatl x)
{
    floatl a = floatl_exp(x);
    return floatl_div(floatl(2.0), floatl_add(a, floatl_div(FLOATL_CONST_1, a)));
}

/**
 * \brief Calculates the hyperbolic cosecant of a 'floatl' number.
 * \param x: The 'floatl' number for which the hyperbolic cosecant is to be calculated.
 * \return A 'floatl' number representing the hyperbolic cosecant of 'x'.
 */
floatl floatl_csch(floatl x)
{
    floatl a = floatl_exp(x);
    return floatl_div(floatl(2.0), floatl_sub(a, floatl_div(FLOATL_CONST_1, a)));
}

/**
 * \brief Calculates the square root of a 'floatl' number using the Newton - Raphson method.
 * \param x: The 'floatl' number for which the square root is to be calculated.
 * \return A 'floatl' number representing the square root of 'x'. Returns -1 if 'x' is less than 0.
 */
floatl floatl_sqrt(floatl x) 
{
    if (floatl_lt(x, FLOATL_CONST_0)) return floatl_neg(FLOATL_CONST_1);
    floatl a = x; // Initial guess
    floatl epsilon = floatl(1e-7); // Precision
    while (floatl_gt(floatl_sub(floatl_mul(a, a), x), epsilon))
    {
        a = floatl_add(a, floatl_div(x, a));
        a = floatl_div(a, floatl(2.0));
    }
    return a;
}

/**
 * \brief Calculates the power of a 'floatl' number with another 'floatl' exponent.
 * \param x: The base 'floatl' number.
 * \param y: The exponent 'floatl' number.
 * \return A 'floatl' number representing 'x' raised to the power of 'y'.
 */
floatl floatl_pow(floatl x, floatl y) 
{
    return floatl_exp(floatl_mul(y, floatl_ln(x)));
}