Merge pull request #8190 from charris/add-fpower-ufunc

ENH: Add a float_power function with at least float64 precision.

4 files changed, 89 insertions, 1 deletions

diff --git a/doc/release/1.12.0-notes.rst b/doc/release/1.12.0-notes.rst
index 4029abc48..473511c3d 100644
--- a/doc/release/1.12.0-notes.rst
+++ b/doc/release/1.12.0-notes.rst
@@ -261,6 +261,15 @@ values. There are two different propagation modes. The default causes masked
 values to contaminate the result with masks, but the other mode only outputs
 masks if there is no alternative.
 
+New ``float_power`` ufunc
+~~~~~~~~~~~~~~~~~~~~~~~~~
+The new ``float_power`` ufunc is like the ``power`` function except all
+computation is done in a minimum precision of float64. There was a long
+discussion on the numpy mailing list of how to treat integers to negative
+integer powers and a popular proposal was that the ``__pow__`` operator should
+always return results of at least float64 precision. The ``float_power``
+function implements that option. Note that it does not support object arrays.
+
 
 Improvements
 ============
diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py
index 8357fc8da..8c3c86ecd 100644
--- a/numpy/core/code_generators/generate_umath.py
+++ b/numpy/core/code_generators/generate_umath.py
@@ -370,6 +370,12 @@ defdict = {
           TD(inexact, f='pow', astype={'e':'f'}),
           TD(O, f='npy_ObjectPower'),
           ),
+'float_power':
+    Ufunc(2, 1, None,
+          docstrings.get('numpy.core.umath.float_power'),
+          None,
+          TD('dgDG', f='pow'),
+          ),
 'absolute':
     Ufunc(1, 1, None,
           docstrings.get('numpy.core.umath.absolute'),
diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py
index 0d506102f..dd4cf1ea8 100644
--- a/numpy/core/code_generators/ufunc_docstrings.py
+++ b/numpy/core/code_generators/ufunc_docstrings.py
@@ -2543,7 +2543,8 @@ add_newdoc('numpy.core.umath', 'power',
     First array elements raised to powers from second array, element-wise.
 
     Raise each base in `x1` to the positionally-corresponding power in
-    `x2`.  `x1` and `x2` must be broadcastable to the same shape.
+    `x2`.  `x1` and `x2` must be broadcastable to the same shape. Note that an
+    integer type raised to a negative integer power will raise a ValueError.
 
     Parameters
     ----------
@@ -2557,6 +2558,10 @@ add_newdoc('numpy.core.umath', 'power',
     y : ndarray
         The bases in `x1` raised to the exponents in `x2`.
 
+    See Also
+    --------
+    float_power : power function that promotes integers to float
+
     Examples
     --------
     Cube each element in a list.
@@ -2585,6 +2590,63 @@ add_newdoc('numpy.core.umath', 'power',
 
     """)
 
+add_newdoc('numpy.core.umath', 'float_power',
+    """
+    First array elements raised to powers from second array, element-wise.
+
+    Raise each base in `x1` to the positionally-corresponding power in `x2`.
+    `x1` and `x2` must be broadcastable to the same shape. This differs from
+    the power function in that integers, float16, and float32  are promoted to
+    floats with a minimum precision of float64 so that the result is always
+    inexact.  The intent is that the function will return a usable result for
+    negative powers and seldom overflow for positive powers.
+
+    .. versionadded:: 1.12.0
+
+    Parameters
+    ----------
+    x1 : array_like
+        The bases.
+    x2 : array_like
+        The exponents.
+
+    Returns
+    -------
+    y : ndarray
+        The bases in `x1` raised to the exponents in `x2`.
+
+    See Also
+    --------
+    power : power function that preserves type
+
+    Examples
+    --------
+    Cube each element in a list.
+
+    >>> x1 = range(6)
+    >>> x1
+    [0, 1, 2, 3, 4, 5]
+    >>> np.float_power(x1, 3)
+    array([   0.,    1.,    8.,   27.,   64.,  125.])
+
+    Raise the bases to different exponents.
+
+    >>> x2 = [1.0, 2.0, 3.0, 3.0, 2.0, 1.0]
+    >>> np.float_power(x1, x2)
+    array([  0.,   1.,   8.,  27.,  16.,   5.])
+
+    The effect of broadcasting.
+
+    >>> x2 = np.array([[1, 2, 3, 3, 2, 1], [1, 2, 3, 3, 2, 1]])
+    >>> x2
+    array([[1, 2, 3, 3, 2, 1],
+           [1, 2, 3, 3, 2, 1]])
+    >>> np.float_power(x1, x2)
+    array([[  0.,   1.,   8.,  27.,  16.,   5.],
+           [  0.,   1.,   8.,  27.,  16.,   5.]])
+
+    """)
+
 add_newdoc('numpy.core.umath', 'radians',
     """
     Convert angles from degrees to radians.
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 4538c53ac..4c0243559 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -465,6 +465,17 @@ class TestPower(TestCase):
             assert_raises(ValueError, np.power, one, minusone)
 
 
+class TestFloat_power(TestCase):
+    def test_type_conversion(self):
+        arg_type = '?bhilBHILefdgFDG'
+        res_type = 'ddddddddddddgDDG'
+        for dtin, dtout in zip(arg_type, res_type):
+            msg = "dtin: %s, dtout: %s" % (dtin, dtout)
+            arg = np.ones(1, dtype=dtin)
+            res = np.float_power(arg, arg)
+            assert_(res.dtype.name == np.dtype(dtout).name, msg)
+
+
 class TestLog2(TestCase):
     def test_log2_values(self):
         x = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]