/* solver.c */

#include "computorv1.h"

static double get_coefficient_of_power(int power, const s_polynom *polynom)
{
    int i;
    double coefficient;

    // if no term found with exponent, default coefficient is 0
    coefficient = 0;

    i = 0;
    while (polynom[i].sign != TERM_SIGN_END)
    {
        if (polynom[i].exponent == power)
        {
            coefficient = polynom[i].coefficient;
            break;
        }
        i++;
    }

    print_debug("- coefficient of power '%i' : %g (polynom[%i])\n", power, coefficient, i);
    return coefficient;
}

static double positiv_zero(double num)
{
    if (num == 0)
        return 0;
    return num;
}

int has_decimal_part(double num)
{
    return (num != (int)num);
}

int any_has_decimal_part(double *num, size_t len)
{
    while (len > 0)
    {
        if (has_decimal_part(num[len - 1]))
            return 1;
        len--;
    }
    return 0;
}

// find GCD of two integers using euclidean algorithm
static int gcd_int(int a, int b)
{
    int tmp;

    while (b != 0)
    {
        tmp = b;
        b = a % b;
        a = tmp;
    }
    return abs(a);
}

static void solve_degree_1(s_solution_degree_1 *solution, double a, double b)
{
    solution->a = a;
    solution->b = b;

    solution->solution = positiv_zero(-b / a); // -b / a
}

static void solve_degree_2(s_solution_degree_2 *solution, double a, double b, double c)
{
    double delta;
    double delta_sqrt;
    int gcd;

    solution->a = a;
    solution->b = b;
    solution->c = c;

    delta = b * b - 4 * a * c;
    solution->delta = delta;                                  // Δ == b² - 4ac
    solution->delta_sign = ft_fsign(delta);                   // DELTA_PLUS or DELTA_MINUS or DELTA_ZERO
    solution->delta_absolute = ft_fabs(delta);                // |Δ|
    solution->delta_sqrt = ft_sqrt(solution->delta_absolute); // √|Δ|
    delta_sqrt = solution->delta_sqrt;

    /**
     * SOLVE LEFT AND RIGHT TERMS AS DOUBLES
     */

    if (b != 0.0)
        solution->left_term = positiv_zero(-b / (2 * a)); // -b / 2a
    if (delta != 0.0)
        solution->right_term = positiv_zero(delta_sqrt / 2 * a); // √|Δ| / 2a

    /**
     * COMPLEX : SOLVE TERMS AS FRACTIONS
     */

    // check
    if (solution->delta_sign != DELTA_MINUS)
        return;
    if (any_has_decimal_part((double[]){b, a, solution->delta_sqrt}, 3))
        return;
    solution->all_int = true;

    // right term
    gcd = gcd_int((int)solution->delta_sqrt, (int)(a * 2));
    solution->right_term_numerator = solution->delta_sqrt / gcd;
    solution->right_term_denominator = (a * 2) / gcd;

    // left term
    if (b == 0.0)
        return;
    gcd = gcd_int((int)b, (int)(a * 2));
    solution->left_term_numerator = -b / gcd;
    solution->left_term_denominator = (a * 2) / gcd;
}

void solve(const s_polynom *polynom, s_solution *solution)
{
    double a;
    double b;
    double c;

    if (solution->degree == 1)
    {
        a = get_coefficient_of_power(1, polynom);
        b = get_coefficient_of_power(0, polynom);

        solve_degree_1(&solution->solution_degree_1, a, b);
    }
    else if (solution->degree == 2)
    {
        a = get_coefficient_of_power(2, polynom);
        b = get_coefficient_of_power(1, polynom);
        c = get_coefficient_of_power(0, polynom);

        solve_degree_2(&solution->solution_degree_2, a, b, c);
    }
}