MAINT: Changed diff to use iterative instead of recursive approach

TST: Added tests for `n` parameter
  Added test for `datetime64` type change
  Added tests for axis normalization
  Added test for subtype handling

DOC: Minor updates to docs:
  Added explanation for `n==0`
  Added documentation describing `datetime64` handling
  Updated formatting
  Added call to normalize_axis_index

1 files changed, 31 insertions, 23 deletions

diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index dda3e021b..f5e9ff2a5 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -15,7 +15,7 @@ from numpy.core.numeric import (
     )
 from numpy.core.umath import (
     pi, multiply, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin,
-    mod, exp, log10
+    mod, exp, log10, not_equal, subtract
     )
 from numpy.core.fromnumeric import (
     ravel, nonzero, sort, partition, mean, any, sum
@@ -1827,7 +1827,7 @@ def gradient(f, *varargs, **kwargs):
 
 def diff(a, n=1, axis=-1):
     """
-    Calculate the n-th discrete difference along given axis.
+    Calculate the n-th discrete difference along the given axis.
 
     The first difference is given by ``out[n] = a[n+1] - a[n]`` along
     the given axis, higher differences are calculated by using `diff`
@@ -1838,16 +1838,21 @@ def diff(a, n=1, axis=-1):
     a : array_like
         Input array
     n : int, optional
-        The number of times values are differenced.
+        The number of times values are differenced. If zero, the input
+        is returned as-is.
     axis : int, optional
-        The axis along which the difference is taken, default is the last axis.
+        The axis along which the difference is taken, default is the
+        last axis.
 
     Returns
     -------
     diff : ndarray
         The n-th differences. The shape of the output is the same as `a`
         except along `axis` where the dimension is smaller by `n`. The
-        type of the output is the same as that of the input.
+        type of the output is the same as the type of the difference
+        between any two elements of `a`. This is the same as the type of
+        `a` in most cases. A notable exception is `datetime64`, which
+        results in a `timedelta64` output array.
 
     See Also
     --------
@@ -1855,13 +1860,13 @@ def diff(a, n=1, axis=-1):
 
     Notes
     -----
-    For boolean arrays, the preservation of type means that the result
-    will contain `False` when consecutive elements are the same and
-    `True` when they differ.
+    Type is preserved for boolean arrays, so the result will contain
+    `False` when consecutive elements are the same and `True` when they
+    differ.
 
-    For unsigned integer arrays, the results will also be unsigned. This should
-    not be surprising, as the result is consistent with calculating the
-    difference directly:
+    For unsigned integer arrays, the results will also be unsigned. This
+    should not be surprising, as the result is consistent with
+    calculating the difference directly:
 
     >>> u8_arr = np.array([1, 0], dtype=np.uint8)
     >>> np.diff(u8_arr)
@@ -1869,8 +1874,8 @@ def diff(a, n=1, axis=-1):
     >>> u8_arr[1,...] - u8_arr[0,...]
     array(255, np.uint8)
 
-    If this is not desirable, then the array should be cast to a larger integer
-    type first:
+    If this is not desirable, then the array should be cast to a larger
+    integer type first:
 
     >>> i16_arr = u8_arr.astype(np.int16)
     >>> np.diff(i16_arr)
@@ -1891,30 +1896,33 @@ def diff(a, n=1, axis=-1):
     >>> np.diff(x, axis=0)
     array([[-1,  2,  0, -2]])
 
+    >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
+    >>> np.diff(x)
+    array([1, 1], dtype='timedelta64[D]')
+
     """
     if n == 0:
         return a
     if n < 0:
         raise ValueError(
             "order must be non-negative but got " + repr(n))
+
     a = asanyarray(a)
     nd = a.ndim
-    slice1 = [slice(None)]*nd
-    slice2 = [slice(None)]*nd
+    axis = normalize_axis_index(axis, nd)
+
+    slice1 = [slice(None)] * nd
+    slice2 = [slice(None)] * nd
     slice1[axis] = slice(1, None)
     slice2[axis] = slice(None, -1)
     slice1 = tuple(slice1)
     slice2 = tuple(slice2)
 
-    if a.dtype == np.bool_:
-        da = a[slice1] != a[slice2]
-    else:
-        da = a[slice1] - a[slice2]
+    op = not_equal if a.dtype == np.bool_ else subtract
+    for _ in range(n):
+        a = op(a[slice1], a[slice2])
 
-    if n > 1:
-        return diff(da, n-1, axis=axis)
-    else:
-        return da
+    return a
 
 
 def interp(x, xp, fp, left=None, right=None, period=None):