#include <stdio.h>
#include <math.h>
#include <float.h>

int newton (double *x, double tolf, double tolx, int maxit,
            void (*fdf)(double,double*,double*,void*), void *prm, int ixrr) {
    double x0 = *x;
    for (int i = 0; i < maxit; i++) {
        double fx = 0.0, dfx = 0.0;
        fdf(x0, &fx, &dfx, prm);

        if (ixrr) fprintf(stderr, "newton: it=%d x=%.16g f(x)=%.16g f'(x)=%.16g\n", i, x0, fx, dfx);

        if (fabs(fx) <= tolf) {
            *x = x0;
            return i + 1;
        }

        *x = x0 - (fx / dfx);
        double dx = fabs(*x - x0);

        if (ixrr) fprintf(stderr, "newton: it=%d x=%.16g |dx|=%.16g\n", i, *x, dx);

        if (dx <= tolx) return i + 1;
        x0 = *x;
    }
    return -1;
}

double neville (int n, double xi[], double yi[], double x) {
    double y[n + 1];
    for (int i = 0; i <= n; i++) {
        y[i] = yi[i];
    }
    for (int i = 1; i <= n; i++) {
        for (int j = 0; j <= n - i; j++) {
            y[j] = ((x - xi[j]) * y[j + 1] - (x - xi[j + i]) * y[j]) / (xi[j + i] - xi[j]);
        }
    }
    double result = y[0];
    return result;
}