import pandas as pd
from cStringIO import StringIO
import random
from decimal import Decimal
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def xstrtod(s):
    return pd.read_csv(StringIO(s), header=None)[0][0]

def randnum():
    s = '1.'  # Start with 1, otherwise you can randomly get very small numbers
    for i in range(20):
        s += random.choice('0123456789')
    return s

def ulp(s):
    num = Decimal(s)
    f = np.float64(s)
    a = f.view((np.uint8, 8))
    # Since this is uint8 make sure the result doesn't accidentally wrap
    if a[0] == 0:
        a[0] = 1
    elif a[0] == 255:
        a[0] = 254
    elif Decimal(f) < num:
        a[0] += 1
    elif Decimal(f) > num:
        a[0] -= 1
    f2 = a.view(np.float64)[0]
    return abs(f2 - f)

n_vals = 10000
y = np.empty(n_vals)
good_vals = 0
ulp_vals = np.empty(n_vals)
diffs = np.empty(n_vals)

for i in range(n_vals):
    val = randnum()
    guess = xstrtod(val)
    decimal_ulp_val = Decimal(ulp(val))
    ulp_vals[i] = float(decimal_ulp_val / Decimal(1e-16))
    decimal_diff = abs(Decimal(val) - Decimal(guess))
    ulp_diff = decimal_diff / decimal_ulp_val
    y[i] = float(ulp_diff)
    diffs[i] = float(decimal_diff)
    if float(ulp_diff) < 1.0:
        good_vals += 1

print('{0}% of values within 1.0 ULP'.format(good_vals * 100.0 / 10000))
plt.hist(y, bins=30, log=True)
plt.show()

# OK, the ULP is always very near 2**-52 which is one bit in the
# 52-bit mantissa (http://kipirvine.com/asm/workbook/floating_tut.htm)
plt.hist(ulp_vals, bins=50)
plt.xlabel('ULP / 1e-16')
plt.grid()
plt.show()
np.min(ulp_vals)

plt.hist(diffs, bins=50, log=True)
plt.xlabel('Diff')
plt.grid()
plt.show()

x = np.float64(1.23123123123123123123)

a = x.view((np.uint8, 8))
a

y0 = a.view(np.float64)[0]
a[0] += 1
y1 = a.view(np.float64)[0]
np.abs(y1 - y0)

2**-52