from random import randint
def bitcoin_hash(block_header,nonce):
    import hashlib
    return int('0x' + hashlib.sha256(hashlib.sha256(str(block_header + nonce)).digest()).hexdigest(), 16)

max_val = 2**256

class SeedObj():
    "A class for managing a random number state, using a linear congruential generator"
    def __init__(self, seed):
        self.seed = seed
        self.a = 3**160
        self.c = 7**80
        self.n = 2**256 - 4294968273 # secp256k1n, why not        

    def rand(self, r):
        self.seed = (self.seed * self.a + self.c) % self.n
        return self.seed % r

def encode_int(x):
    "A helper function to provide a standard 64-bit big-endian encoding for numbers"
    o = ''
    for _ in range(8):
        o = chr(x % 256) + o
        x //= 256 
    return o
    
ops = {
    "plus": lambda x,y: (x + y) % 2**64,
    "times": lambda x,y: (x * y) % 2**64,
    "xor": lambda x,y: x^y,
    "and": lambda x,y: x&y,
    "or": lambda x,y: x|y,
    "not": lambda x,y: 2**64-1-x,
    "nxor": lambda x,y: (2**64-1-x) ^ y,
    "rshift": lambda x,y: x >> (y % 64)
}

def gentape(W, H, SEED):
    "Generate a random tape of binary instructions"
    s = SeedObj(SEED)
    tape = []
    for i in range(H):
        op = ops.keys()[s.rand(len(ops))]
        r = s.rand(100)
        if r < 20 and i > 20:
            x1 = tape[-r]["x1"]
        else:
            x1 = s.rand(W)
        x2 = s.rand(W)
        tape.append({"op": op, "x1": x1, "x2": x2})
    return tape

def runtape(TAPE, SEED, S):
    import hashlib
    s = SeedObj(SEED)
    mem = [s.rand(2**64) for _ in range(S)]
    dir = 1
    for i in range(S // 100):
        for j in (range(100) if dir == 1 else range(99, -1, -1)):
            t = TAPE[i * 100 + j]
            mem[t["x1"]] = ops[t["op"]](mem[t["x1"]], mem[t["x2"]])
        if 2 < mem[t["x1"]] % 37 < 9:
            dir *= -1
    return int('0x' + hashlib.sha256(''.join(encode_int(x) for x in mem)).hexdigest(), 16)

def PoWVerify(block_header, nonce, S, D):
    tape = gentape(S, S, bitcoin_hash(block_header, nonce))
    h = runtape(tape, bitcoin_hash(block_header, nonce), S)
    return h < max_val / D

def rc_mine(block_header, S, D):
    nonce = 1
    while not PoWVerify(block_header, nonce, S, D):
        nonce += 1
    return nonce

# Example Run
rc_mine(randint(0,max_val), 10, 10)

from timeit import timeit

def time_rc_mine(S,D):
    return timeit('rc_mine(randint(0,max_val),  %d, %d)' % (S,D), 
                  'from random import randint ; from __main__ import rc_mine, max_val', 
                  number=1)

rc_S_difficulties = range(10,370,50)
rc_S_samples = [[time_rc_mine(S,2) for _ in range(1000)] for S in rc_S_difficulties]

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import numpy as np
from pylab import polyfit
import matplotlib.pyplot as plt
rc_S_medians = map(np.median, rc_S_samples)
rc_S_medians_A, rc_S_medians_k = \
    polyfit(rc_S_difficulties, rc_S_medians, 1)
plt.plot(rc_S_difficulties, rc_S_medians)
plt.plot(rc_S_difficulties, 
         [rc_S_medians_A*S + rc_S_medians_k 
          for S in rc_S_difficulties], '-r')
plt.title('$S$ vs. Median Mining Time (Linear)')
plt.ylabel('Seconds')
plt.xlabel('$S$')
plt.show()
rc_S_medians_C = - rc_S_medians_A * np.log2(1-1./2)
print "Estimated Ĉ:", rc_S_medians_C

rc_D_difficulties = range(1,502,50)
rc_D_samples = [[time_rc_mine(1,D) for _ in range(1000)] for D in rc_D_difficulties]

rc_D_medians = map(np.median, rc_D_samples)
plt.plot(rc_D_difficulties, rc_D_medians)
rc_D_medians_A, rc_D_medians_k = \
    polyfit(rc_D_difficulties, rc_D_medians, 1)
plt.plot(rc_D_difficulties, 
         [rc_D_medians_A*S + rc_D_medians_k 
          for S in rc_D_difficulties], '-r')
plt.title('$D$ vs. Median Mining Time (Linear)')
plt.ylabel('Seconds')
plt.xlabel('$D$')
plt.show()
rc_D_medians_C = rc_D_medians_A / np.log(2)
print "Estimated Ĉ:"
print "   using S medians:", rc_S_medians_C
print "   using D medians:", rc_D_medians_C

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10,5))
fig.subplots_adjust(wspace=0.3, hspace=0.4)

rc_S_means = map(np.mean, rc_S_samples)
rc_D_means = map(np.mean, rc_D_samples)
ax1.plot(rc_S_difficulties, rc_S_means)
ax2.plot(rc_D_difficulties, rc_D_means)

rc_S_vars = map(np.var, rc_S_samples)
rc_D_vars = map(np.var, rc_D_samples)
ax3.plot(rc_S_difficulties, rc_S_vars)
ax4.plot(rc_D_difficulties, rc_D_vars)

rc_S_means_A, rc_S_means_k = \
    polyfit(rc_S_difficulties, rc_S_means, 1)
ax1.plot(rc_S_difficulties, 
         [rc_S_means_A*S + rc_S_means_k 
          for S in rc_S_difficulties], '-r')
rc_S_means_C = rc_S_means_A / 2  # Recall that we set D=2 for these samples

rc_D_means_C, rc_D_means_k = \
    polyfit(rc_D_difficulties, rc_D_means, 1)
ax2.plot(rc_D_difficulties, 
         [rc_D_means_C*D + rc_D_means_k 
          for D in rc_D_difficulties], '-r')

rc_S_vars_A, rc_S_vars_B, rc_S_vars_k = \
    polyfit(rc_S_difficulties, rc_S_vars, 2)
ax3.plot(rc_S_difficulties, 
         [rc_S_vars_A*S*S + rc_S_vars_B*S + rc_S_vars_k 
          for S in rc_S_difficulties], '-r')
rc_S_vars_C = np.sqrt(rc_S_vars_A / 2)  # Recall that we set D=2 for these samples

rc_D_vars_A, rc_D_vars_B, rc_D_vars_k = \
    polyfit(rc_D_difficulties, rc_D_vars, 2)
ax4.plot(rc_D_difficulties, 
         [rc_D_vars_A*D*D + rc_D_vars_B*D + rc_D_vars_k 
          for D in rc_D_difficulties], '-r')
rc_D_vars_C = np.sqrt(rc_D_vars_A)

ax1.set_title('$S$ vs. Mean Mining Time (Linear)')
ax1.set_ylabel('Seconds')
ax1.set_xlabel('$S$')
ax2.set_title('$D$ vs. Mean Mining Time (Linear)')
ax2.set_ylabel('Seconds')
ax2.set_xlabel('$D$')
ax3.set_title('$S$ vs. Mining Time Variance (Quadratic)')
ax3.set_ylabel('Seconds')
ax3.set_xlabel('$S$')
ax4.set_title('$D$ vs. Mining Time Variance (Quadratic)')
ax4.set_ylabel('Seconds')
ax4.set_xlabel('$D$')
plt.show()

print "Estimated Ĉ:"
print "   using S medians:", rc_S_medians_C
print "   using S means:", rc_S_means_C
print "   using S variance:", rc_S_vars_C
print "   using D medians:", rc_D_medians_C
print "   using D means:", rc_D_means_C
print "   using D variance:", rc_D_vars_C

import scipy.stats
import pandas

C = rc_S_medians_C
S_data = []
for S, samples in zip(rc_S_difficulties, rc_S_samples):
    cdf = lambda t: 1 - (1 - (1. / 2))**(t / (S * C))
    ks_statistic, p_value = scipy.stats.kstest(samples, cdf)
    S_data.append({"Difficulty D": 2,
                   "Difficulty S": S,
                   "KS Statistic": ks_statistic,
                   "P-Value": p_value})


pandas.DataFrame.from_records(S_data)

D_data = []
for D, samples in zip(rc_D_difficulties, rc_D_samples):
    cdf = lambda t: 1 - (1 - (1. / D))**(t / C)
    ks_statistic, p_value = scipy.stats.kstest(samples, cdf)
    D_data.append({"Difficulty D": D,
                   "Difficulty S": 1,
                   "KS Statistic": ks_statistic,
                   "P-Value": p_value})
pandas.DataFrame.from_records(D_data)