From 3397bd982e0690c1d06d43aa54ad86944d693730 Mon Sep 17 00:00:00 2001 From: vub Date: Sun, 23 Apr 2017 08:10:21 -0400 Subject: [PATCH] Added some optimizations --- erasure_code/ec65536/poly_utils.cpp | 39 ++-- erasure_code/ec65536/poly_utils.py | 5 +- .../ec65536/subquadratic_poly_utils.py | 193 ++++++++++++++++++ iceage.py | 8 +- 4 files changed, 217 insertions(+), 28 deletions(-) create mode 100644 erasure_code/ec65536/subquadratic_poly_utils.py diff --git a/erasure_code/ec65536/poly_utils.cpp b/erasure_code/ec65536/poly_utils.cpp index 35b39d5..d0b38b8 100644 --- a/erasure_code/ec65536/poly_utils.cpp +++ b/erasure_code/ec65536/poly_utils.cpp @@ -36,19 +36,21 @@ int eval_poly_at(vector poly, int x) { } vector lagrange_interp(vector ys, vector xs) { - deque root; - root.push_back(1); - for (int i = 0; i < xs.size(); i++) { + int xs_size = xs.size(); + vector root(xs_size + 1); + root[xs_size] = 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]; + int offset = xs_size - 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 > nums; - for (int i = 0; i < xs.size(); i++) { - vector output(root.size() - 1); + vector b(xs_size); + vector 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--) { @@ -57,18 +59,11 @@ vector lagrange_interp(vector ys, vector xs) { else output[j - 1] = root[j]; } - nums.push_back(output); - } - vector denoms; - for (int i = 0; i < xs.size(); i++) { - denoms.push_back(eval_poly_at(nums[i], xs[i])); - } - vector 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]; + 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; diff --git a/erasure_code/ec65536/poly_utils.py b/erasure_code/ec65536/poly_utils.py index ecc7d9d..a02298b 100644 --- a/erasure_code/ec65536/poly_utils.py +++ b/erasure_code/ec65536/poly_utils.py @@ -70,8 +70,9 @@ def lagrange_interp(pieces, xs): for j in range(len(root)-1): if root[j+1] and x: root[j] ^= gexptable[glogtable[root[j+1]] + logx] + #print(root) assert len(root) == len(pieces) + 1 - print(root) + # 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 @@ -86,7 +87,7 @@ def lagrange_interp(pieces, xs): output[j-1] = root[j] assert len(output) == len(pieces) nums.append(output) - print(nums) + #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/subquadratic_poly_utils.py b/erasure_code/ec65536/subquadratic_poly_utils.py new file mode 100644 index 0000000..f9e35ae --- /dev/null +++ b/erasure_code/ec65536/subquadratic_poly_utils.py @@ -0,0 +1,193 @@ +modulus_poly = [1, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 1, 0, 1, 0, 0, 1, + 1] +modulus_poly_as_int = sum([(v << i) for i, v in enumerate(modulus_poly)]) +degree = len(modulus_poly) - 1 + +two_to_the_degree = 2**degree +two_to_the_degree_m1 = 2**degree - 1 + +def galoistpl(a): + # 2 is not a primitive root, so we have to use 3 as our logarithm base + if a * 2 < two_to_the_degree: + return (a * 2) ^ a + else: + return (a * 2) ^ a ^ modulus_poly_as_int + +# Precomputing a log table for increased speed of addition and multiplication +glogtable = [0] * (two_to_the_degree) +gexptable = [] +v = 1 +for i in range(two_to_the_degree_m1): + glogtable[v] = i + gexptable.append(v) + v = galoistpl(v) + +gexptable += gexptable + gexptable + +# Add two values in the Galois field +def galois_add(x, y): + return x ^ y + +# In binary fields, addition and subtraction are the same thing +galois_sub = galois_add + +# Multiply two values in the Galois field +def galois_mul(x, y): + return 0 if x*y == 0 else gexptable[glogtable[x] + glogtable[y]] + +# Divide two values in the Galois field +def galois_div(x, y): + return 0 if x == 0 else gexptable[(glogtable[x] - glogtable[y]) % two_to_the_degree_m1] + +# Evaluate a polynomial at a point +def eval_poly_at(p, x): + if x == 0: + return p[0] + y = 0 + logx = glogtable[x] + for i, p_coeff in enumerate(p): + if p_coeff: + # Add x**i * coeff + y ^= gexptable[(logx * i + glogtable[p_coeff]) % two_to_the_degree_m1] + return y + + +# Given p+1 y values and x values with no errors, recovers the original +# p+1 degree polynomial. +# Lagrange interpolation works roughly in the following way. +# 1. Suppose you have a set of points, eg. x = [1, 2, 3], y = [2, 5, 10] +# 2. For each x, generate a polynomial which equals its corresponding +# y coordinate at that point and 0 at all other points provided. +# 3. Add these polynomials together. + +def lagrange_interp(pieces, xs): + # Generate master numerator polynomial, eg. (x - x1) * (x - x2) * ... * (x - xn) + root = mk_root_2(xs) + #print(root) + assert len(root) == len(pieces) + 1 + # print(root) + # Generate the derivative + d = derivative(root) + # Generate denominators by evaluating numerator polys at each x + denoms = multi_eval_2(d, 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)]) + +def multi_root_derive(xs, muls): + if len(xs) == 1: + return [muls[0]] + R1 = mk_root_2(xs[:len(xs) // 2]) + R2 = mk_root_2(xs[len(xs) // 2:]) + x1 = karatsuba_mul(R1, multi_root_derive(xs[len(xs) // 2:], muls[len(muls) // 2:]) + [0]) + x2 = karatsuba_mul(R2, multi_root_derive(xs[:len(xs) // 2], muls[:len(muls) // 2]) + [0]) + o = [v1 ^ v2 for v1, v2 in zip(x1, x2)][:len(xs)] + # print(len(R1), len(x1), len(xs), len(o)) + return o + +def multi_root_derive_1(xs, muls): + o = [0] * len(xs) + for i in range(len(xs)): + _xs = xs[:i] + xs[(i+1):] + root = mk_root_2(_xs) + for j in range(len(root)): + o[j] ^= galois_mul(root[j], muls[i]) + return o + +a = 124 +b = 8932 +c = 12415 + +assert galois_mul(galois_add(a, b), c) == galois_add(galois_mul(a, c), galois_mul(b, c)) + +def karatsuba_mul(p1, p2): + L = len(p1) + # assert L == len(p2) + if L <= 16: + o = [0] * (L * 2) + for i, v1 in enumerate(p1): + for j, v2 in enumerate(p2): + if v1 and v2: + o[i + j] ^= gexptable[glogtable[v1] + glogtable[v2]] + return o + if L % 2: + p1 = p1 + [0] + p2 = p2 + [0] + L += 1 + halflen = L // 2 + low1 = p1[:halflen] + high1 = p1[halflen:] + sum1 = [l ^ h for l, h in zip(low1, high1)] + low2 = p2[:halflen] + high2 = p2[halflen:] + sum2 = [l ^ h for l, h in zip(low2, high2)] + z2 = karatsuba_mul(high1, high2) + z0 = karatsuba_mul(low1, low2) + z1 = [m ^ _z0 ^ _z2 for m, _z0, _z2 in zip(karatsuba_mul(sum1, sum2), z0, z2)] + o = z0[:halflen] + \ + [a ^ b for a, b in zip(z0[halflen:], z1[:halflen])] + \ + [a ^ b for a, b in zip(z2[:halflen], z1[halflen:])] + \ + z2[halflen:] + return o + +def mk_root_1(xs): + root = [1] + for x in xs: + logx = glogtable[x] + root.insert(0, 0) + for j in range(len(root)-1): + if root[j+1] and x: + root[j] ^= gexptable[glogtable[root[j+1]] + logx] + return root + +def mk_root_2(xs): + if len(xs) >= 128: + return karatsuba_mul(mk_root_2(xs[:len(xs) // 2]), mk_root_2(xs[len(xs) // 2:]))[:len(xs) + 1] + root = [1] + for x in xs: + logx = glogtable[x] + root.insert(0, 0) + for j in range(len(root)-1): + if root[j+1] and x: + root[j] ^= gexptable[glogtable[root[j+1]] + logx] + return root + +def derivative(root): + return [0 if i % 2 else r for i, r in enumerate(root[1:])] + +# Credit to http://people.csail.mit.edu/madhu/ST12/scribe/lect06.pdf for the algorithm +def xn_mod_poly(p): + if len(p) == 1: + return [galois_div(1, p[0])] + halflen = len(p) // 2 + lowinv = xn_mod_poly(p[:halflen]) + submod_high = karatsuba_mul(lowinv, p[:halflen])[halflen:] + med = karatsuba_mul(p[halflen:], lowinv)[:halflen] + med_plus_high = [x ^ y for x, y in zip(med, submod_high)] + highinv = karatsuba_mul(med_plus_high, lowinv) + o = (lowinv + highinv)[:len(p)] + # assert karatsuba_mul(o, p)[:len(p)] == [1] + [0] * (len(p) - 1) + return o + +def mod(a, b): + assert len(a) == 2 * (len(b) - 1) + L = len(b) + inv_rev_b = xn_mod_poly(b[::-1] + [0] * (len(a) - L))[:L] + quot = karatsuba_mul(inv_rev_b, a[::-1][:L])[:L-1][::-1] + subt = karatsuba_mul(b, quot + [0])[:-1] + o = [x ^ y for x, y in zip(a[:L-1], subt[:L-1])] + # assert [x^y for x, y in zip(karatsuba_mul(quot + [0], b), o)] == a + return o + +def multi_eval_1(poly, xs): + return [eval_poly_at(poly, x) for x in xs] + +def multi_eval_2(poly, xs): + if len(xs) <= 1024: + return [eval_poly_at(poly, x) for x in xs] + halflen = len(xs) // 2 + return multi_eval_2(mod(poly, mk_root_2(xs[:halflen])), xs[:halflen]) + \ + multi_eval_2(mod(poly, mk_root_2(xs[halflen:])), xs[halflen:]) + # [eval_poly_at(poly, xs[-2]), eval_poly_at(poly, xs[-1])] diff --git a/iceage.py b/iceage.py index 6c3b5bf..1aed867 100644 --- a/iceage.py +++ b/iceage.py @@ -1,12 +1,12 @@ import random import datetime -diffs = [257.74 * 10**12] -hashpower = diffs[0] / 14.7 -times = [1491920186] +diffs = [286.86 * 10**12] +hashpower = diffs[0] / 15.0 +times = [1492650109] -for i in range(3517029, 6010000): +for i in range(3566076, 6010000): blocktime = random.expovariate(hashpower / diffs[-1]) adjfac = max(1 - int(blocktime / 10), -99) / 2048. newdiff = diffs[-1] * (1 + adjfac)