BF: try other checks for badly-detected longdouble

Numpy gives wrong values from the np.finfo routine for PPC, and I put a
hack in to detect that.  It looks like numpy also gives wrong values on
Debian on the IBM POWER7, and that might have a real binary128
longdouble. Try to make the checks more flexible to recognize both
cases.

Author: matthew-brett

nibabel/casting.py

@@ -244,13 +244,91 @@ def type_info(np_type):
         # at float64 precision. The finfo values give nexp == 15 (as for intel
         # 80) but in calculations nexp in fact appears to be 11 as for float64
         ret.update(dict(width=width))
-    elif vals == (1, 1, 16) and processor() == 'powerpc': # broken PPC
-        ret.update(dict(nmant=106, width=width))
+    # Oh dear, we don't recognize the type information.  Try some known types
+    # and then give up. At this stage we're expecting exotic longdouble or their
+    # complex equivalent.
+    if not np_type in (np.longdouble, np.longcomplex) or width not in (16, 32):
+        raise FloatingError('We had not expected type %s' % np_type)
+    if (vals == (1, 1, 16) and processor() == 'powerpc' and
+        _check_maxexp(np.longdouble, 1024)):
+        # double pair on PPC.  The _check_nmant routine does not work for this
+        # type, hence the processor check
+        ret.update(dict(nmant = 106, width=width))
+    elif (_check_nmant(np.longdouble, 112) and
+          _check_maxexp(np.longdouble, 16384)):
+        # binary 128, but with some busted type information. np.longcomplex
+        # seems to break here too, so we need to use np.longdouble and
+        # complexify
+        two = np.longdouble(2)
+        # See: http://matthew-brett.github.com/pydagogue/floating_point.html
+        max_val = (two ** 113 - 1) / (two ** 112) * two ** 16383
+        if np_type is np.longcomplex:
+            max_val += 0j
+        ret = dict(min = -max_val,
+                   max= max_val,
+                   nmant = 112,
+                   nexp = 15,
+                   minexp = -16382,
+                   maxexp = 16384,
+                   width = width)
     else: # don't recognize the type
         raise FloatingError('We had not expected type %s' % np_type)
     return ret
 
 
+def _check_nmant(np_type, nmant):
+    """ True if fp type `np_type` seems to have `nmant` significand digits
+
+    Note 'digits' does not include implicit digits.  And in fact if there are no
+    implicit digits, the `nmant` number is one less than the actual digits.
+    Assumes base 2 representation.
+
+    Parameters
+    ----------
+    np_type : numpy type specifier
+        Any specifier for a numpy dtype
+    nmant : int
+        Number of digits to test against
+
+    Returns
+    -------
+    tf : bool
+        True if `nmant` is the correct number of significand digits, false
+        otherwise
+    """
+    np_type = np.dtype(np_type).type
+    max_contig = np_type(2 ** (nmant + 1)) # maximum of contiguous integers
+    tests = max_contig + np.array([-2, -1, 0, 1, 2], dtype=np_type)
+    return np.all(tests - max_contig == [-2, -1, 0, 0, 2])
+
+
+def _check_maxexp(np_type, maxexp):
+    """ True if fp type `np_type` seems to have `maxexp` maximum exponent
+
+    We're testing "maxexp" as returned by numpy. This value is set to one
+    greater than the maximum power of 2 that `np_type` can represent.
+
+    Assumes base 2 representation.  Very crude check
+
+    Parameters
+    ----------
+    np_type : numpy type specifier
+        Any specifier for a numpy dtype
+    maxexp : int
+        Maximum exponent to test against
+
+    Returns
+    -------
+    tf : bool
+        True if `maxexp` is the correct maximum exponent, False otherwise.
+    """
+    dt = np.dtype(np_type)
+    np_type = dt.type
+    two = np_type(2).reshape((1,)) # to avoid upcasting
+    return (np.isfinite(two ** (maxexp - 1)) and
+            not np.isfinite(two ** maxexp))
+
+
 def as_int(x, check=True):
     """ Return python integer representation of number
 

nibabel/tests/test_floating.py

@@ -5,10 +5,11 @@
 import numpy as np
 
 from ..casting import (floor_exact, ceil_exact, as_int, FloatingError,
-                       int_to_float, floor_log2, type_info)
+                       int_to_float, floor_log2, type_info, _check_nmant,
+                       _check_maxexp, ok_floats)
 
 from nose import SkipTest
-from nose.tools import assert_equal, assert_raises
+from nose.tools import assert_equal, assert_raises, assert_true, assert_false
 
 IEEE_floats = [np.float32, np.float64]
 try:
@@ -80,6 +81,25 @@ def test_nmant():
         assert_equal(type_info(np.longdouble)['nmant'], 63)
 
 
+def test_check_nmant_nexp():
+    # Routine for checking number of sigificand digits and exponent
+    for t in IEEE_floats:
+        nmant = np.finfo(t).nmant
+        maxexp = np.finfo(t).maxexp
+        assert_true(_check_nmant(t, nmant))
+        assert_false(_check_nmant(t, nmant - 1))
+        assert_false(_check_nmant(t, nmant + 1))
+        assert_true(_check_maxexp(t, maxexp))
+        assert_false(_check_maxexp(t, maxexp - 1))
+        assert_false(_check_maxexp(t, maxexp + 1))
+    # Check against type_info
+    for t in ok_floats():
+        ti = type_info(t)
+        if ti['nmant'] != 106: # This check does not work for PPC double pair
+            assert_true(_check_nmant(t, ti['nmant']))
+        assert_true(_check_maxexp(t, ti['maxexp']))
+
+
 def test_as_int():
     # Integer representation of number
     assert_equal(as_int(2.0), 2)