idk, when i wrote it, but it somehow works

.gitignore

@@ -1,5 +1,4 @@
 a.out
 main
 .vscode
-vgcore*
+vgcore*
-output.txt

README.md

@@ -1,3 +1,3 @@
 # Polynomial Interpolation
 
-ANSI C program which composes polynomial of n - 1 degree that passes through n dots.
+ANSI C program which composes polynomial of n - 1 degree for n dots.

input.py

@@ -1,25 +0,0 @@
-import sys
-
-try:
-  n = int(sys.argv[1])
-except:
-  n = 5
-
-print(n)
-
-def f(x: int) -> int:
-  """
-  f(x) = sum with i from 0 to n-1 (i+1)*x^i
-
-  E.g. f(x) = 5x^4 + 4x^3 + 3x^2 + 2x + 1
-  """
-
-  res: int = 0
-
-  for i in range(n):
-    res += (i+1) * pow(x, i)
-
-  return res
-
-for i in range(n):
-  print(i, f(i))

main.c

@@ -1,196 +1,319 @@
+#include <stdio.h>
+#include <stdlib.h>
+
 #include "./polynominal_interpolation.h"
 
-/* Utils */
+/*
+ Utils
+*/
 
-double fabs(double x)
+int min(int a, int b)
 {
-  return x > 0 ? x : -x;
+    return (a + b - abs(a - b)) / 2;
+}
+
+int max(int a, int b)
+{
+    return (a + b + abs(a - b)) / 2;
 }
 
 /*
- Newton interpolation polynomial
+ Array utils
 */
 
-/* Divided difference is evaluated for:
+arr *init(int n)
-   array y stands for f(x)
-   array x stands for x
-   number i stands for index of evaluated difference (from 0)
-   number d stands for order of difference (from 0)
-   example: https://shorturl.at/tBCPS */
-double div_diff(double *y, double *x, unsigned int i, unsigned int d)
 {
-  return (y[i] - y[i - 1]) / (x[i] - x[i - d]);
+    arr *a = (arr *)malloc(sizeof(arr));
+
+    a->size = n;
+
+    a->p = (double *)malloc(sizeof(double) * n);
+    for (int i = 0; i < n; i++)
+        set(a, i, 0);
+
+    return a;
 }
 
-/* Evaluates divided differences of n values - array of some kind of derivatives with big enough dx
+arr *resize(arr *a, int new_size)
-   Example: https://shorturl.at/tBCPS
-   Warning: result is evaluated in `double *y` array */
-double *div_diff_es(double *x, double *y, unsigned int n)
 {
-  for (int i = 1; i < n; i++)        // first element remains unchanged
+    if (a->size == new_size)
-    for (int j = n - 1; j >= i; j--) // evaluate from the end of array, decreacing number of step every repeation
+        return a;
-      y[j] = div_diff(y, x, j, i);
 
-  return y;
+    double *new_p = (double *)malloc(sizeof(double) * new_size);
+    for (int i = 0; i < min(new_size, a->size); i++)
+        new_p[i] = get(a, i);
+
+    free(a->p);
+
+    for (int i = a->size; i < new_size; i++)
+        new_p = 0;
+
+    a->p = new_p;
+    a->size = new_size;
+
+    return a;
 }
 
-/*
+int convert_addr(arr *a, int pos)
- Coeficients of simplified polynomial computation
-*/
-
-void simplify_polynomial(double *res, double *el_coef, double *x, unsigned int n)
 {
-  double *tmp_polynomial = (double *)malloc(sizeof(double) * n);
+    pos = pos % a->size;
-  tmp_polynomial[0] = 1;
+    if (pos < 0)
+        pos = a->size + pos;
 
-  for (int i = 0; i < n; i++)
+    return pos;
-    if (el_coef[i])
+}
+
+double get(arr *a, int pos)
+{
+    pos = convert_addr(a, pos);
+
+    return a->p[pos];
+}
+
+void set(arr *a, int pos, double val)
+{
+    pos = convert_addr(a, pos);
+
+    a->p[pos] = val;
+
+    // printa(a, 1);
+}
+
+arr *add(arr *a, arr *b)
+{
+    for (int i = 0; i < a->size; i++)
+        set(a, i, a->p[i] + b->p[i]);
+
+    return a;
+}
+
+arr *mult(arr *a, double mul)
+{
+    arr *res = init(a->size);
+
+    for (int i = 0; i < a->size; i++)
+        set(res, i, a->p[i] * mul);
+
+    return res;
+}
+
+void printa(arr *a, int q)
+{
+    if (q)
     {
-      if (i > 0)
+        for (int i = 0; i < a->size; i++)
-        mult_by_root(tmp_polynomial, x[i - 1], i - 1);
+            printf("%f ", get(a, i));
+        printf("\n");
 
-      for (int j = 0; j <= i; j++)
+        return;
-        res[j] += el_coef[i] * tmp_polynomial[j];
     }
+
+    printf("Array of size %d:\n", a->size);
+
+    for (int i = 0; i < a->size; i++)
+        printf("%5d ", i + 1);
+    printf("\n");
+
+    for (int i = 0; i < a->size; i++)
+        printf("%5.2f ", get(a, i));
+    printf("\n");
 }
 
-void mult_by_root(double *res, double root, unsigned int step)
+arr *arr_without_el(arr *a, int ex_pos)
 {
-  for (int j = step + 1; j >= 0; j--)
+    arr *res = init(a->size - 1);
-    res[j] = (j ? res[j - 1] : 0) - (root * res[j]);
+    for (int i = 0, pos = 0; i < a->size; i++)
+    {
+        if (i == ex_pos)
+            continue;
+        set(res, pos, a->p[i]);
+        pos++;
+    }
+
+    return res;
+}
+
+arr *reverse(arr *a)
+{
+    arr *res = init(a->size);
+    for (int i = 0; i < a->size; i++)
+        set(res, i, a->p[a->size - 1 - i]);
+
+    return res;
+}
+
+void free_arr(arr *a)
+{
+    free(a->p);
+    free(a);
 }
 
 /*
- User Interface 
+ Business logic
 */
 
-/* Prints interpolation polynomial in Newton notation */
+int has_comb(int *arr, int n, int k)
-void print_newton_poly(double *f, double *x, unsigned int n)
 {
-  printf("Newton polynomial form:\n");
+    if (n == k)
-  for (int i = 0; i < n; i++)
+        return 0;
-  {
+
-    if (f[i]) // If coefficient != 0
+    int pos = k - 1;
+
+    if (arr[pos] == n - 1)
     {
-      /* Coefficient sign and sum symbol */
+        if (k == 1)
-      if (i > 0 && f[i - 1]) // If it's not the first summond
+            return 0;
-      {
-        if (f[i] > 0)
-          printf("+ ");
-        else
-          printf("- ");
-      }
-      else if (f[i] < 0) // If it is the first summond and coefficient is below zero
-        printf("-");
 
-      printf("%lf", fabs(f[i])); // Print coefficient without sign
+        while ((pos > 0) && arr[pos] == n - 1)
-
-      for (int j = 0; j < i; j++) // For each (x-xi) bracket
-      {
-        if (x[j]) // If summond is not zero, print it
         {
-          if (x[j] > 0)
+            pos--;
-            printf("*(x-%lf)", x[j]);
+            arr[pos]++;
-          else
-            printf("*(x+%lf)", -x[j]);
         }
-        else
-          printf("*x");
-      }
 
-      printf(" ");
+        for (int i = pos + 1; i < k; i++)
+            arr[i] = arr[i - 1] + 1;
+
+        if (arr[0] > n - k)
+            return 0;
     }
-  }
+    else
+        arr[pos]++;
 
-  printf("\n");
+    return 1;
 }
 
-unsigned int insert_n()
+int mult_by_index(arr *a, int *coords, int n)
 {
-  printf("Insert number of dots: ");
+    double res = 1;
-  unsigned int n = 0;
+    for (int i = 0; i < n; i++)
-  scanf("%u", &n);
+        res = res * get(a, coords[i]);
 
-  return n;
+    return res;
 }
 
-void insert_coords(double *xes, double *yes, unsigned int n)
+int sum_of_mult_of_n_combinations(arr *a, int n)
 {
-  printf("Insert dots coordinates in the following format:\n<x> (space) <y>\nEach dot on new line\n");
+    if (n == 0)
+        return 1;
 
-  for (int i = 0; i < n; i++)
+    if (a->size == 1)
-  {
-    double x, y;
-    scanf("%lf %lf", &x, &y);
-
-    xes[i] = x;
-    yes[i] = y;
-  }
-}
-
-void print_array(double *arr, unsigned int n)
-{
-  printf("Simplified coefficients array (starting from 0 upto n-1 power):\n");
-
-  for (int i = 0; i < n; i++)
-    printf("%lf ", arr[i]);
-
-  printf("\n");
-}
-
-void print_poly(double *coef, unsigned int n)
-{
-  printf("Simplified polynom:\n");
-
-  for (int i = 0; i < n; i++)
-  {
-    if (coef[i])
     {
-      if (i > 0 && coef[i - 1])
+        return a->p[0];
-        if (coef[i] > 0)
-          printf("+ ");
-        else
-          printf("- ");
-      else
-        printf("-");
-
-      printf("%lf", fabs(coef[i]));
-      if (i > 0)
-        printf("*x");
-      if (i > 1)
-        printf("^%d ", i);
-      else
-        printf(" ");
     }
-  }
 
-  printf("\n");
+    double acc = 0;
+
+    int coords[n];
+    for (int i = 0; i < n; i++)
+        coords[i] = i;
+
+    acc += mult_by_index(a, coords, n);
+    while (has_comb(coords, a->size, n))
+        acc += mult_by_index(a, coords, n);
+
+    return acc;
 }
 
-/* 
+int compose_denominator(arr *a, int pos)
- Main
-*/
-
-int main()
 {
-  unsigned n = insert_n();
+    double res = 1;
+    for (int i = 0; i < a->size; i++)
+    {
+        if (i == pos)
+            continue;
 
-  double *x = (double *)malloc(sizeof(double) * n),
+        res = res * (get(a, pos) - get(a, i));
-         *y = (double *)malloc(sizeof(double) * n);
+    }
+    return res;
+}
 
-  insert_coords(x, y, n);
+arr *compose_interpolation_polynomial(arr *xes, arr *ys)
+{
+    arr *res = init(xes->size);
 
-  double *f = div_diff_es(x, y, n);
+    arr *jcoef = init(xes->size);
+    for (int j = 0; j < xes->size; j++)
+    {
+        int minus = (!(xes->size % 2) ? -1 : 1);
+        double denominator = compose_denominator(xes, j);
+        double multiplicator = get(ys, j);
 
-  print_newton_poly(f, x, n);
+        arr *xis = arr_without_el(xes, j);
 
-  double *coefficients = (double *)malloc(sizeof(double) * n);
+        for (int i = 0; i < xes->size; i++)
+        {
+            double k_sum = sum_of_mult_of_n_combinations(xis, xes->size - 1 - i);
+            set(jcoef, i, minus * (multiplicator * k_sum) / denominator);
+            minus = -minus;
+        }
 
-  simplify_polynomial(coefficients, f, x, n);
+        res = add(res, jcoef);
 
-  print_array(coefficients, n);
+        free_arr(xis);
+    }
 
-  print_poly(coefficients, n);
+    free_arr(jcoef);
 
-  return 0;
+    return res;
+}
+
+int main(int argc, char *argv[])
+{
+    int quiet_mode = 0;
+    if (argc > 1 && argv[1][0] == '-' && argv[1][1] == 'q')
+        quiet_mode = 1;
+
+    if (!quiet_mode)
+        printf("Insert number of dots: ");
+    int n = 6;
+    scanf("%d", &n);
+
+    if (!quiet_mode)
+        printf("Insert dots coordinates in the following format:\n<x> (space) <y>\nEach dot on new line\n");
+
+    arr *xes = init(n);
+    arr *ys = init(n);
+
+    // set(xes, 0, 1);
+    // set(ys, 0, 1);
+    // set(xes, 1, 2);
+    // set(ys, 1, 2);
+    // set(xes, 2, 3);
+    // set(ys, 2, 3);
+    // set(xes, 3, 4);
+    // set(ys, 3, 4);
+    // set(xes, 4, 5);
+    // set(ys, 4, 5);
+    // set(xes, 5, 6);
+    // set(ys, 5, 6);
+
+    for (int i = 0; i < n; i++)
+    {
+        double x, y;
+        scanf("%lf %lf", &x, &y);
+
+        set(xes, i, x);
+        set(ys, i, y);
+    }
+
+    if (!quiet_mode)
+    {
+        printf("Inserted the following doths:\n");
+        printa(xes, 0);
+        printa(ys, 0);
+    }
+
+    arr *res = compose_interpolation_polynomial(xes, ys);
+
+    if (!quiet_mode)
+        printf("Resulting polynomial will have such coeficients:\n");
+    arr *reversed = reverse(res);
+    printa(reversed, quiet_mode);
+
+    free_arr(reversed);
+    free_arr(res);
+    free_arr(xes);
+    free_arr(ys);
+
+    return 0;
 }

polynominal_interpolation.h

@@ -1,36 +1,42 @@
 #ifndef POLYNOMIAL_INTERPOLATION_H
 #define POLYNOMIAL_INTERPOLATION_H
 
-#include <stdio.h>
-#include <stdlib.h>
-
 /*
  Utils
 */
-double fabs(double x);
+
+int min(int a, int b);
+int max(int a, int b);
+
+/*
+ Array utils
+*/
+
+typedef struct
+{
+    int size;
+    double *p;
+} arr;
+
+arr *init(int n);
+arr *resize(arr *a, int new_size);
+int convert_addr(arr *a, int pos);
+double get(arr *a, int pos);
+void set(arr *a, int pos, double val);
+arr *add(arr *a, arr *b);
+arr *mult(arr *a, double mul);
+void printa(arr *a, int q);
+arr *arr_without_el(arr *a, int ex_pos);
+arr *reverse(arr *a);
 
 /*
  Business logic
 */
 
-double div_diff(double *y, double *x, unsigned int i, unsigned int d);
+int has_comb(int *arr, int n, int k);
-double *div_diff_es(double *x, double *y, unsigned int n);
+int mult_by_index(arr *a, int *coords, int n);
-
+int sum_of_mult_of_n_combinations(arr *a, int n);
-/*
+int compose_denominator(arr *a, int pos);
- User interface
+arr *compose_interpolation_polynomial(arr *xes, arr *ys);
-*/
-
-unsigned int insert_n();
-void print_newton_poly(double *f, double *x, unsigned int n);
-void insert_coords(double *x, double *y, unsigned int n);
-void print_array(double *arr, unsigned int n);
-void print_poly(double *coef, unsigned int n);
-
-/*
- Coeficients of simplified polynomial computation
-*/
-
-void simplify_polynomial(double *res, double *el_coef, double *x, unsigned int n);
-void mult_by_root(double *res, double root, unsigned int step);
 
 #endif

test.sh

@@ -1,4 +0,0 @@
-#!/bin/sh
-
-gcc main.c
-python input.py $1 | tee /dev/fd/2 | ./a.out | tee output.txt