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()
71.2% of values within 1.0 ULP
# 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)
2.2204460492503131
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
array([ 50, 59, 255, 132, 31, 179, 243, 63], dtype=uint8)
y0 = a.view(np.float64)[0]
a[0] += 1
y1 = a.view(np.float64)[0]
np.abs(y1 - y0)
2.2204460492503131e-16
2**-52
2.220446049250313e-16