Changed test generation function

Switch from Lagrange to newton interpolation polynomial

.gitignore

@@ -1,4 +1,5 @@
 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 for n dots.
+ANSI C program which composes polynomial of n - 1 degree that passes through n dots.

input.py

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

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