Нахождение корней уравнений методом половинного деления

Развлечения ради, временами появляется желание покрутить задачки, которые выходят за рамки повседневного программирования. В разделе SICP 1.3.3 есть интересная, но несложная задачка. Более того, решение там приведено, но прежде чем смотреть его я попытался решить задачу самостоятельно и потом сравнить результаты.

Метод половинного деления (half-interval method) — это простой, но мощный способ нахождения корней уравнения f(x)=0, где f — непрерывная функция. Идея состоит в том, что если нам даны такие точки a и b, что f(a)<0<f(b), то функция f должна иметь по крайней мере один ноль на отрезке между a и b. Чтобы найти его, возьмем x, равное среднему между a и b, и вычислим f(x). Если f(x)>0, то f должна иметь ноль на отрезке между a и x. Если f(x)<0, то f должна иметь ноль на отрезке между x и b. Продолжая таким образом, мы сможем находить все более узкие интервалы, на которых f должна иметь ноль. Когда мы дойдем до точки, где этот интервал достаточно мал, процесс останавливается. Поскольку интервал неопределенности уменьшается вдвое на каждом шаге процесса, число требуемых шагов растет как Θ(log(L/T)), где L есть длина исходного интервала, а T есть допуск ошибки (то есть размер интервала, который мы считаем «достаточно малым»).

Касательно реализации. Она достаточно упрощенная(напр. я не осуществляю проверку того, что a меньше b):

(ns sicp.clojure
  (:use [clojure.test :only [is]])
  (:use [clojure.contrib.generic.math-functions :only [sin]])
  (:use [clojure.contrib.math :only [abs]]))

(defn get-middle [left right]
    (/ (+ left right) 2))

(defn get-adge-point [f middle x]
  (let [y1 (f middle) y2 (f x)]
    (cond
      (and (> y1 0) (> y2 0))middle
      (and (< y1 0) (< y2 0))middle
      :else x)))

(defn good-enough? [left right]
  (<  (- right left) 0.0001))

(defn half-interval-method-iter [f left middle right]
  (if
    (good-enough? left right) (get-middle left right)
    (let [left (get-adge-point f middle left)
          middle (get-middle left right)
          right (get-adge-point f middle right)]
      (recur f left middle right))))

(defn half-interval-method [f a b]
  (double (half-interval-method-iter f a (get-middle a b) b)))

;; sin(x)=0
(is (= 3.141571044921875 (half-interval-method sin 2 4)))

;; x^3 − 2x − 3 = 0
(is (= 1.893280029296875 (half-interval-method (fn [x] (- (* x x x) (* 2 x) 3)) 1 2)))

Для получения очередного(половинного) интервала следуем правилу:

 f(x) не меняет свой знак по отношению к f(a), новый интервал -> [x:b]
 f(x) не меняет свой знак по отношению к f(b), новый интервал -> [a:x]  

Функция get-adge-point реализует этот алгоритм. x - среднее значение для a и b, которое вычисляется при помощи get-middle. good-enough? - вспомогательная функция для проверки достигли мы заданной точности при вычислении результата. Для сравнения реализация на java:

public class HalfMidleMethod {

    public static interface Fn {
        double f(double x);
    }

    private boolean isGoodEnough(double left, double right) {
        return (right - left) < 0.0001;
    }

    private double getAdgePoint(Fn y, double middle, double x) {
        double y1 = y.f(middle);
        double y2 = y.f(x);
        if (y1 > 0 && y2 > 0) {
            return middle;
        } else if (y1 < 0 && y2 < 0) {
            return middle;
        } else {
            return x;
        }
    }

    private double getMiddle(double left, double right) {
        return (left + right) / 2;
    }

    private double getRoot(Fn y, double left, double middle, double right) {
        //System.out.println("left=" + left + ", middle=" + middle + ", right=" + right);
        if (isGoodEnough(left, right)) {
            return getMiddle(left, right);
        }

        left = getAdgePoint(y, middle, left);
        right = getAdgePoint(y, middle, right);
        middle = getMiddle(left, right);
        return getRoot(y, left, middle, right);
    }

    public double getRoot(Fn y, double left, double right) {
        return getRoot(y, left, getMiddle(left, right), right);
    }

    public static void main(String[] arvg) {
        Fn sin = new Fn() {
            @Override
            public double f(double x) {
                return Math.sin(x);
            }
        };

        System.out.println("Pi=" + new HalfMidleMethod().getRoot(sin, 4, 7));
    }
}

Т.к. java целиком и полностью объектно ориентированный язык, то как следствие в нем нет функций высших порядков. Для передачи функции в качестве аргумента используется интерфейс HalfMidleMethod.Fn.