11 #include "mathkit.hpp"
16 bool Epsilon::operator()(
double val)
const {
17 if (_eps == 0)
return val == 0;
18 else return abs(val) < _eps;
21 double Epsilon::operator()()
const {
25 double pv(
double amount,
double rate,
double period) {
26 return amount / pow(1 + rate, period);
29 double fv(
double amount,
double rate,
double period) {
30 return amount * pow(1 + rate, period);
33 double pv_coef(
double rate,
double period) {
34 return pow(1 + rate, -period);
37 double fv_coef(
double rate,
double period) {
38 return pow(1 + rate, period);
41 double apv(
double annuity,
double rate,
int period,
bool prepaid) {
42 return annuity *
apv_coef(rate, period, prepaid);
45 double afv(
double annuity,
double rate,
int period,
bool prepaid) {
46 return annuity *
afv_coef(rate, period, prepaid);
49 double apv_coef(
double rate,
int period,
bool prepaid) {
50 if (rate == 0)
return period;
51 double coef = (1 - pow(1 + rate, -period)) / rate;
52 if (prepaid) coef *= (1 + rate);
56 double afv_coef(
double rate,
int period,
bool prepaid) {
57 if (rate == 0)
return period;
58 double coef = (pow(1 + rate, period) - 1) / rate;
59 if (prepaid) coef *= (1 + rate);
63 double spv(
const Vector & amount,
double rate,
bool prepaid) {
66 Vector::const_iterator citer = amount.begin();
71 while (citer != amount.end()) {
72 pvalue +=
pv(*citer, rate, period);
79 double sfv(
const Vector & amount,
double rate,
bool prepaid) {
82 Vector::const_reverse_iterator criter = amount.rbegin();
87 while (criter != amount.rend()) {
88 fvalue +=
fv(*criter, rate, period);
95 double comp_rate(
double pval,
double fval,
double period) {
96 return pow(fval / pval, 1 / period) - 1;
101 Vector::const_iterator citer;
102 for (citer = data.begin(); citer != data.end(); ++citer) sum += *citer;
103 return sum / data.size();
109 sort(temp.begin(), temp.end());
110 return median(temp,
true);
112 if (data.size() % 2 == 0)
return (data[data.size() / 2] + data[data.size() / 2 - 1]) / 2;
113 else return data[data.size() / 2];
117 double avg =
mean(data);
119 Vector::const_iterator citer;
120 for (citer = data.begin(); citer != data.end(); ++citer) dev += (*citer - avg) * (*citer - avg);
121 if (sample)
return dev / (data.size() - 1);
122 else return dev / data.size();
126 return sqrt(
var(data, sample));
130 double avg1 =
mean(data1);
131 double avg2 =
mean(data2);
133 Vector::const_iterator citer1, citer2;
134 citer1 = data1.begin();
135 citer2 = data2.begin();
136 while (citer1 != data1.end() && citer2 != data2.end()) {
137 dev += (*citer1 - avg1) * (*citer2 - avg2);
141 if (sample)
return dev / (data1.size() - 1);
142 else return dev / data1.size();
146 return cov(data1, data2) / (
sd(data1) *
sd(data2));
151 Vector::const_iterator citer;
153 double mu =
mean(data);
154 for (citer = data.begin(); citer != data.end(); ++citer)
155 result += pow(*citer - mu, k);
158 for (citer = data.begin(); citer != data.end(); ++citer)
159 result += pow(*citer, k);
161 return result / data.size();
169 double result =
moment(data, 4) / pow(
moment(data, 2), 2);
170 if (excess) result -= 3;
174 double dexp(
double theta,
double x) {
176 else return 1 / theta * exp(-1 / theta * x);
179 double pexp(
double theta,
double x) {
180 if (x > 0)
return 1 - exp(-1 / theta * x);
185 if (p == 0)
return 0;
186 double lx, rx = theta;
187 while (
pexp(theta, rx) < p) rx += theta;
189 if (p ==
pexp(theta, lx))
return lx;
190 if (p ==
pexp(theta, rx))
return rx;
192 double mid = (lx + rx) / 2;
193 double pmid =
pexp(theta, mid);
194 if (eps(rx - lx) || p == pmid)
return mid;
195 if (pmid > p) rx = mid;
200 double dnorm(
double x,
double mu,
double sigma) {
201 double z = (x - mu) / sigma;
202 return 1 / (sqrt(2 *
pi) * sigma) * exp(-0.5 * z * z);
206 double z = (x - mu) / sigma;
208 double s = z, sn = z;
209 double p = 0.5 + d * s;
214 double pn = 0.5 + d * s;
227 if (peps(p - 0.5))
return mu;
231 while (
pnorm(rx, mu, sigma, peps) < p) rx += sigma;
236 while (
pnorm(lx, mu, sigma, peps) > p) lx -= sigma;
239 if (peps(p -
pnorm(lx, mu, sigma, peps)))
return lx;
240 if (peps(p -
pnorm(rx, mu, sigma, peps)))
return rx;
242 double mid = (lx + rx) / 2;
243 double pmid =
pnorm(mid, mu, sigma, peps);
244 if (eps(rx - lx) || peps(p - pmid))
return mid;
245 if (pmid > p) rx = mid;
250 double dt(
double x,
int n) {
251 return gamma((n + 1) / 2.0) / (sqrt(n *
pi) *
gamma(n / 2.0)) * pow(1 + x * x / n, (n + 1) / -2.0);
255 if (x == 0)
return 0.5;
258 double p =
integrate(bind2nd(ptr_fun<double, int, double>(
dt), n), region, eps);
259 return p + (1 - p) / 2;
263 double p =
integrate(bind2nd(ptr_fun<double, int, double>(
dt), n), region, eps);
270 if (peps(p - 0.5))
return 0;
274 while (
pt(rx, n, peps) < p) rx += 1;
279 while (
pt(lx, n, peps) > p) lx -= 1;
282 if (peps(p -
pt(lx, n, peps)))
return lx;
283 if (peps(p -
pt(rx, n, peps)))
return rx;
285 double mid = (lx + rx) / 2;
286 double pmid =
pt(mid, n, peps);
287 if (eps(rx - lx) || peps(p - pmid))
return mid;
288 if (pmid > p) rx = mid;
294 return x < 0 ? 0 : pow(x, k / 2.0 - 1) * exp(-x / 2) / pow(2, k / 2.0) /
gamma(k / 2.0);
298 if (x <= 0)
return 0;
299 else return igamma(k / 2.0, x / 2) /
gamma(k / 2.0);
303 if (p == 0)
return 0;
306 while (
pchisq(rx, k) < p) rx += k;
308 if (peps(p -
pchisq(lx, k)))
return lx;
309 if (peps(p -
pchisq(rx, k)))
return rx;
311 double mid = (lx + rx) / 2;
312 double pmid =
pchisq(mid, k);
313 if (eps(rx - lx) || peps(p - pmid))
return mid;
314 if (pmid > p) rx = mid;
319 double df(
double x,
int d1,
int d2) {
321 else return pow(
double(d1) / d2, d1 / 2.0) * pow(x, d1 / 2.0 - 1) * pow(1 + d1 * x / d2, -(d1 + d2) / 2.0) /
beta(d1 / 2.0, d2 / 2.0);
324 double pf(
double x,
int d1,
int d2) {
325 if (x <= 0)
return 0;
326 return ibeta(d1 * x / (d1 * x + d2), d1 / 2.0, d2 / 2.0);
329 double qf(
double p,
int d1,
int d2,
const Epsilon & eps) {
330 if (p == 0)
return 0;
333 while (
pf(rx, d1, d2) < p) rx += 1;
335 if (peps(p -
pf(lx, d1, d2)))
return lx;
336 if (peps(p -
pf(rx, d1, d2)))
return rx;
338 double mid = (lx + rx) / 2;
339 double pmid =
pf(mid, d1, d2);
340 if (eps(rx - lx) || peps(p - pmid))
return mid;
341 if (pmid > p) rx = mid;
346 double pois(
double lmd,
int k,
bool cum) {
347 if (!cum)
return pow(lmd, k) * exp(-lmd) /
fac(k);
349 for (
int i = 0; i <= k; ++i) result +=
pois(lmd, i,
false);
353 double binom(
int n,
double p,
int k,
bool cum) {
354 if (!cum)
return comb(n, k) * pow(p, k) * pow(1 - p, n - k);
356 for (
int i = 0; i <= k; ++i) result +=
binom(n, p, i,
false);
360 double nb(
int k,
double p,
int x,
bool cum) {
361 if (!cum)
return comb(k + x - 1, k - 1) * pow(p, k) * pow(1 - p, x);
363 for (
int i = 0; i <= x; ++i) result +=
nb(k, p, i,
false);
367 double hg(
int s,
int k,
int n,
int x,
bool cum) {
368 if (!cum)
return comb(k, x) *
comb(n - k, s - x) /
comb(n, s);
370 for (
int i = 0; i <= x; ++i) result +=
hg(s, k, n, i,
false);
379 if (!(x > 0))
throw invalid_argument(
"Invalid argument.");
380 static double coef[] = {57.1562356658629235, -59.5979603554754912, 14.1360979747417471,
381 -0.491913816097620199, 0.339946499848118887e-4, 0.465236289270485756e-4,
382 -0.983744753048795646e-4, 0.158088703224912494e-3, -0.210264441724104883e-3,
383 0.217439618115212643e-3, -0.164318106536763890e-3, 0.844182239838527433e-4,
384 -0.261908384015814087e-4, 0.368991826595316234e-5};
385 double temp, s, xn = x;
386 temp = x + 671.0 / 128;
387 temp = (x + 0.5) * log(temp) - temp;
388 s = 0.999999999999997092;
389 for (
int i = 0; i < 14; ++i) s += coef[i] / ++xn;
390 return temp + log(2.5066282746310005 * s / x);
394 double g = pow(x, s) / s;
397 double gn = pow(x, s + k) /
fac(k) / (s + k);
398 if (k % 2 != 0) gn = -gn;
409 double beta(
double x,
double y) {
410 return exp(
lbeta(x, y));
418 double f = pow(x, a) * pow(1 - x, b) / a /
beta(a, b);
419 double s = exp(
lbeta(a + 1, 1) -
lbeta(a + b, 1)) * x;
421 double ib = f * (1 + s);
425 sn += exp(
lbeta(a + 1, n) -
lbeta(a + b, n)) * pow(x, n);
426 double ibn = f * (1 + sn);
441 Vector::const_iterator citer;
442 for (citer = xdata.begin(); citer != xdata.end(); ++citer)
444 Vector dev = ydata - yline;
445 Vector::iterator iter;
446 for (iter = dev.begin(); iter != dev.end(); ++iter) *iter *= *iter;
447 double sse =
sum(dev);
448 regline.
variance = sse / (xdata.size() - 2);
449 double sst =
var(ydata) * (ydata.size() - 1);
455 int pn = data.size();
456 if (n < 0 || n > pn)
throw invalid_argument(
"Invalid argument.");
457 if (n == pn)
return data;
460 for (
int i = 0; i < n; ++i) {
461 int order =
randi(0, pn - i);
463 for (
int j = 0; j < pn; ++j) {
464 if (find(mask.begin(), mask.end(), j) == mask.end()) --order;
466 samp.push_back(data[j]);
478 if (from > to)
return vec;
482 }
while (!(from > to));
485 if (from < to)
return vec;
489 }
while (!(from < to));
496 if (count < 2)
return vec;
498 vec.push_back(start);
502 double delta = (end - start) / (count - 1);
503 vec.push_back(start);
504 for (
int i = 0; i < count - 2; ++i) {
506 vec.push_back(start);
514 return *(max_element(data.begin(), data.end()));
518 return *(min_element(data.begin(), data.end()));
522 return accumulate(data.begin(), data.end(), 0.0);
526 return accumulate(data.begin(), data.end(), 1.0, multiplies<double>());
531 Vector::const_iterator citer1, citer2;
532 citer1 = vec1.begin();
533 citer2 = vec2.begin();
534 while (citer1 != vec1.end() && citer2 != vec2.end()) {
535 result.push_back(*citer1 + *citer2);
544 Vector::const_iterator citer;
545 for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer + scalar);
550 return add(vec1, vec2);
554 return add(vec, scalar);
558 return add(vec, scalar);
563 Vector::const_iterator citer1, citer2;
564 citer1 = vec1.begin();
565 citer2 = vec2.begin();
566 while (citer1 != vec1.end() && citer2 != vec2.end()) {
567 result.push_back(*citer1 - *citer2);
576 Vector::const_iterator citer;
577 for (citer = vec.begin(); citer != vec.end(); ++citer) {
578 double element = *citer - scalar;
579 if (!dir) element = -element;
580 result.push_back(element);
586 return sub(vec1, vec2);
590 return sub(vec, scalar);
594 return sub(vec, scalar,
false);
599 Vector::const_iterator citer1, citer2;
600 citer1 = vec1.begin();
601 citer2 = vec2.begin();
602 while (citer1 != vec1.end() && citer2 != vec2.end()) {
603 result.push_back(*citer1 * *citer2);
612 Vector::const_iterator citer;
613 for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer * scalar);
618 return mul(vec1, vec2);
622 return mul(vec, scalar);
626 return mul(vec, scalar);
631 Vector::const_iterator citer1, citer2;
632 citer1 = vec1.begin();
633 citer2 = vec2.begin();
634 while (citer1 != vec1.end() && citer2 != vec2.end()) {
635 result.push_back(*citer1 / *citer2);
644 Vector::const_iterator citer;
645 for (citer = vec.begin(); citer != vec.end(); ++citer) {
646 if (dir) result.push_back(*citer / scalar);
647 else result.push_back(scalar / *citer);
653 return div(vec1, vec2);
657 return div(vec, scalar);
661 return div(vec, scalar,
false);
666 Vector::const_iterator citer1, citer2;
667 citer1 = vec1.begin();
668 citer2 = vec2.begin();
669 while (citer1 != vec1.end() && citer2 != vec2.end()) {
670 result += *citer1 * *citer2;
679 vec.push_back(vec1[1] * vec2[2] - vec1[2] * vec2[1]);
680 vec.push_back(vec1[2] * vec2[0] - vec1[0] * vec2[2]);
681 vec.push_back(vec1[0] * vec2[1] - vec1[1] * vec2[0]);
690 Vector::const_iterator citer;
691 for (citer = vec.begin(); citer != vec.end(); ++citer)
692 if (!eps(*citer))
return false;
699 Vector::size_type n = vec.size();
700 Vector::size_type last = n - 1;
701 for (Vector::size_type i = 0; i < n; ++i) {
703 if (i != last) oss << sep;
705 string result = oss.str();
706 if (delim.size() == 2) result = delim[0] + result + delim[1];
714 for (
int i = 0; i < n; ++i) vec.push_back(va_arg(argptr,
double));
722 if (getline(ins, line)) {
725 while ((i = line.find(delim)) != string::npos) line.replace(i, delim.size(),
" ");
727 istringstream iss(line);
731 data.push_back(item);
737 Table load(istream & ins,
int nrow,
int ncol,
string delim) {
739 for (
int i = 0; i < nrow; ++i) {
741 for (
int j = 0; j < ncol; ++j) data[j].push_back(row[j]);
746 void save(ostream & outs,
const Vector & data,
string delim) {
747 Vector::const_iterator citer, last;
748 last = data.end() - 1;
749 int prec = outs.precision();
751 for (citer = data.begin(); citer != data.end(); ++citer) {
753 if (citer != last) outs << delim;
756 outs.precision(prec);
759 void save(ostream & outs,
const Table & data,
string delim) {
760 int nrow = data[0].size();
761 int ncol = data.size();
762 for (
int i = 0; i < nrow; ++i) {
764 for (
int j = 0; j < ncol; ++j) row.push_back(data[j][i]);
765 save(outs, row, delim);
779 double randf(
double low,
double high) {
780 return rand() / (double) RAND_MAX * (high - low) + low;
785 for (
int i = 0; i < n; ++i) result.push_back(
randf(low, high));
790 return (
int) floor(rand() / ((
double) RAND_MAX + 1) * ((
double) high - low + 1) + low);
795 for (
int i = 0; i < n; ++i) result.push_back(
randi(low, high));
800 return !(
randf(0, 1) > p);
804 srand((
unsigned int) time(NULL));
811 double rnorm(
double mu,
double sigma) {
812 double u =
randf(0, 1);
813 double v =
randf(0, 1);
814 double z = sqrt(-2 * log(u)) * cos(2 *
pi * v);
815 return mu + sigma * z;
820 for (
int i = 0; i < n; ++i) result.push_back(
rnorm(mu, sigma));
825 return -theta * log(
randf(0, 1));
830 for (
int i = 0; i < n; ++i) result.push_back(
rexp(theta));
835 double L = exp(-lmd);
847 for (
int i = 0; i < n; ++i) result.push_back(
rpois(lmd));
853 for (
int i = 0; i < n; ++i)
854 if (
randp(p)) k += 1;
860 for (
int i = 0; i < count; ++i) result.push_back(
rbinom(n, p));
864 int rnb(
int k,
double p) {
867 if (
randp(p)) ++s;
else ++f;
873 for (
int i = 0; i < count; ++i) result.push_back(
rnb(k, p));
877 int rhg(
int s,
int k,
int n) {
879 fill_n(pop.begin(), k, 1.0);
880 random_shuffle(pop.begin(), pop.end());
882 return count(samp.begin(), samp.end(), 1.0);
887 for (
int i = 0; i < count; ++i) result.push_back(
rhg(s, k, n));
891 double prec(
double num,
int ndec) {
894 decpart = modf(num, &intpart);
895 double temp = decpart * 10;
896 while (temp < 1 && temp > -1 && ndec > 0) {
900 if (ndec == 0)
return intpart;
905 return intpart + decpart;
910 Vector::const_iterator citer;
911 for (citer = vec.begin(); citer != vec.end(); ++citer)
912 result.push_back(
prec(*citer, ndec));
917 if (n < 0)
throw invalid_argument(
"Non-negative integer needed.");
918 if (n == 0)
return 1;
925 return prec(result, 0);
930 if (m < 0 || n < 0)
return 0;
937 return prec(result, 0);
941 if (m < 0 || n < 0)
return 0;
943 if (m - n < n)
return comb(m, m - n);
947 unsigned int gcd(
unsigned int a,
unsigned int b) {
956 unsigned int lcm(
unsigned int a,
unsigned int b) {
957 return a * b /
gcd(a, b);