



package jsat.math;

import jsat.distributions.Normal;
import jsat.linear.Vec;
import jsat.math.rootfinding.RiddersMethod;
import static java.lang.Math.*;
import static jsat.math.MathTricks.*;

/**
 * This class provides static methods for computing accurate approximations to 
 * many special functions. <br>
 * <br>
 * All methods should return absolute differences of less than 10<sup>-9</sup> 
 * for all reasonable values. Unreasonable values would be those in areas of 
 * high change (such as those approaching positive / negative infinity). 
 * 
 * @author Edward Raff
 */
public class SpecialMath
{
    public static final double EULER_MASCHERONI =  0.57721566490153286061;
    
    /*
     * TODO this class needs more documentation. All these methods should have 
     * indications of their expected accuracies. 
     */
    
    public static double invXlnX(double y)
    {
        //Method from Numerical Recipies, 3rd edition
        if(y >= 0 || y <= -exp(-1))
            throw new ArithmeticException("Inverse value can not be computed for the range [-e^-1, 0]");
        double u;
        
        if(y < -0.2)
            u = log(exp(-1) - sqrt(2*exp(-1)* (y+exp(-1)) ) );
        else
            u = 10;
        double previousT = 0, t;
        do
        {
            t = (log(y/u)-u)*(u/(1+u));
            u += t;
            if (t < 1e-8 && abs(t + previousT) < 0.01 * abs(t))
                break;
            previousT = t;
        }
        while (abs(t / u) < 1e-15);
        
        return exp(u);

    }

    /**
     * The gamma function is a generalization of the factorial function. 
     * This method provides the gamma function for values from -Infinity to Infinity.
     * <br>Special Values: 
     * <ul>
     * <li>{@link Double#NaN} returned if z = 0 or z = {@link Double#NEGATIVE_INFINITY}</li>
     * </ul>
     * @param z any real value
     * @return Γ(z)
     */
    public static double gamma(double z)
    {
        if(z == 0)//It is actualy infinity*I, where I - sqrt(-1).
            return Double.NaN;
        else if(z == Double.POSITIVE_INFINITY)
            return z;
                
        else if(z < 0)
        {
            /*
             * Using the identity
             *
             *                      __
             *                     -||
             * gamma(-z) = -------------------
             *                            __
             *             z gamma(z) sin || z
             */

            z = -z;
            return -PI/(z*gamma(z)*sin(PI*z));
        }

        /**
         * General case for z > 0, from Numerical Recipes in C (2nd ed. Cambridge University Press, 1992). |error| is <= 2*10^-10 for all z > 0
         *
         *               ____ /       6        \
         *              /  __ |     =====  p   |
         *            \/ 2 || |     \       n  |          z + 0.5  - (z + 5.5)
         * gamma(z) = ------- |p  +  >    -----| (z + 5.5)        e
         *               z    | 0   /     z + n|
         *                    |     =====      |
         *                    \     n = 1      /
         *
         *
         * see http://www.rskey.org/gamma.htm
         * 
         */

        double p[] =
        {
            1.000000000190015,76.18009172947146,-86.50532032941677,
            24.01409824083091,-1.231739572450155,1.208650973866179e-3,-5.395239384953e-6
        };

        double innerLoop = 0;
        for(int n = 1; n < p.length; n++)
            innerLoop += p[n]/(z+n);

        double result = p[0] + innerLoop;

        result *= sqrt(2*PI)/z;

        return result*pow(z+5.5, z+0.5)*exp(-(z+5.5));
    }

    /**
     * Computes the natural logarithm of {@link #gamma(double) }. 
     * This method is more numerically stable than taking the log of the result of Γ(z). 
     * 
     * <br>Special Values: 
     * <ul>
     * <li>{@link Double#POSITIVE_INFINITY} returned if z &le; 0</li>
     * <li>{@link Double#NaN} returned if z = {@link Double#NEGATIVE_INFINITY}</li>
     * </ul>
     * 
     * @param z any real number value
     * @return Log(Γ(z))
     */
    public static double lnGamma(double z)
    {
        if(z == Double.NEGATIVE_INFINITY)
            return Double.NaN;
        else if(z == Double.POSITIVE_INFINITY)
            return z;
        else if(z <= 0)
            return Double.POSITIVE_INFINITY;

        /*
         * Lanczos approximation for the log of the gamma function, with |error| < 10^-15 for all z > 0 (Almost full double precision)
         */

        int j;
        double x, tmp, y, ser = 0.999999999999997092;
        double[] c = new double[]
        {
            57.1562356658629235, -59.5979603554754912,
            14.1360979747417471, -0.491913816097620199, .339946499848118887e-4,
            .465236289270485756e-4, -.983744753048795646e-4, .158088703224912494e-3,
            -.210264441724104883e-3, .217439618115212643e-3, -.164318106536763890e-3,
            .844182239838527433e-4, -.261908384015814087e-4, .368991826595316234e-5
        };

        y = x = z;
        tmp = x+671.0/128.0;
        tmp = (x+0.5)*log(tmp)-tmp;
        for (j = 0; j < 14; j++)
        {
            y++;
            ser += c[j] / y;
        }

        return tmp+log(2.5066282746310005*ser/x);
    }
    
    /**
     * Positive zero of the digamma function
     */
    static final private double digammaPosZero = 1.461632144968362341;
    
    /**
     * Computes the value of the digamma function, &#936;(x), which is the 
     * derivative of {@link #lnGamma(double) }. <br>
     * <br>
     * This method may return:<br>
     * <ul>
     * <li> {@link Double#NaN} for zero and negative integer values (would be complex infinity)</li>
     * </ul>
     * <br>
     * This method should be accurate to an absolute difference of 
     * 10<sup>-14</sup> for all values that are not near an asymptote. 
     * 
     * @param x the value to compute the digamma function at
     * @return the value of &#936;(x)
     */
    public static double digamma(double x)
    {
        if(x == 0)
            return Double.NaN;//complex infinity
        else if(x < 0)//digamma(1-x) == digamma(x)+pi/tan(pi*x), to make x positive
        {
            if(Math.rint(x) == x)
                return Double.NaN;//the zeros are complex infinity
            return digamma(1-x)-PI/tan(PI*x); 
        }
        else if(x < 2)//shift the value into [2, Inf]
            return digamma(x+1) - 1/x;
        double approx= log(x);//rel error of 10^-11 for x >= 100 10^-13 for x >= 250, near machine precision at x >= 500
        double approxS = -1/(2*x);
        double approxSS = -1/(12*x*x);

        if(x <= 7)
            return (x-digammaPosZero)*hornerPolyR(digamma_p_2_7, x)/hornerPolyR(digamma_q_2_7, x);
        else if(x <= 70)
            return approx + (approxS + (approxSS + hornerPolyR(digamma_p_7_70, x)/hornerPolyR(digamma_q_7_70, x)));
        if(x < 500)
            return approx + (approxS + (approxSS + hornerPolyR(digamma_adj_p, x)/hornerPolyR(digamma_adj_q, x)));
        else
            return approx;
    }
    
    /**
     * The upper polynomial for adjustment to the digamma function approximation
     * in the range [5, 500]. Computationally the absolute difference is less 
     * than 1e-14 after 70
     */
    private static final double[] digamma_adj_p = new double[] 
    {
        6.662015538739419312162593679911046099212e-17, 8.191636265707720257639745288823408068504e-14,
        3.594113658754209591796105198667912207196e-11, -6.860373206319231677418396622671250368996e-9,
        0.05718773229071206996201131377297202073545
    };
    /**
     * The lower polynomial for adjustment to the digamma function approximation
     * in the range [5, 500]
     */
    private static final double[] digamma_adj_q = new double[] 
    {
        6.862596697230474833874884670160991920179, 0.002471283005727189026784415544975350355939,
        3.307589698304449008222736059167667287747, 0.2934987180435057544393612182922529976682,
        1.0
    };
    
    /**
     * The upper polynomial for computing the digamma function, fit against 
     * psi(x)/(x-1.46163...)
     * 
     * Absolute difference less than 1e-15 for the range [2,7]
     * 
     * Idea from Chebyshev Approximations for the Psi Function. Mathematics of 
     * Computation, Volume 27, Number 121. By Cody, Strecok, and Thacher, 
     */
    private static final double[] digamma_p_2_7 = new double[] 
    {
        0.00803356767428942100, 8.71902391724677187,
        445.627457353132455, 5678.99715950204957,
        23638.5586690114249, 32569.4896509749708,
        11137.2906503774953
    };
    /**
     * The lower polynomial for computing the digamma function, fit against 
     * psi(x)/(x-1.46163...)
     * 
     * Absolute difference less than 1e-15 for the range [2,7]
     */
    private static final double[] digamma_q_2_7 = new double[] 
    {
        1.63429827694094891, 123.543728215902577,
        2208.64039265513217, 13061.4741667968999,
        27097.1711564673534, 16271.6188328948966, 
        1.0
    };
    
    /**
     * The upper polynomial for approximating the digamma function 
     * in the range [7, 70]. The coefficients are in reverse order 
     * 
     * Absolute difference is less than 1e-14 for the whole range 
     * (1e-15 for most of it!) 
     */
    private static final double[] digamma_p_7_70 = new double[] 
    {
        1.184517214829426509228398e-14, -3.115254117966573324590434e-12,
        3.376965291465973686259273e-10, -1.965130493876512169028824e-8,
        -0.04505398537169693060781875
        
    };
    /**
     * The lower polynomial for approximating the digamma function in the range 
     * [7, 70]. The coefficients are in reverse order
     */
    private static final double[] digamma_q_7_70 = new double[] 
    {
        -5.406558865634368930280970, 0.001668808697639819914012839,
        -2.595288371769254271653798, 0.1463978383403391904410391,
        1
    };

    
    /**
     * Computes the Riemann zeta function &zeta;(x) for some value of x<br>
     * <br>
     * This method may return:<br>
     * <ul>
     * <li> {@link Double#NaN} for x = 1 (would be complex infinity)</li>
     * </ul>
     * <br>
     * <br>
     * <b>NOTE:</b> This method is not yet complete in terms of accuracy. <br>
     * <ul>
     *  <li>For x &lt; -0.5, the values returned will be of the correct 
     *      magnitude - but are not very accurate</li>
     *  <li>For x in [-0.5, 2.5], the values returned will be of reasonable 
     *      accurate (absolute difference around 10<sup>-7</sup>), unless it is very close to 1</li>
     *  <li>For x &gt; 2.5, the result will be very accurate (absolute difference less than 10<sup>-14</sup></li>
     * </ul>
     * @param x a real valued input 
     * @return &zeta;(x)
     */
    public static double zeta(double x)
    {
        if(x == 1)
            return Double.NaN;
        if(x == 0)
            return -0.5;
        if(x < 0 || abs(1-x) <= 0.2)
        {
            if(x <= 0.2 && x > -2.)
            {
                /*
                 * For this specific range we keep our own approximant,
                 * see below comment
                 */
                return hornerPolyR(zeta_p_special, x)/hornerPolyR(zeta_q_special, x);
            }
            /*
             * http://dlmf.nist.gov/25.4#E2
             * 
             * Reflect zeta across 1 if negative to make it positive. 
             * 
             * Reflect zeta across 1 if it is just less than one so that it will
             * be just more than 1 (much easier to compute)
             * 
             * Reflect if just more than 1 to an area that is easier to approximate
             */
            double otherPart = 2*pow(2*PI, x-1)*sin(PI/2*x);
            if(x < 0)
                return otherPart*exp(lnGamma(1-x)*log(zeta(1-x)));
            else//log(zeta(1-x)) would have caused a NaN
                return otherPart*gamma(1-x)*zeta(1-x);
        }
        if(x < 14)
            return hornerPolyR(zeta_p_l14, x)/hornerPolyR(zeta_q_l14, x);
        if(x < 50)
        {
            //use truncated form of http://dlmf.nist.gov/25.2#E3
            double mul = 1/(1-pow(2, 1-x));
            double sumP = 0;
            double sumN = 0;
            for(int i = 11; i >= 1; i-=2)//all odd values are positive
                sumP += pow(i, -x);
            for(int i = 10; i >= 1; i-=2)//all even values are negative
                sumN -= pow(i, -x);
            return mul*(sumP+sumN);
        }
        //else x>=50, 1 is so close we might as well use it
        return 1;
    }
    
    /**
     * Coefficients used for computing the zeta function involving the coef[k_]:=BernoulliB[2 k]/Factorial[2 k] term 
     */
    private static final double[] zetBernCoefs = new double[]
    {
        0.0,//FAKE VALUE, Never used
        0.083333333333333333333,-0.0013888888888888888889,0.000033068783068783068783,
        -8.2671957671957671958e-7,2.0876756987868098979e-8,-5.2841901386874931848e-10,
        1.3382536530684678833e-11,-3.3896802963225828668e-13,8.5860620562778445641e-15,
        -2.1748686985580618730e-16
    };
    
    /**
     * Implements the Hurwitz zeta function &#950;(x, a) <br>
     * The relative error of this implementation is less than 1e-4 for most
     * positive values of x, though is lower for negative values of x. Accuracy
     * improves as a increases.
     *
     * @param x the first argument of the Hurwitz zeta function/
     * @param a the second argument, must be positive in all cases.
     * @return the value of &#950;(x, a), or {@link Double#NaN} if the result
     * would be a Complex Infinity
     */
    public static double zeta(double x, double a)
    {
        if(x == 0)
            return 0.5-a;
        if(a <= 0)
            return Double.NaN;//Results would be complex infinity
        if(x == 1)
            return Double.NaN;//Result would be complex infinity
        if(a > 1e7 || (x < 0 && x >=-100 && a >= 1e3))
        {
            //In the limit this is correct. At 1e6, relative error is <= 10^(-8) for x [1, 10,000]
            //Its also good for  x in [-100, 1] - so use in this case for small q values of 1e3
            return (1/(x - 1) + 1/(2*a))*pow(a, 1 - x);
        }
        if(x < 0)
        {
            if(a <= 1)//Reflection using https://dlmf.nist.gov/25.11#E9
            {
                double s = 1-x;
                double sum = 0;
                for(int n = 1; n <= 20; n++)
                    sum += pow(n, -s) * cos(PI*0.5*s-2*n*PI*a);
                return exp(log(2) + lnGamma(s) - s*log(2*PI))*sum;
            }
            else //reduce a using https://dlmf.nist.gov/25.11.E4
            {
                double m = Math.floor(a);
                if(m == a)//a was an integer, so lets adjust m by 1
                    m--;
                //Now a_new will be in (0, 1]
                double a_new = a-m;
                double sum = 0;
                for(int n = (int) (m-1); n >= 0; n--)
                {
                    double t = pow(n+a_new, -x);
                    sum += t;
                    if (t / sum < 1e-6)//eventually our contributions will be tiny, break when it happens
                        break;
                }
                    
                return zeta(x, a_new) - sum;
            }
        }
        //The error on this seems high... so lets not do that for now
//        if(x <= -1 && Math.rint(x) == x && x >= Integer.MIN_VALUE)//Use https://dlmf.nist.gov/25.11.E13 
//        {
//            int n = (int) -x;
//            return -exp(reLnBn(n+1)+log(a)-log(n+1));
//        }
        
        //General case, done as outlined here https://math.stackexchange.com/questions/917100/numerical-evaluation-of-hurwitz-zeta-function?rq=1 
        
        double part1 = 0;
        double part2 = 0;
        final int n = 9;//Constant chosen by recomendation for double precision
        for(int k = 0; k <= n; k++)
        {
            part1 += 1/pow(a+k, x);
        }
        
        //Second summerant term is zetBernCoefs[k] * ( ((a+n)^(-k-x) Gamma[-2+3 k+x])/Gamma[-2+2 k+x] )
        //simplfy 2nd term with logs as exp( -(k+x) Log[a+n]+Log[Gamma[-2+3 k+x]]-Log[Gamma[-2+2 k+x]] )
        
        for(int k = 1; k < zetBernCoefs.length; k++ )
        {
            part2 += zetBernCoefs[k] * exp(-(k+x)* log(a+n)+lnGamma(-2+3*k+x)-lnGamma(-2+2*k+x));
        }
        
        return part1 + pow(a+n, 1-x)/(x-1) - 1./(2*pow(a+n, x)) + part2 ;
    }
    
    /**
     * Computes the n'th harmonic number <i>H<sub>n</sub></i>. <br>
     * Note that the relative error of this method is less than 10<sup>-6</sup>
     * for the entire range, and decreases as <i>n</i> increases.
     *
     * @param n any non-negative value
     * @return the value of <i>H<sub>n</sub></i>, of {@link Double#NaN} if <i>n</i> &lt; 0
     */
    public static double harmonic(double n)
    {
        if(n == 0)
            return 0;
        else if(n < 0)
            return Double.NaN;
        if(n < 2)//Lets use Economized Rational Approximations 
        {
            final double[] era_up_4 = new double[]
            {
                0.00067701604350, 5483.1926096, 4283.42774099, 602.59193727 
            };
            final double[] era_low_5 = new double[]
            {
                3333.39376261, 5039.82176406, 1856.44874041, 140.548664241, -1.000000000000
            };
            
            return hornerPoly(era_up_4, n)/hornerPoly(era_low_5, n);
        }
        //else, this 
        
        double h_n = log(n) + EULER_MASCHERONI + 1./(2*n);
        
        if(n >= 1000000)
            return h_n;//Close enough with rel error < 5e-15
        //We are going to add at least two terms then
        
        //1/(240 n^8)-1/(252 n^6)+1/(120 n^4)-1/(12 n^2)
        double nSqrd = n*n;
        double nQuad = nSqrd*nSqrd;
        h_n += -1/(12 * nSqrd) + 1/(120 * nQuad);
        
        double n8 = nQuad*nQuad;
        double n6 = nSqrd*nQuad;
        h_n += +(1/(240*n8))-1/(252*n6);
        
        return h_n;
    }
    
    /**
     * Computes the generalized n'th harmonic number of the m'th order
     * <i>H<sub>n</sub><sup>(m)</sup></i>. <br>
     *
     * @param n  any non-negative value
     * @param m the harmonic order
     * @return the value of <i>H<sub>n</sub><sup>(m)</sup></i>
     */
    public static double harmonic(double n, double m)
    {
        if (n == 0)
            return 0;
        else if (n < 0)
            if (m == 0)
                return n;
            else
                return Double.NaN;//truth is complex infinity
        else if(m == 1)
            return harmonic(n);

        return zeta(m)-zeta(m,n+1);
    }
    
    /**
     * upper polynomial approximation of the zeta function between [-0.20, 0.5]
     */
    private static double[] zeta_p_special = 
    {
        -5.276454584406249e-6, 0.00014685004463733906,
        -0.0029925134974932046, -0.03542393126377964,
        -0.2062582384669163, -0.5425801056911627, 
        -0.500000000000001
    };
    /**
     * Lower polynomial approximation of the zeta function between [-0.20, 0.5]
     */
    private static double[] zeta_q_special = 
    {
        -3.9917203368386765e-6, -0.00017067854372054898,
        -0.0029262260848543224, -0.03374331488845255,
        -0.21043893362112356, -0.75271685502711,
        1
    };
    /**
     * Upper polynomial approximation of the zeta function between [1.20, 14]
     */
    private static double[] zeta_p_l14 = 
    {
        -8.60637125178308808e-8, -1.15605577219727645e-6,
        -0.0000431166632155845063, -0.000526489508370994743,
        -0.00681551580000143337, -0.0630725028974146228,
        -0.294618454483912025, -0.634306044238860948,
        -0.500000284011974157
    };
    
    /**
     * Lower polynomial approximation of the zeta function between [1.20, 14]
     */
    private static double[] zeta_q_l14 = 
    {
        -8.61872960765516585e-8, -1.13850129664146526e-6,
        -0.0000442195744615244028, -0.000486888035876849659,
        -0.00768589078947058848, -0.0516310207074782699, 
        -0.370887845301558816, -0.569262911426117336, 
        1
    };
    
    /**
     * Contains the exact values for Re[Log[Bn[x]]] for all even values less than 50
     */
    private static final double[] reLnBn_sub_50 = 
    {
        0, 
        -1.79175946922805500, -3.40119738166215538, -3.73766961828336831, 
        -3.40119738166215538, -2.58021682959232517, -1.37391706441133552, 
        0.154150679827258304, 1.95898950623372649, 4.00680901051075934, 
        6.27122326708809952, 8.73103330984004285, 11.3688270444561602, 
        14.1700452298205194, 17.1223324620051408, 20.2150715380410726, 
        23.4390405118087582, 26.7861544648980854, 30.2492673319812623, 
        33.8220172715222676, 37.4987042389444309, 41.2741917911127339, 
        45.1438274079056960, 49.1033771613463772, 53.1489716411456250
    };
    
    /**
     * Computes the real part of the natural logarithm of the Bernoulli numbers. 
     * <br> The Bernoulli zeros for odd n will return 
     * {@link Double#NEGATIVE_INFINITY} and for any value less than 0 will 
     * return {@link Double#NaN}.  
     * <br><br>
     * <br><br>
     * Currently only accurate to an absolute difference of 10<sup>-11</sup>
     * 
     * @param n the integer Bernoulli value to obtain an approximation of
     * @return <i>Re(Log(B<sub>n</sub>))</i> 
     */
    public static double reLnBn(int n)
    {
        if(n < 0)
            return Double.NaN;
        if(n == 1)
            return -log(2);
        if(n % 2 != 0)
            return Double.NEGATIVE_INFINITY;
        if(n >= 50)//rel err < 1e-14
        {
            //Log[-2^(3/2 - n) ((3 i)/e)^n n^(1/2 + n) ((3 + 40 n^2)/(-1 + 120 n^2))^n Pi^(1/2 - n)]
            double x = 0;
            x += (3.0/2.0-n)*log(2);//ignoring + imaginary pi here
            x += n*(log(3)-1);//ignoring + imaginary pi/2 in the parenthesis 
            x += (n+0.5)*log(n);
            x += n*(log(40*n*n+3)-log(120*n*n-1));
            x += (0.5-n)*log(PI);
            return x;
        }
        return reLnBn_sub_50[n/2];
    }
    
    /**
     * Computes an approximation to the n'th Bernoulli number 
     * <i>B<sub>n</sub></i>. 
     * The Bernoulli numbers grow in value rapidly, and so the accuracy of this
     * method decays quickly. n &gt; 20 should have the correct order of magnitude,
     * but may not have many significant figures. {@link #reLnBn(int) } should
     * be used instead when possible.
     *
     * @param n the bernoulli number to compute
     * @return <i>B<sub>n</sub></i>
     */
    public static double bernoulli(int n)
    {
        if(n < 0)
            return Double.NaN;
        if(n == 0)
            return 1;
        if(n == 1)
            return -0.5;
        if(n % 2 != 0)
            return 0;
        int sign = 1;
        if(n > 2 && n % 4 == 0)
            sign = -1;
        return sign*exp(reLnBn(n));
    }
    
    public static double erf(double x)
    {
       /*
        * erf(x) = 2 * cdf(x sqrt(2)) -1
        *
        * where cdf is the cdf of the normal distribution
        */

        return 2 * Normal.cdf(x * sqrt(2.0), 0, 1)-1; 
    }

    public static double invErf(double x)
    {
        /*
        * inverf(x) = invcdf(x/2+1/2)/sqrt(2)
        *
        * where invcdf is the inverse cdf of the normal distribution
        */

        return Normal.invcdf(x/2+0.5, 0, 1)/sqrt(2.0);
    }

    public static double erfc(double x)
    {
        /*
        * erf(x) = 2 * cdf(-x sqrt(2))
        *
        * where cdf is the cdf of the normal distribution
        */
        return 2 * Normal.cdf(-x * sqrt(2.0), 0, 1);
    }

    public static double invErfc(double x)
    {
       /*
        * inverf(x) = invcdf(x/2)/-sqrt(2)
        *
        * where invcdf is the inverse cdf of the normal distribution
        */

        return Normal.invcdf(x/2, 0, 1)/-sqrt(2.0);
    }

    /**
     * Computes the Beta function B(z,w)
     * @param z
     * @param w
     * @return B(z,w)
     */
    public static double beta(double z, double w)
    {
        return exp(lnBeta(z, w));
    }
    
    public static double lnBeta(double z, double w)
    {
        /*
         * The beta function is defined by
         *
         *           Gamma(z) Gamma(w)
         * B(z, w) = -----------------
         *             Gamma(z + w)
         *
         * However, the definition is numericaly unstable (large value / large value to small result & small input).
         * Taking the log of each size and then exponentiating gives a more stable method of computing the result
         *
         *            lnGamma(z) + lnGamma(w) - lnGamma(z + w)
         * B(z, w) = e
         */
        return lnGamma(z)+lnGamma(w)-lnGamma(z+w);
    }
    
    /**
     * Computes the regularized incomplete beta function, I<sub>x</sub>(a, b). The result of which is always in the range [0, 1]
     * 
     * @param x any value in the range [0, 1]
     * @param a any value &ge; 0
     * @param b any value &ge; 0
     * @return the result in a range of [0,1]
     */
    public static double betaIncReg(double x, double a, double b)
    {
        if(a <= 0 || b <= 0)
            throw new ArithmeticException("a and b must be > 0, not" + a+ ", and " + b);
        if(x == 0 || x == 1)
            return x;
        else if(x < 0 || x > 1)
            throw new ArithmeticException("x must be in the range [0,1], not " + x);
        
        //We use this identity to make sure that our continued fraction is always in a range for which it converges faster
        if(x > (a+1)/(a+b+2) || (1-x) < (b+1)/(a+b+2) )
            return 1-betaIncReg(1-x, b, a);
        
        
        /*
         * All values are from x = 0 to x = 1, in 0.025 increments
         * a = 0.5, b = 0.5: max rel error ~ 2.2e-15
         * a = 0.5, b = 5: max rel error ~ 2e-15
         * a = 5, b = 0.5: max rel error ~ 1.5e-14 @ x ~= 7.75, otherwise rel error ~ 2e-15
         * a = 8, b = 10: max rel error ~ 9e-15, rel error is clearly not uniform but always small
         * a = 80, b = 100: max rel error ~ 1.2e-14, rel error is clearly not uniform but always small
         */
        
        double numer = a*log(x)+b*log(1-x)-(log(a)+lnBeta(a, b));
        
        return exp(numer)/regIncBeta.lentz(x,a,b);
    }
    
    /**
     * Computes the inverse of the incomplete beta function,
     * I<sub>p</sub><sup>-1</sup>(a,b), such that {@link #betaIncReg(double, double, double) I<sub>x</sub>(a, b) } = <tt>p</tt>. 
     * The returned value, x, will always be in the range [0,1].
     * The input <tt>p</tt>, must also be in the range [0,1]. 
     * 
     * @param p any value in the range [0,1]
     * @param a any value &ge; 0
     * @param b any value &ge; 0
     * @return the value x, such that {@link #betaIncReg(double, double, double) I<sub>x</sub>(a, b) }  will return p. 
     */
    public static double invBetaIncReg(double p, double a, double b)
    {
        if(p < 0 || p > 1)
            throw new ArithmeticException("The value p must be in the range [0,1], not" + p);
        return RiddersMethod.root(0, 1, (x) -> betaIncReg(x, a, b)- p);
    }
    
    /**
     * Computes the regularized gamma function Q(a,z) = Γ(a,z)/Γ(a). <br> 
     * This method is more numerically stable and accurate than computing
     * it via the direct method, and is always in the range [0,1]. <br><br>
     * Note: The this method returns {@link Double#NaN} for a&lt;0, though real values of Q(a,z) do exist
     * @param a any value &ge; 0
     * @param z any value &gt; 0
     * @return Q(a,z)
     */
    public static double gammaQ(double a, double z)
    {
        if(z<= 0 || a < 0 )
            return Double.NaN;
        
        if(z < a+1)
            return 1-gammaPSeries(a, z);
        
        /**
         * On the range of x from 0.5 to 50
         * a=0.15, |rel error| is ~ 3e-15 for most values of x, with a bad spot of |rel error| ~ 3e-11 when x ~= 5.75
         * a=0.5, max |rel error| is 3.9e-15
         * a=1, max |rel error| is 4e-15, rel error groining as x -> infinity
         * a=5, max |rel error| is 7e-13, but only near x ~= 0, most is in the range 3e-15
         * a=10, max |rel error| is 4e-6, but only near x ~= 0, most is in the range 3e-15
         */
        return exp(a*log(z)-z-lnGamma(a))/gammaQ.lentz(a, z);
    }
    
    public static double gammaPSeries(double a, double z)
    {
        double ap = a;
        double sum;
        double del = sum = 1.0/a;
        
        do
        {
            ap += 1.0;
            del *= z/ap;
            sum += del;
            if(abs(del) < abs(sum)*1e-15)
                return sum*exp(-z+a*log(z)-lnGamma(a));
        }
        while(true);
        
    }
    
    /**
     * Returns the regularized gamma function P(a,z) = γ(a,z)/Γ(a). <br>
     * This method is more numerically stable and accurate than computing
     * it via the direct method, and is always in the range [0,1]. <br><br>
     * Note: The this method returns {@link Double#NaN} for a&lt;0, though real values of P(a,z) do exist
     * @param a any value &ge; 0
     * @param z any value &gt; 0
     * @return P(a,z)
     */
    public static double gammaP(double a, double z)
    {
        if(z<= 0 || a < 0)
            return Double.NaN;
        
        if(z < a+1)
            return gammaPSeries(a, z);
        
        /*
         * This method is currently usntable for values of z that grow larger, so it is not currently in use
         * return exp(a*log(z)-z-lnGamma(a))/gammaP.lentz(a,z);
         */
        return 1 - gammaQ(a, z);
    }
    
    /**
     * Finds the value <tt>x</tt> such that {@link #gammaP(double,double) P(a,x)} = <tt>p</tt>. 
     * @param a any real value 
     * @param p and value in the range [0, 1]
     * @return the inverse
     */
    public static double invGammaP(double p, double a)
    {
        //Method from Numerical Recipies 3rd edition p 263
        if(p < 0 || p > 1)
            throw new ArithmeticException("Probability p must be in the range [0,1], "+ p + "is not valid");
        //First an estimate must be obtained
        double am1 = a-1;
        double lnGamA = lnGamma(a);
        double x;//the to be returned
        double afac = 1;//ONLY used when a>1
        double lna1 =1;//ONLY used when a>1
        if(a > 1)
        {
            lna1 = log(am1);
            afac = exp(am1*(lna1-1)-lnGamA);
            double pp = (p < 0.5) ? p : 1-p;
            double t = sqrt(-2*log(pp));
            //Now our inital estimate
            x = (2.30753+t*0.27061)/(1.+t*(0.99229+t*0.04481)) - t;
            if(p < 0.5)
                x = -x;
            x = max(1e-3, a * pow(1.0 - (1.0/(9.*a)) - x/(3.*sqrt(a))  , 3) );//if the estimate is too small, increase it
            
        }
        else
        {
            double t = 1.0 - a*(0.253+a*0.12);
            if(p < t)
                x = pow(p/t, 1.0/a);
            else
                x = 1-log(1-(p-t)/(1.0-t));
        }
        
        //Estimate obtained, now refinement
        for(int j = 0; j < 12; j++)
        {
            if (x <= 0)//x is very small, return 0 b/c rounding errors and loss of precision will make accuracy imposible
                return 0;
            double err = gammaP(a, x) - p;
            double t;
            if(a > 1)
                t = afac*exp(-(x-am1)+am1*(log(x)-lna1));
            else
                t = exp(-x+am1*log(x)-lnGamA);
            
            double u = err/t;
            t = u / (1-0.5*min(1 , u*(am1/x - 1)) );//Halley's method
            x -= t;
            if(x <= 0)//bais x from going negative (means we are getting a bad value)
                x = (x+t)/2;
            if(abs(t) < 1e-8*x)
                break;//the error is the (1e-8)^2, if we reach this point we have already converged
        }
        
        return x;
    }
    
    /**
     * Computes the incomplete gamma function, Γ(a,z). 
     * <br>
     * Returns {@link Double#NaN} for z &le; 0
     * @param a any value (-&infin;, &infin;)
     * @param z any value &gt; 0
     * @return Γ(a,z)
     */
    public static double gammaIncUp(double a, double z)
    {
        if(z <= 0)
            return Double.NaN;
        /**
         * On the range of x from 0.5 to 50
         * a=0.15, max |rel error| is ~1.3e-11, but only in a small range x ~= 11. Otherwise rel error ~ 2e-15
         * a=0.5, max |rel error| is ~3.6e-15, less in all places otherwise
         * a=1, max |rel error| is ~4e-15, rel error grows as x-> infinity
         * a=5, max |rel error| is 7e-13, but only near x ~= 0, most is in the range 2e-15
         * a=10, max |rel error| is 4.9e-6, but only near x ~= 0, most is in the range 5e-15
         */
        return exp(a*log(z)-z)/upIncGamma.lentz(a,z);
    }
    
    /**
     * Computes the lower incomplete gamma function, γ(a,z). 
     * <br>
     * Returns {@link Double#NaN} for z &le; 0
     * @param a any value (-&infin;, &infin;)
     * @param z any value &gt; 0
     * @return γ(a,z)
     */
    public static double gammaIncLow(double a, double z)
    {
         if(z <= 0)
            return Double.NaN;
        /**
         * On the range of x from 0.5 to 50
         * a=0.15, see a=1
         * a=0.5, see a=1
         * a=1, max |rel error| is ~1.7e-13, rel error starts at 1e-15 and grows to the max as z increases
         * a=5, max |rel error| is 3e-13, but only near x ~= 0. for x > epsilon error starts at 5e-15 and grows with z up to 1e-13
         * a=10, max |rel error| is 2e-7, but only near x ~= 0, most is in the range 5e-15, grows to 5e-14 as z-> infinity
         */
        return exp(lnLowIncGamma(a, z));
    }
    
    private static double lnLowIncGamma(double a, double x)
    {

        /*
         * We compute the log of the lower incomplete gamma function by taking the log of
         *
         *                   oo
         *                  =====
         *            -x  a \         Gamma(a)      n
         * y(a, x) = e   x   >    ---------------- x
         *                  /     Gamma(a + 1 + n)
         *                  =====
         *                  n = 0
         *
         * Which becomes
         *
         *
         *                                   / oo                      \
         *                                   |=====                    |
         *                                   |\         Gamma(a)      n|
         * log(y(a, x)) = -x + log(x) a + log| >    ---------------- x |
         *                                   |/     Gamma(a + 1 + n)   |
         *                                   |=====                    |
         *                                   \n = 0                    /
         *
         * To reduce over flow of the gammas and exponentation in the summation, we compute the sum as
         *
         *  oo
         * =====
         * \      LogGamma(a) - LogGamma(a + 1 + n) + ln(x) n
         *  >    e
         * /
         * =====
         * n = 0
         *
         * Testin with x=0.5 to 100 (in increments of 0.5)
         * a=0.15, max relative error is ~ 6e-13, with x < 30 having a relative error smaller than 5e-14
         * a=0.5, maximum relative and absolute eror is ~1.5e-12 and ~1.5e-12 repsectivly, the error starts getting above e-15 when x = 25, and x=50 the error is up to 2.2e-13
         * a=10, maximum relative and absolute eror is ~4e-14 and ~4e-13 repsectivly, the error starts getting abovee-15 when x = 49. For x near zero (up to x ~ 2.5) the error is higher before droping, ~10^-14
         * a=25, maximum relative error is ~9.99e-15. From x~ 0 to 14, the error is worse, then droping to ~e-16 .
         * a=50, accuracy starting to degrad badly. From x~ 0 to 18 the error goes from 1.3 to 1e-7, the erro grows at an exponential rate as x -> 0. As x-> Infinity the error gets back down to ~5e-16
         */

        //Sumation first

        double lnGa = lnGamma(a);
        /**
         * This value will be updated by the property Gamma(z+1) = Gamma(z) * z, which - when taken the log of, is <br>
         * LnGamma(z+1) = LnGamma(z) + ln(z)
         */
        double lnGan = lnGa+log(a);
        double n = 0;
        /**
         * this is the n * ln(x)  term. it will be updated by adding the log of x at each step
         */
        double lnXN = 0;
        double lnX = log(x);

        //Set up, now start summing
        double term = exp(lnGa - lnGan + lnXN);
        double sum = term;
        while(term > 1e-15)
        {
            n++;
            lnXN += lnX;
            lnGan += log(a+n);

            term = exp(lnGa - lnGan + lnXN);

            sum += term;
        }

        //now the rest

        return -x + lnX*a +log(sum);
    }
    @SuppressWarnings("unused")
    private static double lnLowIncGamma1(double a, double x)
    {
        double inter = lowIncGamma.lentz(a,x);
        if(inter <= 1e-16)//The result was ~0, in which case Gamma[a,z] ~= Gamma[a]
            return lnGamma(a);
        return a*log(x)-x-log(inter);
    }
    
    private static final ContinuedFraction lowIncGamma = new ContinuedFraction()
    {

        @Override
        public double getA(int pos, double... args)
        {

            if (pos % 2 == 0)
            {
                pos /= 2;//the # of the even term

                return pos * args[1];
            }
            else
            {
                pos = (pos + 0) / 2;

                return -(args[0] + pos) * args[1];
            }
        }

        @Override
        public double getB(int pos, double... args)
        {

            return args[0] + pos;
        }
    };

    /**
     * See http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/
     */
    private static final ContinuedFraction gammaQ = new ContinuedFraction()
    {

        @Override
        public double getA(int pos, double... args)
        {
            return pos*(args[0]-pos);
        }

        @Override
        public double getB(int pos, double... args)
        {

            return (1 +pos*2) - args[0] + args[1];
        }
    };
    
    /**
     * Using the formula from http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0009/
     * 
     * Note the formula is given in terms of gammaQ
     * 
     * {@link #gammaQ} is accurate for P and Q for the range of z > 0, this is used for z <= 0
     * 
     */
    private static final ContinuedFraction gammaP = new ContinuedFraction()
    {

        @Override
        public double getA(int pos, double... args)
        {
            if(pos % 2 == 0)//even step
            {
                pos/=2;
                return args[1]*pos;
            }
            //Else its an odd step
            pos= (pos+1)/2;
            return -args[1]*(args[0]+pos);
        }

        @Override
        public double getB(int pos, double... args)
        {
            return args[0] + pos;
        }
    };
    
    /**
     * continued fraction generated from mathematica
     * 
     * f(a,z) = e^-x*x^a / CF
     * 
     * a_k = (a-k)k
     * b_k = 1+k*2-a+x
     * 
     */
    private static final ContinuedFraction upIncGamma = new ContinuedFraction()
    {

        @Override
        public double getA(int pos, double... args)
        {
            return (args[0]-pos)*pos;
        }

        @Override
        public double getB(int pos, double... args)
        {
            return (1+pos*2)-args[0]+args[1];
        }
    };
    
    /**
     * See http://dlmf.nist.gov/8.17
     */
    private static final ContinuedFraction regIncBeta = new ContinuedFraction() {

        @Override
        public double getA(int pos, double... args)
        {
            if(pos % 2 == 0)
            {
                pos /=2;
                return pos*(args[2]-pos)*args[0]/ ( (args[1] + 2*pos-1)*(args[1] + 2*pos) );
            }
            
            pos = (pos-1)/2;
            
            double numer  = -(args[1] + pos)*(args[1]+args[2]+pos)*args[0];
            double denom = (args[1] + 2*pos)*(args[1]+1+2*pos);
            
            return numer/denom;
        }

        @Override
        public double getB(int pos, double... args)
        {
            return 1.0;
        }
    };
}
