Skip to content

Instantly share code, notes, and snippets.

@al42and
Created January 8, 2014 18:18
Show Gist options
  • Save al42and/8321600 to your computer and use it in GitHub Desktop.
Save al42and/8321600 to your computer and use it in GitHub Desktop.
This code snippet is a part of stackexchange answer http://bitcoin.stackexchange.com/a/20060/2343
#!/usr/bin/env python
#-*- encoding: utf-8 -*-
# This code is a part of stackexchange answer
# http://bitcoin.stackexchange.com/a/20060/2343
# License: cc-wiki with attribution required
import numpy as np
HR = 550 # GHash/s, own hashrate
THR0 = 1.2e7 # GHash/s, network hashrate
DA = 1.2 # Expected difficulty adjustment
N = 26*2016 # For how much block make prediction? ~ 1 year
NB = 2016 # Number of block between difficulty adjustments
T0 = 10 # minutes, desired block time
def THR(n):
return THR0*np.exp((n/NB)*np.log(DA))
p = HR/THR(np.arange(N))
# To caclulate distribution we use algorithm from
# Yili Hong "On Computing the Distribution Function for the Sum of
# Independent and Non-identical Random Indicators", 2011
# http://www.stat.vt.edu/research/Technical_Reports/TechReport11-2.pdf
x = np.zeros(N+1, dtype=np.complex128)
w = 2*np.pi / (N+1)
x[0] = 1.0 + 0.j
for l in xrange(1, int(np.ceil(N/2))+1):
z = 1.0 - p[:] + p[:]*np.cos(w*l) + 1.j * p[:]*np.sin(w*l)
dl = np.exp(np.sum(np.log(np.absolute(z[:]))))
ang = np.sum(np.angle(z[:]))
al = dl*np.cos(ang)
bl = dl*np.sin(ang)
x[l] = al + 1.j*bl
if l%1000 == 0: print l, '/', N/2
for l in xrange(int(np.ceil(N/2))+1, N+1):
x[l] = np.conj(x[N+1-l])
x = x[:] / (N+1)
y = np.fft.fft(x)
print DA
for i in enumerate(np.real(y[:10])):
print i[0], i[1]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment