From 932acdec792df5b98a1f5ff2a651602ff3d25db2 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 18 Sep 2019 07:53:27 -0400
Subject: small changes

---
 n_to_p.cpp      | 10 +++++++---
 percolation.cpp |  3 +--
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/n_to_p.cpp b/n_to_p.cpp
index 04c162f..df954b8 100644
--- a/n_to_p.cpp
+++ b/n_to_p.cpp
@@ -6,6 +6,7 @@
 #include <string>
 #include <vector>
 #include <iostream>
+#include <iomanip>
 
 long double lfactorial(unsigned i) { return lgamma((long double)(i + 1)); }
 
@@ -32,6 +33,7 @@ int main(int argc, char *argv[]) {
 
     unsigned line_num = 0;
     std::string line;
+    std::getline(in_file, line);
     while (std::getline(in_file, line)) {
       dat.resize(line_num + 1);
       std::istringstream iss(line);
@@ -45,7 +47,8 @@ int main(int argc, char *argv[]) {
     in_file.close();
 
     unsigned n_moms = dat.size();
-    unsigned n_ps = 10000;
+    unsigned n_ps = 100;
+    double ε = 1e-5;
 
     std::vector<std::vector<double>> output(n_moms);
 
@@ -56,7 +59,7 @@ int main(int argc, char *argv[]) {
     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 p = 0.5 - 0.5 * expl((j + 1) * logl(ε) / n_ps);
         long double mom_p = pow(1 - p, N);
         long double tot = pow(1 - p, N);
 
@@ -74,7 +77,8 @@ int main(int argc, char *argv[]) {
 
     for (unsigned i = 0; i < n_moms; i++) {
       for (unsigned j = 0; j < output[i].size(); j++) {
-        out_file << output[i][j] << " ";
+        long double p = 0.5 - 0.5 * expl((j + 1) * logl(ε) / n_ps);
+        out_file << std::fixed << std::setprecision(8) << p << " " << output[i][j] << " ";
       }
       out_file << "\n";
     }
diff --git a/percolation.cpp b/percolation.cpp
index 2f17f40..7a6dbef 100644
--- a/percolation.cpp
+++ b/percolation.cpp
@@ -30,8 +30,7 @@ public:
 
   void join(unsigned i, unsigned j) {
     for (unsigned k = 1; k <= m.size(); k++) {
-      m[k - 1] -= pow(o[i], k) / nc;
-      m[k - 1] -= pow(o[j], k) / nc;
+      m[k - 1] -= (pow(o[i], k) + pow(o[j], k)) / nc;
     }
 
     nc--;
-- 
cgit v1.2.3-70-g09d2