cas.c

// anic17 - 2023/12/18 23:46

#include "cas.h"
#include "calculus.c"

const size_t allocsz = 512;
const double epsilon = 1e-13;
const char variable = 'x';

size_t strlen_n(char *s)
{
    if (!s)
        return 0;
    return strlen(s);
}

char *strncpy_n(char *dest, const char *src, size_t count)
{

    if (count)
    {
        memset(dest, 0, count);
        strncat(dest, src, count);
        return dest;
    }
    return NULL;
}

char *substring(const char *str, size_t length)
{
    char *sub = malloc(length + 1);
    if (sub == NULL)
    {
        return NULL;
    }
    strncpy(sub, str, length);
    sub[length] = '\0';

    return sub;
}

int is_number(int c)
{
    char sep[] = ".( +-";
    return isdigit(c) || strchr(sep, c) != NULL;
}

int mcd_euclidean(int a, int b)
{
    int t = 0;
    while (b != 0)
    {
        t = b;
        b = a % b;
        a = t;
    }
    return a;
}

int reduce(int *num1, int *num2)
{
    int divide = mcd_euclidean(*num1, *num2);
    *num1 /= divide;
    *num2 /= divide;
    return 0;
}

double eval_func(term_t **term, double x, size_t termcount)
{
    double y = 0;

    double coef_eval = 1, exp_eval = 1, div_eval = 1, multiply_func = 1;

    for (size_t k = 0; k < termcount; k++)
    {
        coef_eval = 1, exp_eval = 1, div_eval = 1, multiply_func = 1;
        if (term[k]->t_coef != NULL)
            coef_eval = eval_func(&term[k]->t_coef, x, term[k]->t_coef_sz); // funció recursiva que evalua les funcions en el coeficient
        coef_eval *= term[k]->coef;

        if (term[k]->t_exponent != NULL)
            exp_eval = eval_func(&term[k]->t_exponent, x, term[k]->t_exponent_sz);
        exp_eval *= term[k]->exponent;

        if (term[k]->t_division != NULL)
        {
            div_eval = eval_func(&term[k]->t_division, x, term[k]->t_division_sz);
        }
        div_eval *= term[k]->division;

        if (term[k]->t_multiply != NULL)
        {
            for (size_t i = 0; i < term[k]->t_multiply_sz; i++)
            {
                multiply_func *= eval_func(&term[k]->t_multiply[i], x, 1);
            }
        }

        if (term[k]->function != 0)
        {
            y += multiply_func * pow(func_list[term[k]->function].fn(coef_eval * pow(x, exp_eval) / div_eval), term[k]->func_exponent);
        }
        else
        {
            y += multiply_func * coef_eval * pow(x, exp_eval) / div_eval;
        }
    }
    return y;
}

term_t **simp_func(term_t **func, size_t termcount, bool simplify_division)
{
    int coef = 0, division = 0;
    for (size_t i = 0; i < termcount; i++)
    {
        if (simplify_division)
        {
            if (func[i]->t_division != NULL && func[i]->t_division->function == 0 && func[i]->t_division->t_coef == NULL && func[i]->t_division->t_multiply == NULL && func[i]->t_division->t_division == NULL && func[i]->t_division->t_exponent == NULL && func[i]->t_division_sz == 1)
            {
                func[i]->exponent -= func[i]->t_division->exponent;
                func[i]->coef *= func[i]->t_division->division;
                func[i]->division *= func[i]->t_division->coef;
                func[i]->t_division_sz = 0;
                free(func[i]->t_division);
                func[i]->t_division = NULL;
            }
            if (func[i]->function == SQRT)
            {
                if (func[i]->t_division != NULL && func[i]->t_division->t_division == NULL)
                {
                    func[i]->t_division->exponent /= 2;
                    func[i]->t_division->coef = sqrt(func[i]->t_division->coef);
                    func[i]->t_division->division = sqrt(func[i]->t_division->division);
                }
                func[i]->exponent /= 2;
                func[i]->coef = sqrt(func[i]->coef);
                func[i]->division = sqrt(func[i]->division);
                func[i]->function = 0;
            }
        }

        if (func[i]->coef == floor(func[i]->coef) && func[i]->division == floor(func[i]->division))
        {
            coef = (int)func[i]->coef;
            division = (int)func[i]->division;
            reduce(&coef, &division);
            func[i]->coef = (double)coef;
            func[i]->division = (double)division;
        }
        else if (modf(func[i]->coef / func[i]->division, &(double){0}) == 0.0)
        {
            func[i]->coef /= func[i]->division;
            func[i]->division = 1;
        }
    }
    return func;
}

double bisection_method(term_t **func, size_t count, double a, double b)
{
    size_t iterations = 100;
    double x = 0;
    size_t n = 0;
    if (a > b)
    {
        double tmp = a;
        a = b;
        b = tmp;
    }
    while (sgn_fn(eval_func(func, a, count)) == sgn_fn(eval_func(func, b, count)) && n < 50)
    {

        if (!isnan(a * 2 + 5))
        {
            a = a * 2 + 5;
        }
        if (!isnan(b / 3 - 12))
        {
            b = b / 3 - 12;
        }

        if (a != 0.0)
        {
            b /= a;
        }

        n++;
    }
    if (n >= 50)
    {
        return NAN; // l'algoritme no troba un interval amb dos signes oposats
    }
    double midpoint = (a + b) / 2;
    double ev_a = 0, ev_b = 0;
    for (size_t i = 0; i < iterations; i++)
    {
        x = eval_func(func, midpoint, count);
        ev_a = eval_func(func, a, count);
        ev_b = eval_func(func, b, count);
        if (isnan(x) || isinf(x) || isnan(ev_a) || isinf(ev_a) || isnan(ev_b) || isinf(ev_b)) // si la funció no es continua en l'interval el mètode no funciona
        {
            return NAN;
        }
        if (signbit(eval_func(func, a, count)) == signbit(eval_func(func, b, count)))
        {
            return NAN;
        }
        if (x < 0.0)
        {
            a = midpoint;
        }
        else
        {
            b = midpoint;
        }
        midpoint = (a + b) / 2;
    }
    if (fabs(eval_func(func, midpoint, count)) <= epsilon)
    {
        return midpoint;
    }
    return NAN; // no troba cap solució real
}

int is_polynomial(term_t **func, size_t termcount, bool allow_func)
{
    for (size_t i = 0; i < termcount; i++)
    {
        if (func[i]->t_coef != NULL || func[i]->t_division != NULL || func[i]->t_exponent != NULL || func[i]->t_multiply != NULL || (func[i]->function != 0 && !allow_func))
        {
            return false;
        }
    }
    return true;
}

char *print_num(double num, char *buf, size_t allocsz)
{
    typedef struct number_codes_t
    {
        double num;
        char name[10];
    } number_codes_t;
    if (num == 0.0)
    {
        snprintf(buf, allocsz, "0");
        return buf;
    }

    if (isnan(num))
    {
        snprintf(buf, allocsz, "NaN");
        return buf;
    }
    else if (isinf(num) || fabs(num) > 1e300)
    {
        snprintf(buf, allocsz, signbit(num) ? "-oo" : "+oo");
        return buf;
    }

#ifdef WIN32
    number_codes_t s_nums[3] = {{M_PI, "pi"}, {M_E, "e"}, {1.61803398874989484820, "phi"}};
#else
    number_codes_t s_nums[3] = {{M_PI, "π"}, {M_E, "e"}, {1.61803398874989484820, "φ"}};
#endif

    for (size_t i = 0; i < 2; i++)
    {
        if (pow(s_nums[i].num, 2) == num)
        {
            snprintf(buf, allocsz, "%s^2", s_nums->name);
            return buf;
        }
        else if (s_nums[i].num == pow(num, 2))
        {

            snprintf(buf, allocsz, "sqrt(%s)", s_nums->name);
            return buf;
        }
        else if (fabs(num) == s_nums[i].num)
        {
            snprintf(buf, allocsz, num < 0 ? "-%s" : "%s", s_nums->name);
            return buf;
        }
        else if (fabs(num) > 1e-17 && num != 0.0 && floor(100 * num / s_nums[i].num) == 100 * num / s_nums[i].num && num != s_nums[i].num)
        {
            if ((fabs(num) / s_nums[i].num) - floor(fabs(num) / s_nums[i].num) < epsilon)
            {
                snprintf(buf, allocsz, num < 0 ? "-%.08g%s" : "%.08g%s", fabs(num) / s_nums[i].num, s_nums[i].name);
            }
            else
            {
                snprintf(buf, allocsz, num < 0 ? "-%s/%.08g" : "%s/%.08g", s_nums[i].name, s_nums[i].num / fabs(num));
            }

            return buf;
        }
    }

    double intp = 0, fracp = 0;
    fracp = modf(pow(num, 2), &intp);
    if (fracp <= 1e-15 && floor(num) != num && intp != 0.0)
    {
        snprintf(buf, allocsz, "sqrt(%.08g)", intp);
    }
    else
    {
        if (fabs(num - floor(num)) < 1e-15) // la precisió del double és de 14 digits
        {
            snprintf(buf, allocsz, "%.08g", floor(num));
        }
        else
        {
            snprintf(buf, allocsz, "%.08g", num);
        }
    }
    return buf;
}

int print_func(term_t **func, size_t termcount)
{
    char *buf = calloc(allocsz + 1, 1);
    for (size_t i = 0; i < termcount; i++)
    {
        if (func[i] == NULL)
        {
            continue;
        }

        if (func[i]->coef != 0.0)
        {
            if (func[i]->function != 0 && func[i]->function != FACT)
            {
                if (func[i]->function == ABS)
                {
                    printf("|");
                }
                else
                {
                    printf("%s", func_list[func[i]->function].name);
                    if (func[i]->function != 0 && func[i]->func_exponent != 1)
                    {
                        print_num(func[i]->func_exponent, buf, allocsz);
                        printf("^%s(", buf);
                    }
                    else
                    {
                        printf("(");
                    }
                }
            }

            if (func[i]->coef != 1.0)
            {
                print_num(func[i]->coef, buf, allocsz);
                printf("%s", buf);
            }
            if (func[i]->t_multiply != NULL)
            {
                for (size_t k = 0; k < func[i]->t_multiply_sz; k++)
                {
                    print_func(&func[i]->t_multiply[k], 1);
                }
            }
            if (func[i]->t_coef != NULL)
            {
                print_func(&func[i]->t_coef, func[i]->t_coef_sz);
            }

            if (!(func[i]->exponent == 0.0 && func[i]->t_exponent == NULL))
            {
                printf("%c", func[i]->variable);
            }
            else if (func[i]->t_coef == NULL && func[i]->t_multiply == NULL && func[i]->coef == 1.0)
            {
                printf("1");
            }

            if (func[i]->t_exponent == NULL)
            {
                if (func[i]->exponent != 1.0 && func[i]->exponent != 0.0)
                {
                    print_num(func[i]->exponent, buf, allocsz);
                    printf("^%s", buf);
                }
            }
            else
            {
                printf("^(");
                print_func(&func[i]->t_exponent, func[i]->t_exponent_sz);
                printf(")");
            }

            if (func[i]->division != 1.0)
            {
                print_num(func[i]->division, buf, allocsz);
                if (func[i]->t_division == NULL)
                {
                    printf("/%s", buf);
                }
                else
                {
                    printf("/(%s", buf);
                }
            }
            if (func[i]->t_division != NULL)
            {
                if (func[i]->division == 1.0)
                {
                    printf("/(");
                }
                print_func(&func[i]->t_division, 1);
                printf(")");
            }
            switch (func[i]->function)
            {
            case ABS:
                printf("|");
                break;
            case FACT:
                printf("!");
                break;
            case 0:
                break;
            default:
                printf(")");
                break;
            }
        }

        if (i < termcount - 1 && func[i]->coef != 0.0)
        {
            printf(" + ");
        }
    }
    free(buf);
    return 0;
}

double newton_method(term_t **func, size_t count, double x0)
{
    term_t **dx = malloc(sizeof(term_t) * count);
    memcpy(dx, func, count * sizeof(dx));

    bool derivative_valid = false;

    double x = x0, eval_dx = 0, ev_func = 0;

    if (x0 == 0.0 && eval_func(func, count, 0.0) < epsilon)
    {
        return 0.0;
    }
    /*
        Algoritme del mètode de Newton:
        x_n+1 = x_n - f(x_n)/f'(x_n)
    */

    size_t iterations = 100;

    for (size_t k = 0; k < iterations; k++)
    {
        if (derivative_valid)
        {
            eval_dx = eval_func(func, x, count);
        }
        else
        {
            eval_dx = diff_point(func, x, count);
        }
        ev_func = eval_func(func, x, count);

        x = x - ev_func / eval_dx;
        if (eval_dx == 0.0 || isnan(ev_func) || isnan(eval_dx))
        {
            return NAN;
        }
    }
    if (fabs(x) < 1e-8)
    {
        if (fabs(eval_func(func, 0, count)) < epsilon)
            return 0;
    }
    if (fabs(eval_func(func, x, count)) < epsilon)
        return x;

    return NAN;
}

double secant_method(term_t **func, size_t count, double a, double b)
{
    term_t **dx = malloc(sizeof(term_t) * count);
    memcpy(dx, func, count * sizeof(dx));

    if (eval_func(func, 0.0, count) == 0.0)
    {
        return 0;
    }

    double x = 0, x1 = a, x2 = b;

    double ev_x1 = 0, ev_x2 = 0;
    const size_t iterations = 07;

    for (size_t k = 0; k < iterations; k++)
    {
        ev_x1 = eval_func(func, x1, count);
        ev_x2 = eval_func(func, x2, count);

        x = x1 - ev_x1 * (x1 - x2) / (ev_x1 - ev_x2);
        x2 = x1;
        x1 = x;
    }
    if (fabs(eval_func(func, x, count)) < epsilon)
        return x;
    return NAN;
}

int cmpfunc(const void *a, const void *b)
{
    return (*(double *)a > *(double *)b) - (*(double *)a < *(double *)b);
}

double solve_equation(term_t **func, size_t termcount, double x0, double bs_a, double bs_b)
{
    double result = 0, old_result = 0;
    char *fx = calloc(allocsz, sizeof(char));
    double found[162] = {0};
    size_t found_index = 0;
    bool found_loop = false;
    for (int i = -40; i < 40; i++)
    {
        for (int k = -1; k <= 1; k += 2)
        {
            if (i == 0)
                result = newton_method(func, termcount, 0);
            else
                result = newton_method(func, termcount, sgn_fn(i) * pow(1.5, k * abs(i)));

            if (!isnan(result) && old_result != result)
            {
                found_loop = false;
                for (size_t k = 0; k <= found_index; k++)
                {
                    if (fabs(found[k] - result) < epsilon)
                    {
                        found_loop = true;
                    }
                }
                if (!found_loop)
                {
                    if (found_index == 0 || (found_index > 0 && fabs(found[found_index - 1] - result) > 1e-8) || i == 0)
                    {
                        found[found_index] = result;
                        found_index++;
                    }
                }
            }
            old_result = result;
        }
    }
    if (found_index > 0)
    {
        qsort(found, found_index, sizeof(found[0]), cmpfunc);
        printf("Arrels trobades de ");
        print_func(func, termcount);
        printf(" :\n");
        for (size_t i = 0; i < found_index; i++)
        {
            print_num(found[i], fx, 120);
            printf("x = %s\n", fx);
        }

        return found[0];
    }
    result = bisection_method(func, termcount, bs_a, bs_b);
    if (!isnan(result))
    {
        printf("Arrel trobada de ");
        print_func(func, termcount);
        print_num(result, fx, 120);
        printf(" :\nx = %s\n", fx);
        return result;
    }
    result = secant_method(func, termcount, x0, bs_a);
    if (!isnan(result))
    {
        printf("Arrel trobada de ");
        print_func(func, termcount);
        print_num(result, fx, 120);
        printf(" :\nx = %s\n", fx);
    }
    return result;
}

char *remove_whitespace(char *s)
{

    char *s2 = calloc(strlen_n(s) + 1, sizeof(char)); // variable temporal
    size_t count = 0;
    for (size_t i = 0; i < strlen_n(s); i++)
    {
        if (!isspace(s[i]))
            s2[count++] = s[i];
    }
    strncpy_n(s, s2, strlen_n(s2));
    free(s2);
    return s;
}

term_t *create_term()
{
    term_t *term = calloc(sizeof(term_t), 1);
    if (term == NULL)
    {
        fprintf(stderr, "Error: No hi ha prou espai en memòria per crear un tipus 'term_t'.\n");
        exit(EXIT_FAILURE);
    }
    term->coef = 0;
    term->exponent = 0;
    term->division = 1;
    term->function = 0;
    term->variable = variable;
    term->t_coef = NULL;
    term->t_exponent = NULL;
    term->t_division = NULL;
    term->t_multiply = NULL;
    term->t_coef_sz = 0;
    term->t_exponent_sz = 0;
    term->t_division_sz = 0;
    term->t_multiply_sz = 0;

    return term;
}

term_t **multiply_terms(term_t *f1, term_t *f2)
{

    term_t **mult_func = calloc(sizeof(term_t **), 1);

    mult_func[0] = create_term();
    mult_func[0]->coef = 1;

    if (f2 == NULL)
    {
        memcpy(mult_func[0], f1, sizeof(term_t));
        return mult_func;
    }
    if (f1 == NULL)
    {
        memcpy(mult_func[0], f2, sizeof(term_t));
        return mult_func;
    }

    size_t m_count = f1->t_multiply_sz + f2->t_multiply_sz;

    if (f1->function != 0)
        m_count++;
    if (f2->function != 0)
        m_count++;

    mult_func[0]->coef = 1;
    mult_func[0]->division = 1;
    mult_func[0]->exponent = 0;
    mult_func[0]->t_multiply_sz = 0;

    if (m_count)
    {
        mult_func[0]->t_multiply = calloc(sizeof(term_t **), m_count);
    }
    for (size_t i = 0; i < m_count; i++)
    {
        mult_func[0]->t_multiply[i] = create_term();
    }

    for (size_t k = 0; k < f1->t_multiply_sz; k++)
    {
        memcpy(mult_func[0]->t_multiply[k], f1->t_multiply[k], sizeof(term_t));
    }
    mult_func[0]->t_multiply_sz += f1->t_multiply_sz;

    for (size_t k = 0; k < f2->t_multiply_sz; k++)
    {
        memcpy(mult_func[0]->t_multiply[mult_func[0]->t_multiply_sz], f2->t_multiply[k], sizeof(term_t));
    }

    mult_func[0]->t_multiply_sz += f2->t_multiply_sz;

    if (f1->function != 0)
    {
        memcpy(mult_func[0]->t_multiply[mult_func[0]->t_multiply_sz++], f1, sizeof(term_t));
    }
    else
    {
        mult_func[0]->coef *= f1->coef;
        mult_func[0]->division *= f1->division;
        mult_func[0]->exponent += f1->exponent;
    }
    if (f2->function != 0)
    {
        memcpy(mult_func[0]->t_multiply[mult_func[0]->t_multiply_sz++], f2, sizeof(term_t));
    }
    else
    {
        mult_func[0]->coef *= f2->coef;
        mult_func[0]->division *= f2->division;
        mult_func[0]->exponent += f2->exponent;
    }
    return mult_func;
}

term_t **internal_add_term(term_t **sum_func, term_t *f1, term_t *f2, size_t *finalsize)
{
    if (is_polynomial(&f1, 1, false) && is_polynomial(&f2, 1, false) && f1->variable == f2->variable && f1->exponent == f2->exponent)
    {
        sum_func[0]->coef = f1->coef * f2->division + f2->coef * f1->division;
        sum_func[0]->division = f1->division * f2->division;
    }
    else
    {
        memcpy(sum_func[1], f2, sizeof(*f2));
        *finalsize = 2;
    }
    return sum_func;
}

term_t **add_terms(term_t *f1, term_t *f2, size_t *finalsize) // suma f1 i f2 a f1
{
    term_t **sum_func = calloc(sizeof(f1), 2);
    sum_func[0] = calloc(sizeof(term_t), 1);
    sum_func[1] = calloc(sizeof(term_t), 1);
    memcpy(sum_func[0], f1, sizeof(*f1));
    *finalsize = 1;
    internal_add_term(sum_func, f1, f2, finalsize);

    if (*finalsize == 1)
    {
        free(sum_func[1]);
    }

    return sum_func;
}

void sigint_h()
{
    printf("\nSenyal SIGINT: Control-C. Finalitzant programa...\n");
    exit(0);
}

void version_info()
{
    printf("Computer Algebra System (CAS) v%s (%s)\nAutor: anic17\n\n", cas_version, date);
}

domain_t *alloc_dom(domain_t *dom, size_t domsize)
{
    dom->lowbound = calloc(sizeof(double), domsize);
    dom->highbound = calloc(sizeof(double), domsize);
    dom->is_open_low = calloc(sizeof(int), domsize);
    dom->is_open_high = calloc(sizeof(int), domsize);
    dom->is_periodic = calloc(sizeof(int), domsize);
    dom->period = calloc(sizeof(double), domsize);
    dom->x0 = calloc(sizeof(double), domsize);
    dom->holes = calloc(sizeof(double), domsize);
    dom->holecount = 0;

    return dom;
}

domain_t *set_domain(domain_t *dom, size_t index, double lowbound, double highbound, int is_open_low, int is_open_high, int is_periodic, double period, double x0)
{
    if (is_periodic)
    {
        dom->is_periodic[index] = is_periodic;
        dom->period[index] = period;
        dom->x0[index] = x0;
    }
    dom->lowbound[index] = lowbound;
    dom->highbound[index] = highbound;
    dom->is_open_low[index] = is_open_low;
    dom->is_open_high[index] = is_open_high;
    dom->domcount++;
    return dom;
}

domain_t *domain_func(term_t *func, size_t count, domain_t *dom, size_t *index)
{
    if (func->t_coef != NULL)
        domain_func(func->t_coef, func->t_coef_sz, dom, index);
    if (func->t_exponent != NULL)
        domain_func(func->t_exponent, func->t_exponent_sz, dom, index);
    if (func->t_division != NULL)
    {

        double root1;
        bool foundHole = false;
        for (int i = -40; i < 40; i++)
        {
            root1 = newton_method(&func->t_division, func->t_division_sz, pow(1.5, i));

            if (!isnan(root1))
            {
                foundHole = false;
                for (size_t k = 0; k < dom->holecount; k++)
                {
                    if (dom->holes[k] == root1)
                    {
                        foundHole = true;
                    }
                }
                if (!foundHole)
                {
                    if (dom->holecount == 0 || (dom->holecount > 0 && fabs(dom->holes[dom->holecount] - root1) > 1e-8))
                        dom->holes[dom->holecount++] = root1;
                }
            }
        }

        domain_func(func->t_division, func->t_division_sz, dom, index);
    }
    if (func->function != 0)
    {
        switch (func->function)
        {
        case SIN:
        case COS:
        case ABS:
        case DELTA:
        case SGN:
        case ATAN:
        case ACOT:
        case SINH:
        case COSH:
        case TANH:
        case SECH:
            set_domain(dom, (*index)++, -INF, +INF, true, true, false, 0, 0);
            break;
        case INV:
            set_domain(dom, (*index)++, -INF, +INF, true, true, false, 0, 0);
            dom->holes[dom->holecount++] = 0.0;
            break;

        case SEC:
        case TAN:
            set_domain(dom, (*index)++, -INF, +INF, true, true, true, M_PI, -M_PI / 2);
            break;

        case COTH:
            set_domain(dom, (*index)++, -INF, +INF, true, true, false, 0, 0);
            dom->holes[dom->holecount++] = 0.0;
            break;

        case CSC:
        case COT:
            set_domain(dom, (*index)++, -INF, +INF, true, true, true, M_PI, 0);
            break;
        case ASIN:
        case ACOS:
            set_domain(dom, (*index)++, -1, 1, false, false, false, 0, 0);
            break;
        case ASECH:
            set_domain(dom, (*index)++, 0, 1, true, false, false, 0, 0);
            break;
        case CSCH:
        case ACSCH:

        case LN:
        case LOG:
        case LOG2:
        case LOG10:
        case LB:
            set_domain(dom, (*index)++, 0, +INF, true, true, false, 0, 0);
            break;
        case SQRT:
            set_domain(dom, (*index)++, 0, +INF, false, true, false, 0, 0);
            break;
        case GAMMA:
        case FACT:
            set_domain(dom, (*index)++, -1, +INF, false, true, false, 0, 0);
            break;
        default:
            break;
        }
        if (((func->t_coef != NULL || func->t_multiply != NULL || func->t_division != NULL || func->t_exponent != NULL) && !(isinf(dom->lowbound[*index - 1]) && isinf(dom->highbound[*index - 1]))) || (dom->is_periodic[*index - 1] && func->exponent != 1 && func->exponent != 0))
        {
            puts("\nError: No es pot obtenir el domini d'una funció composta");
            set_domain(dom, (*index)++, 0, 0, true, true, false, 0, 0);
            return NULL;
        }

        dom->lowbound[*index - 1] /= (func->coef / func->division);
        dom->highbound[*index - 1] /= (func->coef / func->division);
        if (dom->is_periodic[*index - 1])
        {
            dom->period[*index - 1] /= (func->coef / func->division);
            dom->x0[*index - 1] /= (func->coef / func->division);
        }
    }
    else
    {
        set_domain(dom, (*index)++, -INF, +INF, true, true, false, 0, 0);
    }

    return dom;
}

domain_t *get_domain(term_t **func, size_t termcount)
{
    size_t domsize = 1000;
    domain_t *dom = calloc(sizeof(domain_t), domsize);

    alloc_dom(dom, domsize);

    double dswitch = 0;
    int openswitch = 0;
    size_t index = 0;
    for (size_t i = 0; i < termcount; i++)
    {
        if (!domain_func(func[i], 1, dom, &index))
            return NULL;
    }
    char *buf1 = calloc(allocsz + 1, 1), *buf2 = calloc(allocsz + 1, 1);

    for (size_t k = 0; k < index; k++)
    {

        if (dom->highbound[k] < dom->lowbound[k])
        { // inverteix els intervals en cas que siguin contràris
            dswitch = dom->lowbound[k];
            dom->lowbound[k] = dom->highbound[k];
            dom->highbound[k] = dswitch;

            openswitch = dom->is_open_low[k];
            dom->is_open_low[k] = dom->is_open_high[k];
            dom->is_open_high[k] = openswitch;
        }

        if ((dom->lowbound[k] >= dom->highbound[k] && (dom->is_open_high[k] || dom->is_open_low)) || isnan(dom->lowbound[k]) || isnan(dom->highbound[k]))
        {
            fprintf(stderr, "No s'ha pogut obtenir el domini de la funció o bé la funció no té cap imatge real");
            return NULL;
        }
    }
    domain_t secdom = {0};
    secdom.lowbound = calloc(sizeof(double), 1);
    secdom.highbound = calloc(sizeof(double), 1);
    secdom.is_open_low = calloc(sizeof(int), 1);
    secdom.is_open_high = calloc(sizeof(int), 1);

    secdom.lowbound[0] = dom->lowbound[0];
    secdom.highbound[0] = dom->highbound[0];
    secdom.is_open_low[0] = dom->is_open_low[0];
    secdom.is_open_high[0] = dom->is_open_high[0];

    double *x0_arr = calloc(sizeof(double), domsize);
    double *period_arr = calloc(sizeof(double), domsize);
    size_t period_count = 0;

    for (size_t k = 0; k < index; k++)
    {
        if (dom->lowbound[k] > secdom.lowbound[0])
        {
            secdom.is_open_low[0] = dom->is_open_low[k];
            secdom.lowbound[0] = dom->lowbound[k];
        }

        if (dom->highbound[k] < secdom.highbound[0])
        {
            secdom.is_open_high[0] = dom->is_open_high[k];
            secdom.highbound[0] = dom->highbound[k];
        }
        if (dom->is_periodic[k])
        {
            x0_arr[period_count] = dom->x0[k];
            period_arr[period_count] = dom->period[k];

            period_count++;
        }
    }

    if (*secdom.lowbound == *secdom.highbound || (check_if_no_func(func, termcount) && isnan(eval_func(func, (*secdom.lowbound), termcount))))
    {
        fprintf(stderr, "La funció no té cap imatge real");
        return NULL;
    }

    print_num(*secdom.lowbound, buf1, allocsz);
    print_num(*secdom.highbound, buf2, allocsz);
    printf("%c%s, %s%c", *(secdom.is_open_low) ? '(' : '[', buf1, buf2, *(secdom.is_open_high) ? ')' : ']');

    if (period_count)
    {
        if (!check_if_no_func(func, termcount))
        {
            printf(" - {");
            for (size_t k = 0; k < period_count; k++)
            {

                memset(buf1, 0, allocsz);
                memset(buf2, 0, allocsz);

                print_num(period_arr[k], buf1, allocsz);
                print_num(x0_arr[k], buf2, allocsz);
                printf("%s + %s*k, k = Z", buf2, buf1);
                if (k + 1 < period_count)
                {
                    printf(" | ");
                }
            }
            putchar('}');
        }
    }

    if (dom->holecount > 0)
    {
        char *holeprint = calloc(allocsz + 1, sizeof(char));
        printf(" - {");
        for (size_t i = 0; i < dom->holecount; i++)
        {
            print_num(dom->holes[i], holeprint, allocsz);
            printf("%s", holeprint);

            if (i < dom->holecount - 1)
            {
                printf(", ");
            }
        }
        printf("}");
        free(holeprint);
    }

    free(buf1);
    free(buf2);
    free(x0_arr);
    free(period_arr);
    free(dom->holes);

    free(dom);
    return NULL;
}

int get_func_num(char *func_name, bool checkwhole)
{
    for (size_t i = 1; i <= FUNC_COUNT; i++)
    {
        if (!strcmp(func_list[i].name, func_name))
        {
            return i;
        }
    }

    return -1;
}

int check_if_no_func(term_t **func, size_t termcount)
{
    for (size_t i = 0; i < termcount; i++)
    {
        if (func[i]->exponent != 0.0 || func[i]->t_multiply != NULL || func[i]->t_exponent != NULL || func[i]->t_coef != NULL || (func[i]->t_division != NULL && func[i]->t_division->exponent != 0.0))
        {
            return false;
        }
    }
    return true;
}

double get_num(char *s)
{
    size_t inc = 0;
    int sign = 1;
    double num = 0;
    while (*s != '\0' && *s == '-')
    {
        sign *= -1;
        s++;
    }
    while (*(s + inc) != '\0' && is_number(*(s + inc)))
    {
        inc++;
    }
    char *s2;
    if (inc == 0)
    {
        num = sign;
    }
    else
    {
        s2 = substring(s, inc);
        num = atof(s2);
        free(s2);
    }
    return num;
}

term_t *get_expression_inside(char *ptr, term_t *term)
{
    char *s_num;
    size_t num_index = 0;
    while (ptr[num_index] != '\0' && is_number(ptr[num_index]))
    {
        num_index++;
    }
    if (num_index)
    {
        s_num = substring(ptr, num_index);
        term->division = atof(s_num);
        free(s_num);
        ptr += num_index;
    }
    if (*ptr == '\0' || *ptr == ')')
        return term;

    term->t_division = create_term();

    for (size_t k = 0; k < strlen_n(ptr); k++)
    {
        if (ptr[k] == variable)
        {
            term->t_division->coef = 1;
            break;
        }
    }

    if (term->t_division->coef == 1)
    {
        term->t_division->coef = get_num(ptr);
        term->t_division->exponent = 1;
        term->t_division_sz = 1;
        parse_input(ptr, term->t_division);
    }
    else
    {
        free(term->t_division);
    }
    return term;
}

term_t *parse_input(char *func, term_t *term)
{
    char *fx = strdup(func);
    int sign = 1;

    term->coef = 1.0;
    term->exponent = 0.0;
    term->division = 1.0;
    term->func_exponent = 0.0;
    term->variable = variable;

    term->t_coef = NULL;
    term->t_exponent = NULL;
    term->t_division = NULL;
    term->t_multiply = NULL;

    size_t n_index = 0, func_index = 0;

    while (fx[func_index] != '\0' && fx[func_index] != '(' && fx[func_index] != ')')
    {
        func_index++;
    }
    char *s = substring(fx, func_index);
    int func_code = get_func_num(s, false);
    if (func_code > 0)
    {
        term->function = func_code;
        term->func_exponent = 1;
        fx += func_index;
    }
    free(s);
    while (isblank(*fx) || *fx == '(' || *fx == ')')
    {
        fx++;
    }
    if (!strncmp(fx, "PI", 2) || !strncmp(fx, "pi", 2))
    {
        term->coef = M_PI;
        fx += 2;
    }
    else
    {
        while (*fx != '\0' && *fx == '-')
        {
            sign *= -1;

            fx++;
        }
        while (is_number(fx[n_index]))
        {
            n_index++;
        }

        s = substring(fx, n_index);
        if (s[0] == '\0')
        {
            term->coef = sign;
        }
        else
        {
            term->coef = sign * atof(s);
        }
        fx += n_index;
        free(s);
    }

    int mult_count = 0;

    while (*fx == variable)
    {
        mult_count++;
        fx++;
    }
    if (*fx && *fx == '^')
    {
        fx++;

        char *s2 = strdup(fx);
        size_t charcount = 0;

        while (*s2 && *s2 != '/')
        {
            charcount++;
            s2++;
        }

        if (charcount)
        {
            char *sbs = substring(fx, charcount);
            term->exponent = atof(sbs) - 1;
            fx += strlen_n(sbs);
            free(sbs);
        }
    }

    if (mult_count == 0)
    {
        if (term->exponent != 0.0)
        {
            term->coef = pow(term->coef, term->exponent);
            term->exponent = 0;
        }
    }
    else
    {
        term->exponent += (double)mult_count;
    }
    while (*fx == '/')
    {
        fx++;
       // char *ptr = strtok(fx, "/");
        
            get_expression_inside(fx, term);
            //ptr = strtok(NULL, "/");
        
    }
    return term;
}

size_t count_terms_input(const char *expression, size_t *termcount)
{
    size_t len = strlen(expression);
    int parenthesis_cnt = 0;
    *termcount = 1;

    for (size_t i = 0; i < len; ++i)
    {
        if (expression[i] == '(')
        {
            parenthesis_cnt++;
        }
        else if (expression[i] == ')')
        {
            if (parenthesis_cnt > 0)
            {
                parenthesis_cnt--;
            }
        }
        else if (expression[i] == '+' && parenthesis_cnt == 0 && i > 0)
        {

            (*termcount)++;
        }
    }
    return *termcount;
}

int transform_input(char *str, size_t size)
{
    str[strcspn(str, "\r\n")] = '\0';

    remove_whitespace(str);

    size_t len = strlen_n(str);
    int parenthesis_level = 0;
    for (size_t i = 0; i < len; ++i)
    {
        if (str[i] == '(')
        {
            parenthesis_level++;
        }
        else if (str[i] == ')')
        {
            if (parenthesis_level <= 0)
            {
                fprintf(stderr, "Error: No hi ha cap parèntesi obert.\n");
                return false;
            }
            parenthesis_level--;
        }
        if (str[i] == '-' && i < len - 1 && str[i + 1] == '-')
        {
            memmove(&str[i], &str[i + 2], len - i);
            len -= 2;
            continue;
        }
        if ((i >= len - 1 && str[i] == '-') || (str[i] == '-' && i > 0 && !is_number(str[i])))
        {
            continue;
        }
        if (str[i] == '-' && parenthesis_level == 0)
        {
            if (len + 1 < size)
            {
                memmove(&str[i + 1], &str[i], len - i);
                str[i] = '+';
                len++;
                i++;
            }
            else
            {
                fprintf(stderr, "Error: No hi ha prou espai en memòria per processar l'entrada\nAbortant procés...\n");
                return false;
            }
        }
    }
    if (parenthesis_level > 0)
    {
        fprintf(stderr, "Error: Hi ha parèntesis oberts.\n");
        return false;
    }
    return true;
}

term_t **build_func(size_t *termcount, char *fname, bool emptyfunc)
{

    char *func_name = calloc(allocsz + 1, sizeof(char));
    printf("f(x)=");

    if (fname != NULL)
    {
        memcpy(func_name, fname, allocsz);
        puts(fname);
    }
    else
    {
        if (fgets(func_name, allocsz, stdin) == NULL)
        {
            fprintf(stderr, "Error: No s'ha pogut llegir l'entrada\n");
            free(func_name);
            return NULL;
        }
    }
    term_t **function = NULL;
    function = calloc(sizeof(term_t *), *termcount);

    if (function == NULL)
    {
        fprintf(stderr, "Error: No hi ha prou espai en memòria.\n");
        free(func_name);
        return NULL;
    }
    if (!transform_input(func_name, allocsz))
    {
        exit(0);
    }
    if (func_name[0] == '\0')
    {
        if (emptyfunc)
        {
            fprintf(stderr, "Error: Funció buida.\n");
            exit(0);
        }
        return NULL;
    }

    *termcount = count_terms_input(func_name, termcount);

    for (size_t i = 0; i < *termcount; i++)
    {
        function[i] = create_term();
    }
    size_t termindex = 0;

    char *ptr = strtok(func_name, "+");
    while (ptr != NULL && termindex < *termcount)
    {
        parse_input(ptr, function[termindex]);

        ptr = strtok(NULL, "+");
        termindex++;
    }

    printf("\nExpressió interpretada com a:\n");
    print_func(function, termindex);

    free(func_name);

    return function;
}

int eat_stdin_newline()
{
    int ch = getchar();
    if (ch != EOF && ch != '\n')
    {
        ungetc(ch, stdin);
    }
    return ch;
}

double initial_value(char *str, double defaultnum)
{

    char *buffer = calloc(sizeof(char), allocsz + 1);
    double initial_val = 0;
    double mult = 1;
    printf("%s", str);
    if (fgets(buffer, allocsz, stdin) != NULL)
    {
        buffer[strcspn(buffer, "\r\n")] = '\0';

        if (buffer[0] != '\n')
        {
            while (!strncmp(buffer, "PI", 2) || !strncmp(buffer, "pi", 2))
            {
                mult *= M_PI;
                buffer += 2;
            }

            if (!strncmp(buffer, "-oo", 3))
            {
                return -pow(2, 1023.99);
            }
            if (!strncmp(buffer, "+oo", 3) || !strncmp(buffer, "oo", 2))
            {
                return pow(2, 1023.99);
            }
            if (buffer[0] != '\0' && sscanf(buffer, "%lf", &initial_val) != 1 && mult == 1)
            {
                printf("Error: Entrada invàlida! Utilitzant el valor per defecte (%.08g).\n", defaultnum);
                return defaultnum;
            }

            if (initial_val == 0 && mult != 1)
            {
                return mult;
            }
            return initial_val * mult;
        }
    }
    return defaultnum;
}

void domain_program(term_t **func, size_t termcount)
{

    printf("Dom(");
    print_func(func, termcount);
    printf(") = ");
    get_domain(func, termcount);
}

void root_program(term_t **func, size_t termcount)
{

    double bs_a = initial_value("1r valor mètode bisecció: ", 1);
    double bs_b = initial_value("2n valor mètode bisecció: ", 8);
    printf("\n");
    double eqsolve = solve_equation(func, termcount, 1, bs_a, bs_b);

    if (isnan(eqsolve))
    {
        printf("No s'han trobat solucions per l'equació ");
        print_func(func, termcount);
    }
}

void integral_program(term_t **func, size_t termcount)
{
    double int_lowbound = initial_value("Límit inferior: ", 0);
    double int_highbound = initial_value("Límit superior: ", 0);

    double integral_value = riemann_sum(func, termcount, int_lowbound, int_highbound);
    if (isnan(integral_value))
    {
        printf("No s'ha pogut calcular l'integral definida\n");
    }
    else
    {
        char *fx1 = calloc(allocsz + 1, 1), *fx2 = calloc(allocsz + 1, 1), *print_intg = calloc(allocsz + 1, 1);
        print_num(int_lowbound, fx1, allocsz);
        print_num(int_highbound, fx2, allocsz);
        print_num(integral_value, print_intg, allocsz);

        printf("%s\n", fx2);
        printf("∫ ");
        print_func(func, termcount);
        printf("dx = %s\n", print_intg);
        printf("%s\n", fx1);
        free(fx1);
        free(fx2);
        free(print_intg);
    }
}

void derivative_program(term_t **func, size_t termcount)
{
    size_t dx_termcount = 0, tmp_termcount = 0;
    simp_func(func, tmp_termcount, true);
    term_t **function;
    printf("Segona funció per multiplicar la derivada (opcional)\n");
    function = build_func(&termcount, NULL, false);
    term_t **dx_func;
    if (function != NULL)
    {
        simp_func(function, termcount, true);
        dx_func = derivative(multiply_terms(func[0], function[0]), termcount, &dx_termcount);
    }
    else
    {
        dx_func = derivative(func, termcount, &dx_termcount);
    }
    printf("\nDerivada: \n");
    print_func(dx_func, dx_termcount);
}

bool eval_constant_expression(term_t **function, size_t termcount)
{
    size_t strallocsz = 512;
    char *fx;

    if (check_if_no_func(function, termcount))
    {
        fx = calloc(strallocsz, sizeof(char));
        print_func(function, termcount);
        double eval_point = eval_func(function, 0, termcount);
        print_num(eval_point, fx, strallocsz);
        printf(" = %s", fx);
        free(fx);
        return true;
    }
    return false;
}

term_t **menu(size_t *termcount, int argc, char *argv[])
{

    printf("Operació a realitzar:\n\n");
    printf("1. Cerca arrels\n2. Domini\n3. Derivada\n4. Integral definida\n\nOpció: ");
    char input[100];
    if (argc <= 1)
    {
        fgets(input, sizeof input, stdin);
        if (tolower(input[0]) == 'q')
        {
            printf("Sortint del programa...\n");
            exit(0);
        }
        if (input[0] < '1' || input[0] > '4')
        {
            printf("Opció inexistent!\n");
            return NULL;
        }
        // eat_stdin_newline();
    }

    term_t **function = build_func(termcount, NULL, true);

    if (!function)
    {
        fprintf(stderr, "Error: No s'ha pogut construïr la funció");
        exit(0);
    }

    puts("\n");
    switch (input[0])
    {
    case '1':
        if (eval_constant_expression(function, *termcount))
        {
            break;
        }
        simp_func(function, *termcount, true);
        root_program(function, *termcount);
        break;
    case '2':
        simp_func(function, *termcount, false);
        domain_program(function, *termcount);
        break;
    case '3':
        simp_func(function, *termcount, true);
        derivative_program(function, *termcount);
        break;
    case '4':
        simp_func(function, *termcount, true);
        integral_program(function, *termcount);
        break;
    default:
        fprintf(stderr, "Error: Selecció incorrecta.");
        exit(0);
    }

    return function;
}

int main(int argc, char *argv[])
{
    // inicialització

    signal(SIGINT, sigint_h);
    size_t termcount = 0;

    while (true)
    {
        version_info();
        menu(&termcount, argc, argv);
        printf("\n\nPrém Retorn per tornar al menú principal.");
        getchar();
#ifdef WIN32
        system("cls");
#else
        system("clear");
#endif

        printf("\n");
    }
    return 0;
}