/*  davidson.c -- Davidson's Method to Find an Extremum
    Written Oct 29, 2007 by Eric Olson for Math/CS 466/666

    For more information see Chapter 7.2 in Numerical Analysis and
    Scientific Computation by Jeffery J. Leader */

#include <stdio.h>
#include <math.h>
#include <bits/nan.h>

typedef double function(double x);

double findmax(function * f, function * df, double a, double b)
{
    double ba, fb, dfb, xa;
    double b0, b1, b2, b3;
    int i, findmax = 0;
    fb = f(b);
    dfb = df(b);
    for (i = 0; i < 100; i++) {
/*      printf("a=%g, b=%g\n",a,b); */
        if (fabs(a - b) < 1e-8 * (fabs(a) + fabs(b) + 1))
            return a;
        b0 = f(a);
        b1 = df(a);
        ba = b - a;
        {
            double t1 = ba * ba;
            double t5 = b1 * ba;
            double t7 = ba * dfb;
            b3 = 1 / t1 / ba * (2.0 * b0 + t5 - 2.0 * fb + t7);
            b2 = -(2.0 * t5 + 3.0 * b0 - 3.0 * fb + t7) / t1;
        }
/*      printf("a=%g, b0=%g, b1=%g, b2=%g, b3=%g\n", a, b0, b1, b2, b3); */
        {
            double t1 = b2 * b2;
            double t5 = sqrt(t1 - 3.0 * b3 * b1);
            if (!findmax) {
                if (b2 < 0) {   /* Minimum */
                    xa = -1 / b3 * (b2 - t5) / 3.0;
                } else {
                    xa = 1 / (-b2 - t5) * b1;
                }
            } else {
                if (b2 > 0) {   /* Maximum */
                    xa = -1 / b3 * (b2 + t5) / 3.0;
                } else {
                    xa = 1 / (-b2 + t5) * b1;
                }
            }
        }
        b = a;
        fb = b0;
        dfb = b1;
        a += xa;
    }
    printf("Failed to converge!\n");
    return NAN;
}

double f(double x)
{
    return cos(x);
}
double df(double x)
{
    return -sin(x);
}

int main()
{
    double x = findmax(f, df, 0.1, 3.8);
    printf("Minimum is at %20.16g.\n", x);
    return 0;
}