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

namespace StatisticalDistribution.SpecialFunctions
{
    class Functions
    {
        #region Erf1
        public static double Erf1(double x)
        {
            if (double.IsNegativeInfinity(x)) { return -1.0; }
            if (double.IsPositiveInfinity(x)) { return 1.0; }
            if (x < 0.0) { return -Functions.GammaRegularized(0.5, x * x); }
            return Functions.GammaRegularized(0.5, x * x);
        } 

        public static double Erf1Inverse(double x)
        {
            //if(x < -1.0 || x > 1.0)
            //{
            //    throw new ArgumentOutOfRangeException("x", x, Properties.LocalStrings.ArgumentInIntervalXYInclusive(-1, 1));
            //}

            x = 0.5 * (x + 1.0);

            // Define break-points.
            const double Plow = 0.02425;
            const double Phigh = 1 - Plow;

            double q;

            // Rational approximation for lower region:
            if (x < Plow)
            {
                q = Math.Sqrt(-2 * Math.Log(x));
                return ((((((((((ErfInvC[0] * q) + ErfInvC[1]) * q) + ErfInvC[2]) * q) + ErfInvC[3]) * q) + ErfInvC[4]) * q) + ErfInvC[5]) /
                    ((((((((ErfInvD[0] * q) + ErfInvD[1]) * q) + ErfInvD[2]) * q) + ErfInvD[3]) * q) + 1)
                    * Constants.Sqrt1_2;
            }

            // Rational approximation for upper region:
            if (Phigh < x)
            {
                q = Math.Sqrt(-2 * Math.Log(1 - x));
                return -((((((((((ErfInvC[0] * q) + ErfInvC[1]) * q) + ErfInvC[2]) * q) + ErfInvC[3]) * q) + ErfInvC[4]) * q) + ErfInvC[5]) /
                    ((((((((ErfInvD[0] * q) + ErfInvD[1]) * q) + ErfInvD[2]) * q) + ErfInvD[3]) * q) + 1)
                    * Constants.Sqrt1_2;
            }

            // Rational approximation for central region:
            q = x - 0.5;
            double r = q * q;
            return ((((((((((ErfInvA[0] * r) + ErfInvA[1]) * r) + ErfInvA[2]) * r) + ErfInvA[3]) * r) + ErfInvA[4]) * r) + ErfInvA[5]) * q /
                ((((((((((ErfInvB[0] * r) + ErfInvB[1]) * r) + ErfInvB[2]) * r) + ErfInvB[3]) * r) + ErfInvB[4]) * r) + 1)
                * Constants.Sqrt1_2;
        }
        private static readonly double[] ErfInvA = {
            -3.969683028665376e+01, 2.209460984245205e+02,
            -2.759285104469687e+02, 1.383577518672690e+02,
            -3.066479806614716e+01, 2.506628277459239e+00
            };

        private static readonly double[] ErfInvB = {
            -5.447609879822406e+01, 1.615858368580409e+02,
            -1.556989798598866e+02, 6.680131188771972e+01,
            -1.328068155288572e+01
            };

        private static readonly double[] ErfInvC = {
            -7.784894002430293e-03, -3.223964580411365e-01,
            -2.400758277161838e+00, -2.549732539343734e+00,
            4.374664141464968e+00, 2.938163982698783e+00
            };

        private static readonly double[] ErfInvD = {
            7.784695709041462e-03, 3.224671290700398e-01,
            2.445134137142996e+00, 3.754408661907416e+00
            };
        #endregion

        #region Erf2
        // ERF-2a (accurate in a neighborhood of 0 and a neighborhood of infinity, and the error is less than 0.00035 for all x)
        private static double erf2a_a = 0.140012;
        public static double Erf2(double x)
        {
            double xSq = Math.Pow(x, 2.0);
            return sgn(x) * Math.Sqrt(1.0 - Math.Exp(-xSq * (4 / Math.PI + erf2a_a * xSq) / (1 + erf2a_a * xSq)));
        }
        public static double Erf2Inverse(double y)
        {
            double ln_1_ySq = Math.Log(1 - Math.Pow(y, 2.0));
            return sgn(y) * Math.Sqrt(
                Math.Sqrt(Math.Pow(2 / Math.PI / erf2a_a + ln_1_ySq / 2, 2.0) - ln_1_ySq / erf2a_a)
                - 2 / Math.PI / erf2a_a - ln_1_ySq / 2);
        }
        private static double sgn(double x) { if (x < 0) { return -1.0; } if (x > 0) { return 1.0; } return 0.0; } 
        #endregion

        #region Erf3
        // ERF-2b (Using the alternate value a ≈ 0.147 reduces the maximum error to about 0.00012)
        private static double erf2b_a = 0.147;
        public static double Erf3(double x)
        {
            double xSq = Math.Pow(x, 2.0);
            return sgn(x) * Math.Sqrt(1.0 - Math.Exp(-xSq * (4 / Math.PI + erf2b_a * xSq) / (1 + erf2b_a * xSq)));
        }
        public static double Erf3Inverse(double y)
        {
            double ln_1_ySq = Math.Log(1 - Math.Pow(y, 2.0));
            return sgn(y) * Math.Sqrt(
                Math.Sqrt(Math.Pow(2 / Math.PI / erf2b_a + ln_1_ySq / 2, 2.0) - ln_1_ySq / erf2b_a)
                - 2 / Math.PI / erf2b_a - ln_1_ySq / 2);
        } 
        #endregion



        #region GammaLn
        public static double GammaLn(double x)
        {
            if (x < -34.0)
            {
                double y = -x;
                double decimalPart = y - Math.Floor(y);
                double z = y * Math.Sin(Math.PI * decimalPart);
                return Constants.LnPi - Math.Log(z) - GammaLnLargePositive(y);
            }
            if (x < 13) { return GammaLnSmall(x); }
            return GammaLnLargePositive(x);
        }

        static double GammaLnLargePositive(double x)
        {
            // Stirling
            double q = ((x - 0.5) * Math.Log(x)) - x + Constants.Ln2Pi_2;
            if (x > 100000000)
            {
                return q;
            }

            double p = 1 / (x * x);
            if (x >= 1000.0)
            {
                double a = 7.9365079365079365079365e-4;
                a = -2.7777777777777777777778e-3 + (p * a);
                a = 0.0833333333333333333333 + (p * a);

                return q + (a / x);
            }

            double b = 8.11614167470508450300e-4;
            b = -5.95061904284301438324e-4 + (p * b);
            b = 7.93650340457716943945e-4 + (p * b);
            b = -2.77777777730099687205e-3 + (p * b);
            b = 8.33333333333331927722e-2 + (p * b);

            return q + (b / x);
        }
        static double GammaLnSmall(double x)
        {
            double z = 1;

            // normalize to interval [2..3) by incr/decr, build z
            double normalized = x;
            double offset = 0;

            while (normalized >= 3)
            {
                offset--;
                normalized = x + offset;
                z *= normalized;
            }

            while (normalized < 2)
            {
                z /= normalized;
                offset++;
                normalized = x + offset;
            }

            if (z < 0)
            {
                z = -z;
            }

            // integer case
            if (normalized == 2)
            {
                return Math.Log(z);
            }

            // normalize x to range [0..1) by incr/decr
            x = x + offset - 2;

            double b = -1378.25152569120859100;
            b = -38801.6315134637840924 + (x * b);
            b = -331612.992738871184744 + (x * b);
            b = -1162370.97492762307383 + (x * b);
            b = -1721737.00820839662146 + (x * b);
            b = -853555.664245765465627 + (x * b);

            double c = 1;
            c = -351.815701436523470549 + (x * c);
            c = -17064.2106651881159223 + (x * c);
            c = -220528.590553854454839 + (x * c);
            c = -1139334.44367982507207 + (x * c);
            c = -2532523.07177582951285 + (x * c);
            c = -2018891.41433532773231 + (x * c);

            return Math.Log(z) + (x * b / c);
        }
        #endregion

        #region GammaRegularized
        public static double GammaRegularized(double a, double x)
        {
            if (a < 0.0) { throw new ArgumentOutOfRangeException(); }
            if (x < 0.0) { throw new ArgumentOutOfRangeException(); }

            const double Epsilon = 0.000000000000001;
            const double BigNumber = 4503599627370496.0;
            const double BigNumberInverse = 2.22044604925031308085e-16;

            if (Number.AlmostZero(a))
            {
                if (Number.AlmostZero(x)) { return Double.NaN; }
                return 1d;
            }
            if (Number.AlmostZero(x)) { return 0d; }
            double ax = (a * Math.Log(x)) - x - GammaLn(a);
            if (ax < -709.78271289338399) { return 1d; }
            if (x <= 1 || x <= a)
            {
                double r2 = a;
                double c2 = 1;
                double ans2 = 1;
                do
                {
                    r2 = r2 + 1;
                    c2 = c2 * x / r2;
                    ans2 += c2;
                }
                while ((c2 / ans2) > Epsilon);
                return Math.Exp(ax) * ans2 / a;
            }
            int c = 0;
            double y = 1 - a;
            double z = x + y + 1;

            double p3 = 1;
            double q3 = x;
            double p2 = x + 1;
            double q2 = z * x;
            double ans = p2 / q2;

            double error;

            do
            {
                c++;
                y += 1;
                z += 2;
                double yc = y * c;

                double p = (p2 * z) - (p3 * yc);
                double q = (q2 * z) - (q3 * yc);

                if (q != 0)
                {
                    double nextans = p / q;
                    error = Math.Abs((ans - nextans) / nextans);
                    ans = nextans;
                }
                else
                {
                    // zero div, skip
                    error = 1;
                }

                // shift
                p3 = p2;
                p2 = p;
                q3 = q2;
                q2 = q;

                // normalize fraction when the numerator becomes large
                if (Math.Abs(p) > BigNumber)
                {
                    p3 *= BigNumberInverse;
                    p2 *= BigNumberInverse;
                    q3 *= BigNumberInverse;
                    q2 *= BigNumberInverse;
                }
            }
            while (error > Epsilon);

            return 1d - (Math.Exp(ax) * ans);
        }
        #endregion

        #region GammaRegularizedInverse
        public static double GammaRegularizedInverse(double a, double y0)
        {
            const double Epsilon = 0.000000000000001;
            const double BigNumber = 4503599627370496.0;
            const double Threshold = 5 * Epsilon;

            // TODO: Consider to throw an out-of-range exception instead of NaN
            if (a < 0 || Number.AlmostZero(a) || y0 < 0 || y0 > 1)
            {
                return Double.NaN;
            }

            if (Number.AlmostZero(y0))
            {
                return 0d;
            }

            if (Number.AlmostEqual(y0, 1))
            {
                return Double.PositiveInfinity;
            }

            y0 = 1 - y0;

            double xUpper = BigNumber;
            double xLower = 0;
            double yUpper = 1;
            double yLower = 0;

            // Initial Guess
            double d = 1 / (9 * a);
            double y = 1 - d - (0.98 * Constants.Sqrt2 * Functions.Erf1Inverse((2.0 * y0) - 1.0) * Math.Sqrt(d));
            double x = a * y * y * y;
            double lgm = GammaLn(a);

            for (int i = 0; i < 10; i++)
            {
                if (x < xLower || x > xUpper)
                {
                    d = 0.0625;
                    break;
                }

                y = 1 - GammaRegularized(a, x);
                if (y < yLower || y > yUpper)
                {
                    d = 0.0625;
                    break;
                }

                if (y < y0)
                {
                    xUpper = x;
                    yLower = y;
                }
                else
                {
                    xLower = x;
                    yUpper = y;
                }

                d = ((a - 1) * Math.Log(x)) - x - lgm;
                if (d < -709.78271289338399)
                {
                    d = 0.0625;
                    break;
                }

                d = -Math.Exp(d);
                d = (y - y0) / d;
                if (Math.Abs(d / x) < Epsilon)
                {
                    return x;
                }

                if ((d > (x / 4)) && (y0 < 0.05))
                {
                    // Naive heuristics for cases near the singularity
                    d = x / 10;
                }

                x -= d;
            }

            if (xUpper == BigNumber)
            {
                if (x <= 0)
                {
                    x = 1;
                }

                while (xUpper == BigNumber)
                {
                    x = (1 + d) * x;
                    y = 1 - GammaRegularized(a, x);
                    if (y < y0)
                    {
                        xUpper = x;
                        yLower = y;
                        break;
                    }

                    d = d + d;
                }
            }

            int dir = 0;
            d = 0.5;
            for (int i = 0; i < 400; i++)
            {
                x = xLower + (d * (xUpper - xLower));
                y = 1 - GammaRegularized(a, x);
                lgm = (xUpper - xLower) / (xLower + xUpper);
                if (Math.Abs(lgm) < Threshold)
                {
                    return x;
                }

                lgm = (y - y0) / y0;
                if (Math.Abs(lgm) < Threshold)
                {
                    return x;
                }

                if (x <= 0d)
                {
                    return 0d;
                }

                if (y >= y0)
                {
                    xLower = x;
                    yUpper = y;
                    if (dir < 0)
                    {
                        dir = 0;
                        d = 0.5;
                    }
                    else
                    {
                        if (dir > 1)
                        {
                            d = (0.5 * d) + 0.5;
                        }
                        else
                        {
                            d = (y0 - yLower) / (yUpper - yLower);
                        }
                    }

                    dir = dir + 1;
                }
                else
                {
                    xUpper = x;
                    yLower = y;
                    if (dir > 0)
                    {
                        dir = 0;
                        d = 0.5;
                    }
                    else
                    {
                        if (dir < -1)
                        {
                            d = 0.5 * d;
                        }
                        else
                        {
                            d = (y0 - yLower) / (yUpper - yLower);
                        }
                    }

                    dir = dir - 1;
                }
            }

            return x;
        }
        #endregion


        


    }
}
