diff --git a/erasure_code/ec65536/poly_utils.cpp b/erasure_code/ec65536/poly_utils.cpp
index d0b38b8..8501d41 100644
--- a/erasure_code/ec65536/poly_utils.cpp
+++ b/erasure_code/ec65536/poly_utils.cpp
@@ -9,6 +9,8 @@ using namespace std;
 vector<int> glogtable(65536, 0);
 vector<int> gexptable(196608, 0);
 
+const int ROOT_CUTOFF = 32;
+
 void initialize_tables() {
     int v = 1;
     for (int i = 0; i < 65536; i++) {
@@ -35,57 +37,189 @@ int eval_poly_at(vector<int> poly, int x) {
     return y;
 }
 
-vector<int> lagrange_interp(vector<int> ys, vector<int> xs) {
-    int xs_size = xs.size();
-    vector<int> root(xs_size + 1);
-    root[xs_size] = 1;
-    for (int i = 0; i < xs_size; i++) {
+int eval_log_poly_at(vector<int> poly, int x) {
+    if (x == 0)
+        return poly[0] == 65537 ? 0 : gexptable[poly[0]];
+    int logx = glogtable[x];
+    int y = 0;
+    for (int i = 0; i < poly.size(); i++) {
+        if (poly[i] != 65537)
+            y ^= gexptable[(logx * i + poly[i]) % 65535];
+    }
+    return y;
+}
+
+// Compute the product of two (equal length) polynomials. Takes ~O(N ** 1.59) time.
+vector<int> karatsuba_mul(vector<int> p, vector<int> q) {
+    int L = p.size();
+    if (L <= 64) {
+        vector<int> o(L * 2);
+        for (int i = 0; i < L; i++) {
+            for (int j = 0; j < L; j++) {
+                if (p[i] && q[j])
+                    o[i + j] ^= gexptable[glogtable[p[i]] + glogtable[q[j]]];
+            }
+        }
+        return o;
+    }
+    if (L % 2) {
+        L += 1;
+        p.push_back(0);
+        q.push_back(0);
+    }
+    int halflen = L / 2;
+    vector<int> low1 = vector<int>(p.begin(), p.begin() + halflen);
+    vector<int> low2 = vector<int>(q.begin(), q.begin() + halflen);
+    vector<int> high1 = vector<int>(p.begin() + halflen, p.end());
+    vector<int> high2 = vector<int>(q.begin() + halflen, q.end());
+    vector<int> sum1(halflen);
+    vector<int> sum2(halflen);
+    for (int i = 0; i < halflen; i++) {
+        sum1[i] = low1[i] ^ high1[i];
+        sum2[i] = low2[i] ^ high2[i];
+    }
+    vector<int> z0 = karatsuba_mul(low1, low2);
+    vector<int> z2 = karatsuba_mul(high1, high2);
+    vector<int> m = karatsuba_mul(sum1, sum2);
+    vector<int> o(L * 2);
+    for (int i = 0; i < L; i++) {
+        o[i] ^= z0[i];
+        o[i + halflen] ^= (m[i] ^ z0[i] ^ z2[i]);
+        o[i + L] ^= z2[i];
+    }
+    return o;
+}
+
+vector<int> mk_root(vector<int> xs) {
+    int L = xs.size();
+    if (L >= ROOT_CUTOFF) {
+        int halflen = L / 2;
+        vector<int> left = vector<int>(xs.begin(), xs.begin() + halflen);
+        vector<int> right = vector<int>(xs.begin() + halflen, xs.end());
+        vector<int> o = karatsuba_mul(mk_root(left), mk_root(right));
+        o.resize(L + 1);
+        return o;
+    }
+    vector<int> root(L + 1);
+    root[L] = 1;
+    for (int i = 0; i < L; i++) {
         int logx = glogtable[xs[i]];
-        int offset = xs_size - i - 1;
+        int offset = L - i - 1;
         root[offset] = 0;
         for (int j = 0; j < i + 1; j++) {
             if (root[j + 1 + offset] and xs[i])
                 root[j + offset] ^= gexptable[glogtable[root[j+1 + offset]] + logx];
         }
     }
-    vector<int> b(xs_size);
-    vector<int> output(root.size() - 1);
-    for (int i = 0; i < xs_size; i++) {
-        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];
-        }
-        int denom = eval_poly_at(output, xs[i]);
-        int log_yslice = glogtable[ys[i]] - glogtable[denom] + 65535;
-        for (int j = 0; j < xs_size; j++) {
-            if(output[j] and ys[i])
-                b[j] ^= gexptable[glogtable[output[j]] + log_yslice];
-        }
-    }
-    return b;
+    return root;
 }
 
+vector<int> subroot_linear_combination(vector<int> xs, vector<int> factors) {
+    int L = xs.size();
+    /*if (L <= ROOT_CUTOFF) {
+        vector<int> out(L + 1);
+        vector<int> root = mk_root(xs);
+        for (int i = 0; i < L; i++) {
+            vector<int> output(L + 1);
+            output[root.size() - 2] = 1;
+            int logx = glogtable[xs[i]];
+            if (factors[i]) {
+                int log_fac = glogtable[factors[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];
+                    out[j] ^= gexptable[glogtable[output[j]] + log_fac];
+                }
+                out[0] ^= gexptable[glogtable[output[0]] + log_fac];
+            }
+        }
+        return out;
+    }*/
+    if (L == 1) {
+        vector<int> o(2);
+        o[0] = factors[0];
+        return o;
+    }
+    int halflen = L / 2;
+    vector<int> xs_left = vector<int>(xs.begin(), xs.begin() + halflen);
+    vector<int> xs_right = vector<int>(xs.begin() + halflen, xs.end());
+    vector<int> factors_left = vector<int>(factors.begin(), factors.begin() + halflen);
+    vector<int> factors_right = vector<int>(factors.begin() + halflen, factors.end());
+    vector<int> R1 = mk_root(xs_left);
+    vector<int> R2 = mk_root(xs_right);
+    vector<int> o1 = karatsuba_mul(R1, subroot_linear_combination(xs_right, factors_right));
+    vector<int> o2 = karatsuba_mul(R2, subroot_linear_combination(xs_left, factors_left));
+    vector<int> o(L + 1);
+    for (int i = 0; i < L; i++) {
+        o[i] = o1[i] ^ o2[i];
+    }
+    return o;
+}
+
+
+vector<int> derivative(vector<int> p) {
+    vector<int> o(p.size() - 1);
+    for (int i = 0; i < o.size(); i+= 2) {
+        o[i] = p[i + 1];
+    }
+    return o;
+}
+
+vector<int> poly_to_logs(vector<int> p) {
+    vector<int> o(p.size());
+    for (int i = 0; i < p.size(); i++) {
+        if (p[i])
+            o[i] = glogtable[p[i]];
+        else
+            o[i] = 65537;
+    }
+    return o;
+}
+
+vector<int> lagrange_interp(vector<int> ys, vector<int> xs) {
+    int xs_size = xs.size();
+    vector<int> root = mk_root(xs);
+    vector<int> log_rootprime = poly_to_logs(derivative(root));
+    vector<int> factors(xs_size);
+    for (int i = 0; i < xs_size; i++) {
+        int denom = eval_log_poly_at(log_rootprime, xs[i]);
+        if (ys[i])
+            factors[i] = gexptable[glogtable[ys[i]] + 65535 - glogtable[denom]];
+    }
+    return subroot_linear_combination(xs, factors);
+}
+
+const int SIZE = 1024;
+
 
 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> xs(SIZE);
+    vector<int> ys(SIZE);
+    for (int v = 0; v < SIZE; v++) {
+        ys[v] = v * 3;
+        xs[v] = 1000 + v * 7;
     }
+    //vector<int> d = derivative(mk_root(xs));
+    //for (int i = 0; i < d.size(); i++) cout << d[i] << " ";
+    //cout << "\n";
+    /*vector<int> prod = mk_root(xs);
+    vector<int> prod = karatsuba_mul(xs, ys);
+    for (int i = 0; i < SIZE + 1; i++)
+        cout << prod[i] << " ";
+    cout << "\n";
+    cout << eval_poly_at(prod, 189) << " " << gexptable[glogtable[eval_poly_at(xs, 189)] + glogtable[eval_poly_at(ys, 189)]] << "\n";*/
     vector<int> poly = lagrange_interp(ys, xs);
+    cout << eval_poly_at(poly, 1700) << "\n";
     unsigned int o = 0;
-    for (int i = 4096; i < 8192; i++) {
+    for (int i = SIZE; i < SIZE * 2; i++) {
         o += eval_poly_at(poly, i);
     }
-    cout << o;
+    cout << o << "\n";
     //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";
diff --git a/erasure_code/ec65536/subquadratic_poly_utils.py b/erasure_code/ec65536/subquadratic_poly_utils.py
index f9e35ae..1aa0838 100644
--- a/erasure_code/ec65536/subquadratic_poly_utils.py
+++ b/erasure_code/ec65536/subquadratic_poly_utils.py
@@ -74,7 +74,10 @@ def lagrange_interp(pieces, xs):
     # denoms = [eval_poly_at(d, xs[i]) for i in range(len(xs))]
     # Generate output polynomial, which is the sum of the per-value numerator
     # polynomials rescaled to have the right y values
-    return multi_root_derive(xs, [galois_div(p, d) for p, d in zip(pieces, denoms)])
+    factors = [galois_div(p, d) for p, d in zip(pieces, denoms)]
+    o = multi_root_derive(xs, factors)
+    print(o)
+    return o
 
 def multi_root_derive(xs, muls):
     if len(xs) == 1: