From 9cf2a9f889afb3ac26d2b3466dc7a3adef96bbdc Mon Sep 17 00:00:00 2001
From: vub <Vitalik Buterin>
Date: Tue, 18 Apr 2017 07:26:05 -0400
Subject: [PATCH] Added c++ implementation

---
 erasure_code/ec65536/poly_utils.cpp  | 101 +++++++++++++++++++++++++++
 erasure_code/ec65536/poly_utils.py   |   2 +
 erasure_code/ec65536/test_ec65536.py |  10 ++-
 3 files changed, 112 insertions(+), 1 deletion(-)
 create mode 100644 erasure_code/ec65536/poly_utils.cpp

diff --git a/erasure_code/ec65536/poly_utils.cpp b/erasure_code/ec65536/poly_utils.cpp
new file mode 100644
index 0000000..35b39d5
--- /dev/null
+++ b/erasure_code/ec65536/poly_utils.cpp
@@ -0,0 +1,101 @@
+#include "stdlib.h"
+#include "stdio.h"
+#include <iostream>
+#include <vector>
+#include <deque>
+
+using namespace std;
+
+vector<int> glogtable(65536, 0);
+vector<int> gexptable(196608, 0);
+
+void initialize_tables() {
+    int v = 1;
+    for (int i = 0; i < 65536; i++) {
+        glogtable[v] = i;
+        gexptable[i] = v;
+        gexptable[i + 65535] = v;
+        gexptable[i + 131070] = v;
+        if (v & 32768)
+            v = (v * 2) ^ v ^ 103425;
+        else
+            v = (v * 2) ^ v;
+    }
+}
+
+int eval_poly_at(vector<int> poly, int x) {
+    if (x == 0)
+        return poly[0];
+    int logx = glogtable[x];
+    int y = 0;
+    for (int i = 0; i < poly.size(); i++) {
+        if (poly[i])
+            y ^= gexptable[(logx * i + glogtable[poly[i]]) % 65535];
+    }
+    return y;
+}
+
+vector<int> lagrange_interp(vector<int> ys, vector<int> xs) {
+    deque<int> root;
+    root.push_back(1);
+    for (int i = 0; i < xs.size(); i++) {
+        int logx = glogtable[xs[i]];
+        root.push_front(0);
+        for (int j = 0; j < root.size() - 1; j++) {
+            if (root[j + 1] and xs[i])
+                root[j] ^= gexptable[glogtable[root[j+1]] + logx];
+        }
+    }
+    vector<vector<int> > nums;
+    for (int i = 0; i < xs.size(); i++) {
+        vector<int> output(root.size() - 1);
+        output[root.size() - 2] = 1;
+        int logx = glogtable[xs[i]];
+        for (int j = root.size() - 2; j > 0; j--) { 
+            if (output[j] and xs[i])
+                output[j - 1] = root[j] ^ gexptable[glogtable[output[j]] + logx];
+            else
+                output[j - 1] = root[j];
+        }
+        nums.push_back(output);
+    }
+    vector<int> denoms;
+    for (int i = 0; i < xs.size(); i++) {
+        denoms.push_back(eval_poly_at(nums[i], xs[i]));
+    }
+    vector<int> b(xs.size());
+    for (int i = 0; i < xs.size(); i++) {
+        int log_yslice = glogtable[ys[i]] - glogtable[denoms[i]] + 65535;
+        for (int j = 0; j < xs.size(); j++) {
+            if(nums[i][j] and ys[i])
+                b[j] ^= gexptable[glogtable[nums[i][j]] + log_yslice];
+        }
+    }
+    return b;
+}
+
+
+int main() {
+    initialize_tables();
+    //int myxs[] = {1, 2, 3, 4};
+    //std::vector<int> xs (myxs, myxs + sizeof(myxs) / sizeof(int) );
+    vector<int> xs(4096);
+    vector<int> ys(4096);
+    for (int v = 0; v < 4096; v++) {
+        xs[v] = v;
+        ys[v] = (v * v) % 65536;
+    }
+    vector<int> poly = lagrange_interp(ys, xs);
+    unsigned int o = 0;
+    for (int i = 4096; i < 8192; i++) {
+        o += eval_poly_at(poly, i);
+    }
+    cout << o;
+    //cout << eval_poly_at(poly, 0) << " " << ys[0] << "\n";
+    //cout << eval_poly_at(poly, 134) << " " << ys[134] << "\n";
+    //cout << eval_poly_at(poly, 375) << " " << ys[375] << "\n";
+    //int o;
+    //for (int i = 0; i < 524288; i ++)
+    //    o += eval_poly_at(poly, i % 65536);
+    //std::cout << o;
+}
diff --git a/erasure_code/ec65536/poly_utils.py b/erasure_code/ec65536/poly_utils.py
index 8469ad9..ecc7d9d 100644
--- a/erasure_code/ec65536/poly_utils.py
+++ b/erasure_code/ec65536/poly_utils.py
@@ -71,6 +71,7 @@ def lagrange_interp(pieces, xs):
             if root[j+1] and x:
                 root[j] ^= gexptable[glogtable[root[j+1]] + logx]
     assert len(root) == len(pieces) + 1
+    print(root)
     # Generate per-value numerator polynomials, eg. for x=x2,
     # (x - x1) * (x - x3) * ... * (x - xn), by dividing the master
     # polynomial back by each x coordinate
@@ -85,6 +86,7 @@ def lagrange_interp(pieces, xs):
                 output[j-1] = root[j]
         assert len(output) == len(pieces)
         nums.append(output)
+    print(nums)
     # Generate denominators by evaluating numerator polys at each x
     denoms = [eval_poly_at(nums[i], xs[i]) for i in range(len(xs))]
     # Generate output polynomial, which is the sum of the per-value numerator
diff --git a/erasure_code/ec65536/test_ec65536.py b/erasure_code/ec65536/test_ec65536.py
index 0c74729..b8a70a6 100644
--- a/erasure_code/ec65536/test_ec65536.py
+++ b/erasure_code/ec65536/test_ec65536.py
@@ -1,19 +1,27 @@
 import ec65536
 import rlp
+import time
 
 # 12.8 kilobyte test string
 testdata = 'the cow jumped over the moon!!! ' * 400
 
+t1 = time.time()
 prover = ec65536.Prover(testdata)
-print("Created prover")
+t2 = time.time()
+print("Created prover in %.2f sec" % (t2 - t1))
 
 assert ec65536.verify_proof(prover.merkle_root, prover.prove(13), 13)
+t3 = time.time()
+print("Created and verified a proof in %.2f sec" % (t3 - t2))
 
 proofs = [prover.prove(i) for i in range(0, prover.length, 2)]
 print("Created merkle proofs")
 
+t4 = time.time()
 print("Starting to attempt fill")
 response = ec65536.fill(prover.merkle_root, prover.length // 2, proofs, list(range(0, prover.length, 2)))
+t5 = time.time()
+print("Completed fill in %.2f sec" % (t5 - t4))
 assert response is not False
 assert b''.join(response)[:len(rlp.encode(testdata))] == rlp.encode(testdata)
 print("Fill successful")