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 that passes through 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

@@ -16,16 +16,16 @@ double fabs(double 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 */
+   example: https://en.wikipedia.org/wiki/Newton_polynomial#Examples */
-double div_diff(double *y, double *x, unsigned int i, unsigned int d)
+double div_diff(double *y, double *x, unsigned i, unsigned d)
 {
   return (y[i] - y[i - 1]) / (x[i] - x[i - d]);
 }
 
 /* Evaluates divided differences of n values - array of some kind of derivatives with big enough dx
-   Example: https://shorturl.at/tBCPS
+   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 int n)
+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
@@ -38,66 +38,35 @@ double *div_diff_es(double *x, double *y, unsigned int n)
  Coeficients of simplified polynomial computation
 */
 
-void simplify_polynomial(double *res, double *rev_el_coef, double *x, unsigned int n)
+/* 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)
 {
-  for (int i = 0; i < n; i++)
+  double *tmp_polynomial // Temporary array for storage of coefficients of multiplication of (x-x_i) polynomial
-    if (rev_el_coef[i])
+      = (double *)malloc(sizeof(double) * n);
-      for (int j = 0; j <= i; j++)
+  for (int i = 1; i < n; i++)
-        res[i - j] += (j % 2 ? -1 : 1) * rev_el_coef[i] * compute_sum_of_multiplications_of_k(x, j, i);
+    tmp_polynomial[i] = 0;
-}
+  tmp_polynomial[0] = 1; // Set polynomial to 1 to start multiplication with it
 
-double compute_sum_of_multiplications_of_k(double *arr, unsigned int k, unsigned int n)
+  for (int i = 0; i < n; i++) // For each elemnt of sum
-{
-  if (k == 0)
-    return 1;
-
-  if (k == 1 && n == 1)
-    return arr[0];
-
-  unsigned int *selected = (unsigned int *)malloc(sizeof(unsigned int) * k); // Indexes of selected for multiplication elements
-
-  int i = 0, // index of `arr` array
-      j = 0; // index of `selected` array
-
-  double sum = 0;
-
-  while (j >= 0)
   {
-    if (i <= (n + (j - k)))
+    if (i > 0) // Start multiplication from second element of sum
-    {
+      mult_by_root(tmp_polynomial, x[i - 1], i - 1);
-      selected[j] = i;
 
-      if (j == k - 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
-        sum += mult_by_indexes(arr, selected, k);
-        i++;
-      }
-      else
-      {
-        i = selected[j] + 1;
-        j++;
-      }
-    }
-    else
-    {
-      j--;
-      if (j >= 0)
-        i = selected[j] + 1;
-    }
   }
 
-  free(selected);
+  free(tmp_polynomial);
-
-  return sum;
 }
 
-double mult_by_indexes(double *arr, unsigned int *indexes, unsigned int size)
+/* `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)
 {
-  double res = 1;
+  for (int j = power + 1; j >= 0; j--)
-  for (int i = 0; i < size; i++)
+    res[j] = (j ? res[j - 1] : 0) - (root * res[j]); // coefficient is k_i-1 - root * k_i
-    res *= arr[indexes[i]];
-
-  return res;
 }
 
 /*
@@ -105,7 +74,7 @@ double mult_by_indexes(double *arr, unsigned int *indexes, unsigned int size)
 */
 
 /* Prints interpolation polynomial in Newton notation */
-void print_newton_poly(double *f, double *x, unsigned int n)
+void print_newton_poly(double *f, double *x, unsigned n)
 {
   printf("Newton polynomial form:\n");
   for (int i = 0; i < n; i++)
@@ -123,16 +92,16 @@ void print_newton_poly(double *f, double *x, unsigned int n)
       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
+      printf("%g", 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
         {
           if (x[j] > 0)
-            printf("*(x-%lf)", x[j]);
+            printf("*(x-%g)", x[j]);
           else
-            printf("*(x+%lf)", -x[j]);
+            printf("*(x+%g)", -x[j]);
         }
         else
           printf("*x");
@@ -145,16 +114,18 @@ void print_newton_poly(double *f, double *x, unsigned int n)
   printf("\n");
 }
 
-unsigned int insert_n()
+/* Returns inputed by user number of dots */
+unsigned insert_n()
 {
   printf("Insert number of dots: ");
-  unsigned int n = 0;
+  unsigned n = 0;
   scanf("%u", &n);
 
   return n;
 }
 
-void insert_coords(double *xes, double *yes, unsigned int 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");
 
@@ -168,19 +139,22 @@ void insert_coords(double *xes, double *yes, unsigned int n)
   }
 }
 
-void print_array(double *arr, unsigned int n)
+/* 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("%lf ", arr[i]);
+    printf("%g ", arr[i]);
 
   printf("\n");
 }
 
-void print_poly(double *coef, unsigned int n)
+/* Prints interpolation polynomial in standart form
+   e.g. a*x^2 + b*x + c */
+void print_poly(double *coef, unsigned n)
 {
-  printf("Simplified polynom:\n");
+  printf("Polynomial in standart form:\n");
 
   for (int i = 0; i < n; i++)
   {
@@ -191,15 +165,16 @@ void print_poly(double *coef, unsigned int n)
           printf("+ ");
         else
           printf("- ");
-      else
+      else if (coef[i] < 0)
         printf("-");
 
-      printf("%lf", fabs(coef[i]));
+      printf("%g", fabs(coef[i]));
       if (i > 0)
         printf("*x");
       if (i > 1)
         printf("^%d ", i);
-      else printf(" ");
+      else
+        printf(" ");
     }
   }
 
@@ -224,6 +199,8 @@ int main()
   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);
 
@@ -231,5 +208,9 @@ int main()
 
   print_poly(coefficients, n);
 
+  free(x);
+  free(y);
+  free(coefficients);
+
   return 0;
 }

polynominal_interpolation.h

@@ -13,25 +13,24 @@ double fabs(double x);
  Business logic
 */
 
-double div_diff(double *y, double *x, unsigned int i, unsigned int d);
+double div_diff(double *y, double *x, unsigned i, unsigned d);
-double *div_diff_es(double *x, double *y, unsigned int n);
+double *div_diff_es(double *x, double *y, unsigned n);
 
 /*
  User interface
 */
 
-unsigned int insert_n();
+unsigned insert_n();
-void print_newton_poly(double *f, double *x, unsigned int n);
+void print_newton_poly(double *f, double *x, unsigned n);
-void insert_coords(double *x, double *y, unsigned int n);
+void insert_coords(double *x, double *y, unsigned n);
-void print_array(double *arr, unsigned int n);
+void print_array(double *arr, unsigned n);
-void print_poly(double *coef, unsigned int n);
+void print_poly(double *coef, unsigned n);
 
 /*
  Coeficients of simplified polynomial computation
 */
 
-void simplify_polynomial(double *res, double *rev_el_coef, double *x, unsigned int n);
+void simplify_polynomial(double *res, double *el_coef, double *x, unsigned n);
-double compute_sum_of_multiplications_of_k(double *x, unsigned int k, unsigned int n);
+void mult_by_root(double *res, double root, unsigned step);
-double mult_by_indexes(double *arr, unsigned int *indexes, unsigned int size);
 
 #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