using System;

namespace Psychlops
{

	public static class Math
	{
		public const double PI = 3.14159265, E = 2.718281828459045, LOG2E = 1.44269504088896340736;
		public static Random random_generator;
		static Math()
		{
			random_generator = new Random();
		}

		public static T max<T>(T val1, T val2) where T : IComparable
		{
			return val1.CompareTo(val2) > 0 ? val1 : val2;
		}
		public static T min<T>(T val1, T val2) where T : IComparable
		{
			return val1.CompareTo(val2) < 0 ? val1 : val2;
		}
		public static void shuffle<X>(X[] array, int n)
		{
			int a;
			X tmp;
			for(int i = 1; i < n; i++){
				a = random(i + 1);
				tmp = array[i];
				array[i] = array[a];
				array[a] = tmp;
			}
		}


		public static double mod(double lhs, double rhs)
		{
			return lhs - System.Math.Floor(lhs/rhs)*rhs;
		}
		public static double abs(double x)
		{
			return System.Math.Abs(x);
		}
		public static double floor(double x)
		{
			return System.Math.Floor(x);
		}
		public static double ceil(double x)
		{
			return System.Math.Ceiling(x);
		}
		public static double sin(double x)
		{
			return System.Math.Sin(x);
		}
		public static double cos(double x)
		{
			return System.Math.Cos(x);
		}
		public static double tan(double x)
		{
			return System.Math.Tan(x);
		}
		public static double asin(double x)
		{
			return System.Math.Asin(x);
		}
		public static double acos(double x)
		{
			return System.Math.Acos(x);
		}
		public static double atan(double x)
		{
			return System.Math.Atan(x);
		}
		public static double atan(double y, double x)
		{
			return System.Math.Atan2(y, x);
		}
		public static double atan2(double y, double x)
		{
			return System.Math.Atan2(y, x);
		}
		public static double sqrt(double x)
		{
			return System.Math.Sqrt(x);
		}
		public static double pow(double x, double y)
		{
			return System.Math.Pow(x,y);
		}
		public static double exp(double x)
		{
			return System.Math.Exp(x);
		}
		public static double log(double x)
		{
			return System.Math.Log(x);
		}
		public static double log2(double val)
		{
			return log(val) * LOG2E;
		}
		public static double round(double val)
		{
			return System.Math.Round(val);
		}

		public static double radius(double x, double y)
		{
			return System.Math.Sqrt(x * x + y * y);
		}

		public static double random()
		{
			return (random_generator.NextDouble());
		}
		public static int random(int x)
		{
			return (int)((random_generator.NextDouble()) * x);
		}
		public static double random(double x)
		{
			return (random_generator.NextDouble()) * x;
		}
		public static double random(double x, double y)
		{
			return (random_generator.NextDouble()) * (y-x) + x;
		}


		public static double gaussian(double x, double sigma)
		{
			return exp(- (x*x) / (2*sigma*sigma));
		}

		public static double normalDistibution(double x, double mu, double sigma)
		{
			return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);
		}

		public static double cumulativeNormalDistibution(double x, double mu, double sigma)
		{
			return .5 + .5*Internal.GammaFunction.erf( (x-mu)/(sigma*sqrt(2.0) ) );
		}

	}

	namespace Internal
	{
		public static class GammaFunction
		{
			/************ loggamma(x) -- gamma.c より再掲 *************/

			static readonly double PI      = 3.14159265358979324;  /* $\pi$ */
			static readonly double LOG_2PI = 1.83787706640934548;  /* $\log 2\pi$ */
			static readonly double N       = 8;

			static readonly double B0  = 1            ;     /* 以下はBernoulli数 */
			static readonly double B1  = (-1.0 / 2.0);
			static readonly double B2  = ( 1.0 / 6.0);
			static readonly double B4  = (-1.0 / 30.0);
			static readonly double B6  = ( 1.0 / 42.0);
			static readonly double B8  = (-1.0 / 30.0);
			static readonly double B10 = ( 5.0 / 66.0);
			static readonly double B12 = (-691.0 / 2730.0);
			static readonly double B14 = ( 7.0 / 6.0);
			static readonly double B16 = (-3617.0 / 510.0);

			public static double loggamma(double x)  /* ガンマ関数の対数 */
			{
				double v, w;

				v = 1;
				while (x < N) {  v *= x;  x++;  }
				w = 1 / (x * x);
				return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w
							+ (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w
							+ (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w
							+ (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x
							+ 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);
			}

			public static double p_gamma(double a, double x, double loggamma_a)  /* 本文参照 */
			{
				int k;
				double result, term, previous;

				if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
				if (x == 0)     return 0;
				result = term = Math.exp(a * Math.log(x) - x - loggamma_a) / a;
				for (k = 1; k < 1000; k++) {
					term *= x / (a + k);
					previous = result;  result += term;
					if (result == previous) return result;
				}
				//throw new Exception("p_gamma(): the sequence is not convergent.");
				return result;
			}

			public static double q_gamma(double a, double x, double loggamma_a)  /* 本文参照 */
			{
				int k;
				double result, w, temp, previous;
				double la = 1, lb = 1 + x - a;  /* Laguerreの多項式 */

				if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
				w = Math.exp(a * Math.log(x) - x - loggamma_a);
				result = w / lb;
				for (k = 2; k < 1000; k++) {
					temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
					la = lb;  lb = temp;
					w *= (k - 1 - a) / k;
					temp = w / (la * lb);
					previous = result;  result += temp;
					if (result == previous) return result;
				}
				//throw new Exception("q_gamma(): the sequence is not convergent.");
				return result;
			}

			public static double p_chisq(double chisq, int df)  /* カイ2乗分布の下側確率 */
			{
				return p_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));
			}

			public static double q_chisq(double chisq, int df)  /* カイ2乗分布の上側確率 */
			{
				return q_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));
			}

			static readonly double LOG_PI = 1.14472988584940017;  /* $\log_e \pi$ */

			public static double erf(double x)  /* Gaussの誤差関数 ${\rm erf}(x)$ */
			{
				if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI / 2);
				else        return - p_gamma(0.5, x * x, LOG_PI / 2);
			}

			public static double erfc(double x)  /* $1 - {\rm erf}(x)$ */
			{
				if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI / 2);
				else        return  1 + p_gamma(0.5, x * x, LOG_PI / 2);
			}

			public static double p_normal(double x)  /* 標準正規分布の下側確率 */
			{
				if (x >= 0) return
					0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));
				else return
					0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);
			}

			public static double q_normal(double x)  /* 標準正規分布の上側確率 */
			{
				if (x >= 0) return
					0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);
				else return
					0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));
			}
		}
	}
}