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

namespace Science.Statistics.BasicStatistics
{
    public class NormalCurve
    {
        public NormalCurve()
        {
        }
        private double func(double x)
        {
            return 100.0 / Math.Sqrt(2.0 * Math.PI) * Math.Exp(- x * x / 2.0);
        }
        public double ValueAt(double x)
        {
            return func(x);
        }
        public double AreaBetween(double from, double to)
        {
            Function.ToLastType<double, double> del 
                = new Function.ToLastType<double,double>(func);
            IntegrationMidpoint obj = new IntegrationMidpoint(del, from, to);
            obj.Compute();
            return obj.Result;
        }
        public double ToStandardUnit(double symmetricalArea)
        {
            area = symmetricalArea;
            Function.ToLastType<double, double> ff 
                = new Function.ToLastType<double, double>(whatIsSU);
            return zbrent(ff, 0.0, 10.0, 0.0000000001);
        }

        private double area;
        private double whatIsSU(double x)
        {
            return (area - AreaBetween(-x, x));
        }
        private double zbrent(Function.ToLastType<double, double> func, double x1, double x2, double tol)
        {
            int ITMAX = 100;
            double EPS = 0.00000000001;
            double a = x1, b = x2, c = x2, d = 0.0, e = 0.0, fa = func(a), fb = func(b), fc, p, q, r, s, tol1, xm;
            
            fc = fb;
            for (int iter = 0; iter < ITMAX; iter++)
            {
                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;
                }
                tol1 = 2.0 * EPS * Math.Abs(b) + 0.5 * tol;
                xm = 0.5 * (c - b);
                if (Math.Abs(xm) <= tol1 || fb == 0.0) return b;
                if (Math.Abs(e) >= tol1 && Math.Abs(fa) > Math.Abs(fb))
                {
                    s = fb / fa;
                    if (a == c)
                    {
                        p = 2.0 * xm * s;
                        q = 1.0 - s;
                    }
                    else
                    {
                        q = fa / fc;
                        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(tol1 * 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) > tol1)
                    b += d;
                else
                    b += SIGN(tol1, xm);
                fb = func(b);
            }
            return 0.0;
        }
        private double SIGN(double a, double b) { return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); }
	}


    class Integration
    {
        private Function.ToLastType<double, double> del;
        private double from1, to1, result;
        public Integration()
        {
        }
        public Integration(Function.ToLastType<double, double> f)
        {
            del = f;
        }
        public Integration(Function.ToLastType<double, double> f, double from, double to)
        {
            del = f;
            from1 = from;
            to1 = to;
        }
        public void Compute()
        {
            Qromb qr = new Qromb();
            result = qr.qromb(del, from1, to1, 10e-15);
        }

        public Function.ToLastType<double, double> Function
        {
            set { del = value; }
        }
        public double From
        {
            set { from1 = value; }
        }
        public double To
        {
            set { to1 = value; }
        }
        public double Result
        {
            get { return result; }
        }
    }

    class Qromb
    {
        public Qromb()
        {
        }
        public double qromb(Function.ToLastType<double, double> func, double a, double b, double eps)
        {
            int JMAX = 20, JMAXP = JMAX + 1, K = 5;
            double[] s = new double[JMAX];
            double[] h = new double[JMAXP];
            Poly_interp polint = new Poly_interp(h, s, K);
            h[0] = 1.0;
            Trapzd t = new Trapzd(func, a, b);
            for (int j = 1; j <= JMAX; j++)
            {
                s[j - 1] = t.next();
                if (j >= K)
                {
                    double ss = polint.rawinterp(j - K, 0.0);
                    if (Math.Abs(polint.dy) <= eps * Math.Abs(ss)) return ss;
                }
                h[j] = 0.25 * h[j - 1];
            }
            try { throw new Science.Error(); }
            catch (Science.Error e)
            {
                e.Write("Too many steps in routine qromb");
            }
            return 0.0;
        }
    }

    class Base_interp
    {
        public int n, mm, jsav, cor, dj;
        public double[] xx, yy;
        public double[] Y
        {
            set { yy = value; }
        }
        public Base_interp(double[] x, double[] y, int m)
        {
            n = x.Length;
            mm = m;
            jsav = 0;
            cor = 0;
            xx = x;
            yy = y;
            dj = Math.Min(1, (int)Math.Pow((double)n, 0.25));
        }

        public virtual double interp(double x)
        {
            int jlo = (cor != 0) ? hunt(x) : locate(x);
            return rawinterp(jlo, x);
        }

        public int locate(double x)
        {
            int ju, jm, jl;
            if (n < 2 || mm < 2 || mm > n) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("locate size error");
                }
            bool ascnd = (xx[n - 1] >= xx[0]);
            jl = 0;
            ju = n - 1;
            while (ju - jl > 1)
            {
                jm = (ju + jl) >> 1;
                if (x >= xx[jm] == ascnd)
                    jl = jm;
                else
                    ju = jm;
            }
            cor = Math.Abs(jl - jsav) > dj ? 0 : 1;
            jsav = jl;
            return Math.Max(0, Math.Min(n - mm, jl - ((mm - 2) >> 1)));
        }
        public int hunt(double x)
        {
            int jl = jsav, jm, ju, inc = 1;
            if (n < 2 || mm < 2 || mm > n) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("hunt size error");
                }
            bool ascnd = (xx[n - 1] >= xx[0]);
            if (jl < 0 || jl > n - 1)
            {
                jl = 0;
                ju = n - 1;
            }
            else
            {
                if (x >= xx[jl] == ascnd)
                {
                    for (; ; )
                    {
                        ju = jl + inc;
                        if (ju >= n - 1) { ju = n - 1; break; }
                        else if (x < xx[ju] == ascnd) break;
                        else
                        {
                            jl = ju;
                            inc += inc;
                        }
                    }
                }
                else
                {
                    ju = jl;
                    for (; ; )
                    {
                        jl = jl - inc;
                        if (jl <= 0) { jl = 0; break; }
                        else if (x >= xx[jl] == ascnd) break;
                        else
                        {
                            ju = jl;
                            inc += inc;
                        }
                    }
                }
            }
            while (ju - jl > 1)
            {
                jm = (ju + jl) >> 1;
                if (x >= xx[jm] == ascnd)
                    jl = jm;
                else
                    ju = jm;
            }
            cor = Math.Abs(jl - jsav) > dj ? 0 : 1;
            jsav = jl;
            return Math.Max(0, Math.Min(n - mm, jl - ((mm - 2) >> 1)));
        }

        public virtual double rawinterp(int jlo, double x)
        {
            return 0.0;
        }
    }

    class Poly_interp : Base_interp
    {
        public double dy;
        public Poly_interp(double[] xv, double[] yv, int m)
            : base(xv, yv, m)
        {
            dy = 0.0;
        }
        public override double rawinterp(int jl, double x)
        {
            int i, m, ns = 0;
            double y, den, dif, dift, ho, hp, w;
            double[] xa = xx;
            double[] ya = yy;
            double[] c = new double[mm];
            double[] d = new double[mm];
            dif = Math.Abs(x - xa[jl + 0]);
            for (i = 0; i < mm; i++)
            {
                if ((dift = Math.Abs(x - xa[jl + i])) < dif)
                {
                    ns = i;
                    dif = dift;
                }
                c[i] = ya[jl + i];
                d[i] = ya[jl + i];
            }
            y = ya[ns--];
            for (m = 1; m < mm; m++)
            {
                for (i = 0; i < mm - m; i++)
                {
                    ho = xa[jl + i] - x;
                    hp = xa[jl + i + m] - x;
                    w = c[i + 1] - d[i];
                    if ((den = ho - hp) == 0.0) try { throw new Science.Error(); }
                        catch (Science.Error e)
                        {
                            e.Write("Poly_interp error");
                        }
                    den = w / den;
                    d[i] = hp * den;
                    c[i] = ho * den;
                }
                y += (dy = (2 * (ns + 1) < (mm - m) ? c[ns + 1] : d[ns--]));
            }
            return y;
        }
    }

    class Quadrature
    {
        public Quadrature()
        {
        }
        public int n;
        public virtual double next()
        {
            return 0.0;
        }
    }

    class Trapzd : Quadrature
    {
        double a, b, s;
        Function.ToLastType<double, double> func;
        public Trapzd(Function.ToLastType<double, double> funcc, double aa, double bb)
        {
            func = funcc;
            a = aa;
            b = bb;
            n = 0;
        }
        public override double next()
        {
            double x, tnm, sum, del;
            int it, j;
            n++;
            if (n == 1)
            {
                return (s = 0.5 * (b - a) * (func(a) + func(b)));
            }
            else
            {
                for (it = 1, j = 1; j < n - 1; j++) it <<= 1;
                tnm = it;
                del = (b - a) / tnm;
                x = a + 0.5 * del;
                for (sum = 0.0, j = 0; j < it; j++, x += del) sum += func(x);
                s = 0.5 * (s + (b - a) * sum / tnm);
                return s;
            }
        }
    }

    class IntegrationMidpoint
    {
        private Function.ToLastType<double, double> del;
        private double from1, to1, result;
        public IntegrationMidpoint()
        {
        }
        public IntegrationMidpoint(Function.ToLastType<double, double> f)
        {
            del = f;
        }
        public IntegrationMidpoint(Function.ToLastType<double, double> f, double from, double to)
        {
            del = f;
            from1 = from;
            to1 = to;
        }
        public void Compute()
        {
            Midpnt obj = new Midpnt(del, from1, to1);
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            obj.next();
            result = obj.next();
        }

        public Function.ToLastType<double, double> Function
        {
            set { del = value; }
        }
        public double From
        {
            set { from1 = value; }
        }
        public double To
        {
            set { to1 = value; }
        }
        public double Result
        {
            get { return result; }
        }
    }

    class Midpnt : Quadrature
    {
        protected double a, b, s = 0.0;
        protected Function.ToLastType<double, double> funk;
        public virtual double func(double x)
        {
            return funk(x);
        }
        public Midpnt(Function.ToLastType<double, double> funcc, double aa, double bb)
        {
            funk = funcc;
            a = aa;
            b = bb;
            n = 0;
        }
        public override double next()
        {
            int it, j;
            double x, tnm, sum, del, ddel;
            n++;
            if (n == 1)
            {
                return (s = (b - a) * func(0.5 * (a + b)));
            }
            else
            {
                for (it = 1, j = 1; j < n - 1; j++) it *= 3;
                tnm = it;
                del = (b - a) / (3.0 * tnm);
                ddel = del + del;
                x = a + 0.5 * del;
                sum = 0.0;
                for (j = 0; j < it; j++)
                {
                    sum += func(x);
                    x += ddel;
                    sum += func(x);
                    x += del;
                }
                s = (s + (b - a) * sum / tnm) / 3.0;
                return s;
            }
        }
    }
}