using System;

namespace Science.Mathematics.SpecialFunction
{
	/// <summary>
	/// Factorial Function, Gamma Function.
	/// </summary>
	public class GammaRelated
	{
		public GammaRelated()
		{
		}
		public static double LogGamma(double xx)
		{
			double x,y,tmp,ser;
			double[]  cof = new Double[6];
			cof[0] = 76.18009172947146;
			cof[1] = -86.50532032941677;
			cof[2] = 24.01409824083091;
			cof[3] = -1.231739572450155; 
			cof[4] = 0.1208650973866179e-2;
			cof[5] = -0.5395239384953e-5;
			int j;
			y=x=xx;
			tmp=x+5.5;
			tmp -= (x+0.5)*Math.Log(tmp);
			ser=1.000000000190015;
			for (j=0;j<=5;j++) ser += cof[j]/++y;
			return -tmp+Math.Log(2.5066282746310005*ser/x);
		}
		public static double Gamma(double x)
		{
			return Math.Exp(LogGamma(x));
		}
		public static double Factorial1(int n)
		{
			int ntop = 4;
			double[] a = new Double[33];
			a[0] = 1.0;
			a[1] = 1.0;
			a[2] = 2.0;
			a[3] = 6.0;
			a[4] = 24.0;
			int j;
            if (n < 0) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Negative factorial in routine factrl");
                }
			if (n > 32) return Math.Exp(LogGamma(n+1.0));
			while (ntop < n) 
			{
				j=ntop++;
				a[ntop]=a[j]*ntop;
			}
			return a[n];
		}
		public static double Factorial2(int n)
		{
			if(n==0 || n==1) return 1.0;
			else return n*Factorial2(n-2);			
		}
		public static double LogFactorial(int n)
		{
			double[] a = new Double[101];
            if (n < 0) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Negative factorial in routine factln");
                }
			if (n <= 1) return 0.0;
			if (n <= 100) return a[n]!=0.0 ? a[n] : (a[n]=LogGamma(n+1.0));
			else return LogGamma(n+1.0);
		}
		public static double Binomial(int n, int k)
		{
			return Math.Floor(0.5+Math.Exp(LogFactorial(n)-LogFactorial(k)
				-LogFactorial(n-k)));
		}
		public static double Beta(double z, double w)
		{
			return Math.Exp(LogGamma(z)+LogGamma(w)-LogGamma(z+w));
		}
		public static double GammaRegularized(double a, double x)
		{
            if (x < 0.0 || a <= 0.0) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Invalid arguments in routine gammp");
                } 
			if (x < (a+1.0)) 
			{
				return 1.0-gser(a,x);
			} 
			else 
			{
				return gcf(a,x); 
			}
		}		
		private static double gser(double a, double x)
		{
			double gamser=0.0, gln;
		    int ITMAX = 100;
		    double EPS = 3.0e-7;

			int n;
			double sum,del,ap;

			gln=LogGamma(a);
			if (x <= 0.0) 
			{
                if (x < 0.0) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("x less than 0 in routine gser");
                    }
				gamser=0.0;
				return gamser;
			} 
			else 
			{
				ap=a;
				del=sum=1.0/a;
				for (n=1;n<=ITMAX;n++) 
				{
					++ap;
					del *= x/ap;
					sum += del;
					if (Math.Abs(del) < Math.Abs(sum)*EPS) 
					{
						gamser=sum*Math.Exp(-x+a*Math.Log(x)-gln);
						return gamser;
					}
				}
                try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("a too large, ITMAX too small in routine gser");
                }
				return gamser;
			}
		}		
		private static double gcf(double a, double x)
		{
			double gammcf, gln;
			int ITMAX = 100;
			double EPS = 3.0e-7;
		    double FPMIN = 1.0e-30;

			int i;
			double an,b,c,d,del,h;

			gln=LogGamma(a);
			b=x+1.0-a;
			c=1.0/FPMIN;
			d=1.0/b;
			h=d;
			for (i=1;i<=ITMAX;i++) 
			{
				an = -i*(i-a);
				b += 2.0;
				d=an*d+b;
				if (Math.Abs(d) < FPMIN) d=FPMIN;
				c=b+an/c;
				if (Math.Abs(c) < FPMIN) c=FPMIN;
				d=1.0/d;
				del=d*c;
				h *= del;
				if (Math.Abs(del-1.0) < EPS) break;
			}
            if (i > ITMAX) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("a too large, ITMAX too small in gcf");
                }
			gammcf=Math.Exp(-x+a*Math.Log(x)-gln)*h;
			return gammcf;
		}
		public static double Error(double x)
		{
			return x < 0.0 ? GammaRegularized(0.5,x*x)-1.0 : 1.0-GammaRegularized(0.5,x*x);
		}
		public static double ErrorComplementary(double x)
		{
			return x < 0.0 ? 2.0-GammaRegularized(0.5,x*x) : GammaRegularized(0.5,x*x);
		}
		public static double ExponentialIntegral(int n, double x)
		{
			int MAXIT = 100;
			double EULER = 0.5772156649;
			double FPMIN = 1.0e-30;
			double EPS = 1.0e-7;

			int i,ii,nm1;
			double a,b,c,d,del,fact,h,psi,ans = 0.0;

			nm1=n-1;
			if (n < 0 || x < 0.0 || (x==0.0 && (n==0 || n==1)))
				try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("bad arguments in expint");
                } 
			else 
			{
				if (n == 0) ans=Math.Exp(-x)/x;
				else 
				{
					if (x == 0.0) ans=1.0/nm1;

					else 
					{
						if (x > 1.0) 
						{
							b=x+n;
							c=1.0/FPMIN;
							d=1.0/b;
							h=d;
							for (i=1;i<=MAXIT;i++) 
							{
								a = -i*(nm1+i);
								b += 2.0;
								d=1.0/(a*d+b);
								c=b+a/c;
								del=c*d;
								h *= del;
								if (Math.Abs(del-1.0) < EPS) 
								{
									ans=h*Math.Exp(-x);
									return ans;
								}
							}
							try { throw new Science.Error(); }
                            catch (Science.Error e)
                            {
                                 e.Write("continued fraction failed in expint");
                            } 
						} 
						else 
						{
							ans = (nm1!=0 ? 1.0/nm1 : -Math.Log(x)-EULER);
							fact=1.0;
							for (i=1;i<=MAXIT;i++) 
							{
								fact *= -x/i;
								if (i != nm1) del = -fact/(i-nm1);
								else 
								{
									psi = -EULER;
									for (ii=1;ii<=nm1;ii++) psi += 1.0/ii;
									del=fact*(-Math.Log(x)+psi);
								}
								ans += del;
								if (Math.Abs(del) < Math.Abs(ans)*EPS) return ans;
							}
							try { throw new Science.Error(); }
                            catch (Science.Error e)
                            {
                                e.Write("series failed in expint");
                            } 
						}
					}
				}
			}
			return ans;
		}
		public static double ExponentialIntegralPrincipalValue(double x)
		{
			double EULER = 0.57721566;
			int MAXIT = 100;
			double FPMIN = 1.0e-30;
			double EPS = 6.0e-8;

			int k;
			double fact,prev,sum,term;

            if (x <= 0.0) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad argument in ei");
                } 
			if (x < FPMIN) return Math.Log(x)+EULER;
			if (x <= -Math.Log(EPS)) 
			{
				sum=0.0;
				fact=1.0;
				for (k=1;k<=MAXIT;k++) 
				{
					fact *= x/k;
					term=fact/k;
					sum += term;
					if (term < EPS*sum) break;
				}
                if (k > MAXIT) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Series failed in ei");
                    } 
				return sum+Math.Log(x)+EULER;
			} 
			else 
			{
				sum=0.0;
				term=1.0;
				for (k=1;k<=MAXIT;k++) 
				{
					prev=term;
					term *= k/x;
					if (term < EPS) break;
					if (term < prev) sum += term;
					else 
					{
						sum -= prev;
						break;
					}
				}
				return Math.Exp(x)*(1.0+sum)/x;
			}
		}
		public static double BetaRegularized(double a, double b, double x)
		{
			double bt;

            if (x < 0.0 || x > 1.0) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad x in routine betai");
                } 
			if (x == 0.0 || x == 1.0) bt=0.0;
			else
				bt=Math.Exp(LogGamma(a+b)-LogGamma(a)-LogGamma(b)+a*Math.Log(x)+b*Math.Log(1.0-x));
			if (x < (a+1.0)/(a+b+2.0))
				return bt*betacf(a,b,x)/a;
			else
				return 1.0-bt*betacf(b,a,1.0-x)/b;
		}		
		private static double betacf(double a, double b, double x)
		{
			int MAXIT = 100;
			double EPS = 3.0e-7;
			double FPMIN = 1.0e-30;

			int m,m2;
			double aa,c,d,del,h,qab,qam,qap;

			qab=a+b;
			qap=a+1.0;
			qam=a-1.0;
			c=1.0;
			d=1.0-qab*x/qap;
			if (Math.Abs(d) < FPMIN) d=FPMIN;
			d=1.0/d;
			h=d;
			for (m=1;m<=MAXIT;m++) 
			{
				m2=2*m;
				aa=m*(b-m)*x/((qam+m2)*(a+m2));
				d=1.0+aa*d;
				if (Math.Abs(d) < FPMIN) d=FPMIN;
				c=1.0+aa/c;
				if (Math.Abs(c) < FPMIN) c=FPMIN;
				d=1.0/d;
				h *= d*c;
				aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
				d=1.0+aa*d;
				if (Math.Abs(d) < FPMIN) d=FPMIN;
				c=1.0+aa/c;
				if (Math.Abs(c) < FPMIN) c=FPMIN;
				d=1.0/d;
				del=d*c;
				h *= del;
				if (Math.Abs(del-1.0) < EPS) break;
			}
            if (m > MAXIT) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("a or b too big, or MAXIT too small in betacf");
                } 
			return h;
		}

	/*	public static double InverseBetaRegularized(double x)
		{
			return x;
		}
		public static double InverseGammaRegularized(double x)
		{
			return x;
		}
		public static double Multinomial(double x)
		{
			return x;
		}
		public static double Pochhammer(double x)
		{
			return x;
		}
		public static double PolyGamma(double x)
		{
			return x;
		}  */
	}
}