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 pv(
double amount,
double rate,
double period) {
22 return amount / pow(1 + rate, period);
25 double fv(
double amount,
double rate,
double period) {
26 return amount * pow(1 + rate, period);
29 double pv_coef(
double rate,
double period) {
30 return pow(1 + rate, -period);
33 double fv_coef(
double rate,
double period) {
34 return pow(1 + rate, period);
37 double apv(
double annuity,
double rate,
int period,
bool prepaid) {
38 return annuity *
apv_coef(rate, period, prepaid);
41 double afv(
double annuity,
double rate,
int period,
bool prepaid) {
42 return annuity *
afv_coef(rate, period, prepaid);
45 double apv_coef(
double rate,
int period,
bool prepaid) {
46 if (rate == 0)
return period;
47 double coef = (1 - pow(1 + rate, -period)) / rate;
48 if (prepaid) coef *= (1 + rate);
52 double afv_coef(
double rate,
int period,
bool prepaid) {
53 if (rate == 0)
return period;
54 double coef = (pow(1 + rate, period) - 1) / rate;
55 if (prepaid) coef *= (1 + rate);
59 double spv(
const Vector & amount,
double rate,
bool prepaid) {
62 Vector::const_iterator citer = amount.begin();
67 while (citer != amount.end()) {
68 pvalue +=
pv(*citer, rate, period);
75 double sfv(
const Vector & amount,
double rate,
bool prepaid) {
78 Vector::const_reverse_iterator criter = amount.rbegin();
83 while (criter != amount.rend()) {
84 fvalue +=
fv(*criter, rate, period);
91 double comp_rate(
double pval,
double fval,
double period) {
92 return pow(fval / pval, 1 / period) - 1;
97 Vector::const_iterator citer;
98 for (citer = data.begin(); citer != data.end(); ++citer) sum += *citer;
99 return sum / data.size();
105 sort(temp.begin(), temp.end());
106 return median(temp,
true);
108 if (data.size() % 2 == 0)
return (data[data.size() / 2] + data[data.size() / 2 - 1]) / 2;
109 else return data[data.size() / 2];
113 double avg =
mean(data);
115 Vector::const_iterator citer;
116 for (citer = data.begin(); citer != data.end(); ++citer) dev += (*citer - avg) * (*citer - avg);
117 if (sample)
return dev / (data.size() - 1);
118 else return dev / data.size();
122 return sqrt(
var(data, sample));
126 double avg1 =
mean(data1);
127 double avg2 =
mean(data2);
129 Vector::const_iterator citer1, citer2;
130 citer1 = data1.begin();
131 citer2 = data2.begin();
132 while (citer1 != data1.end() && citer2 != data2.end()) {
133 dev += (*citer1 - avg1) * (*citer2 - avg2);
137 if (sample)
return dev / (data1.size() - 1);
138 else return dev / data1.size();
142 return cov(data1, data2) / (
sd(data1) *
sd(data2));
147 Vector::const_iterator citer;
149 double mu =
mean(data);
150 for (citer = data.begin(); citer != data.end(); ++citer)
151 result += pow(*citer - mu, k);
154 for (citer = data.begin(); citer != data.end(); ++citer)
155 result += pow(*citer, k);
157 return result / data.size();
165 double result =
moment(data, 4) / pow(
moment(data, 2), 2);
166 if (excess) result -= 3;
170 double dnorm(
double x,
double mu,
double sigma) {
171 double z = (x - mu) / sigma;
172 return 1 / (sqrt(2 *
pi) * sigma) * exp(-0.5 * z * z);
176 double z = (x - mu) / sigma;
178 double s = z, sn = z;
179 double p = 0.5 + d * s;
184 double pn = 0.5 + d * s;
197 if (peps(p - 0.5))
return mu;
201 while (
pnorm(rx, mu, sigma, peps) < p) rx += sigma;
206 while (
pnorm(lx, mu, sigma, peps) > p) lx -= sigma;
209 if (peps(p -
pnorm(lx, mu, sigma, peps)))
return lx;
210 if (peps(p -
pnorm(rx, mu, sigma, peps)))
return rx;
212 double mid = (lx + rx) / 2;
213 double pmid =
pnorm(mid, mu, sigma, peps);
214 if (eps(rx - lx) || peps(p - pmid))
return mid;
215 if (pmid > p) rx = mid;
220 double dt(
double x,
int n) {
221 return gamma((n + 1) / 2.0) / (sqrt(n *
pi) *
gamma(n / 2.0)) * pow(1 + x * x / n, (n + 1) / -2.0);
225 if (x == 0)
return 0.5;
228 double p =
integrate(bind2nd(ptr_fun<double, int, double>(
dt), n), region, eps);
229 return p + (1 - p) / 2;
233 double p =
integrate(bind2nd(ptr_fun<double, int, double>(
dt), n), region, eps);
240 if (peps(p - 0.5))
return 0;
244 while (
pt(rx, n, peps) < p) rx += 1;
249 while (
pt(lx, n, peps) > p) lx -= 1;
252 if (peps(p -
pt(lx, n, peps)))
return lx;
253 if (peps(p -
pt(rx, n, peps)))
return rx;
255 double mid = (lx + rx) / 2;
256 double pmid =
pt(mid, n, peps);
257 if (eps(rx - lx) || peps(p - pmid))
return mid;
258 if (pmid > p) rx = mid;
263 double pois(
double lmd,
int k,
bool cum) {
264 if (!cum)
return pow(lmd, k) * exp(-lmd) /
fac(k);
266 for (
int i = 0; i <= k; ++i) result +=
pois(lmd, i,
false);
270 double binom(
int n,
double p,
int k,
bool cum) {
271 if (!cum)
return comb(n, k) * pow(p, k) * pow(1 - p, n - k);
273 for (
int i = 0; i <= k; ++i) result +=
binom(n, p, i,
false);
277 double expo(
double theta,
double x) {
278 if (x > 0)
return 1 - exp(-1 / theta * x);
283 if (!(x > 0))
throw invalid_argument(
"Invalid argument.");
286 if (modf(n, &ipart) == 0 && n <= UINT_MAX)
return fac((
unsigned int) n);
287 double p[] = {1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 1.208650973866179e-3, -5.395239384953e-6};
289 for (
int i = 1; i <= 6; ++i) s += p[i] / (x + i);
290 return sqrt(2 *
pi) / x * s * pow(x + 5.5, x + 0.5) * exp(-x - 5.5);
303 if (from > to)
return vec;
307 }
while (!(from > to));
310 if (from < to)
return vec;
314 }
while (!(from < to));
321 if (count < 2)
return vec;
323 vec.push_back(start);
327 double delta = (end - start) / (count - 1);
328 vec.push_back(start);
329 for (
int i = 0; i < count - 2; ++i) {
331 vec.push_back(start);
339 return *(max_element(data.begin(), data.end()));
343 return *(min_element(data.begin(), data.end()));
347 return accumulate(data.begin(), data.end(), 0.0);
351 return accumulate(data.begin(), data.end(), 1.0, multiplies<double>());
356 Vector::const_iterator citer1, citer2;
357 citer1 = vec1.begin();
358 citer2 = vec2.begin();
359 while (citer1 != vec1.end() && citer2 != vec2.end()) {
360 result.push_back(*citer1 + *citer2);
369 Vector::const_iterator citer;
370 for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer + scalar);
375 return add(vec1, vec2);
379 return add(vec, scalar);
383 return add(vec, scalar);
388 Vector::const_iterator citer1, citer2;
389 citer1 = vec1.begin();
390 citer2 = vec2.begin();
391 while (citer1 != vec1.end() && citer2 != vec2.end()) {
392 result.push_back(*citer1 - *citer2);
401 Vector::const_iterator citer;
402 for (citer = vec.begin(); citer != vec.end(); ++citer) {
403 double element = *citer - scalar;
404 if (!dir) element = -element;
405 result.push_back(element);
411 return sub(vec1, vec2);
415 return sub(vec, scalar);
419 return sub(vec, scalar,
false);
424 Vector::const_iterator citer1, citer2;
425 citer1 = vec1.begin();
426 citer2 = vec2.begin();
427 while (citer1 != vec1.end() && citer2 != vec2.end()) {
428 result.push_back(*citer1 * *citer2);
437 Vector::const_iterator citer;
438 for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer * scalar);
443 return mul(vec1, vec2);
447 return mul(vec, scalar);
451 return mul(vec, scalar);
456 Vector::const_iterator citer1, citer2;
457 citer1 = vec1.begin();
458 citer2 = vec2.begin();
459 while (citer1 != vec1.end() && citer2 != vec2.end()) {
460 result.push_back(*citer1 / *citer2);
469 Vector::const_iterator citer;
470 for (citer = vec.begin(); citer != vec.end(); ++citer) {
471 if (dir) result.push_back(*citer / scalar);
472 else result.push_back(scalar / *citer);
478 return div(vec1, vec2);
482 return div(vec, scalar);
486 return div(vec, scalar,
false);
491 Vector::const_iterator citer1, citer2;
492 citer1 = vec1.begin();
493 citer2 = vec2.begin();
494 while (citer1 != vec1.end() && citer2 != vec2.end()) {
495 result += *citer1 * *citer2;
504 vec.push_back(vec1[1] * vec2[2] - vec1[2] * vec2[1]);
505 vec.push_back(vec1[2] * vec2[0] - vec1[0] * vec2[2]);
506 vec.push_back(vec1[0] * vec2[1] - vec1[1] * vec2[0]);
516 Vector::size_type n = vec.size();
517 Vector::size_type last = n - 1;
518 for (Vector::size_type i = 0; i < n; ++i) {
520 if (i != last) oss << sep;
522 string result = oss.str();
523 if (delim.size() == 2) result = delim[0] + result + delim[1];
531 for (
int i = 0; i < n; ++i) vec.push_back(va_arg(argptr,
double));
539 if (getline(ins, line)) {
542 while ((i = line.find(delim)) != string::npos) line.replace(i, delim.size(),
" ");
544 istringstream iss(line);
548 data.push_back(item);
554 Table load(istream & ins,
int nrow,
int ncol,
string delim) {
556 for (
int i = 0; i < nrow; ++i) {
558 for (
int j = 0; j < ncol; ++j) data[j].push_back(row[j]);
563 void save(ostream & outs,
const Vector & data,
string delim) {
564 Vector::const_iterator citer, last;
565 last = data.end() - 1;
566 int prec = outs.precision();
568 for (citer = data.begin(); citer != data.end(); ++citer) {
570 if (citer != last) outs << delim;
573 outs.precision(prec);
576 void save(ostream & outs,
const Table & data,
string delim) {
577 int nrow = data[0].size();
578 int ncol = data.size();
579 for (
int i = 0; i < nrow; ++i) {
581 for (
int j = 0; j < ncol; ++j) row.push_back(data[j][i]);
582 save(outs, row, delim);
596 double randf(
double low,
double high) {
597 return rand() / (double) RAND_MAX * (high - low) + low;
602 for (
int i = 0; i < n; ++i) result.push_back(
randf(low, high));
607 return (
int) floor(rand() / ((
double) RAND_MAX + 1) * ((
double) high - low + 1) + low);
612 for (
int i = 0; i < n; ++i) result.push_back(
randi(low, high));
617 return !(
randf(0, 1) > p);
621 srand((
unsigned int) time(NULL));
628 double rnorm(
double mu,
double sigma) {
629 double u =
randf(0, 1);
630 double v =
randf(0, 1);
631 double z = sqrt(-2 * log(u)) * cos(2 *
pi * v);
632 return mu + sigma * z;
637 for (
int i = 0; i < n; ++i) result.push_back(
rnorm(mu, sigma));
642 double L = exp(-lmd);
654 for (
int i = 0; i < n; ++i) result.push_back(
rpois(lmd));
660 for (
int i = 0; i < n; ++i)
661 if (
randp(p)) k += 1;
667 for (
int i = 0; i < count; ++i) result.push_back(
rbinom(n, p));
672 return -theta * log(
randf(0, 1));
677 for (
int i = 0; i < n; ++i) result.push_back(
rexp(theta));
681 double prec(
double num,
int ndec) {
684 decpart = modf(num, &intpart);
686 ss.precision(ndec + 1);
689 return intpart + decpart - 1;
694 Vector::const_iterator citer;
695 for (citer = vec.begin(); citer != vec.end(); ++citer)
696 result.push_back(
prec(*citer, ndec));
700 double fac(
unsigned int n) {
701 if (n == 0)
return 1;
708 return prec(result, 0);
712 double perm(
unsigned int m,
unsigned int n) {
719 return prec(result, 0);
722 double comb(
unsigned int m,
unsigned int n) {
723 if (m - n < n)
return comb(m, m - n);
727 unsigned int gcd(
unsigned int a,
unsigned int b) {
736 unsigned int lcm(
unsigned int a,
unsigned int b) {
737 return a * b /
gcd(a, b);