6 public static class Math
\r
8 public const double PI = 3.14159265, E = 2.718281828459045, LOG2E = 1.44269504088896340736;
\r
9 public static Random random_generator;
\r
12 random_generator = new Random();
\r
15 public static T max<T>(T val1, T val2) where T : IComparable
\r
17 return val1.CompareTo(val2) > 0 ? val1 : val2;
\r
19 public static T min<T>(T val1, T val2) where T : IComparable
\r
21 return val1.CompareTo(val2) < 0 ? val1 : val2;
\r
23 public static void shuffle<X>(X[] array, int n)
\r
27 for(int i = 1; i < n; i++){
\r
30 array[i] = array[a];
\r
36 public static double mod(double lhs, double rhs)
\r
38 return lhs - System.Math.Floor(lhs/rhs)*rhs;
\r
40 public static double abs(double x)
\r
42 return System.Math.Abs(x);
\r
44 public static double floor(double x)
\r
46 return System.Math.Floor(x);
\r
48 public static double ceil(double x)
\r
50 return System.Math.Ceiling(x);
\r
52 public static double sin(double x)
\r
54 return System.Math.Sin(x);
\r
56 public static double cos(double x)
\r
58 return System.Math.Cos(x);
\r
60 public static double tan(double x)
\r
62 return System.Math.Tan(x);
\r
64 public static double asin(double x)
\r
66 return System.Math.Asin(x);
\r
68 public static double acos(double x)
\r
70 return System.Math.Acos(x);
\r
72 public static double atan(double x)
\r
74 return System.Math.Atan(x);
\r
76 public static double atan(double y, double x)
\r
78 return System.Math.Atan2(y, x);
\r
80 public static double atan2(double y, double x)
\r
82 return System.Math.Atan2(y, x);
\r
84 public static double sqrt(double x)
\r
86 return System.Math.Sqrt(x);
\r
88 public static double pow(double x, double y)
\r
90 return System.Math.Pow(x,y);
\r
92 public static double exp(double x)
\r
94 return System.Math.Exp(x);
\r
96 public static double log(double x)
\r
98 return System.Math.Log(x);
\r
100 public static double log2(double val)
\r
102 return log(val) * LOG2E;
\r
104 public static double round(double val)
\r
106 return System.Math.Round(val);
\r
109 public static double radius(double x, double y)
\r
111 return System.Math.Sqrt(x * x + y * y);
\r
114 public static double random()
\r
116 return (random_generator.NextDouble());
\r
118 public static int random(int x)
\r
120 return (int)((random_generator.NextDouble()) * x);
\r
122 public static double random(double x)
\r
124 return (random_generator.NextDouble()) * x;
\r
126 public static double random(double x, double y)
\r
128 return (random_generator.NextDouble()) * (y-x) + x;
\r
132 public static double gaussian(double x, double sigma)
\r
134 return exp(- (x*x) / (2*sigma*sigma));
\r
137 public static double normalDistibution(double x, double mu, double sigma)
\r
139 return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);
\r
142 public static double cumulativeNormalDistibution(double x, double mu, double sigma)
\r
144 return .5 + .5*Internal.GammaFunction.erf( (x-mu)/(sigma*sqrt(2.0) ) );
\r
151 public static class GammaFunction
\r
153 /************ loggamma(x) -- gamma.c より再掲 *************/
\r
155 static readonly double PI = 3.14159265358979324; /* $\pi$ */
\r
156 static readonly double LOG_2PI = 1.83787706640934548; /* $\log 2\pi$ */
\r
157 static readonly double N = 8;
\r
159 static readonly double B0 = 1 ; /* 以下はBernoulli数 */
\r
160 static readonly double B1 = (-1.0 / 2.0);
\r
161 static readonly double B2 = ( 1.0 / 6.0);
\r
162 static readonly double B4 = (-1.0 / 30.0);
\r
163 static readonly double B6 = ( 1.0 / 42.0);
\r
164 static readonly double B8 = (-1.0 / 30.0);
\r
165 static readonly double B10 = ( 5.0 / 66.0);
\r
166 static readonly double B12 = (-691.0 / 2730.0);
\r
167 static readonly double B14 = ( 7.0 / 6.0);
\r
168 static readonly double B16 = (-3617.0 / 510.0);
\r
170 public static double loggamma(double x) /* ガンマ関数の対数 */
\r
175 while (x < N) { v *= x; x++; }
\r
177 return ((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w
\r
178 + (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w
\r
179 + (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w
\r
180 + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x
\r
181 + 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);
\r
184 public static double p_gamma(double a, double x, double loggamma_a) /* 本文参照 */
\r
187 double result, term, previous;
\r
189 if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
\r
190 if (x == 0) return 0;
\r
191 result = term = Math.exp(a * Math.log(x) - x - loggamma_a) / a;
\r
192 for (k = 1; k < 1000; k++) {
\r
193 term *= x / (a + k);
\r
194 previous = result; result += term;
\r
195 if (result == previous) return result;
\r
197 //throw new Exception("p_gamma(): the sequence is not convergent.");
\r
201 public static double q_gamma(double a, double x, double loggamma_a) /* 本文参照 */
\r
204 double result, w, temp, previous;
\r
205 double la = 1, lb = 1 + x - a; /* Laguerreの多項式 */
\r
207 if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
\r
208 w = Math.exp(a * Math.log(x) - x - loggamma_a);
\r
210 for (k = 2; k < 1000; k++) {
\r
211 temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
\r
212 la = lb; lb = temp;
\r
213 w *= (k - 1 - a) / k;
\r
214 temp = w / (la * lb);
\r
215 previous = result; result += temp;
\r
216 if (result == previous) return result;
\r
218 //throw new Exception("q_gamma(): the sequence is not convergent.");
\r
222 public static double p_chisq(double chisq, int df) /* カイ2乗分布の下側確率 */
\r
224 return p_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));
\r
227 public static double q_chisq(double chisq, int df) /* カイ2乗分布の上側確率 */
\r
229 return q_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));
\r
232 static readonly double LOG_PI = 1.14472988584940017; /* $\log_e \pi$ */
\r
234 public static double erf(double x) /* Gaussの誤差関数 ${\rm erf}(x)$ */
\r
236 if (x >= 0) return p_gamma(0.5, x * x, LOG_PI / 2);
\r
237 else return - p_gamma(0.5, x * x, LOG_PI / 2);
\r
240 public static double erfc(double x) /* $1 - {\rm erf}(x)$ */
\r
242 if (x >= 0) return q_gamma(0.5, x * x, LOG_PI / 2);
\r
243 else return 1 + p_gamma(0.5, x * x, LOG_PI / 2);
\r
246 public static double p_normal(double x) /* 標準正規分布の下側確率 */
\r
249 0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));
\r
251 0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);
\r
254 public static double q_normal(double x) /* 標準正規分布の上側確率 */
\r
257 0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);
\r
259 0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));
\r