Compute local maxima using Newton's method

maxima.c

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

#define N 3

typedef struct {
    double x[N];
} vector;

void print_vector(vector v)
{
    putchar('(');
    int i;
    for (i = 0; i < N; i++)
        printf("%f%s", v.x[i], i == N - 1 ? ")\n" : ", ");
}

/* Return true if all of v's elements are 0. */
int zero(vector v)
{
    int i;
    for (i = 0; i < N; i++) {
        if (v.x[i] >= DBL_EPSILON)
            return 0;
    }
    return 1;
}

/* Central finite difference coefficients. */
const double coef[][9] = {
    { 1./280, -4./105,  1./5, -4./5,        0, 4./5, -1./5, 4./105, -1./280},
    {-1./560,  8./315, -1./5,  8./5, -205./72, 8./5, -1./5, 8./315, -1./560},
};

/* Compute a reasonable h for floating point precision. */
#define compute_h(x) (x) != 0 ? sqrt(fabs(x) * FLT_EPSILON) : FLT_EPSILON

/* Compute the nth derivatives of f() at v. */
vector df(int n, double (*f)(vector), vector v)
{
    vector result = {{0}};
    int i, d;
    for (d = 0; d < N; d++) {
        vector vh = v;
        double h = compute_h(v.x[d]);
        for (i = -4; i <= 4; i++) {
            vh.x[d] = v.x[d] + h * i;
            result.x[d] += coef[n - 1][i + 4] * f(vh);
        }
        result.x[d] /= pow(h, n);
    }
    return result;
}

/* Find the local minima/maxima via Newton's method and gradient descent. */
vector fmin(double (*f)(vector), vector v)
{
    while (!zero(df(1, f, v))) {
        vector d1 = df(1, f, v), d2 = df(2, f, v);
        int i;
        for (i = 0; i < N; i++)
            v.x[i] -= d1.x[i] / d2.x[i];
    }
    return v;
}

/* Example function (higher-order paraboloid). */
double f(vector v)
{
    return pow(v.x[0], 2) + pow(v.x[1], 2) + pow(v.x[2], 2);
}

int main()
{
    vector v = {{1, 1, 1}};
    print_vector(fmin(f, v));
    return 0;
}