1 files changed, 87 insertions, 0 deletions
diff --git a/n_to_p.cpp b/n_to_p.cpp
new file mode 100644
index 0000000..04c162f
--- /dev/null
+++ b/n_to_p.cpp
@@ -0,0 +1,87 @@
+
+#include <cinttypes>
+#include <cmath>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <iostream>
+
+long double lfactorial(unsigned i) { return lgamma((long double)(i + 1)); }
+
+long double lbinom(unsigned N, unsigned n) {
+  return lfactorial(N) - lfactorial(n) - lfactorial(N-n);
+}
+
+long double B(unsigned N, unsigned n, long double p) {
+  return expl(lbinom(N, n) + log(p) * n + log(1-p) * (N - n));
+}
+
+int main(int argc, char *argv[]) {
+  std::vector<std::string> all_args;
+
+  if (argc > 1) {
+    all_args.assign(argv + 1, argv + argc);
+  }
+
+  for (std::string filename : all_args) {
+    std::ifstream in_file;
+    in_file.open(filename);
+
+    std::vector<std::vector<double>> dat;
+
+    unsigned line_num = 0;
+    std::string line;
+    while (std::getline(in_file, line)) {
+      dat.resize(line_num + 1);
+      std::istringstream iss(line);
+      double a;
+      while (iss >> a) {
+        dat[line_num].push_back(a);
+      }
+      line_num++;
+    }
+
+    in_file.close();
+
+    unsigned n_moms = dat.size();
+    unsigned n_ps = 10000;
+
+    std::vector<std::vector<double>> output(n_moms);
+
+    for (unsigned i = 0; i < n_moms; i++) {
+      output[i].resize(n_ps);
+    }
+
+    for (unsigned i = 0; i < n_moms; i++) {
+      unsigned N = dat[i].size();
+      for (unsigned j = 0; j < n_ps; j++) {
+        long double p = (1.0 / (n_ps + 1)) * (j + 1);
+        long double mom_p = pow(1 - p, N);
+        long double tot = pow(1 - p, N);
+
+        for (unsigned n = 1; n <= N; n++) {
+          mom_p += B(N, n, p) * dat[i][n - 1];
+          tot += B(N, n, p);
+        }
+
+        output[i][j] = mom_p;
+      }
+    }
+
+    std::ofstream out_file;
+    out_file.open(filename + ".ps");
+
+    for (unsigned i = 0; i < n_moms; i++) {
+      for (unsigned j = 0; j < output[i].size(); j++) {
+        out_file << output[i][j] << " ";
+      }
+      out_file << "\n";
+    }
+
+    out_file.close();
+  }
+
+  return 0;
+}
+