BUG: don't use pow for integer power ufunc loops.

Fixes gh-7405.
numpy · Apr 1, 2016 · e7ddb39 · e7ddb39
1 parent 2af06c8
commit e7ddb39
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 3 deletions.
diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
@@ -913,9 +913,29 @@ NPY_NO_EXPORT void
 @TYPE@_power(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
     BINARY_LOOP {
-        const @ftype@ in1 = (@ftype@)*(@type@ *)ip1;
-        const @ftype@ in2 = (@ftype@)*(@type@ *)ip2;
-        *((@type@ *)op1) = (@type@) pow(in1, in2);
+        @type@ in1 = *(@type@ *)ip1;
+        @type@ in2 = *(@type@ *)ip2;
+        @type@ out;
+
+        if (in2 < 0 || in1 == 0) {
+            *((@type@ *)op1) = 0;
+            continue;
+        }
+        if (in2 == 0) {
+            *((@type@ *)op1) = 1;
+            continue;
+        }
+
+        out = in2 & 1 ? in1 : 1;
+        in2 >>= 1;
+        while (in2 > 0) {
+            in1 *= in1;
+            if (in2 & 1) {
+                out *= in1;
+            }
+            in2 >>= 1;
+        }
+        *((@type@ *) op1) = out;
     }
 }
 

diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
@@ -425,6 +425,11 @@ def test_fast_power(self):
         res = x ** np.array([[[2]]])
         assert_equal(res.shape, (1, 1, 3))
 
+    def test_integer_power(self):
+        a = np.array([15, 15], 'i8')
+        b = a ** a
+        assert_equal(b, [437893890380859375, 437893890380859375])
+
 
 class TestLog2(TestCase):
     def test_log2_values(self):