On various websites, I've come across a lot of realizations of Newton polynomial interpolation that use separate functions for forward and backward interpolation, respectively. Such solutions seem ugly to me because the functions look pretty much the same.
For a single Newton interpolation function, I treat the backward interpolation case as forward interpolation with the data points written in reverse order. Usually the data points come already sorted in ascending order by x
values, so the reverse order means simply going through the data points from right to left over an interval.
With all these in mind, my Newton interpolation function in C is:
/*
Error codes:
-2 - malloc() error
-1 - range check error
0 - no error (success)
*/
// The necessary typedefs
typedef struct {
int points;
double *xList, *yList;
} DataTable;
typedef struct {
double *coeff;
int degree;
double *xList;
} NewtonPoly;
// I've excluded the code for checking if the data table is filled properly.
int getNewtonPoly(DataTable *table, int nodes, NewtonPoly *poly)
{
int n = abs(nodes);
int points = table->points;
if (n > points || n == 0)
return -1;
double *result = (double *)malloc(n * sizeof(double));
if (result == NULL)
return -2;
double *tmp = (double *)malloc(n * sizeof(double));
if (tmp == NULL)
return -2;
double *xList = (double *)malloc(n * sizeof(double));
if (xList == NULL)
return -2;
if (nodes > 0)
for (int i = 0; i < n; ++i) {
xList[i] = table->xList[i];
tmp[i] = table->yList[i];
}
else
for (int i = points - 1, j = 0; i >= points - n; --i, ++j) {
xList[j] = table->xList[i];
tmp[j] = table->yList[i];
}
result[0] = tmp[0];
for (int i = 1; i < n; ++i) {
for (int j = 0; j <= n - i - 1; ++j)
tmp[j] = (tmp[j + 1] - tmp[j]) / (xList[j + i] - xList[j]);
result[i] = tmp[0];
}
poly->degree = n - 1;
poly->coeff = result;
poly->xList = xList;
free(tmp);
return 0;
}
Are there some better approaches to this problem?
1 Answer 1
Don't cast the result from malloc()
Instead, include <stdlib.h>
so that the compiler knows it returns a void*
. Also, it's a good practice to use the object name as argument to sizeof
, so its type can be changed in one place (e.g. if you wanted to use float
or long double
instead). It's idiomatic to test a pointer isn't null by using the !
operator on it.
double *result = malloc(sizeof result * n);
if (!result)
return -2;
Be careful about memory leaks
I don't see a corresponding free()
of result
or xList
- we should be clear that this is the caller's responsibility. It's also easily missed that the function won't free an existing poly->coeff
or poly->xList
before overwriting - I think that's a big magnet for inviting leaks.