mathkit  2.0
 All Classes Namespaces Functions Variables Typedefs
mathkit.cpp
1 #include <ctime>
2 #include <cstdlib>
3 #include <cstdarg>
4 #include <climits>
5 #include <sstream>
6 #include <numeric>
7 #include <stdexcept>
8 #include <algorithm>
9 #include <functional>
10 
11 #include "mathkit.hpp"
12 
13 using namespace std;
14 
15 namespace mathkit {
16  bool Epsilon::operator()(double val) const {
17  if (_eps == 0) return val == 0;
18  else return abs(val) < _eps;
19  }
20 
21  double Epsilon::operator()() const {
22  return _eps;
23  }
24 
25  double pv(double amount, double rate, double period) {
26  return amount / pow(1 + rate, period);
27  }
28 
29  double fv(double amount, double rate, double period) {
30  return amount * pow(1 + rate, period);
31  }
32 
33  double pv_coef(double rate, double period) {
34  return pow(1 + rate, -period);
35  }
36 
37  double fv_coef(double rate, double period) {
38  return pow(1 + rate, period);
39  }
40 
41  double apv(double annuity, double rate, int period, bool prepaid) {
42  return annuity * apv_coef(rate, period, prepaid);
43  }
44 
45  double afv(double annuity, double rate, int period, bool prepaid) {
46  return annuity * afv_coef(rate, period, prepaid);
47  }
48 
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);
53  return coef;
54  }
55 
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);
60  return coef;
61  }
62 
63  double spv(const Vector & amount, double rate, bool prepaid) {
64  double pvalue = 0;
65  int period = 1;
66  Vector::const_iterator citer = amount.begin();
67  if (prepaid) {
68  pvalue += *citer;
69  ++citer;
70  }
71  while (citer != amount.end()) {
72  pvalue += pv(*citer, rate, period);
73  ++period;
74  ++citer;
75  }
76  return pvalue;
77  }
78 
79  double sfv(const Vector & amount, double rate, bool prepaid) {
80  double fvalue = 0;
81  int period = 1;
82  Vector::const_reverse_iterator criter = amount.rbegin();
83  if (!prepaid) {
84  fvalue += *criter;
85  ++criter;
86  }
87  while (criter != amount.rend()) {
88  fvalue += fv(*criter, rate, period);
89  ++period;
90  ++criter;
91  }
92  return fvalue;
93  }
94 
95  double comp_rate(double pval, double fval, double period) {
96  return pow(fval / pval, 1 / period) - 1;
97  }
98 
99  double mean(const Vector & data) {
100  double sum = 0;
101  Vector::const_iterator citer;
102  for (citer = data.begin(); citer != data.end(); ++citer) sum += *citer;
103  return sum / data.size();
104  }
105 
106  double median(const Vector & data, bool sorted) {
107  if (!sorted) {
108  Vector temp = data;
109  sort(temp.begin(), temp.end());
110  return median(temp, true);
111  }
112  if (data.size() % 2 == 0) return (data[data.size() / 2] + data[data.size() / 2 - 1]) / 2;
113  else return data[data.size() / 2];
114  }
115 
116  double var(const Vector & data, bool sample) {
117  double avg = mean(data);
118  double dev = 0;
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();
123  }
124 
125  double sd(const Vector & data, bool sample) {
126  return sqrt(var(data, sample));
127  }
128 
129  double cov(const Vector & data1, const Vector & data2, bool sample) {
130  double avg1 = mean(data1);
131  double avg2 = mean(data2);
132  double dev = 0;
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);
138  ++citer1;
139  ++citer2;
140  }
141  if (sample) return dev / (data1.size() - 1);
142  else return dev / data1.size();
143  }
144 
145  double cor(const Vector & data1, const Vector & data2) {
146  return cov(data1, data2) / (sd(data1) * sd(data2));
147  }
148 
149  double moment(const Vector & data, int k, bool central) {
150  double result = 0;
151  Vector::const_iterator citer;
152  if (central) {
153  double mu = mean(data);
154  for (citer = data.begin(); citer != data.end(); ++citer)
155  result += pow(*citer - mu, k);
156  }
157  else {
158  for (citer = data.begin(); citer != data.end(); ++citer)
159  result += pow(*citer, k);
160  }
161  return result / data.size();
162  }
163 
164  double skew(const Vector & data) {
165  return moment(data, 3) / pow(moment(data, 2), 1.5);
166  }
167 
168  double kurt(const Vector & data, bool excess) {
169  double result = moment(data, 4) / pow(moment(data, 2), 2);
170  if (excess) result -= 3;
171  return result;
172  }
173 
174  double dexp(double theta, double x) {
175  if (x < 0) return 0;
176  else return 1 / theta * exp(-1 / theta * x);
177  }
178 
179  double pexp(double theta, double x) {
180  if (x > 0) return 1 - exp(-1 / theta * x);
181  else return 0;
182  }
183 
184  double qexp(double theta, double p, const Epsilon & eps) {
185  if (p == 0) return 0;
186  double lx, rx = theta;
187  while (pexp(theta, rx) < p) rx += theta;
188  lx = rx - theta;
189  if (p == pexp(theta, lx)) return lx;
190  if (p == pexp(theta, rx)) return rx;
191  while (true) {
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;
196  else lx = mid;
197  }
198  }
199 
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);
203  }
204 
205  double pnorm(double x, double mu, double sigma, const Epsilon & eps) {
206  double z = (x - mu) / sigma;
207  double d = dnorm(z);
208  double s = z, sn = z;
209  double p = 0.5 + d * s;
210  int i = 3;
211  while (true) {
212  sn *= z * z / i;
213  s += sn;
214  double pn = 0.5 + d * s;
215  if (eps(pn - p)) {
216  p = pn;
217  break;
218  }
219  p = pn;
220  i += 2;
221  }
222  return p;
223  }
224 
225  double qnorm(double p, double mu, double sigma, const Epsilon & eps) {
226  Epsilon peps(1e-15);
227  if (peps(p - 0.5)) return mu;
228  double lx, rx;
229  if (p > 0.5) {
230  rx = sigma;
231  while (pnorm(rx, mu, sigma, peps) < p) rx += sigma;
232  lx = rx - sigma;
233  }
234  else {
235  lx = -sigma;
236  while (pnorm(lx, mu, sigma, peps) > p) lx -= sigma;
237  rx = lx + sigma;
238  }
239  if (peps(p - pnorm(lx, mu, sigma, peps))) return lx;
240  if (peps(p - pnorm(rx, mu, sigma, peps))) return rx;
241  while (true) {
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;
246  else lx = mid;
247  }
248  }
249 
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);
252  }
253 
254  double pt(double x, int n, const Epsilon & eps) {
255  if (x == 0) return 0.5;
256  if (x > 0) {
257  Pair region(-x, x);
258  double p = integrate(bind2nd(ptr_fun<double, int, double>(dt), n), region, eps);
259  return p + (1 - p) / 2;
260  }
261  else {
262  Pair region(x, -x);
263  double p = integrate(bind2nd(ptr_fun<double, int, double>(dt), n), region, eps);
264  return (1 - p) / 2;
265  }
266  }
267 
268  double qt(double p, int n, const Epsilon & eps) {
269  Epsilon peps(1e-10);
270  if (peps(p - 0.5)) return 0;
271  double lx, rx;
272  if (p > 0.5) {
273  rx = 1;
274  while (pt(rx, n, peps) < p) rx += 1;
275  lx = rx - 1;
276  }
277  else {
278  lx = -1;
279  while (pt(lx, n, peps) > p) lx -= 1;
280  rx = lx + 1;
281  }
282  if (peps(p - pt(lx, n, peps))) return lx;
283  if (peps(p - pt(rx, n, peps))) return rx;
284  while (true) {
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;
289  else lx = mid;
290  }
291  }
292 
293  double dchisq(double x, int k) {
294  return x < 0 ? 0 : pow(x, k / 2.0 - 1) * exp(-x / 2) / pow(2, k / 2.0) / gamma(k / 2.0);
295  }
296 
297  double pchisq(double x, int k) {
298  if (x <= 0) return 0;
299  else return igamma(k / 2.0, x / 2) / gamma(k / 2.0);
300  }
301 
302  double qchisq(double p, int k, const Epsilon & eps) {
303  if (p == 0) return 0;
304  Epsilon peps(1e-10);
305  double lx, rx = k;
306  while (pchisq(rx, k) < p) rx += k;
307  lx = rx - k;
308  if (peps(p - pchisq(lx, k))) return lx;
309  if (peps(p - pchisq(rx, k))) return rx;
310  while (true) {
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;
315  else lx = mid;
316  }
317  }
318 
319  double df(double x, int d1, int d2) {
320  if (x < 0) return 0;
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);
322  }
323 
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);
327  }
328 
329  double qf(double p, int d1, int d2, const Epsilon & eps) {
330  if (p == 0) return 0;
331  Epsilon peps(1e-10);
332  double lx, rx = 1;
333  while (pf(rx, d1, d2) < p) rx += 1;
334  lx = rx - 1;
335  if (peps(p - pf(lx, d1, d2))) return lx;
336  if (peps(p - pf(rx, d1, d2))) return rx;
337  while (true) {
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;
342  else lx = mid;
343  }
344  }
345 
346  double pois(double lmd, int k, bool cum) {
347  if (!cum) return pow(lmd, k) * exp(-lmd) / fac(k);
348  double result = 0;
349  for (int i = 0; i <= k; ++i) result += pois(lmd, i, false);
350  return result;
351  }
352 
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);
355  double result = 0;
356  for (int i = 0; i <= k; ++i) result += binom(n, p, i, false);
357  return result;
358  }
359 
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);
362  double result = 0;
363  for (int i = 0; i <= x; ++i) result += nb(k, p, i, false);
364  return result;
365  }
366 
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);
369  double result = 0;
370  for (int i = 0; i <= x; ++i) result += hg(s, k, n, i, false);
371  return result;
372  }
373 
374  double gamma(double x) {
375  return exp(lgamma(x));
376  }
377 
378  double lgamma(double x) {
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);
391  }
392 
393  double igamma(double s, double x, const Epsilon & eps) {
394  double g = pow(x, s) / s;
395  int k = 1;
396  while (true) {
397  double gn = pow(x, s + k) / fac(k) / (s + k);
398  if (k % 2 != 0) gn = -gn;
399  if (eps(gn)) {
400  g += gn;
401  break;
402  }
403  g += gn;
404  k += 1;
405  }
406  return g;
407  }
408 
409  double beta(double x, double y) {
410  return exp(lbeta(x, y));
411  }
412 
413  double lbeta(double x, double y) {
414  return lgamma(x) + lgamma(y) - lgamma(x + y);
415  }
416 
417  double ibeta(double x, double a, double b, const Epsilon & eps) {
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;
420  double sn = s;
421  double ib = f * (1 + s);
422  int n = 1;
423  while (true) {
424  n += 1;
425  sn += exp(lbeta(a + 1, n) - lbeta(a + b, n)) * pow(x, n);
426  double ibn = f * (1 + sn);
427  if (eps(ibn - ib)) {
428  ib = ibn;
429  break;
430  }
431  ib = ibn;
432  }
433  return ib;
434  }
435 
436  LineParam linregress(const Vector & xdata, const Vector & ydata) {
437  LineParam regline;
438  regline.slope = cov(xdata, ydata) / var(xdata);
439  regline.intercept = mean(ydata) - regline.slope * mean(xdata);
440  Vector yline;
441  Vector::const_iterator citer;
442  for (citer = xdata.begin(); citer != xdata.end(); ++citer)
443  yline.push_back(regline.intercept + regline.slope * *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);
450  regline.det_coef = 1 - sse / sst;
451  return regline;
452  }
453 
454  Vector sample(const Vector & data, int n) {
455  int pn = data.size();
456  if (n < 0 || n > pn) throw invalid_argument("Invalid argument.");
457  if (n == pn) return data;
458  Vector samp;
459  vector<int> mask;
460  for (int i = 0; i < n; ++i) {
461  int order = randi(0, pn - i);
462  order += 1;
463  for (int j = 0; j < pn; ++j) {
464  if (find(mask.begin(), mask.end(), j) == mask.end()) --order;
465  if (order == 0) {
466  samp.push_back(data[j]);
467  mask.push_back(j);
468  break;
469  }
470  }
471  }
472  return samp;
473  }
474 
475  Vector seq(double from, double to, double step) {
476  Vector vec;
477  if (step > 0) {
478  if (from > to) return vec;
479  do {
480  vec.push_back(from);
481  from += step;
482  } while (!(from > to));
483  }
484  else if (step < 0) {
485  if (from < to) return vec;
486  do {
487  vec.push_back(from);
488  from += step;
489  } while (!(from < to));
490  }
491  return vec;
492  }
493 
494  Vector linspace(double start, double end, int count) {
495  Vector vec;
496  if (count < 2) return vec;
497  if (count == 2) {
498  vec.push_back(start);
499  vec.push_back(end);
500  }
501  else {
502  double delta = (end - start) / (count - 1);
503  vec.push_back(start);
504  for (int i = 0; i < count - 2; ++i) {
505  start += delta;
506  vec.push_back(start);
507  }
508  vec.push_back(end);
509  }
510  return vec;
511  }
512 
513  double max(const Vector & data) {
514  return *(max_element(data.begin(), data.end()));
515  }
516 
517  double min(const Vector & data) {
518  return *(min_element(data.begin(), data.end()));
519  }
520 
521  double sum(const Vector & data) {
522  return accumulate(data.begin(), data.end(), 0.0);
523  }
524 
525  double prod(const Vector & data) {
526  return accumulate(data.begin(), data.end(), 1.0, multiplies<double>());
527  }
528 
529  Vector add(const Vector & vec1, const Vector & vec2) {
530  Vector result;
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);
536  ++citer1;
537  ++citer2;
538  }
539  return result;
540  }
541 
542  Vector add(const Vector & vec, double scalar) {
543  Vector result;
544  Vector::const_iterator citer;
545  for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer + scalar);
546  return result;
547  }
548 
549  Vector operator+(const Vector & vec1, const Vector & vec2) {
550  return add(vec1, vec2);
551  }
552 
553  Vector operator+(const Vector & vec, double scalar) {
554  return add(vec, scalar);
555  }
556 
557  Vector operator+(double scalar, const Vector & vec) {
558  return add(vec, scalar);
559  }
560 
561  Vector sub(const Vector & vec1, const Vector & vec2) {
562  Vector result;
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);
568  ++citer1;
569  ++citer2;
570  }
571  return result;
572  }
573 
574  Vector sub(const Vector & vec, double scalar, bool dir) {
575  Vector result;
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);
581  }
582  return result;
583  }
584 
585  Vector operator-(const Vector & vec1, const Vector & vec2) {
586  return sub(vec1, vec2);
587  }
588 
589  Vector operator-(const Vector & vec, double scalar) {
590  return sub(vec, scalar);
591  }
592 
593  Vector operator-(double scalar, const Vector & vec) {
594  return sub(vec, scalar, false);
595  }
596 
597  Vector mul(const Vector & vec1, const Vector & vec2) {
598  Vector result;
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);
604  ++citer1;
605  ++citer2;
606  }
607  return result;
608  }
609 
610  Vector mul(const Vector & vec, double scalar) {
611  Vector result;
612  Vector::const_iterator citer;
613  for (citer = vec.begin(); citer != vec.end(); ++citer) result.push_back(*citer * scalar);
614  return result;
615  }
616 
617  Vector operator*(const Vector & vec1, const Vector & vec2) {
618  return mul(vec1, vec2);
619  }
620 
621  Vector operator*(const Vector & vec, double scalar) {
622  return mul(vec, scalar);
623  }
624 
625  Vector operator*(double scalar, const Vector & vec) {
626  return mul(vec, scalar);
627  }
628 
629  Vector div(const Vector & vec1, const Vector & vec2) {
630  Vector result;
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);
636  ++citer1;
637  ++citer2;
638  }
639  return result;
640  }
641 
642  Vector div(const Vector & vec, double scalar, bool dir) {
643  Vector result;
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);
648  }
649  return result;
650  }
651 
652  Vector operator/(const Vector & vec1, const Vector & vec2) {
653  return div(vec1, vec2);
654  }
655 
656  Vector operator/(const Vector & vec, double scalar) {
657  return div(vec, scalar);
658  }
659 
660  Vector operator/(double scalar, const Vector & vec) {
661  return div(vec, scalar, false);
662  }
663 
664  double dot_prod(const Vector & vec1, const Vector & vec2) {
665  double result = 0;
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;
671  ++citer1;
672  ++citer2;
673  }
674  return result;
675  }
676 
677  Vector cross_prod(const Vector & vec1, const Vector & vec2) {
678  Vector vec;
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]);
682  return vec;
683  }
684 
685  double norm(const Vector & vec) {
686  return sqrt(dot_prod(vec, vec));
687  }
688 
689  bool zero(const Vector & vec, const Epsilon & eps) {
690  Vector::const_iterator citer;
691  for (citer = vec.begin(); citer != vec.end(); ++citer)
692  if (!eps(*citer)) return false;
693  return true;
694  }
695 
696  string to_str(const Vector & vec, string sep, string delim) {
697  ostringstream oss;
698  oss.precision(15);
699  Vector::size_type n = vec.size();
700  Vector::size_type last = n - 1;
701  for (Vector::size_type i = 0; i < n; ++i) {
702  oss << vec[i];
703  if (i != last) oss << sep;
704  }
705  string result = oss.str();
706  if (delim.size() == 2) result = delim[0] + result + delim[1];
707  return result;
708  }
709 
710  Vector make_vec(int n, ...) {
711  Vector vec;
712  va_list argptr;
713  va_start(argptr, n);
714  for (int i = 0; i < n; ++i) vec.push_back(va_arg(argptr, double));
715  va_end(argptr);
716  return vec;
717  }
718 
719  Vector load(istream & ins, string delim) {
720  Vector data;
721  string line;
722  if (getline(ins, line)) {
723  if (delim != " ") {
724  string::size_type i;
725  while ((i = line.find(delim)) != string::npos) line.replace(i, delim.size(), " ");
726  }
727  istringstream iss(line);
728  while (!iss.eof()) {
729  double item;
730  iss >> item;
731  data.push_back(item);
732  }
733  }
734  return data;
735  }
736 
737  Table load(istream & ins, int nrow, int ncol, string delim) {
738  Table data(ncol, Vector());
739  for (int i = 0; i < nrow; ++i) {
740  Vector row = load(ins, delim);
741  for (int j = 0; j < ncol; ++j) data[j].push_back(row[j]);
742  }
743  return data;
744  }
745 
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();
750  outs.precision(15);
751  for (citer = data.begin(); citer != data.end(); ++citer) {
752  outs << *citer;
753  if (citer != last) outs << delim;
754  }
755  outs << endl;
756  outs.precision(prec);
757  }
758 
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) {
763  Vector row;
764  for (int j = 0; j < ncol; ++j) row.push_back(data[j][i]);
765  save(outs, row, delim);
766  }
767  }
768 
769  ostream & operator<<(ostream & outs, const Vector & vec) {
770  save(outs, vec);
771  return outs;
772  }
773 
774  istream & operator>>(istream & ins, Vector & vec) {
775  vec = load(ins);
776  return ins;
777  }
778 
779  double randf(double low, double high) {
780  return rand() / (double) RAND_MAX * (high - low) + low;
781  }
782 
783  Vector randf(double low, double high, int n) {
784  Vector result;
785  for (int i = 0; i < n; ++i) result.push_back(randf(low, high));
786  return result;
787  }
788 
789  int randi(int low, int high) {
790  return (int) floor(rand() / ((double) RAND_MAX + 1) * ((double) high - low + 1) + low);
791  }
792 
793  Vector randi(int low, int high, int n) {
794  Vector result;
795  for (int i = 0; i < n; ++i) result.push_back(randi(low, high));
796  return result;
797  }
798 
799  bool randp(double p) {
800  return !(randf(0, 1) > p);
801  }
802 
803  void randseed() {
804  srand((unsigned int) time(NULL));
805  }
806 
807  void randseed(unsigned int s) {
808  srand(s);
809  }
810 
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;
816  }
817 
818  Vector rnorm(int n, double mu, double sigma) {
819  Vector result;
820  for (int i = 0; i < n; ++i) result.push_back(rnorm(mu, sigma));
821  return result;
822  }
823 
824  double rexp(double theta) {
825  return -theta * log(randf(0, 1));
826  }
827 
828  Vector rexp(double theta, int n) {
829  Vector result;
830  for (int i = 0; i < n; ++i) result.push_back(rexp(theta));
831  return result;
832  }
833 
834  int rpois(double lmd) {
835  double L = exp(-lmd);
836  double p = 1;
837  int k = 0;
838  do {
839  k += 1;
840  p *= randf(0, 1);
841  } while (p > L);
842  return k - 1;
843  }
844 
845  Vector rpois(double lmd, int n) {
846  Vector result;
847  for (int i = 0; i < n; ++i) result.push_back(rpois(lmd));
848  return result;
849  }
850 
851  int rbinom(int n, double p) {
852  int k = 0;
853  for (int i = 0; i < n; ++i)
854  if (randp(p)) k += 1;
855  return k;
856  }
857 
858  Vector rbinom(int n, double p, int count) {
859  Vector result;
860  for (int i = 0; i < count; ++i) result.push_back(rbinom(n, p));
861  return result;
862  }
863 
864  int rnb(int k, double p) {
865  int s = 0, f = 0;
866  while (s < k)
867  if (randp(p)) ++s; else ++f;
868  return f;
869  }
870 
871  Vector rnb(int k, double p, int count) {
872  Vector result;
873  for (int i = 0; i < count; ++i) result.push_back(rnb(k, p));
874  return result;
875  }
876 
877  int rhg(int s, int k, int n) {
878  Vector pop(n, 0);
879  fill_n(pop.begin(), k, 1.0);
880  random_shuffle(pop.begin(), pop.end());
881  Vector samp = sample(pop, s);
882  return count(samp.begin(), samp.end(), 1.0);
883  }
884 
885  Vector rhg(int s, int k, int n, int count) {
886  Vector result;
887  for (int i = 0; i < count; ++i) result.push_back(rhg(s, k, n));
888  return result;
889  }
890 
891  double prec(double num, int ndec) {
892  double intpart;
893  double decpart;
894  decpart = modf(num, &intpart);
895  double temp = decpart * 10;
896  while (temp < 1 && temp > -1 && ndec > 0) {
897  ndec -= 1;
898  temp *= 10;
899  }
900  if (ndec == 0) return intpart;
901  stringstream ss;
902  ss.precision(ndec);
903  ss << decpart;
904  ss >> decpart;
905  return intpart + decpart;
906  }
907 
908  Vector prec(const Vector & vec, int ndec) {
909  Vector result;
910  Vector::const_iterator citer;
911  for (citer = vec.begin(); citer != vec.end(); ++citer)
912  result.push_back(prec(*citer, ndec));
913  return result;
914  }
915 
916  double fac(int n) {
917  if (n < 0) throw invalid_argument("Non-negative integer needed.");
918  if (n == 0) return 1;
919  else {
920  double result = 1;
921  while (n > 0) {
922  result *= n;
923  --n;
924  }
925  return prec(result, 0);
926  }
927  }
928 
929  double perm(int m, int n) {
930  if (m < 0 || n < 0) return 0;
931  double result = 1;
932  while (n > 0) {
933  result *= m;
934  --m;
935  --n;
936  }
937  return prec(result, 0);
938  }
939 
940  double comb(int m, int n) {
941  if (m < 0 || n < 0) return 0;
942  if (n > m) return 0;
943  if (m - n < n) return comb(m, m - n);
944  return prec(perm(m, n) / fac(n), 0);
945  }
946 
947  unsigned int gcd(unsigned int a, unsigned int b) {
948  while (b != 0) {
949  int temp = b;
950  b = a % b;
951  a = temp;
952  }
953  return a;
954  }
955 
956  unsigned int lcm(unsigned int a, unsigned int b) {
957  return a * b / gcd(a, b);
958  }
959 }