diff --git a/segpy/caching.py b/segpy/caching.py
new file mode 100644
index 0000000..ae95d6f
--- /dev/null
+++ b/segpy/caching.py
@@ -0,0 +1,48 @@
+def lfu_cache(maxsize=100):
+    '''Least-frequently-used cache decorator.
+
+    Arguments to the cached function must be hashable.
+    Cache performance statistics stored in f.hits and f.misses.
+    Clear the cache with f.clear().
+    http://en.wikipedia.org/wiki/Least_Frequently_Used
+
+    '''
+    def decorating_function(user_function):
+        cache = {}                      # mapping of args to results
+        use_count = Counter()           # times each key has been accessed
+        kwd_mark = object()             # separate positional and keyword args
+
+        @functools.wraps(user_function)
+        def wrapper(*args, **kwds):
+            key = args
+            if kwds:
+                key += (kwd_mark,) + tuple(sorted(kwds.items()))
+            use_count[key] += 1
+
+            # get cache entry or compute if not found
+            try:
+                result = cache[key]
+                wrapper.hits += 1
+            except KeyError:
+                result = user_function(*args, **kwds)
+                cache[key] = result
+                wrapper.misses += 1
+
+                # purge least frequently used cache entry
+                if len(cache) > maxsize:
+                    for key, _ in nsmallest(maxsize // 10,
+                                            use_count.iteritems(),
+                                            key=itemgetter(1)):
+                        del cache[key], use_count[key]
+
+            return result
+
+        def clear():
+            cache.clear()
+            use_count.clear()
+            wrapper.hits = wrapper.misses = 0
+
+        wrapper.hits = wrapper.misses = 0
+        wrapper.clear = clear
+        return wrapper
+    return decorating_function
diff --git a/segpy/ibm_float.py b/segpy/ibm_float.py
index 1a53701..076867e 100644
--- a/segpy/ibm_float.py
+++ b/segpy/ibm_float.py
@@ -1,15 +1,31 @@
-from math import frexp, isnan, isinf
+from math import frexp, isnan, isinf, ceil, floor, trunc
+from numbers import Real
 
 from segpy.portability import long_int, byte_string, four_bytes
 
+
+
+IBM_ZERO_BYTES = b'\x00\x00\x00\x00'
+IBM_NEGATIVE_ONE_BYTES = b'\xc1\x10\x00\x00'
+IBM_POSITIVE_ONE_BYTES = b'A\x10\x00\x00'
+
 MIN_IBM_FLOAT = -7.2370051459731155e+75
 LARGEST_NEGATIVE_NORMAL_IBM_FLOAT = -5.397605346934028e-79
 SMALLEST_POSITIVE_NORMAL_IBM_FLOAT = 5.397605346934028e-79
 MAX_IBM_FLOAT = 7.2370051459731155e+75
 
-_IBM_FLOAT32_BITS_PRECISION = 24
-_L24 = long_int(2) ** _IBM_FLOAT32_BITS_PRECISION
-_F24 = float(pow(2, _IBM_FLOAT32_BITS_PRECISION))
+MAX_BITS_PRECISION_IBM_FLOAT = 24
+MIN_BITS_PRECISION_IBM_FLOAT = 21  # The first 3 bits of the mantissa may be zero
+EPSILON_IBM_FLOAT = pow(2.0, -(MIN_BITS_PRECISION_IBM_FLOAT - 1))
+_L24 = long_int(2) ** MAX_BITS_PRECISION_IBM_FLOAT
+_F24 = float(pow(2, MAX_BITS_PRECISION_IBM_FLOAT))
+
+_L21 = long_int(2) ** MIN_BITS_PRECISION_IBM_FLOAT
+
+EXPONENT_BIAS = 64
+
+MIN_EXACT_INTEGER_IBM_FLOAT = -2**MAX_BITS_PRECISION_IBM_FLOAT
+MAX_EXACT_INTEGER_IBM_FLOAT = 2**MIN_BITS_PRECISION_IBM_FLOAT
 
 
 def ibm2ieee(big_endian_bytes):
@@ -27,15 +43,55 @@ def ibm2ieee(big_endian_bytes):
         return 0.0
 
     sign = -1 if (a & 0x80) else 1
-    exponent = a & 0x7f
+    exponent_16_biased = a & 0x7f
     mantissa = ((b << 16) | (c << 8) | d) / _F24
 
-    value = sign * mantissa * pow(16, exponent - 64)
+    value = sign * mantissa * pow(16, exponent_16_biased - EXPONENT_BIAS)
     return value
 
+BITS_PER_NYBBLE = 4
+
+
+def truncate(big_endian_bytes):
+    a, b, c, d = four_bytes(big_endian_bytes)
+
+    sign = -1 if (a & 0x80) else 1
+    exponent_16_biased = a & 0x7f
+    exponent_16 = exponent_16_biased - EXPONENT_BIAS
+    mantissa = ((b << 16) | (c << 8) | d)
+
+    print("sign =", sign)
+    print("exponent_16", exponent_16, hex(exponent_16), bin(exponent_16))
+    print("mantissa", mantissa, hex(mantissa), bin(mantissa))
+
+    num_nybbles_to_preserve = min(exponent_16, MAX_BITS_PRECISION_IBM_FLOAT // BITS_PER_NYBBLE)
+    num_bits_to_clear = MAX_BITS_PRECISION_IBM_FLOAT - num_nybbles_to_preserve * BITS_PER_NYBBLE
+    clear_mask = 2**num_bits_to_clear - 1
+    preserve_mask = (2**MAX_BITS_PRECISION_IBM_FLOAT - 1) & ~clear_mask
+
+    print("num_nybbles_to_preserve", num_nybbles_to_preserve)
+    print("num_bits_to_clear", num_bits_to_clear)
+    print("clear_mask   ", bin(clear_mask))
+    print("preserve_mask", bin(preserve_mask))
+
+    truncated_mantissa = mantissa & preserve_mask
+
+    value = truncated_mantissa * pow(16, exponent_16)
+
+    scaled_value = value >> MAX_BITS_PRECISION_IBM_FLOAT
+
+    return scaled_value
+
+    #
+    # tb = (preserve_mask >> 16) & b
+    # tc = (preserve_mask >> 8) & c
+    # td = preserve_mask & d
+    #
+    # return byte_string((a, tb, tc, td))
+
 
 def ieee2ibm(f):
-    """Covert a float to four big-endian bytes representing an IBM float.
+    """Convert a float to four big-endian bytes representing an IBM float.
 
     Args:
         f (float): The value to be converted.
@@ -112,3 +168,324 @@ def ieee2ibm(f):
     d = mantissa & 0xff
 
     return byte_string((a, b, c, d))
+
+
+class IBMFloat(Real):
+
+    __slots__ = ['_data']
+
+    _INTERNED = {IBM_ZERO_BYTES: None,
+                 IBM_NEGATIVE_ONE_BYTES: None,
+                 IBM_POSITIVE_ONE_BYTES: None}
+
+    # noinspection PyUnresolvedReferences
+    def __new__(cls, b):
+        obj = object.__new__(cls)
+
+        data = bytes(b)
+        num_bytes = len(data)
+        if num_bytes != 4:
+            raise ValueError("{} cannot be constructed from {} values".format(cls.__name__, num_bytes))
+        obj._data = data
+
+        # Intern common values
+        if data in cls._INTERNED:
+            if cls._INTERNED[data] is None:
+                cls._INTERNED[data] = obj
+            return cls._INTERNED[data]
+
+        return obj
+
+
+    @classmethod
+    def from_float(cls, f):
+        """Construct an IBMFloat from an IEEE float.
+
+        Args:
+            f (float): The value to be converted.
+
+        Returns:
+            An IBMFloat.
+
+        Raises:
+            OverflowError: If f is outside the representable range.
+            ValueError: If f is NaN or infinite.
+            FloatingPointError: If f cannot be represented without total loss of precision.
+        """
+        return cls(ieee2ibm(f))
+
+    @classmethod
+    def from_real(cls, f):
+        if isinstance(f, IBMFloat):
+            return f
+        return cls.from_float(f)
+
+    @classmethod
+    def from_bytes(cls, b):
+        return cls(b)
+
+    @classmethod
+    def ldexp(cls, fraction, exponent):
+        """Make an IBMFloat from fraction and exponent.
+
+        The is the inverse function of IBMFloat.frexp()
+
+        Args:
+            fraction: A Real in the range -1.0 to 1.0.
+            exponent: An integer in the range -256 to 255 inclusive.
+        """
+        if not (-1.0 <= fraction <= 1.0):
+            raise ValueError("ldexp fraction {!r} out of range -1.0 to +1.0")
+
+        if not (-256 <= exponent < 256):
+            raise ValueError("ldexp exponent {!r} out of range -256 to 256")
+
+        ieee = fraction * 2**exponent
+        return IBMFloat.from_float(ieee)
+
+    @property
+    def signbit(self):
+        """True if the value is negative, otherwise False."""
+        return bool(self._data[0] & 0x80)
+
+    def __float__(self):
+        return ibm2ieee(self._data)
+
+    def __bytes__(self):
+        return self._data
+
+    def __repr__(self):
+        return "{}.from_float({!r}) ~{!r}".format(self.__class__.__name__, self._data, float(self))
+
+    def __str__(self):
+        return str(float(self))
+
+    def __bool__(self):
+        return not self.is_zero()
+
+    def is_zero(self):
+        return self.int_mantissa == 0
+
+    def __nonzero__(self):
+        return not self.is_zero()
+
+    def is_subnormal(self):
+        if self.is_zero():
+            # Only one of the many possible representations of zero is considered 'normal' - all the zeros
+            return not all(b == 0 for b in self._data)
+
+        return self._data[1] < 16  # TODO: Replace magic number with constant
+
+    def zero_subnormal(self):
+        return IBM_FLOAT_ZERO if self.is_subnormal() else self
+
+    def frexp(self):
+        """Obtain the fraction and exponent.
+
+        Returns:
+            A pair where the first item is the fraction in the range -1.0 and +1.0 and the
+            exponent is an integer such that f = fraction * 2**exponent
+        """
+        sign = -1 if self.signbit else 1
+        mantissa = sign * self.int_mantissa / _F24
+        exp_2 = self.exp16 * 4
+        return mantissa, exp_2
+
+    def __pos__(self):
+        return self
+
+    def __neg__(self):
+        if self.is_zero():
+            return IBM_FLOAT_ZERO
+
+        data = self._data
+        return IBMFloat((data[0] ^ 0b10000000,
+                         data[1],
+                         data[2],
+                         data[3]))
+
+    def __abs__(self):
+        if self.is_zero():
+            return IBM_FLOAT_ZERO
+
+        data = self._data
+        return IBMFloat((data[0] & 0b01111111,
+                         data[1],
+                         data[2],
+                         data[3]))
+
+    def __eq__(self, rhs):
+        lhs = self
+
+        if not isinstance(rhs, IBMFloat):
+            return float(lhs) == float(rhs)
+
+        lhs_sign = lhs.signbit
+        rhs_sign = rhs.signbit
+
+        if lhs_sign != rhs_sign:
+            return False
+
+        nlhs = lhs.normalize()
+        nrhs = rhs.normalize()
+
+        if not (nlhs.is_subnormal() or nrhs.is_subnormal()):
+            # Both of the numbers are normalised
+            return nlhs._data == nrhs._data
+
+        # Either or both of the numbers are subnormal
+        lhs_exp16 = nlhs.exp16
+        rhs_exp16 = nrhs.exp16
+
+        lhs_mantissa = nlhs.int_mantissa
+        rhs_mantissa = nrhs.int_mantissa
+
+        if lhs_exp16 < rhs_exp16:
+            delta_exp16 = rhs_exp16 - lhs_exp16
+            lhs_mantissa >>= 4 * delta_exp16
+            lhs_exp16 += delta_exp16
+
+        if lhs_exp16 > rhs_exp16:
+            delta_exp16 = lhs_exp16 - rhs_exp16
+            rhs_mantissa >>= 4 * delta_exp16
+            rhs_exp16 += delta_exp16
+
+        assert lhs_exp16 == rhs_exp16
+        return lhs_mantissa == rhs_mantissa
+
+    def __floordiv__(self, rhs):
+        return float(self) // float(rhs)
+
+    def __rfloordiv__(self, lhs):
+        return float(lhs) // float(self)
+
+    def __rtruediv__(self, lhs):
+        q = float(lhs) / float(self)
+        return IBMFloat.from_float(q) if isinstance(lhs, float) else q
+
+    def __pow__(self, exponent):
+        p = pow(float(self), float(exponent))
+        return IBMFloat.from_float(p) if isinstance(exponent, IBMFloat) else p
+
+    def __rpow__(self, base):
+        return IBMFloat.from_float(pow(float(base), float(self)))
+
+    def __mod__(self, rhs):
+        m = float(self) % float(rhs)
+        return IBMFloat.from_float(m) if isinstance(rhs, IBMFloat) else m
+
+    def __rmod__(self, lhs):
+        m = float(lhs) % float(self)
+        return IBMFloat.from_float(m) if isinstance(lhs, IBMFloat) else m
+
+    def __rmul__(self, lhs):
+        p = float(lhs) * float(self)
+        return IBMFloat.from_float(p) if isinstance(lhs, IBMFloat) else p
+
+    def __radd__(self, lhs):
+        s = float(lhs) + float(self)
+        return IBMFloat.from_float(s) if isinstance(lhs, IBMFloat) else s
+
+    def __lt__(self, rhs):
+        return float(self) < float(rhs)
+
+    def __le__(self, rhs):
+        return float(self) <= float(rhs)
+
+    def __gt__(self, rhs):
+        return float(self) > float(rhs)
+
+    def __ge__(self, rhs):
+        return float(self) >= float(rhs)
+
+    def __ceil__(self):
+        t = trunc(self)
+        return t if self.signbit else t + 1
+
+    def __floor__(self):
+        t = trunc(self)
+        return t - 1 if self.signbit else t
+
+    @property
+    def exp16(self):
+        """The base 16 exponent."""
+        exponent_16_biased = self._data[0] & 0x7f
+        exponent_16 = exponent_16_biased - EXPONENT_BIAS
+        return exponent_16
+
+    @property
+    def int_mantissa(self):
+        data = self._data
+        return (data[1] << 16) | (data[2] << 8) | data[3]
+
+    def __trunc__(self):
+        sign = -1 if self.signbit else 1
+        exponent_16 = self.exp16
+        mantissa = self.int_mantissa
+
+        num_nybbles_to_preserve = min(exponent_16, MAX_BITS_PRECISION_IBM_FLOAT // BITS_PER_NYBBLE)
+        num_bits_to_clear = MAX_BITS_PRECISION_IBM_FLOAT - num_nybbles_to_preserve * BITS_PER_NYBBLE
+        clear_mask = 2**num_bits_to_clear - 1
+        preserve_mask = (2**MAX_BITS_PRECISION_IBM_FLOAT - 1) & ~clear_mask
+
+        truncated_mantissa = mantissa & preserve_mask
+        magnitude = truncated_mantissa * pow(16, exponent_16) >> MAX_BITS_PRECISION_IBM_FLOAT
+        return sign * magnitude
+
+    def normalize(self):
+        """Normalize the floating point value.
+
+        Returns:
+            A normalized IBMFloat equal in value to this object.
+
+        Raises:
+            FloatingPointError: If the number could not be normalized.
+        """
+        if self.is_zero():
+            return IBM_FLOAT_ZERO
+
+        exponent_16 = self.exp16
+        mantissa = self.int_mantissa
+
+        while mantissa < (1 << 20):
+            new_exponent_16 = exponent_16 - 1
+            if not (-64 <= new_exponent_16 < 64):
+                raise FloatingPointError("Could not normalize {!r} without causing exponent overflow.".format(self))
+
+            mantissa <<= 4
+            exponent_16 = new_exponent_16
+
+        exponent_16_biased = exponent_16 + EXPONENT_BIAS
+
+        sign = int(self.signbit) << 7
+
+        a = sign | exponent_16_biased
+        b = (mantissa >> 16) & 0xff
+        c = (mantissa >> 8) & 0xff
+        d = mantissa & 0xff
+
+        return IBMFloat.from_bytes((a, b, c, d))
+
+    def __round__(self, ndigits=None):
+        return IBMFloat.from_float(round(float(self), ndigits))
+
+    def __truediv__(self, rhs):
+        q = float(self) / float(rhs)
+        return IBMFloat.from_float(q) if isinstance(rhs, IBMFloat) else q
+
+    def __mul__(self, rhs):
+        p = float(self) * float(rhs)
+        return IBMFloat.from_float(p) if isinstance(rhs, IBMFloat) else p
+
+    def __add__(self, rhs):
+        p = float(self) + float(rhs)
+        return IBMFloat.from_float(p) if isinstance(rhs, IBMFloat) else p
+
+    def __int__(self):
+        return trunc(self)
+
+
+IBM_FLOAT_ZERO = IBMFloat.from_bytes(IBM_ZERO_BYTES)
+
+
+
diff --git a/segpy/toolkit.py b/segpy/toolkit.py
index 293903c..1769c8a 100644
--- a/segpy/toolkit.py
+++ b/segpy/toolkit.py
@@ -15,7 +15,7 @@ from segpy.catalog import CatalogBuilder
 from segpy.datatypes import CTYPES, size_in_bytes
 from segpy.encoding import guess_encoding, is_supported_encoding, UnsupportedEncodingError
 from segpy.binary_reel_header_definition import HEADER_DEF
-from segpy.ibm_float import ibm2ieee, ieee2ibm
+from segpy.ibm_float import IBMFloat
 from segpy.revisions import canonicalize_revision
 from segpy.trace_header_definition import TRACE_HEADER_DEF
 from segpy.util import file_length, batched, pad, complementary_intervals, NATIVE_ENDIANNESS
@@ -461,7 +461,7 @@ def unpack_ibm_floats(data, count):
     Returns:
         A sequence of floats.
     """
-    return array('d', (ibm2ieee(data[i: i+4]) for i in range(0, count * 4, 4)))
+    return [IBMFloat.from_bytes(data[i: i+4]) for i in range(0, count * 4, 4)]
 
 
 def unpack_values(buf, count, fmt, endian='>'):
@@ -817,7 +817,7 @@ def pack_ibm_floats(values):
         A sequence of bytes. (Python 2 - a str object, Python 3 - a bytes
             object)
     """
-    return EMPTY_BYTE_STRING.join(ieee2ibm(value) for value in values)
+    return EMPTY_BYTE_STRING.join(bytes(IBMFloat.from_real(value)) for value in values)
 
 
 def pack_values(values, fmt, endian='>'):
diff --git a/segpy/util.py b/segpy/util.py
index 20d00e2..9ac147e 100644
--- a/segpy/util.py
+++ b/segpy/util.py
@@ -240,4 +240,11 @@ def first_sentence(s):
 
 def lower_first(s):
     """Lower case the first character of a string."""
-    return s[:1].lower() + s[1:]
\ No newline at end of file
+    return s[:1].lower() + s[1:]
+
+
+def almost_equal(x, y, epsilon):
+    max_xy_one = max(1.0, abs(x), abs(y))
+    e = epsilon * max_xy_one
+    delta = abs(x - y)
+    return delta <= e
\ No newline at end of file
diff --git a/test/strategies.py b/test/strategies.py
new file mode 100644
index 0000000..c56725d
--- /dev/null
+++ b/test/strategies.py
@@ -0,0 +1,21 @@
+from hypothesis import strategy
+from hypothesis.specifiers import integers_in_range
+
+PRINTABLE_ASCII_RANGE = (32, 127)
+
+
+def multiline_ascii_encodable_text(min_num_lines, max_num_lines):
+    """A Hypothesis strategy to produce a multiline Unicode string.
+
+    Args:
+        min_num_lines: The minimum number of lines in the produced strings.
+        max_num_lines: The maximum number of lines in the produced strings.
+
+    Returns:
+        A strategy for generating Unicode strings containing only newlines
+        and characters which are encodable as printable 7-bit ASCII characters.
+    """
+
+    return strategy(integers_in_range(0, 10))                                      \
+           .flatmap(lambda n: ([integers_in_range(*PRINTABLE_ASCII_RANGE)],) * n)  \
+           .map(lambda xs: '\n'.join(bytes(x).decode('ascii') for x in xs))
diff --git a/test/test-requirements.txt b/test/test-requirements.txt
index fd4f18c..bcf2f58 100644
--- a/test/test-requirements.txt
+++ b/test/test-requirements.txt
@@ -1 +1 @@
-hypothesis==0.4.0
+hypothesis>=1.2
diff --git a/test/test_extended_textual_header.py b/test/test_extended_textual_header.py
index 21bd2b3..81adb86 100644
--- a/test/test_extended_textual_header.py
+++ b/test/test_extended_textual_header.py
@@ -1,57 +1,37 @@
 import unittest
+
 from hypothesis import given
-from hypothesis.descriptors import one_of, SampledFrom, Just, sampled_from, just
-from hypothesis.searchstrategy import MappedSearchStrategy, StringStrategy
-from hypothesis.strategytable import StrategyTable
+from hypothesis.specifiers import sampled_from, just
+
 from segpy.encoding import EBCDIC, ASCII
 from segpy.toolkit import format_extended_textual_header, CARDS_PER_HEADER, END_TEXT_STANZA, CARD_LENGTH
-from segpy.portability import unicode
-
-
-class MultiLineString(unicode):
-    pass
-
-
-class MultiLineStringStrategy(MappedSearchStrategy):
-
-    def pack(self, x):
-        return MultiLineString(unicode('\n').join(x))
-
-    def unpack(self, x):
-        return x.splitlines()
-
-
-StrategyTable.default().define_specification_for(
-    MultiLineString,
-    lambda s, d: MultiLineStringStrategy(
-        strategy=s.strategy([unicode]),
-        descriptor=MultiLineString))
+from test.strategies import multiline_ascii_encodable_text
 
 
 class TestFormatExtendedTextualHeader(unittest.TestCase):
 
-    @given(MultiLineString,
+    @given(multiline_ascii_encodable_text(0, 100),
            sampled_from([ASCII, EBCDIC]),
            bool)
     def test_forty_lines_per_page(self, text, encoding, include_text_stop):
         pages = format_extended_textual_header(text, encoding, include_text_stop)
         self.assertTrue(all(len(page) == CARDS_PER_HEADER for page in pages))
 
-    @given(MultiLineString,
+    @given(multiline_ascii_encodable_text(0, 100),
            sampled_from([ASCII, EBCDIC]),
            bool)
     def test_eighty_bytes_per_encoded_line(self, text, encoding, include_text_stop):
         pages = format_extended_textual_header(text, encoding, include_text_stop)
         self.assertTrue(all([len(line.encode(encoding)) == CARD_LENGTH for page in pages for line in page]))
 
-    @given(MultiLineString,
+    @given(multiline_ascii_encodable_text(0, 100),
            sampled_from([ASCII, EBCDIC]),
            bool)
     def test_lines_end_with_cr_lf(self, text, encoding, include_text_stop):
         pages = format_extended_textual_header(text, encoding, include_text_stop)
         self.assertTrue(all([line.endswith('\r\n') for page in pages for line in page]))
 
-    @given(MultiLineString,
+    @given(multiline_ascii_encodable_text(0, 100),
            sampled_from([ASCII, EBCDIC]),
            just(True))
     def test_end_text_stanza_present(self, text, encoding, include_text_stop):
diff --git a/test/test_float.py b/test/test_float.py
index 2388dab..e345b87 100644
--- a/test/test_float.py
+++ b/test/test_float.py
@@ -1,8 +1,15 @@
+from math import trunc, floor
 import unittest
 
+from hypothesis import given, assume
+import math
+from hypothesis.specifiers import integers_in_range, floats_in_range
+
 from segpy.portability import byte_string
 from segpy.ibm_float import (ieee2ibm, ibm2ieee, MAX_IBM_FLOAT, SMALLEST_POSITIVE_NORMAL_IBM_FLOAT,
-                             LARGEST_NEGATIVE_NORMAL_IBM_FLOAT, MIN_IBM_FLOAT)
+                             LARGEST_NEGATIVE_NORMAL_IBM_FLOAT, MIN_IBM_FLOAT, IBMFloat, EPSILON_IBM_FLOAT,
+                             MAX_EXACT_INTEGER_IBM_FLOAT, MIN_EXACT_INTEGER_IBM_FLOAT, EXPONENT_BIAS)
+from segpy.util import almost_equal
 
 
 class Ibm2Ieee(unittest.TestCase):
@@ -134,5 +141,224 @@ class Ibm2IeeeRoundtrip(unittest.TestCase):
         self.assertEqual(ibm_start, ibm_result)
 
 
+class TestIBMFloat(unittest.TestCase):
+
+    def test_zero_from_float(self):
+        zero = IBMFloat.from_float(0.0)
+        self.assertTrue(zero.is_zero())
+
+    def test_zero_from_bytes(self):
+        zero = IBMFloat.from_bytes(b'\x00\x00\x00\x00')
+        self.assertTrue(zero.is_zero())
+
+    def test_subnormal(self):
+        ibm = IBMFloat.from_float(1.6472184286297693e-83)
+        self.assertTrue(ibm.is_subnormal())
+
+    def test_smallest_subnormal(self):
+        ibm = IBMFloat.from_float(5.147557589468029e-85)
+        self.assertEqual(bytes(ibm), byte_string((0x00, 0x00, 0x00, 0x01)))
+
+    def test_too_small_subnormal(self):
+        with self.assertRaises(FloatingPointError):
+            IBMFloat.from_float(1e-86)
+
+    def test_nan(self):
+        with self.assertRaises(ValueError):
+            IBMFloat.from_float(float('nan'))
+
+    def test_inf(self):
+        with self.assertRaises(ValueError):
+            IBMFloat.from_float(float('inf'))
+
+    def test_too_large(self):
+        with self.assertRaises(OverflowError):
+            IBMFloat.from_float(MAX_IBM_FLOAT * 10)
+
+    def test_too_small(self):
+        with self.assertRaises(OverflowError):
+            IBMFloat.from_float(MIN_IBM_FLOAT * 10)
+
+    @given(floats_in_range(MIN_IBM_FLOAT, MAX_IBM_FLOAT))
+    def test_bool(self, f):
+        self.assertEqual(bool(IBMFloat.from_float(f)), bool(f))
+
+    @given(integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255))
+    def test_bytes_roundtrip(self, a, b, c, d):
+        b = byte_string((a, b, c, d))
+        ibm = IBMFloat.from_bytes(b)
+        self.assertEqual(bytes(ibm), b)
+
+    @given(floats_in_range(MIN_IBM_FLOAT, MAX_IBM_FLOAT))
+    def test_floats_roundtrip(self, f):
+        ibm = IBMFloat.from_float(f)
+        self.assertTrue(almost_equal(f, float(ibm), epsilon=EPSILON_IBM_FLOAT))
+
+    @given(integers_in_range(0, MAX_EXACT_INTEGER_IBM_FLOAT - 1),
+           floats_in_range(0.0, 1.0))
+    def test_trunc_above_zero(self, i, f):
+        assume(f != 1.0)
+        ieee = i + f
+        ibm = IBMFloat.from_float(ieee)
+        self.assertEqual(trunc(ibm), i)
+
+    @given(integers_in_range(MIN_EXACT_INTEGER_IBM_FLOAT + 1, 0),
+           floats_in_range(0.0, 1.0))
+    def test_trunc_below_zero(self, i, f):
+        assume(f != 1.0)
+        ieee = i - f
+        ibm = IBMFloat.from_float(ieee)
+        self.assertEqual(trunc(ibm), i)
+
+    @given(integers_in_range(MIN_EXACT_INTEGER_IBM_FLOAT, MAX_EXACT_INTEGER_IBM_FLOAT - 1),
+           floats_in_range(0.0, 1.0))
+    def test_ceil(self, i, f):
+        assume(f != 1.0)
+        ieee = i + f
+        ibm = IBMFloat.from_float(ieee)
+        self.assertEqual(math.ceil(ibm), i + 1)
+
+    @given(integers_in_range(MIN_EXACT_INTEGER_IBM_FLOAT, MAX_EXACT_INTEGER_IBM_FLOAT - 1),
+           floats_in_range(0.0, 1.0))
+    def test_floor(self, i, f):
+        assume(f != 1.0)
+        ieee = i + f
+        ibm = IBMFloat.from_float(ieee)
+        self.assertEqual(math.floor(ibm), i)
+
+    def test_normalise_subnormal_expect_failure(self):
+        # This float has an base-16 exponent of -64 (the minimum) and cannot be normalised
+        ibm = IBMFloat.from_float(1.6472184286297693e-83)
+        assert ibm.is_subnormal()
+        with self.assertRaises(FloatingPointError):
+            ibm.normalize()
+
+    def test_normalise_subnormal1(self):
+        ibm = IBMFloat.from_bytes((0b01000000, 0b00000000, 0b11111111, 0b00000000))
+        assert ibm.is_subnormal()
+        normalized = ibm.normalize()
+        self.assertFalse(normalized.is_subnormal())
+
+    def test_normalise_subnormal2(self):
+        ibm = IBMFloat.from_bytes((64, 1, 0, 0))
+        assert ibm.is_subnormal()
+        normalized = ibm.normalize()
+        self.assertFalse(normalized.is_subnormal())
+
+    @given(integers_in_range(128, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(4, 23))
+    def test_normalise_subnormal(self, b, c, d, shift):
+        mantissa = (b << 16) | (c << 8) | d
+        assume(mantissa != 0)
+        mantissa >>= shift
+        assert mantissa != 0
+
+        sa = EXPONENT_BIAS
+        sb = (mantissa >> 16) & 0xff
+        sc = (mantissa >> 8) & 0xff
+        sd = mantissa & 0xff
+
+        ibm = IBMFloat.from_bytes((sa, sb, sc, sd))
+        assert ibm.is_subnormal()
+        normalized = ibm.normalize()
+        self.assertFalse(normalized.is_subnormal())
+
+    @given(integers_in_range(128, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(4, 23))
+    def test_zero_subnormal(self, b, c, d, shift):
+        mantissa = (b << 16) | (c << 8) | d
+        assume(mantissa != 0)
+        mantissa >>= shift
+        assert mantissa != 0
+
+        sa = EXPONENT_BIAS
+        sb = (mantissa >> 16) & 0xff
+        sc = (mantissa >> 8) & 0xff
+        sd = mantissa & 0xff
+
+        ibm = IBMFloat.from_bytes((sa, sb, sc, sd))
+        assert ibm.is_subnormal()
+        z = ibm.zero_subnormal()
+        self.assertTrue(z.is_zero())
+
+    @given(integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255))
+    def test_abs(self, a, b, c, d):
+        ibm = IBMFloat.from_bytes((a, b, c, d))
+        abs_ibm = abs(ibm)
+        self.assertGreaterEqual(abs_ibm.signbit, 0)
+
+    @given(integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255),
+           integers_in_range(0, 255))
+    def test_negate_non_zero(self, a, b, c, d):
+        ibm = IBMFloat.from_bytes((a, b, c, d))
+        assume(not ibm.is_zero())
+        negated = -ibm
+        self.assertNotEqual(ibm.signbit, negated.signbit)
+
+    def test_negate_zero(self):
+        zero = IBMFloat.from_float(0.0)
+        negated = -zero
+        self.assertTrue(negated.is_zero())
+
+    @given(floats_in_range(MIN_IBM_FLOAT, MAX_IBM_FLOAT))
+    def test_signbit(self, f):
+        ltz = f < 0
+        ibm = IBMFloat.from_float(f)
+        self.assertEqual(ltz, ibm.signbit)
+
+    @given(floats_in_range(-1.0, +1.0),
+           integers_in_range(-256, 255))
+    def test_ldexp_frexp(self, fraction, exponent):
+        try:
+            ibm = IBMFloat.ldexp(fraction, exponent)
+        except OverflowError:
+            assume(False)
+        else:
+            f, e = ibm.frexp()
+            self.assertTrue(almost_equal(fraction * 2**exponent, f * 2**e, epsilon=EPSILON_IBM_FLOAT))
+
+    @given(floats_in_range(MIN_IBM_FLOAT, MAX_IBM_FLOAT),
+           floats_in_range(0.0, 1.0))
+    def test_add(self, f, p):
+        a = f * p
+        b = f - a
+
+        ibm_a = IBMFloat.from_float(a)
+        ibm_b = IBMFloat.from_float(b)
+        ibm_c = ibm_a + ibm_b
+
+        ieee_a = float(ibm_a)
+        ieee_b = float(ibm_b)
+        ieee_c = ieee_a + ieee_b
+
+        self.assertTrue(almost_equal(ieee_c, ibm_c, epsilon=EPSILON_IBM_FLOAT * 4))
+
+    @given(floats_in_range(0, MAX_IBM_FLOAT),
+           floats_in_range(0, MAX_IBM_FLOAT))
+    def test_sub(self, a, b):
+        ibm_a = IBMFloat.from_float(a)
+        ibm_b = IBMFloat.from_float(b)
+        ibm_c = ibm_a - ibm_b
+
+        ieee_a = float(ibm_a)
+        ieee_b = float(ibm_b)
+        ieee_c = ieee_a - ieee_b
+
+        self.assertTrue(almost_equal(ieee_c, ibm_c, epsilon=EPSILON_IBM_FLOAT))
+
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/test/test_util.py b/test/test_util.py
index 1b71b94..47d19af 100644
--- a/test/test_util.py
+++ b/test/test_util.py
@@ -1,7 +1,7 @@
 import unittest
 
 from hypothesis import given, assume
-from hypothesis.descriptors import integers_in_range
+from hypothesis.specifiers import integers_in_range
 from segpy.util import batched