#!/usr/bin/python
"""Using the method of secants to find the zeroes of an N-d continuous function.
"""
from __future__ import division
import itertools
import numpy


def sec(f, x0, x1, e):
    y = f(x0), f(x1)

    while abs(y[1]) > e:
        x0, x1 = x1, x1 - y[1]*(x1 - x0)/(y[1] - y[0])
        y = y[1], f(x1)

    return x1

def sec_seq(f, x0, x1):
    y = f(x0), f(x1)

    while True:
        yield x1
        x0, x1 = x1, x1 - y[1]*(x1 - x0)/(y[1] - y[0])
        y = y[1], f(x1)

def secnd_seq(f, x):
    x = list(x)
    y = [f(xi) for xi in x]

    while True:
        yield x[-1], y[-1]

        div = y[-1] - y[0]
        if not div:
            return

        x.append(x[-1] - y[-1] * (x[-1] - x[0]) / div)
        y.append(f(x[-1]))
        x.pop(0)
        y.pop(0)
        
unit_circle = lambda (x, y): x**2 + y**2 - 1

def circle_intersection(x0, y0, r0, x1, y1, r1):
    def f(xv):
        x, y = xv
        return (abs((x - x0)**2 + (y - y0)**2 - r0**2) +
                abs((x - x1)**2 + (y - y1)**2 - r1**2))

    return f

def circle_intersection_L2(x0, y0, r0, x1, y1, r1):
    def f(xv):
        x, y = xv
        return (((x - x0)**2 + (y - y0)**2 - r0**2)**2 +
                ((x - x1)**2 + (y - y1)**2 - r1**2)**2)

    return f

def circle_intersection_L8(x0, y0, r0, x1, y1, r1):
    def f(xv):
        x, y = xv
        return max(abs((x - x0)**2 + (y - y0)**2 - r0**2)**2,
                   abs((x - x1)**2 + (y - y1)**2 - r1**2)**2)

    return f

def test_nd(d=3, n=1000, maxiter=100, eps=1e-15, f=unit_circle):
    ok = not_ok = 0

    for i in range(n):
        x = [numpy.random.random(2) * 10 - 5 for j in range(d)]
        items = list(itertools.islice(secnd_seq(f, x), maxiter))
        if abs(items[-1][-1]) < eps:
            print "ok:", x, items[-1][0], len(items)
            ok += 1
        else:
            print "not ok", x, items[-d:]
            not_ok += 1

    return ok, not_ok

if __name__ == '__main__':
    test_nd()