Switched from broken generated math to LaTex

Switched to "recursive" polynomial simplification

.gitignore

@@ -2,3 +2,4 @@ a.out
 main
 .vscode
 vgcore*
+output.txt

README.md

@@ -1,3 +1,89 @@
 # 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. It presents it in Newton interpolation polynomial and monic form.
+
+## Interface
+
+Application accepts as standart input decimal below 2147483647 `n` as number of dots, followed by n dots in format: `<x> (space) <y>` on each line, where `x` is an abscisse and `y` is an ordinate of single dot. Dot coordinates must fit [2.22507e-308;1.79769e+308] range by modulo.
+
+Result will be printed to standart output in the following format:
+
+Newton polynomial form:
+
+$$f_0 - f_1*(x-x_0) + ... + f_n(x-x_0)*(x-x_1)*...*(x-x_{n-1})$$
+
+Simplified coefficients array (starting from 0 upto n-1 power):
+
+$$a_0 a_1 ... a_{n-1} a_n$$
+
+Polynomial in monic form:
+
+$$a_0 - a_1*x + ... + a_{n-1}*x^(n-2) + a_n*x^(n-1)$$
+
+Where $f_i$ is a divided difference of $y_1,...,y_i$, $a_i$ are coefficients of resulting monic polynomial
+
+## Data structure
+
+- `n` is an `unsigned int` variable, that is used to input and store number of dots
+
+- `x` is a pointer to array of `n` `double`s, that is used to store abscisses of dots
+
+- `y` is a pointer to array of `n` `double`s, that is used to store ordinates of dots
+
+- `coefficients` is a pointer to array of `n` `double`s, that is used to store coefficients of monic interpolation polynomial
+
+- `i`, `j` are `int` variables, those are used in loops as iterators
+
+- `tmp_polynomial` is a pointer to array of `n` `double`s, that is used to store coefficients of polynomial, resulting during simplification of interpolation polynomial summands.
+
+## Example
+
+Build and run application:
+
+```bash
+gcc main.c
+./a.out
+```
+
+### Input/output
+
+For input n = 3 and the following dots
+
+```plain
+1 5
+2 3
+4 8
+```
+
+Output is
+
+```plain
+Newton polynomial form:
+5 - 2*(x-1) + 1.5*(x-1)*(x-2)
+Simplified coefficients array (starting from 0 upto n-1 power):
+10 -6.5 1.5
+Polynomial in standart form:
+10 - 6.5*x + 1.5*x^2
+```
+
+### Illustrations
+
+#### Example 1
+
+<img src="./img/console.png" />
+
+<img src="./img/wolfram.png" />
+
+<img src="./img/plot.png" />
+
+#### Example 2
+
+<img src="./img/console2.png" />
+
+or
+
+<img src="./img/console3.png" />
+
+<img src="./img/wolfram2.png" />
+
+<img src="./img/plot2.png" />

input.py

@@ -0,0 +1,26 @@
+import sys
+import math
+
+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, math.sin(i))

main.c

@@ -1,319 +1,216 @@
-#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://en.wikipedia.org/wiki/Newton_polynomial#Examples */
+double div_diff(double *y, double *x, unsigned i, unsigned d)
 {
-    arr *a = (arr *)malloc(sizeof(arr));
+  return (y[i] - y[i - 1]) / (x[i] - x[i - d]);
+}
 
-    a->size = n;
+/* Evaluates divided differences of n values - array of some kind of derivatives with big enough dx
+   Example: https://en.wikipedia.org/wiki/Newton_polynomial#Examples
+   Warning: result is evaluated in `double *y` array */
+double *div_diff_es(double *x, double *y, unsigned n)
+{
+  for (int i = 1; i < n; i++)        // first element remains unchanged
+    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);
 
-    a->p = (double *)malloc(sizeof(double) * n);
+  return y;
+}
+
+/*
+ Coeficients of simplified polynomial computation
+*/
+
+/* Simplifies Newton polynomial with `el_coef` array of divided differences,
+   and `x` as array of x coordinates of dots,
+   and `n` is number of elements of this sum */
+void simplify_polynomial(double *res, double *el_coef, double *x, unsigned n)
+{
+  double *tmp_polynomial // Temporary array for storage of coefficients of multiplication of (x-x_i) polynomial
+      = (double *)malloc(sizeof(double) * n);
+  for (int i = 1; i < n; i++)
+    tmp_polynomial[i] = 0;
+  tmp_polynomial[0] = 1; // Set polynomial to 1 to start multiplication with it
+
+  for (int i = 0; i < n; i++) // For each elemnt of sum
+  {
+    if (i > 0) // Start multiplication from second element of sum
+      mult_by_root(tmp_polynomial, x[i - 1], i - 1);
+
+    for (int j = 0; j <= i; j++)                // For each cumputed coefficient of i'th polynomial of sum
+      res[j] += el_coef[i] * tmp_polynomial[j]; // Add it, multiplied with divided difference, to sum
+  }
+
+  free(tmp_polynomial);
+}
+
+/* `res` is an array of coefficients of polynomial, which is multiplied with (x - `root`) polynomial.
+   `power` is the power of `res` polynomial */
+void mult_by_root(double *res, double root, unsigned power)
+{
+  for (int j = power + 1; j >= 0; j--)
+    res[j] = (j ? res[j - 1] : 0) - (root * res[j]); // coefficient is k_i-1 - root * k_i
+}
+
+/*
+ User Interface 
+*/
+
+/* Prints interpolation polynomial in Newton notation */
+void print_newton_poly(double *f, double *x, unsigned n)
+{
+  printf("Newton polynomial form:\n");
   for (int i = 0; i < n; i++)
-        set(a, i, 0);
-
-    return a;
-}
-
-arr *resize(arr *a, int new_size)
-{
-    if (a->size == new_size)
-        return a;
-
-    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)
-{
-    pos = pos % a->size;
-    if (pos < 0)
-        pos = a->size + pos;
-
-    return pos;
-}
-
-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 (f[i]) // If coefficient != 0
-            printf("%f ", get(a, i));
+    {
-        printf("\n");
+      /* Coefficient sign and sum symbol */
-
+      if (i > 0 && f[i - 1]) // If it's not the first summond
-        return;
+      {
+        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("Array of size %d:\n", a->size);
+      printf("%g", fabs(f[i])); // Print coefficient without sign
 
-    for (int i = 0; i < a->size; i++)
+      for (int j = 0; j < i; j++) // For each (x-xi) bracket
-        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)
-{
-    arr *res = init(a->size - 1);
-    for (int i = 0, pos = 0; i < a->size; i++)
       {
-        if (i == ex_pos)
+        if (x[j]) // If summond is not zero, print it
-            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
-*/
-
-int has_comb(int *arr, int n, int k)
-{
-    if (n == k)
-        return 0;
-
-    int pos = k - 1;
-
-    if (arr[pos] == n - 1)
         {
-        if (k == 1)
+          if (x[j] > 0)
-            return 0;
+            printf("*(x-%g)", x[j]);
-
+          else
-        while ((pos > 0) && arr[pos] == n - 1)
+            printf("*(x+%g)", -x[j]);
-        {
-            pos--;
-            arr[pos]++;
-        }
-
-        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("*x");
-
-    return 1;
-}
-
-int mult_by_index(arr *a, int *coords, int n)
-{
-    double res = 1;
-    for (int i = 0; i < n; i++)
-        res = res * get(a, coords[i]);
-
-    return res;
-}
-
-int sum_of_mult_of_n_combinations(arr *a, int n)
-{
-    if (n == 0)
-        return 1;
-
-    if (a->size == 1)
-    {
-        return a->p[0];
       }
 
-    double acc = 0;
+      printf(" ");
-
-    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)
-{
-    double res = 1;
-    for (int i = 0; i < a->size; i++)
-    {
-        if (i == pos)
-            continue;
-
-        res = res * (get(a, pos) - get(a, i));
     }
-    return res;
-}
-
-arr *compose_interpolation_polynomial(arr *xes, arr *ys)
-{
-    arr *res = init(xes->size);
-
-    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);
-
-        arr *xis = arr_without_el(xes, j);
-
-        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;
   }
 
-        res = add(res, jcoef);
+  printf("\n");
-
-        free_arr(xis);
-    }
-
-    free_arr(jcoef);
-
-    return res;
 }
 
-int main(int argc, char *argv[])
+/* Returns inputed by user number of dots */
+unsigned insert_n()
 {
-    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;
+  unsigned n = 0;
-    scanf("%d", &n);
+  scanf("%u", &n);
 
-    if (!quiet_mode)
+  return n;
+}
+
+/* Reads pairs of x'es and y'es of n dots to corresponding array */
+void insert_coords(double *xes, double *yes, unsigned n)
+{
   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);
+    xes[i] = x;
-        set(ys, i, y);
+    yes[i] = y;
   }
+}
 
-    if (!quiet_mode)
+/* Prints array of n doubles */
+void print_array(double *arr, unsigned n)
+{
+  printf("Simplified coefficients array (starting from 0 upto n-1 power):\n");
+
+  for (int i = 0; i < n; i++)
+    printf("%g ", arr[i]);
+
+  printf("\n");
+}
+
+/* Prints interpolation polynomial in standart form
+   e.g. a*x^2 + b*x + c */
+void print_poly(double *coef, unsigned n)
+{
+  printf("Polynomial in standart form:\n");
+
+  for (int i = 0; i < n; i++)
   {
-        printf("Inserted the following doths:\n");
+    if (coef[i])
-        printa(xes, 0);
+    {
-        printa(ys, 0);
+      if (i > 0 && coef[i - 1])
+        if (coef[i] > 0)
+          printf("+ ");
+        else
+          printf("- ");
+      else if (coef[i] < 0)
+        printf("-");
+
+      printf("%g", fabs(coef[i]));
+      if (i > 0)
+        printf("*x");
+      if (i > 1)
+        printf("^%d ", i);
+      else
+        printf(" ");
+    }
   }
 
-    arr *res = compose_interpolation_polynomial(xes, ys);
+  printf("\n");
+}
 
-    if (!quiet_mode)
+/* 
-        printf("Resulting polynomial will have such coeficients:\n");
+ Main
-    arr *reversed = reverse(res);
+*/
-    printa(reversed, quiet_mode);
 
-    free_arr(reversed);
+int main()
-    free_arr(res);
+{
-    free_arr(xes);
+  unsigned n = insert_n();
-    free_arr(ys);
+
+  double *x = (double *)malloc(sizeof(double) * n),
+         *y = (double *)malloc(sizeof(double) * n);
+
+  insert_coords(x, y, n);
+
+  double *f = div_diff_es(x, y, n);
+
+  print_newton_poly(f, x, n);
+
+  double *coefficients = (double *)malloc(sizeof(double) * n);
+  for (unsigned i = 0; i < n; i++)
+    coefficients[i] = 0;
+
+  simplify_polynomial(coefficients, f, x, n);
+
+  print_array(coefficients, n);
+
+  print_poly(coefficients, n);
+
+  free(x);
+  free(y);
+  free(coefficients);
 
   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 i, unsigned d);
-int mult_by_index(arr *a, int *coords, int n);
+double *div_diff_es(double *x, double *y, unsigned 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 insert_n();
+void print_newton_poly(double *f, double *x, unsigned n);
+void insert_coords(double *x, double *y, unsigned n);
+void print_array(double *arr, unsigned n);
+void print_poly(double *coef, unsigned n);
+
+/*
+ Coeficients of simplified polynomial computation
+*/
+
+void simplify_polynomial(double *res, double *el_coef, double *x, unsigned n);
+void mult_by_root(double *res, double root, unsigned 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