summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--doc/release/1.13.0-notes.rst16
-rw-r--r--numpy/core/getlimits.py281
-rw-r--r--numpy/core/tests/test_getlimits.py37
3 files changed, 308 insertions, 26 deletions
diff --git a/doc/release/1.13.0-notes.rst b/doc/release/1.13.0-notes.rst
index 69e2a6b5b..ac4d4b2c1 100644
--- a/doc/release/1.13.0-notes.rst
+++ b/doc/release/1.13.0-notes.rst
@@ -105,6 +105,22 @@ Performance improvements for ``packbits`` and ``unpackbits``
The functions ``numpy.packbits`` with boolean input and ``numpy.unpackbits`` have
been optimized to be a significantly faster for contiguous data.
+Fix for PPC long double floating point information
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+In previous versions of numpy, the ``finfo`` function returned invalid
+information about the `double double`_ format of the ``longdouble`` float type
+on Power PC (PPC). The invalid values resulted from the failure of the numpy
+algorithm to deal with the `variable number of digits in the significand
+<https://www.ibm.com/support/knowledgecenter/en/ssw_aix_71/com.ibm.aix.genprogc/128bit_long_double_floating-point_datatype.htm>`_
+that are a feature of PPC long doubles. This release by-passes the failing
+algorithm by using heuristics to detect the presence of the PPC double double
+format. A side-effect of using these heuristics is that the ``finfo``
+function is faster than previous releases.
+
+.. _issues: https://github.com/numpy/numpy/issues/2669
+
+.. _double double: https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic
+
Changes
=======
diff --git a/numpy/core/getlimits.py b/numpy/core/getlimits.py
index d4025cb3b..5294f4682 100644
--- a/numpy/core/getlimits.py
+++ b/numpy/core/getlimits.py
@@ -5,10 +5,15 @@ from __future__ import division, absolute_import, print_function
__all__ = ['finfo', 'iinfo']
+import warnings
+
from .machar import MachAr
from . import numeric
from . import numerictypes as ntypes
-from .numeric import array
+from .numeric import array, inf
+from .umath import log10, exp2
+from . import umath
+
def _frz(a):
"""fix rank-0 --> rank-1"""
@@ -16,12 +21,261 @@ def _frz(a):
a.shape = (1,)
return a
+
_convert_to_float = {
ntypes.csingle: ntypes.single,
ntypes.complex_: ntypes.float_,
ntypes.clongfloat: ntypes.longfloat
}
+
+# Parameters for creating MachAr / MachAr-like objects
+_MACHAR_PARAMS = {
+ ntypes.double: dict(
+ itype = ntypes.int64,
+ fmt = '%24.16e',
+ precname = 'double'),
+ ntypes.single: dict(
+ itype = ntypes.int32,
+ fmt = '%15.7e',
+ precname = 'single'),
+ ntypes.longdouble: dict(
+ itype = ntypes.longlong,
+ fmt = '%s',
+ precname = 'long double'),
+ ntypes.half: dict(
+ itype = ntypes.int16,
+ fmt = '%12.5e',
+ precname = 'half')}
+
+
+class MachArLike(object):
+ """ Object to simulate MachAr instance """
+
+ # These attributes should be of characteristic dtype
+ _native_attrs = ('tiny', 'huge', 'epsneg', 'eps')
+
+ def __init__(self,
+ ftype,
+ **kwargs):
+ params = _MACHAR_PARAMS[ftype]
+ float_conv = lambda v: array([v], ftype)
+ float_to_str = lambda v: params['fmt'] % array(_frz(v)[0], ftype)
+ self.title = 'numpy {} precision floating point number'.format(
+ params['precname'])
+ for key, value in kwargs.items():
+ if key in self._native_attrs:
+ value = float_conv(value)
+ self.__dict__[key] = value
+ self.epsilon = self.eps
+ self.xmax = self.huge
+ self.xmin = self.tiny
+ self.precision = int(-log10(self.eps))
+ self.resolution = float_conv(10) ** (-self.precision)
+ self._str_eps = float_to_str(self.eps)
+ self._str_epsneg = float_to_str(self.epsneg)
+ self._str_xmin = float_to_str(self.xmin)
+ self._str_xmax = float_to_str(self.xmax)
+ self._str_resolution = float_to_str(self.resolution)
+
+
+# Known parameters for float16
+# See docstring of MachAr class for description of parameters.
+_f16 = ntypes.float16
+_float16_ma = MachArLike(_f16,
+ machep=-10,
+ negep=-11,
+ minexp=-14,
+ maxexp=16,
+ it=10,
+ iexp=5,
+ ibeta=2,
+ irnd=5,
+ ngrd=0,
+ eps=exp2(_f16(-10)),
+ epsneg=exp2(_f16(-11)),
+ huge=_f16(65504),
+ tiny=_f16(2 ** -14))
+
+# Known parameters for float32
+_f32 = ntypes.float32
+_float32_ma = MachArLike(_f32,
+ machep=-23,
+ negep=-24,
+ minexp=-126,
+ maxexp=128,
+ it=23,
+ iexp=8,
+ ibeta=2,
+ irnd=5,
+ ngrd=0,
+ eps=exp2(_f32(-23)),
+ epsneg=exp2(_f32(-24)),
+ huge=_f32((1 - 2 ** -24) * 2**128),
+ tiny=exp2(_f32(-126)))
+
+# Known parameters for float64
+_f64 = ntypes.float64
+_epsneg_f64 = 2.0 ** -53.0
+_tiny_f64 = 2.0 ** -1022.0
+_float64_ma = MachArLike(_f64,
+ machep=-52,
+ negep=-53,
+ minexp=-1022,
+ maxexp=1024,
+ it=52,
+ iexp=11,
+ ibeta=2,
+ irnd=5,
+ ngrd=0,
+ eps=2.0 ** -52.0,
+ epsneg=_epsneg_f64,
+ huge=(1.0 - _epsneg_f64) / _tiny_f64 * _f64(4),
+ tiny=_tiny_f64)
+
+# Known parameters for IEEE 754 128-bit binary float
+_ld = ntypes.longdouble
+_epsneg_f128 = exp2(_ld(-113))
+_tiny_f128 = exp2(_ld(-16382))
+# Ignore runtime error when this is not f128
+with numeric.errstate(all='ignore'):
+ _huge_f128 = (_ld(1) - _epsneg_f128) / _tiny_f128 * _ld(4)
+_float128_ma = MachArLike(_ld,
+ machep=-112,
+ negep=-113,
+ minexp=-16382,
+ maxexp=16384,
+ it=112,
+ iexp=15,
+ ibeta=2,
+ irnd=5,
+ ngrd=0,
+ eps=exp2(_ld(-112)),
+ epsneg=_epsneg_f128,
+ huge=_huge_f128,
+ tiny=_tiny_f128)
+
+# Known parameters for float80 (Intel 80-bit extended precision)
+_epsneg_f80 = exp2(_ld(-64))
+_tiny_f80 = exp2(_ld(-16382))
+# Ignore runtime error when this is not f80
+with numeric.errstate(all='ignore'):
+ _huge_f80 = (_ld(1) - _epsneg_f80) / _tiny_f80 * _ld(4)
+_float80_ma = MachArLike(_ld,
+ machep=-63,
+ negep=-64,
+ minexp=-16382,
+ maxexp=16384,
+ it=63,
+ iexp=15,
+ ibeta=2,
+ irnd=5,
+ ngrd=0,
+ eps=exp2(_ld(-63)),
+ epsneg=_epsneg_f80,
+ huge=_huge_f80,
+ tiny=_tiny_f80)
+
+# Guessed / known parameters for double double; see:
+# https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic
+# These numbers have the same exponent range as float64, but extended number of
+# digits in the significand.
+_huge_dd = (umath.nextafter(_ld(inf), _ld(0))
+ if hasattr(umath, 'nextafter') # Missing on some platforms?
+ else _float64_ma.huge)
+_float_dd_ma = MachArLike(_ld,
+ machep=-105,
+ negep=-106,
+ minexp=-1022,
+ maxexp=1024,
+ it=105,
+ iexp=11,
+ ibeta=2,
+ irnd=5,
+ ngrd=0,
+ eps=exp2(_ld(-105)),
+ epsneg= exp2(_ld(-106)),
+ huge=_huge_dd,
+ tiny=exp2(_ld(-1022)))
+
+
+# Key to identify the floating point type. Key is result of
+# ftype('-0.1').newbyteorder('<').tobytes()
+# See:
+# https://perl5.git.perl.org/perl.git/blob/3118d7d684b56cbeb702af874f4326683c45f045:/Configure
+_KNOWN_TYPES = {
+ b'\x9a\x99\x99\x99\x99\x99\xb9\xbf' : _float64_ma,
+ b'\xcd\xcc\xcc\xbd' : _float32_ma,
+ b'f\xae' : _float16_ma,
+ # float80, first 10 bytes containing actual storage
+ b'\xcd\xcc\xcc\xcc\xcc\xcc\xcc\xcc\xfb\xbf' : _float80_ma,
+ # double double; low, high order (e.g. PPC 64)
+ b'\x9a\x99\x99\x99\x99\x99Y<\x9a\x99\x99\x99\x99\x99\xb9\xbf' :
+ _float_dd_ma,
+ # double double; high, low order (e.g. PPC 64 le)
+ b'\x9a\x99\x99\x99\x99\x99\xb9\xbf\x9a\x99\x99\x99\x99\x99Y<' :
+ _float_dd_ma,
+ # IEEE 754 128-bit binary float
+ b'\x9a\x99\x99\x99\x99\x99\x99\x99\x99\x99\x99\x99\x99\x99\xfb\xbf' :
+ _float128_ma,
+}
+
+
+def _get_machar(ftype):
+ """ Get MachAr instance or MachAr-like instance
+
+ Get parameters for floating point type, by first trying signatures of
+ various known floating point types, then, if none match, attempting to
+ identify parameters by analysis.
+
+ Parameters
+ ----------
+ ftype : class
+ Numpy floating point type class (e.g. ``np.float64``)
+
+ Returns
+ -------
+ ma_like : instance of :class:`MachAr` or :class:`MachArLike`
+ Object giving floating point parameters for `ftype`.
+
+ Warns
+ -----
+ UserWarning
+ If the binary signature of the float type is not in the dictionary of
+ known float types.
+ """
+ params = _MACHAR_PARAMS.get(ftype)
+ if params is None:
+ raise ValueError(repr(ftype))
+ # Detect known / suspected types
+ key = ftype('-0.1').newbyteorder('<').tobytes()
+ ma_like = _KNOWN_TYPES.get(key)
+ # Could be 80 bit == 10 byte extended precision, where last bytes can be
+ # random garbage. Try comparing first 10 bytes to pattern.
+ if ma_like is None and ftype == ntypes.longdouble:
+ ma_like = _KNOWN_TYPES.get(key[:10])
+ if ma_like is not None:
+ return ma_like
+ # Fall back to parameter discovery
+ warnings.warn(
+ 'Signature {} for {} does not match any known type: '
+ 'falling back to type probe function'.format(key, ftype),
+ UserWarning, stacklevel=2)
+ return _discovered_machar(ftype)
+
+
+def _discovered_machar(ftype):
+ """ Create MachAr instance with found information on float types
+ """
+ params = _MACHAR_PARAMS[ftype]
+ title = 'numpy %s precision floating point number' % params['precname']
+ return MachAr(lambda v: array([v], ftype),
+ lambda v:_frz(v.astype(params['itype']))[0],
+ lambda v:array(_frz(v)[0], ftype),
+ lambda v: params['fmt'] % array(_frz(v)[0], ftype),
+ title)
+
+
class finfo(object):
"""
finfo(dtype)
@@ -128,30 +382,7 @@ class finfo(object):
def _init(self, dtype):
self.dtype = numeric.dtype(dtype)
- if dtype is ntypes.double:
- itype = ntypes.int64
- fmt = '%24.16e'
- precname = 'double'
- elif dtype is ntypes.single:
- itype = ntypes.int32
- fmt = '%15.7e'
- precname = 'single'
- elif dtype is ntypes.longdouble:
- itype = ntypes.longlong
- fmt = '%s'
- precname = 'long double'
- elif dtype is ntypes.half:
- itype = ntypes.int16
- fmt = '%12.5e'
- precname = 'half'
- else:
- raise ValueError(repr(dtype))
-
- machar = MachAr(lambda v:array([v], dtype),
- lambda v:_frz(v.astype(itype))[0],
- lambda v:array(_frz(v)[0], dtype),
- lambda v: fmt % array(_frz(v)[0], dtype),
- 'numpy %s precision floating point number' % precname)
+ machar = _get_machar(dtype)
for word in ['precision', 'iexp',
'maxexp', 'minexp', 'negep',
diff --git a/numpy/core/tests/test_getlimits.py b/numpy/core/tests/test_getlimits.py
index 600f8f52c..adc7aac20 100644
--- a/numpy/core/tests/test_getlimits.py
+++ b/numpy/core/tests/test_getlimits.py
@@ -7,8 +7,10 @@ import numpy as np
from numpy.core import finfo, iinfo
from numpy import half, single, double, longdouble
from numpy.testing import (
- TestCase, run_module_suite, assert_equal
+ TestCase, run_module_suite, assert_equal, assert_
)
+from numpy.core.getlimits import (_discovered_machar, _float16_ma, _float32_ma,
+ _float64_ma, _float128_ma, _float80_ma)
##################################################
@@ -87,5 +89,38 @@ def test_instances():
iinfo(10)
finfo(3.0)
+
+def assert_ma_equal(discovered, ma_like):
+ for key, value in discovered.__dict__.items():
+ assert_equal(value, getattr(ma_like, key))
+
+
+def test_known_types():
+ # Test we are correctly compiling parameters for known types
+ for ftype, ma_like in ((np.float16, _float16_ma),
+ (np.float32, _float32_ma),
+ (np.float64, _float64_ma)):
+ assert_ma_equal(_discovered_machar(ftype), ma_like)
+ # Suppress warning for broken discovery of double double on PPC
+ with np.errstate(all='ignore'):
+ ld_ma = _discovered_machar(np.longdouble)
+ bytes = np.dtype(np.longdouble).itemsize
+ if (ld_ma.it, ld_ma.maxexp) == (63, 16384) and bytes in (12, 16):
+ # 80-bit extended precision
+ assert_ma_equal(ld_ma, _float80_ma)
+ elif (ld_ma.it, ld_ma.maxexp) == (112, 16384) and bytes == 16:
+ # IEE 754 128-bit
+ assert_ma_equal(ld_ma, _float128_ma)
+
+
+def test_plausible_finfo():
+ # Assert that finfo returns reasonable results for all types
+ for ftype in np.sctypes['float'] + np.sctypes['complex']:
+ info = np.finfo(ftype)
+ assert_(info.nmant > 1)
+ assert_(info.minexp < -1)
+ assert_(info.maxexp > 1)
+
+
if __name__ == "__main__":
run_module_suite()