Subversion Repositories mkgmap

/*
 * Copyright (C) 2009 Christian Gawron
 * 
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License version 2 as
 *  published by the Free Software Foundation.
 * 
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 * 
 * 
 * Author: Christian Gawron
 * Create date: 03-Jul-2009
 */
package uk.me.parabola.mkgmap.reader.dem;

/**
 * Find zero of a function using Brent's method.
 */
public class Brent {
    
        public interface Function {
                public double eval(double x);
        }

    private static final double epsilon = 3.0e-10;

    public static double zero(Function f, double x1, double x2) {
                return zero(f, x1, x2, 1e-8, 100);
    }

        public static double zero(Function f, double x1, double x2, double tol, int maxit) {
                double a=x1, b=x2, c=x2;
                double fa=f.eval(a), fb=f.eval(b);

                if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
                        throw new ArithmeticException("Root must be bracketed");

                double fc = fb;
                double d = 0;
                double e = 0;
                for (int iterations=0; iterations < maxit; iterations++) {
                        if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
                                c=a;
                                fc=fa;
                                e=d=b-a;
                        }
                        if (Math.abs(fc) < Math.abs(fb)) {
                                a=b;
                                b=c;
                                c=a;
                                fa=fb;
                                fb=fc;
                                fc=fa;
                        }
                        double tolerance = 2.0 * epsilon * Math.abs(b) + 0.5 * tol;
                        double xm = 0.5 * (c - b);
                        if (Math.abs(xm) <= tolerance || fb == 0.0) return b;
                        if (Math.abs(e) >= tolerance && Math.abs(fa) > Math.abs(fb)) {
                                double s = fb / fa;
                                double p;
                                double q;
                                if (a == c) {
                                        p=2.0*xm*s;
                                        q=1.0-s;
                                } else {
                                        q=fa/fc;
                                        double r = fb / fc;
                                        p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
                                        q=(q-1.0)*(r-1.0)*(s-1.0);
                                }
                                if (p > 0.0) q = -q;
                                p=Math.abs(p);
                                double min1 = 3.0 * xm * q - Math.abs(tolerance * q);
                                double min2 = Math.abs(e * q);
                                if (2.0*p < (min1 < min2 ? min1 : min2)) {
                                        e=d;
                                        d=p/q;
                                } else {
                                        d=xm;
                                        e=d;
                                }
                        } else {
                                d=xm;
                                e=d;
                        }
                        a=b;
                        fa=fb;
                        if (Math.abs(d) > tolerance)
                                b += d;
                        else
                                b += xm >= 0 ? tolerance : -tolerance;
                        fb=f.eval(b);
                }
                throw new ArithmeticException("Maximum number of iterations exceeded");
        }
}