#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define OK 0
#define ALLOC_ERROR 1
#define ERROR 2
#define NOT_FOUND -1
#define EPS 1e-5
#define INTERPOLATION 0
#define INVERSE_INTERPOLATION 1

/* Дисциплина: Вычислительные алгоритмы, ИУ7-43Б */
/* Интерполя́ция, интерполи́рование — в вычислительной математике способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений */
/* Задание: 1 - Построить интерполяционный полином Ньютона с помощью таблицы; 2 - Найти корень табличной функции методом ПОЛОВИННОГО ДЕЛЕНИЯ и методом ОБРАТНОЙ ИНТЕРПОЛЯЦИИ */

typedef struct
{                   //Структура для хранения таблицы. Таблица хранится в виде матрицы, строки
    double **mtr;   //которой считываются из файла, а число столбцов равно (n + 2), где n - это
    int row;        //степень полинома
    int start;
    int end;
} table;

void menu()
{
    printf("1 - Calculate Newton interpolation polynomial\n2 - Find the radix of a table function using the dichotomy method\n3 - Find the radix of a table function by inverse interpolation\n0 - Exit\n**Warning! the column of x should be sorted in increasing or descending order**\n");
}

int fill_matrix(FILE *f, table *tbl, int n, int indicate)
{
    fseek(f, 0, 0);
    fscanf(f, "%d", &tbl->row);

    if (!(tbl->mtr = (double**)malloc(tbl->row * sizeof(double*))))
        return ALLOC_ERROR;
    for (int i = 0; i < tbl->row; i++)
        if (!(tbl->mtr[i] = (double*)malloc((n + 2) * sizeof(double))))
            return ALLOC_ERROR;

    for (int i = 0; i < tbl->row; i++)
        for (int j = 0; j < n + 2; j++)
            tbl->mtr[i][j] = 0;

    if (indicate == INTERPOLATION)
        for (int i = 0; i < tbl->row + 1; i++)
            for (int j = 0; j < 2; j++)
                fscanf(f, "%lf", &tbl->mtr[i][j]);
    else
        for (int i = 0; i < tbl->row + 1; i++)
            for (int j = 1; j >= 0; j--)
                fscanf(f, "%lf", &tbl->mtr[i][j]);

    return OK;
}

void print_matrix(table *tbl, int n, int m)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < m; j++)
        {
            printf("%lf ", tbl->mtr[i][j]);
        }
        printf("\n");
    }
}

int find_insert(table *tbl, double x)
{
    for (int i = 0; i < tbl->row - 1; i++)
    {
        if ((x >= tbl->mtr[i][0] && x <= tbl->mtr[i + 1][0]) || (x <= tbl->mtr[i][0] && x >= tbl->mtr[i + 1][0]))
        {
            tbl->start = i;
            tbl->end = i + 1;
            printf("%lf %lf\n", tbl->mtr[i][0], tbl->mtr[i + 1][0]);
            return OK;
        }
    }
    return ERROR; //Extrapolation
}

int find_interval(table *tbl, int n)
{
    int h = n - 1;

    while ((int)(h / 2) > 0)
    {
        if (tbl->end == tbl->row - 1 || tbl->start == 0)
            break;
        tbl->end++;
        tbl->start--;
        if  (tbl->start < 0 && tbl->end == tbl->row)
        {
            tbl->start++;
            tbl->end--;
            return ERROR;
        }
        else
        {
            if (tbl->end == tbl->row)
            {
                tbl->end--;
                h++;
                break;
            }
            if (tbl->start < 0)
            {
                tbl->start++;
                h++;
                break;
            }
        }
        h -= 2;
    }
    if (h != 0)
    {
        printf("\nAsymmetric\n");
        if (tbl->start == 0)
        {
            tbl->end += h;
            if (tbl->end >= tbl->row)
                return ERROR;
        }
        else if (tbl->end == tbl->row - 1)
        {
            tbl->start -= h;
            if (tbl->start <= 0)
                return ERROR;
        }
        else
        {
            if (tbl->start >= h)
                tbl->start -= h;
            else if (tbl->row - tbl->end - 1 >= h)
                tbl->end += h;
            else
                return ERROR;
        }
    }
    return OK;
}

int find_difference(table *tbl, int n)
{
    int start = tbl->start, end = tbl->end;
    int interval = n + 1;
    int j = 0, q = 1;

    for (int k = 2; k < n + 2; k++)
    {
        j = 0;
        for (int i = tbl->start; i < end; i++)
        {
            tbl->mtr[j][k] = (tbl->mtr[start][k - 1] - tbl->mtr[start + 1][k - 1]) / (tbl->mtr[i][0] - tbl->mtr[i + q][0]);
            j++;
            start++;
        }
        start = 0;
        end--;
        q++;
    }

    print_matrix(tbl, tbl->row, n + 2);

    return OK;
}

double polynom(table tbl, int n, double x)
{
    double p = tbl.mtr[tbl.start][1];
    double mul = 1;
    int start = tbl.start;

    for (int i = 2; i < n + 2; i++)
    {
        mul *= x - tbl.mtr[start][0];
        start++;
        p += mul * tbl.mtr[0][i];
    }
    return p;
}

int dichotomy(table tbl, int n)
{
    int flag = NOT_FOUND;
    double a = tbl.mtr[tbl.start][0], b = tbl.mtr[tbl.end][0], c = 0;

    if (fabs(polynom(tbl, n, a)) < EPS)
    {
        printf("\nX = %lf\n", a);
        return OK;
    }
    if (fabs(polynom(tbl, n, b)) < EPS)
    {
        printf("\nX = %lf\n", b);
        return OK;
    }
    if (b < a)
    {
        b = tbl.mtr[tbl.start][0];
        a = tbl.mtr[tbl.end][0];
    }
    c = (a + b) / 2;
    while (fabs(b - a) >= fabs(c) * EPS + EPS)
    {
        if (polynom(tbl, n, b) * polynom(tbl, n, c) < 0)
        {
            flag = 1;
            a = c;
        }
        else
            b = c;
        c = (a + b) / 2;
    }
    if (flag == NOT_FOUND)
        return NOT_FOUND;

    printf("\nX = %lf\n", c);

    return OK;
}

void del_tbl(table *tbl)
{
    if (tbl->mtr)
    {
        for (int i = 0; i < tbl->row; i++)
            if (tbl->mtr[i])
                free(tbl->mtr[i]);
        free(tbl->mtr);
    }
}

int main()
{
    FILE *f = fopen("table_x.txt", "r");
    if (!f)
    {
        printf("\nFile was not found\n");
        return -1;
    }
    int tmp;
    double rc;

    fscanf(f, "%d", &tmp);
    printf("\n");
    for (int i = 0; i < tmp; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            fscanf(f, "%lf", &rc);
            printf("%lf ", rc);
        }
        printf("\n");
    }
    printf("\n");

    table tbl;
    int n = 0;
    double x = 0;

    int ch;
    menu();
    printf("Input your choice: ");
    scanf("%d", &ch);

    while (ch != 0)
    {
        switch (ch)
        {
            case 1:
                printf("Input the power of polynom (n): ");
                scanf("%d", &n);
                printf("Input the value (x): ");
                scanf("%lf", &x);

                fill_matrix(f, &tbl, n, INTERPOLATION);
                if (find_insert(&tbl, x) != OK)
                {
                    printf("Extrapolation! (the input valuse is out of table)\n");
                    return ERROR;
                }
                printf("\n%d %d\n", tbl.start, tbl.end);

                if (find_interval(&tbl, n) != OK)
                {
                    printf("\nCannot take an interval: the table is too small\n");
                    return ERROR;
                }
                printf("\n%d %d\n", tbl.start, tbl.end);

                find_difference(&tbl, n);
                double p = polynom(tbl, n, x);
                printf("\nValue = %lf\n", p);

                del_tbl(&tbl);
                break;

            case 2:
                printf("Input the power of polynom (n): ");
                scanf("%d", &n);
                printf("Input the value (x): ");
                scanf("%lf", &x);

                fill_matrix(f, &tbl, n, INTERPOLATION);

                if (find_insert(&tbl, x) != OK)
                {
                    printf("Extrapolation! (the input valuse is out of table)\n");
                    return ERROR;
                }
                printf("\n%d %d\n", tbl.start, tbl.end);

                if (find_interval(&tbl, n) != OK)
                {
                    printf("\nCannot take an interval: the table is too small\n");
                    return ERROR;
                }
                printf("\n%d %d\n", tbl.start, tbl.end);
                find_difference(&tbl, n);
                printf("-------------\n| DICHOTOMY |\n-------------\n");
                if (dichotomy(tbl, n) != OK)
                    printf("The radix was not found (the interval does not contain it)");

                del_tbl(&tbl);
                break;

            case 3:
                printf("Input the power of polynom (n): ");
                scanf("%d", &n);

                fill_matrix(f, &tbl, n, INVERSE_INTERPOLATION);

                if (find_insert(&tbl, 0) != OK)
                {
                    printf("Extrapolation! (the input valuse is out of table)\n");
                    return ERROR;
                }
                //printf("\n%d %d\n", tbl.start, tbl.end);
                if (find_interval(&tbl, n) != OK)
                {
                    printf("\nCannot take an interval: the table is too small\n");
                    return ERROR;
                }
                printf("\n%d %d\n", tbl.start, tbl.end);
                find_difference(&tbl, n);
                printf("-------------------------\n| INVERSE_INTERPOLATION |\n-------------------------\n");
                double radix = polynom(tbl, n, 0);
                printf("\nX = %lf\n", radix);

                del_tbl(&tbl);
                break;

            case 0:

                del_tbl(&tbl);
                break;

            default:

                del_tbl(&tbl);
                printf("\nWrong input\n");
                break;
        }
        printf("\nInput your choice: ");
        scanf("%d", &ch);
    }

    fclose(f);
    return 0;
}