Fast vectorized arithmetic with ~32 significant digits under Numpy

dd.py

# 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)