Skip to content

Instantly share code, notes, and snippets.

@baruchel
Last active Dec 11, 2015
Embed
What would you like to do?
Fast vectorized arithmetic with ~32 significant digits under Numpy
# Adapted from:
# "A Floating-Point Technique for Extending the vailable Precision"
# by T.J. Dekker in Numer. Math. 18, 224-242, 1971.
# (adapted to Numpy from the Algol 60 code in the paper)
import numpy as np
import decimal
_npDecimal = np.vectorize(decimal.Decimal, otypes=[object])
try:
import mpmath
_npMp = np.vectorize(mpmath.mpf, otypes=[object])
except:
pass
try:
import gmpy
_npGmpy = np.frompyfunc(gmpy.mpf, 1, 1)
except:
pass
constant = 2**(53-53/float(2)) + 1
class DD:
def __init__(self, d, copy=True):
if d==None:
return
if isinstance(d, np.ndarray):
if copy:
self.major = np.array(d)
else:
self.major = d
self.minor = np.zeros(d.shape)
return
if isinstance(d, tuple) or isinstance(d, int):
self.major = np.zeros(d)
self.minor = np.zeros(d)
return
def as_floats(self):
return self.major
def as_decimals(self):
return _npDecimal(self.major) + _npDecimal(self.minor)
def as_mpmath(self):
return _npMp(self.major) + _npMp(self.minor)
def as_gmpy(self):
return _npGmpy(self.major) + _npGmpy(self.minor)
@staticmethod
def from_decimals(a):
z = DD(None)
z.major = np.array(a, dtype=np.float64)
z.minor = np.array(a - _npDecimal(z.major), dtype=np.float64)
return z
@staticmethod
def from_mpf(a):
z = DD(None)
z.major = np.array(a, dtype=np.float64)
z.minor = np.array(a - z.major, dtype=np.float64)
return z
def add(self, b):
# solution 1
r = self.major + b.major
test = abs(self.major) > abs(b.major)
ntest = np.logical_not(test)
s = (
test*(self.major - r + b.major + b.minor + self.minor)
+ ntest*(b.major - r + self.major + self.minor + b.minor)
)
z = DD(None)
z.major = r + s
z.minor = r - z.major + s
return z
# # solution 2
# r = self.major + b.major
# s = (abs(self.major) > abs(b.major)).choose(
# np.array([
# b.major - r + self.major + self.minor + b.minor,
# self.major - r + b.major + b.minor + self.minor ]) )
# z = DD(None)
# z.major = r + s
# z.minor = r - z.major + s
# return z
# # solution 3
# r = self.major + b.major
# test = abs(self.major) > abs(b.major)
# s = b.major - r + self.major + self.minor + b.minor
# s[test] = (self.major - r + b.major + b.minor + self.minor)[test]
# z = DD(None)
# z.major = r + s
# z.minor = r - z.major + s
# return z
def addf(self, b):
r = self.major + b
test = abs(self.major) > abs(b)
ntest = np.logical_not(test)
s = (
test*(self.major - r + b + self.minor)
+ ntest*(b - r + self.major + self.minor)
)
z = DD(None)
z.major = r + s
z.minor = r - z.major + s
return z
@staticmethod
def from_sum_floats(a, b):
r = a + b
test = abs(a) > abs(b)
ntest = np.logical_not(test)
s = (
test*(a - r + b)
+ ntest*(b - r + a)
)
z = DD(None)
z.major = r + s
z.minor = r - z.major + s
return z
def sub(self, b):
r = self.major - b.major
test = abs(self.major) > abs(b.major)
ntest = np.logical_not(test)
s = (
test*(self.major - r - b.major - b.minor + self.minor)
+ ntest*(-b.major - r + self.major + self.minor - b.minor)
)
z = DD(None)
z.major = r + s
z.minor = r - z.major + s
return z
def subf(self, b):
r = self.major - b
test = abs(self.major) > abs(b)
ntest = np.logical_not(test)
s = (
test*(self.major - r - b + self.minor)
+ ntest*(-b - r + self.major + self.minor)
)
z = DD(None)
z.major = r + s
z.minor = r - z.major + s
return z
@staticmethod
def from_product_floats(x,y):
global constant
p = x * constant
hx = x - p + p
tx = x - hx
p = y * constant
hy = y - p + p
ty = y - hy
p = hx * hy
q = hx * ty + tx * hy
z = DD(None)
z.major = p + q
z.minor = p - z.major + q + tx * ty
return z
def mul(self,b):
global constant
c = DD.from_product_floats(self.major, b.major)
c.minor = self.major * b.minor + self.minor * b.major + c.minor
z = DD(None)
z.major = c.major + c.minor
z.minor = c.major - z.major + c.minor
return z
def mulf(self,b):
global constant
c = DD.from_product_floats(self.major, b)
c.minor = self.minor * b + c.minor
z = DD(None)
z.major = c.major + c.minor
z.minor = c.major - z.major + c.minor
return z
def div(self,b):
c = self.major / b.major
u = DD.from_product_floats(c, b.major)
cc = (self.major - u.major - u.minor + self.minor
- c * b.minor) / b.major
z = DD(None)
z.major = c + cc
z.minor = c - z.major + cc
return z
def divf(self,b):
c = self.major / b
u = DD.from_product_floats(c, b)
cc = (self.major - u.major - u.minor + self.minor) / b
z = DD(None)
z.major = c + cc
z.minor = c - z.major + cc
return z
def sqrt(self):
c = np.sqrt(self.major)
u = DD.from_product_floats(c, c)
cc = (self.major - u.major - u.minor + self.minor) * 0.5 / c
z = DD(None)
z.major = c + cc
z.minor = c - z.major + cc
return z
# basic test
import gmpy
import mpmath
from time import time
a = np.array([ gmpy.mpf(k) for k in range(1,100000) ], dtype=object)
b = np.array([ mpmath.mpf(k) for k in range(1,100000) ], dtype=object)
c = DD(np.array([ k for k in range(1,100000) ]))
print("\nWith GMPY (dtype=object) ...")
t = time(); d = 1/a; ta = time() - t
print(d); print(" --> time =", ta)
print("\nWith Mpmath (dtype=object) ...")
t = time(); e = 1/b; tb = time() - t
print(e); print(" --> time =", tb)
print("\nWith extended-double ...")
t = time(); f = DD(np.ones(99999)).div(c); tc = time() - t
print(f.as_gmpy()); print(" --> time =", tc)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment