#include <stdio.h>
#include <malloc.h>
#include <string.h>
#include <conio.h>
#include <locale.h>


//Функция считает биномиальный коэффициент, используя треугольник паскаля
long long binom(int n, int k){
    if(k > n) return 0;
//нам незачем хранить весь треугольник, т.к. для подсчета одного уровня,
//нам необходимо знать только значение предыдущего, поэтому, чтобы не тратить
//попусту память, мы заводим только два массива под текущий и предыдущий уровень треугольника:
    long long** pascals_triangle_lines = (long long**)malloc(sizeof(long long*)*2);
    pascals_triangle_lines[0] = (long long*)malloc(sizeof(long long)*1);
    pascals_triangle_lines[0][0] = 1;
    pascals_triangle_lines[1] = (long long*)malloc(sizeof(long long)*2);
    pascals_triangle_lines[1][0] = pascals_triangle_lines[1][1] = 1;

//индекс принимает значение 0 или 1 в зависимости от того, какой массив из двух считается текущим, а какой - предыдущим
    int index = 1;

//доходим до нужного нам уровня по треугольнику:
    for(int i = 2; i <= n; ++i){
        index = 1-index;
        free(pascals_triangle_lines[index]);
        pascals_triangle_lines[index] = (long long*)malloc(sizeof(long long)*(i+1));
        pascals_triangle_lines[index][0] = pascals_triangle_lines[index][i] = 1;
        for(int l = 1; l < i; ++l)
            pascals_triangle_lines[index][l] = pascals_triangle_lines[1-index][l] + pascals_triangle_lines[1-index][l-1];
    }

//очищаем память
    long long ans = pascals_triangle_lines[index][k];
    free(pascals_triangle_lines[0]);
    free(pascals_triangle_lines[1]);
    free(pascals_triangle_lines);
//возвращаем результат
    return ans;
}

//тривиальная функция, считающая факториал числа
long long factorial(int n){
    long long answer = 1;
    for(int i = 1; i <= n; ++i)
        answer *= i;
    return answer;
}

//функция считает конечную разность n-го порядка, принимая на вход массив точек
double finite_difference(int n, double* points){
    double answer = 0;
    int sign = 1;
    for(int i = 0; i <= n; ++i){
        answer += sign*(points[n-i]*binom(n,i));
        sign *= -1;
    }
    return answer;
}

//функция считает значение полинома в заданной точке:
double solve_interpolation(double* factors, int n, double q){
    double answer = factors[0];
    for(int i = 1; i < n; ++i){
        double arg = 1;
        for(int j = 0; j < i; ++j)
            arg *= q-j;
        answer += factors[i]*arg;
    }
    return answer;
}

//функция печатает таблицу конечных разностей:
void print_difference_table(int n, double* points, double x0, double h){
    printf("Таблица конечных разностей:\n");
    printf("X\tY\t");
    for(int i = 0; i < n; ++i)
        printf("d%dy\t", i+1);
    printf("\n");
    for(int i = 0; i <= n; ++i){
        printf("%.4f\t%.4f\t", x0+i*h, points[i]);
        for(int j = 0; j < n-i; ++j)
            printf("%.4f\t", finite_difference(j+1, points+i));
        printf("\n");
    }
}

int main(){
    setlocale(0,"russian");
    int n;
    double x0, h;
    double* points = NULL;
    int correct_input;

    do{
        printf("Введите степень интерполяционного полинома: ");
        scanf("%d", &n);

        printf("Введите x_0: ");
        scanf("%lf", &x0);

        printf("Введите шаг интерполяции: ");
        scanf("%lf", &h);

        if(points) free(points);
        points = (double*)malloc(sizeof(double)*(n+1));

        printf("Последовательно введите %d значений y_i:\n", n+1);
        for(int i = 0; i <= n; ++i)
            scanf("%lf", points + i);

        printf("Строим интерполяционный полином степени %d на следующих точках:\n", n);
        for(int i = 0; i <= n; ++i)
            printf("%.3lf ", points[i]);
        printf("\nНа интервале [%.3lf,%.3lf] с шагом %.3lf\n", x0, x0+n*h, h);
        printf("Подтвердите корректность введенных данных. Введите 0 для продолжения или 1 для ввода данных повторно: ");
        scanf("%d", &correct_input);

    }while(correct_input);

    double* factors = (double*)malloc(sizeof(double)*(n+1));

    print_difference_table(n, points, x0, h);
    printf("Интерполяционный полином степени %d на точках:\n", n);
    for(int i = 0; i <= n; ++i)
        printf("%.3lf ", points[i]);
    printf("\nНа интервале [%.3lf,%.3lf] с шагом %.3lf выглядит следующим образом:\nP_%d(x) = %.3lf ", x0, x0+n*h, h, n, points[0]);

    factors[0] = points[0];
    for(int i = 1; i <= n; ++i){
        char str[64];
        strcpy(str, "q");
        for(int j = 1; j < i; ++j){
            char buf[10];
            sprintf(buf, "(q-%d)", j);
            strcat(str, buf);
        }
        factors[i] = finite_difference(i, points)/factorial(i);
        printf("+ %.5lf%s ", factors[i], str);
    }
    printf("\nгде для удобства представления мы ввели новую переменную q = (x-%.3lf)/%.3lf\n", x0, h);

    int check = 0;
    do{
        printf("Желаете проверить результат интерполяции? (1 - да, 0 - нет): ");
        scanf("%d", &check);
        if(!check) break;
        double point;
        printf("Введите точку, в какой вы хотите вычислить значение полинома: ");
        scanf("%lf", &point);
        printf("Значение полинома в точке %.3lf равно %.3lf\n", point, solve_interpolation(factors, n+1, (point - x0)/h));
    }while(true);

    printf("Для выхода нажмите любую клавишу...");
    free(points);
    free(factors);
    getch();
}