// C# で解く「なぜ関数プログラミングは重要か」
// H27.1/15 SUZUKI Hisao
// cf. http://www.sampou.org/haskell/article/whyfp.html

using System;
using System.Linq;
using System.Collections;
using System.Collections.Generic;

// 4. プログラムの貼り合せ
static class WhyFp4
{
    // 4.1 ニュートン-ラプソン法による平方根
    static IEnumerable<T> repeat<T>(Func<T, T> f, T a) {
        for (;;) {
            yield return a;
            a = f(a);
        }
    }

    // seq から pred(A, B) が成り立つ隣り合った項 A, B を見つける。
    static Tuple<T, T> FindPair<T>(this IEnumerable<T> seq,
                                 Func<T, T, bool> pred) {
        using (var e = seq.GetEnumerator()) {
            if (! e.MoveNext())
                return null;
            T a = e.Current;
            if (! e.MoveNext())
                return null;
            T b = e.Current;
            for (;;) {
                if (pred(a, b))
                    return new Tuple<T, T>(a, b);
                a = b;
                if (! e.MoveNext())
                    return null;
                b = e.Current;
            }
        }
    }

    static double within(double eps, IEnumerable<double> seq) {
        return seq.FindPair((a, b)=> Math.Abs(a - b) <= eps).Item2;
    }

    static double sqrt_v1(double a0, double eps, double N) {
        return within(eps, repeat((x)=> (x + N/x) / 2.0, a0));
    }

    static double relative(double eps, IEnumerable<double> seq) {
        var pair= seq.FindPair((a, b)=> Math.Abs(a - b) <= eps * Math.Abs(b));
        return pair.Item2;
    }

    static double sqrt_v2(double a0, double eps, double N) {
        return relative(eps, repeat((x)=> (x + N/x) / 2.0, a0));
    }

    // 4.2 数値微分
    static double easydiff(Func<double, double> f, double x, double h) {
        double diff = (f(x + h) - f(x)) / h;
        Console.WriteLine(" -- h={0}, diff={1}", h, diff); // DEBUG
        return diff;
    }

    static IEnumerable<double>
    differentiate(double h0, Func<double, double> f, double x) {
        return repeat(h => h / 2.0, h0).Select(h => easydiff(f, x, h));
    }

    // seq の各要素 a(i) を f(a(i), a(i+1)) に換えた列を返す。
    static IEnumerable<R> SelectWithNeighbor<T, R>(this IEnumerable<T> seq,
                                                   Func<T, T, R> f) {
        using (var e = seq.GetEnumerator()) {
            if (! e.MoveNext())
                yield break;
            T a = e.Current;
            if (! e.MoveNext())
                yield break;
            T b = e.Current;
            for (;;) {
                yield return f(a, b);
                a = b;
                if (! e.MoveNext())
                    yield break;
                b = e.Current;
            }
        }
    }

    static IEnumerable<double> elimerror(double n, IEnumerable<double> seq) {
        double p = Math.Pow(2, n);
        if (Double.IsNaN(p) || p == 1) {
            Console.WriteLine(" -- p={0}", p);
            return seq;
        } else {
            return seq.SelectWithNeighbor((double a, double b) => {
                    double r =  (b * p - a) / (p - 1);
                    Console.WriteLine(" -- r={0}, b={1}, n={2}", r, b, n);
                    return r;
                });
        }
    }


    // 一度評価した値をキャッシュする列 (≒ lazy list)
    class Cached<T>: IEnumerable<T> {
        List<T> cache = new List<T>();
        IEnumerator<T> generator;
        int index = -1;         // generator で生成した値の添え字位置

        public Cached(IEnumerable<T> sequence) {
            generator = sequence.GetEnumerator();
        }

        IEnumerator IEnumerable.GetEnumerator() {
            return GetEnumerator();
        }

        // 注意: 要素の評価時に this を lock する。
        public IEnumerator<T> GetEnumerator() {
            for (int i = 0;; i++) {
                T e;
                lock (this) {
                    if (i <= index) {
                        e = cache[i];
                    } else {
                        if (! generator.MoveNext())
                            break;
                        e = generator.Current;
                        cache.Add(e);
                        index++;
                    }
                }
                yield return e;
            }
        }
    }


    static double order(IEnumerable<double> seq) {
        var a = new double[3];
        int i = 0;
        foreach (double e in seq) {
            a[i] = e;
            if (++i == 3)
                break;
        }
        double x = (a[0] - a[2])/(a[1] - a[2]) - 1;
        double n = Math.Round(Math.Log(x, 2.0));
        return n;
    }

    static IEnumerable<double> improve(IEnumerable<double> seq) {
        var s = new Cached<double>(seq);
        return elimerror(order(s), s);
    }

    static IEnumerable<double> super(IEnumerable<double> seq) {
        return repeat(s => new Cached<double>(improve(s)),
                      new Cached<double>(seq))
            .Select(s => s.Skip(1).First());
    }

    // 4.3 数値積分
    static double easyintegrate(Func<double, double> f, double a, double b) {
        double integ = (f(a) + f(b)) * (b - a) / 2;
        Console.WriteLine(" -- a={0}, b={1}, integ={2}", a, b, integ);
        return integ;
    }

    static IEnumerable<double> integrate_v1
    (Func<double, double> f, double a, double b) {
        yield return easyintegrate(f, a, b);
        double mid = (a + b) / 2;
        var xx = integrate_v1(f, a, mid);
        var yy = integrate_v1(f, mid, b);
        foreach (double e in xx.Zip(yy, (x, y)=> x + y)) {
            yield return e;
        }
    }

    static IEnumerable<double> integrate
    (Func<double, double> f, double a, double b) {
        return integ(f, a, b, f(a), f(b));
    }

    static IEnumerable<double> integ
    (Func<double, double> f, double a, double b, double fa, double fb) {
        yield return (fa + fb) * (b - a) / 2;
        double m = (a + b) / 2;
        double fm = f(m);
        var xx = integ(f, a, m, fa, fm);
        var yy = integ(f, m, b, fm, fb);
        foreach (double e in xx.Zip(yy, (x, y)=> x + y)) {
            yield return e;
        }
    }

    static void Main() {
        Console.WriteLine(sqrt_v1(1, 0.1, 0.3)); // => 0.555769230769231
        Console.WriteLine(sqrt_v2(1, 0.1, 0.3)); // => 0.547780809156242

        Console.WriteLine
            (within(10, differentiate(10, x => x*x*x, 100)));
        // => 30005.8597564697

        Console.WriteLine
            (within(10, improve(differentiate(10, x => x*x*x, 100))));
        // => 29996.875

        Console.WriteLine
            (within(10, improve(improve(differentiate(10, x => x*x*x, 100)))));
        // => 30000

       Console.WriteLine
           (within
            (10, improve(improve(improve(differentiate(10, x=>x*x*x, 100))))));
        // => 30000

       Console.WriteLine
           (within(10, super(differentiate(10, x => x*x*x, 100))));
       // => 30000

       Console.WriteLine
           (string.Join(", ", integrate_v1(x => x*x*x, 1, 2).Take(4)));
       // => 4.5, 3.9375, 3.796875, 3.76171875

       Console.WriteLine
           (string.Join(", ", integrate(x => {
                       Console.WriteLine(" -- x={0}", x);
                       return x*x*x;
                   }, 1, 2).Take(4)));
       // => 4.5, 3.9375, 3.796875, 3.76171875

       Console.WriteLine
           (within(0.005, integrate(x => {
                       Console.WriteLine(" -- x={0}", x);
                       return x*x*x;
                   }, 1, 2)));
       // => 3.750732421875

       Console.WriteLine
           (within(0.005, improve(integrate(x => {
                           Console.WriteLine(" -- x={0}", x);
                           return x*x*x;
                       }, 1, 2))));
       // => 3.75

       Console.WriteLine
           (string.Join(", ", improve(integrate(x => {
                           Console.WriteLine(" -- x={0}", x);
                           return 1.0 / (1 + x * x);
                       }, 0, 1)).Take(2)));
       // => 0.783333333333333, 0.785392156862745
    }
}