1 Commits

Author SHA1 Message Date
75dce20390 idk, when i wrote it, but it somehow works 2022-05-24 23:27:04 +03:00
6 changed files with 289 additions and 190 deletions

3
.gitignore vendored
View File

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

View File

@ -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.

View File

@ -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))

393
main.c
View File

@ -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:
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 *init(int n)
{
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
Example: https://shorturl.at/tBCPS
Warning: result is evaluated in `double *y` array */
double *div_diff_es(double *x, double *y, unsigned int n)
arr *resize(arr *a, int new_size)
{
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);
if (a->size == new_size)
return a;
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;
}
/*
Coeficients of simplified polynomial computation
*/
void simplify_polynomial(double *res, double *el_coef, double *x, unsigned int n)
int convert_addr(arr *a, int pos)
{
double *tmp_polynomial = (double *)malloc(sizeof(double) * n);
tmp_polynomial[0] = 1;
pos = pos % a->size;
if (pos < 0)
pos = a->size + pos;
for (int i = 0; i < n; i++)
if (el_coef[i])
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)
{
if (i > 0)
mult_by_root(tmp_polynomial, x[i - 1], i - 1);
for (int i = 0; i < a->size; i++)
printf("%f ", get(a, i));
printf("\n");
for (int j = 0; j <= i; j++)
res[j] += el_coef[i] * tmp_polynomial[j];
return;
}
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--)
res[j] = (j ? res[j - 1] : 0) - (root * res[j]);
arr *res = init(a->size - 1);
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 */
void print_newton_poly(double *f, double *x, unsigned int n)
int has_comb(int *arr, int n, int k)
{
printf("Newton polynomial form:\n");
for (int i = 0; i < n; i++)
{
if (f[i]) // If coefficient != 0
if (n == k)
return 0;
int pos = k - 1;
if (arr[pos] == n - 1)
{
/* Coefficient sign and sum symbol */
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("-");
if (k == 1)
return 0;
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
while ((pos > 0) && arr[pos] == n - 1)
{
if (x[j] > 0)
printf("*(x-%lf)", x[j]);
else
printf("*(x+%lf)", -x[j]);
pos--;
arr[pos]++;
}
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: ");
unsigned int n = 0;
scanf("%u", &n);
double res = 1;
for (int i = 0; i < n; i++)
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++)
{
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 (a->size == 1)
{
if (i > 0 && coef[i - 1])
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(" ");
return a->p[0];
}
}
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;
}
/*
Main
*/
int main()
int compose_denominator(arr *a, int pos)
{
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),
*y = (double *)malloc(sizeof(double) * n);
res = res * (get(a, pos) - get(a, i));
}
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;
}

View File

@ -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);
double *div_diff_es(double *x, double *y, unsigned int n);
/*
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);
int has_comb(int *arr, int n, int k);
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);
arr *compose_interpolation_polynomial(arr *xes, arr *ys);
#endif

View File

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