Merge branch 'uarray' into uarray-me

74 files changed, 1085 insertions, 710 deletions

diff --git a/numpy/__init__.py b/numpy/__init__.py
index ba88c733f..ae297597e 100644
--- a/numpy/__init__.py
+++ b/numpy/__init__.py
@@ -166,6 +166,8 @@ else:
     # now that numpy modules are imported, can initialize limits
     core.getlimits._register_known_types()
 
+    __all__.extend(['bool', 'int', 'float', 'complex', 'object', 'unicode',
+                    'str'])
     __all__.extend(['__version__', 'show_config'])
     __all__.extend(core.__all__)
     __all__.extend(_mat.__all__)
@@ -182,9 +184,35 @@ else:
     oldnumeric = 'removed'
     numarray = 'removed'
 
-    # We don't actually use this ourselves anymore, but I'm not 100% sure that
-    # no-one else in the world is using it (though I hope not)
-    from .testing import Tester
+    if sys.version_info[:2] >= (3, 7):
+        # Importing Tester requires importing all of UnitTest which is not a
+        # cheap import Since it is mainly used in test suits, we lazy import it
+        # here to save on the order of 10 ms of import time for most users
+        #
+        # The previous way Tester was imported also had a side effect of adding
+        # the full `numpy.testing` namespace
+        #
+        # module level getattr is only supported in 3.7 onwards
+        # https://www.python.org/dev/peps/pep-0562/
+        def __getattr__(attr):
+            if attr == 'testing':
+                import numpy.testing as testing
+                return testing
+            elif attr == 'Tester':
+                from .testing import Tester
+                return Tester
+            else:
+                raise AttributeError(
+                    "module %s has no attribute $s".format(__name__, attr))
+
+
+        def __dir__():
+            return __all__ + ['Tester', 'testing']
+
+    else:
+        # We don't actually use this ourselves anymore, but I'm not 100% sure that
+        # no-one else in the world is using it (though I hope not)
+        from .testing import Tester
 
     # Pytest testing
     from numpy._pytesttester import PytestTester
diff --git a/numpy/core/_dtype.py b/numpy/core/_dtype.py
index 3a12c8fad..092b848dc 100644
--- a/numpy/core/_dtype.py
+++ b/numpy/core/_dtype.py
@@ -252,7 +252,7 @@ def _is_packed(dtype):
     from a list of the field names and dtypes with no additional
     dtype parameters.
 
-    Duplicates the C `is_dtype_struct_simple_unaligned_layout` functio.
+    Duplicates the C `is_dtype_struct_simple_unaligned_layout` function.
     """
     total_offset = 0
     for name in dtype.names:
diff --git a/numpy/core/_internal.py b/numpy/core/_internal.py
index c70718cb6..b0ea603e1 100644
--- a/numpy/core/_internal.py
+++ b/numpy/core/_internal.py
@@ -459,7 +459,7 @@ def _getfield_is_safe(oldtype, newtype, offset):
     if newtype.hasobject or oldtype.hasobject:
         if offset == 0 and newtype == oldtype:
             return
-        if oldtype.names:
+        if oldtype.names is not None:
             for name in oldtype.names:
                 if (oldtype.fields[name][1] == offset and
                         oldtype.fields[name][0] == newtype):
diff --git a/numpy/core/arrayprint.py b/numpy/core/arrayprint.py
index ecd05d3ac..5761c4875 100644
--- a/numpy/core/arrayprint.py
+++ b/numpy/core/arrayprint.py
@@ -685,7 +685,7 @@ def array2string(a, max_line_width=None, precision=None,
         if style is np._NoValue:
             style = repr
 
-        if a.shape == () and not a.dtype.names:
+        if a.shape == () and a.dtype.names is None:
             return style(a.item())
     elif style is not np._NoValue:
         # Deprecation 11-9-2017  v1.14
@@ -984,20 +984,6 @@ class FloatingFormat(object):
                                       pad_left=self.pad_left,
                                       pad_right=self.pad_right)
 
-# for back-compatibility, we keep the classes for each float type too
-class FloatFormat(FloatingFormat):
-    def __init__(self, *args, **kwargs):
-        warnings.warn("FloatFormat has been replaced by FloatingFormat",
-                      DeprecationWarning, stacklevel=2)
-        super(FloatFormat, self).__init__(*args, **kwargs)
-
-
-class LongFloatFormat(FloatingFormat):
-    def __init__(self, *args, **kwargs):
-        warnings.warn("LongFloatFormat has been replaced by FloatingFormat",
-                      DeprecationWarning, stacklevel=2)
-        super(LongFloatFormat, self).__init__(*args, **kwargs)
-
 
 @set_module('numpy')
 def format_float_scientific(x, precision=None, unique=True, trim='k',
@@ -1196,21 +1182,6 @@ class ComplexFloatingFormat(object):
 
         return r + i
 
-# for back-compatibility, we keep the classes for each complex type too
-class ComplexFormat(ComplexFloatingFormat):
-    def __init__(self, *args, **kwargs):
-        warnings.warn(
-            "ComplexFormat has been replaced by ComplexFloatingFormat",
-            DeprecationWarning, stacklevel=2)
-        super(ComplexFormat, self).__init__(*args, **kwargs)
-
-class LongComplexFormat(ComplexFloatingFormat):
-    def __init__(self, *args, **kwargs):
-        warnings.warn(
-            "LongComplexFormat has been replaced by ComplexFloatingFormat",
-            DeprecationWarning, stacklevel=2)
-        super(LongComplexFormat, self).__init__(*args, **kwargs)
-
 
 class _TimelikeFormat(object):
     def __init__(self, data):
@@ -1321,16 +1292,6 @@ class StructuredVoidFormat(object):
             return "({})".format(", ".join(str_fields))
 
 
-# for backwards compatibility
-class StructureFormat(StructuredVoidFormat):
-    def __init__(self, *args, **kwargs):
-        # NumPy 1.14, 2018-02-14
-        warnings.warn(
-            "StructureFormat has been replaced by StructuredVoidFormat",
-            DeprecationWarning, stacklevel=2)
-        super(StructureFormat, self).__init__(*args, **kwargs)
-
-
 def _void_scalar_repr(x):
     """
     Implements the repr for structured-void scalars. It is called from the
diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py
index bf1747272..ae871ea6f 100644
--- a/numpy/core/code_generators/generate_umath.py
+++ b/numpy/core/code_generators/generate_umath.py
@@ -662,6 +662,8 @@ defdict = {
     Ufunc(1, 1, None,
           docstrings.get('numpy.core.umath.cos'),
           None,
+          TD('e', f='cos', astype={'e':'f'}),
+          TD('f', simd=[('fma', 'f'), ('avx512f', 'f')]),
           TD(inexact, f='cos', astype={'e':'f'}),
           TD(P, f='cos'),
           ),
@@ -669,6 +671,8 @@ defdict = {
     Ufunc(1, 1, None,
           docstrings.get('numpy.core.umath.sin'),
           None,
+          TD('e', f='sin', astype={'e':'f'}),
+          TD('f', simd=[('fma', 'f'), ('avx512f', 'f')]),
           TD(inexact, f='sin', astype={'e':'f'}),
           TD(P, f='sin'),
           ),
@@ -705,7 +709,7 @@ defdict = {
           docstrings.get('numpy.core.umath.exp'),
           None,
           TD('e', f='exp', astype={'e':'f'}),
-          TD('f', simd=[('avx2', 'f'), ('avx512f', 'f')]),
+          TD('f', simd=[('fma', 'f'), ('avx512f', 'f')]),
           TD(inexact, f='exp', astype={'e':'f'}),
           TD(P, f='exp'),
           ),
@@ -728,7 +732,7 @@ defdict = {
           docstrings.get('numpy.core.umath.log'),
           None,
           TD('e', f='log', astype={'e':'f'}),
-          TD('f', simd=[('avx2', 'f'), ('avx512f', 'f')]),
+          TD('f', simd=[('fma', 'f'), ('avx512f', 'f')]),
           TD(inexact, f='log', astype={'e':'f'}),
           TD(P, f='log'),
           ),
diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py
index 3389e7d66..3314e516e 100644
--- a/numpy/core/fromnumeric.py
+++ b/numpy/core/fromnumeric.py
@@ -25,7 +25,7 @@ __all__ = [
     'argmin', 'argpartition', 'argsort', 'around', 'choose', 'clip',
     'compress', 'cumprod', 'cumproduct', 'cumsum', 'diagonal', 'mean',
     'ndim', 'nonzero', 'partition', 'prod', 'product', 'ptp', 'put',
-    'rank', 'ravel', 'repeat', 'reshape', 'resize', 'round_',
+    'ravel', 'repeat', 'reshape', 'resize', 'round_',
     'searchsorted', 'shape', 'size', 'sometrue', 'sort', 'squeeze',
     'std', 'sum', 'swapaxes', 'take', 'trace', 'transpose', 'var',
 ]
@@ -380,6 +380,7 @@ def choose(a, choices, out=None, mode='raise'):
     See Also
     --------
     ndarray.choose : equivalent method
+    numpy.take_along_axis : Preferable if `choices` is an array
 
     Notes
     -----
@@ -908,17 +909,17 @@ def sort(a, axis=-1, kind=None, order=None):
 
     .. versionadded:: 1.12.0
 
-    quicksort has been changed to `introsort <https://en.wikipedia.org/wiki/Introsort>`_. 
+    quicksort has been changed to `introsort <https://en.wikipedia.org/wiki/Introsort>`_.
     When sorting does not make enough progress it switches to
-    `heapsort <https://en.wikipedia.org/wiki/Heapsort>`_. 
+    `heapsort <https://en.wikipedia.org/wiki/Heapsort>`_.
     This implementation makes quicksort O(n*log(n)) in the worst case.
 
     'stable' automatically chooses the best stable sorting algorithm
-    for the data type being sorted. 
-    It, along with 'mergesort' is currently mapped to 
-    `timsort <https://en.wikipedia.org/wiki/Timsort>`_ 
-    or `radix sort <https://en.wikipedia.org/wiki/Radix_sort>`_ 
-    depending on the data type. 
+    for the data type being sorted.
+    It, along with 'mergesort' is currently mapped to
+    `timsort <https://en.wikipedia.org/wiki/Timsort>`_
+    or `radix sort <https://en.wikipedia.org/wiki/Radix_sort>`_
+    depending on the data type.
     API forward compatibility currently limits the
     ability to select the implementation and it is hardwired for the different
     data types.
@@ -2782,6 +2783,10 @@ def alen(a):
     7
 
     """
+    # NumPy 1.18.0, 2019-08-02
+    warnings.warn(
+        "`np.alen` is deprecated, use `len` instead",
+        DeprecationWarning, stacklevel=2)
     try:
         return len(a)
     except TypeError:
@@ -3414,7 +3419,7 @@ def var(a, axis=None, dtype=None, out=None, ddof=0, keepdims=np._NoValue):
         instead of a single axis or all the axes as before.
     dtype : data-type, optional
         Type to use in computing the variance.  For arrays of integer type
-        the default is `float32`; for arrays of float types it is the same as
+        the default is `float64`; for arrays of float types it is the same as
         the array type.
     out : ndarray, optional
         Alternate output array in which to place the result.  It must have
@@ -3573,30 +3578,3 @@ def alltrue(*args, **kwargs):
     numpy.all : Equivalent function; see for details.
     """
     return all(*args, **kwargs)
-
-
-@array_function_dispatch(_ndim_dispatcher)
-def rank(a):
-    """
-    Return the number of dimensions of an array.
-
-    .. note::
-        This function is deprecated in NumPy 1.9 to avoid confusion with
-        `numpy.linalg.matrix_rank`. The ``ndim`` attribute or function
-        should be used instead.
-
-    See Also
-    --------
-    ndim : equivalent non-deprecated function
-
-    Notes
-    -----
-    In the old Numeric package, `rank` was the term used for the number of
-    dimensions, but in NumPy `ndim` is used instead.
-    """
-    # 2014-04-12, 1.9
-    warnings.warn(
-        "`rank` is deprecated; use the `ndim` attribute or function instead. "
-        "To find the rank of a matrix see `numpy.linalg.matrix_rank`.",
-        VisibleDeprecationWarning, stacklevel=3)
-    return ndim(a)
diff --git a/numpy/core/include/numpy/npy_common.h b/numpy/core/include/numpy/npy_common.h
index 108c0a202..27b83f7b5 100644
--- a/numpy/core/include/numpy/npy_common.h
+++ b/numpy/core/include/numpy/npy_common.h
@@ -44,10 +44,14 @@
 #else
 #define NPY_GCC_TARGET_AVX
 #endif
+
+#if defined HAVE_ATTRIBUTE_TARGET_AVX2_WITH_INTRINSICS
+#define HAVE_ATTRIBUTE_TARGET_FMA
+#define NPY_GCC_TARGET_FMA __attribute__((target("avx2,fma")))
+#endif
+
 #if defined HAVE_ATTRIBUTE_TARGET_AVX2 && defined HAVE_LINK_AVX2
 #define NPY_GCC_TARGET_AVX2 __attribute__((target("avx2")))
-#elif defined HAVE_ATTRIBUTE_TARGET_AVX2_WITH_INTRINSICS
-#define NPY_GCC_TARGET_AVX2 __attribute__((target("avx2")))
 #else
 #define NPY_GCC_TARGET_AVX2
 #endif
diff --git a/numpy/core/include/numpy/npy_math.h b/numpy/core/include/numpy/npy_math.h
index 6a78ff3c2..126b861bf 100644
--- a/numpy/core/include/numpy/npy_math.h
+++ b/numpy/core/include/numpy/npy_math.h
@@ -144,7 +144,22 @@ NPY_INLINE static float __npy_nzerof(void)
 #define NPY_COEFF_Q3_LOGf 9.864942958519418960339e-01f
 #define NPY_COEFF_Q4_LOGf 1.546476374983906719538e-01f
 #define NPY_COEFF_Q5_LOGf 5.875095403124574342950e-03f
-
+/*
+ * Constants used in vector implementation of sinf/cosf(x)
+ */
+#define NPY_TWO_O_PIf 0x1.45f306p-1f
+#define NPY_CODY_WAITE_PI_O_2_HIGHf -0x1.921fb0p+00f
+#define NPY_CODY_WAITE_PI_O_2_MEDf -0x1.5110b4p-22f
+#define NPY_CODY_WAITE_PI_O_2_LOWf -0x1.846988p-48f
+#define NPY_COEFF_INVF0_COSINEf 0x1.000000p+00f
+#define NPY_COEFF_INVF2_COSINEf -0x1.000000p-01f
+#define NPY_COEFF_INVF4_COSINEf 0x1.55553cp-05f
+#define NPY_COEFF_INVF6_COSINEf -0x1.6c06dcp-10f
+#define NPY_COEFF_INVF8_COSINEf 0x1.98e616p-16f
+#define NPY_COEFF_INVF3_SINEf -0x1.555556p-03f
+#define NPY_COEFF_INVF5_SINEf 0x1.11119ap-07f
+#define NPY_COEFF_INVF7_SINEf -0x1.a06bbap-13f
+#define NPY_COEFF_INVF9_SINEf 0x1.7d3bbcp-19f
 /*
  * Integer functions.
  */
@@ -163,6 +178,15 @@ NPY_INPLACE npy_longlong npy_gcdll(npy_longlong a, npy_longlong b);
 NPY_INPLACE npy_longlong npy_lcmll(npy_longlong a, npy_longlong b);
 
 /*
+ * avx function has a common API for both sin & cos. This enum is used to
+ * distinguish between the two
+ */
+typedef enum {
+    npy_compute_sin,
+    npy_compute_cos
+} NPY_TRIG_OP;
+
+/*
  * C99 double math funcs
  */
 NPY_INPLACE double npy_sin(double x);
diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index ea2ef900e..8ada87b9f 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -48,14 +48,6 @@ array_function_dispatch = functools.partial(
     overrides.array_function_dispatch, module='numpy')
 
 
-def loads(*args, **kwargs):
-    # NumPy 1.15.0, 2017-12-10
-    warnings.warn(
-        "np.core.numeric.loads is deprecated, use pickle.loads instead",
-        DeprecationWarning, stacklevel=2)
-    return pickle.loads(*args, **kwargs)
-
-
 __all__ = [
     'newaxis', 'ndarray', 'flatiter', 'nditer', 'nested_iters', 'ufunc',
     'arange', 'array', 'zeros', 'count_nonzero', 'empty', 'broadcast', 'dtype',
@@ -66,7 +58,7 @@ __all__ = [
     'correlate', 'convolve', 'inner', 'dot', 'outer', 'vdot', 'roll',
     'rollaxis', 'moveaxis', 'cross', 'tensordot', 'little_endian',
     'fromiter', 'array_equal', 'array_equiv', 'indices', 'fromfunction',
-    'isclose', 'load', 'loads', 'isscalar', 'binary_repr', 'base_repr', 'ones',
+    'isclose', 'isscalar', 'binary_repr', 'base_repr', 'ones',
     'identity', 'allclose', 'compare_chararrays', 'putmask',
     'flatnonzero', 'Inf', 'inf', 'infty', 'Infinity', 'nan', 'NaN',
     'False_', 'True_', 'bitwise_not', 'CLIP', 'RAISE', 'WRAP', 'MAXDIMS',
@@ -1935,6 +1927,10 @@ def binary_repr(num, width=None):
                 "will raise an error in the future.", DeprecationWarning,
                 stacklevel=3)
 
+    # Ensure that num is a Python integer to avoid overflow or unwanted
+    # casts to floating point.
+    num = operator.index(num)
+
     if num == 0:
         return '0' * (width or 1)
 
@@ -2024,30 +2020,6 @@ def base_repr(number, base=2, padding=0):
     return ''.join(reversed(res or '0'))
 
 
-def load(file):
-    """
-    Wrapper around cPickle.load which accepts either a file-like object or
-    a filename.
-
-    Note that the NumPy binary format is not based on pickle/cPickle anymore.
-    For details on the preferred way of loading and saving files, see `load`
-    and `save`.
-
-    See Also
-    --------
-    load, save
-
-    """
-    # NumPy 1.15.0, 2017-12-10
-    warnings.warn(
-        "np.core.numeric.load is deprecated, use pickle.load instead",
-        DeprecationWarning, stacklevel=2)
-    if isinstance(file, type("")):
-        with open(file, "rb") as file_pointer:
-            return pickle.load(file_pointer)
-    return pickle.load(file)
-
-
 # These are all essentially abbreviations
 # These might wind up in a special abbreviations module
 
@@ -2124,7 +2096,7 @@ def allclose(a, b, rtol=1.e-5, atol=1.e-8, equal_nan=False):
         The absolute tolerance parameter (see Notes).
     equal_nan : bool
         Whether to compare NaN's as equal.  If True, NaN's in `a` will be
-        considered equal to NaN's in `b` in the output array.
+        considered equal to NaN's in `b`.
 
         .. versionadded:: 1.10.0
 
diff --git a/numpy/core/records.py b/numpy/core/records.py
index 0576005e7..a1439f9df 100644
--- a/numpy/core/records.py
+++ b/numpy/core/records.py
@@ -268,8 +268,8 @@ class record(nt.void):
             except AttributeError:
                 #happens if field is Object type
                 return obj
-            if dt.fields:
-                return obj.view((self.__class__, obj.dtype.fields))
+            if dt.names is not None:
+                return obj.view((self.__class__, obj.dtype))
             return obj
         else:
             raise AttributeError("'record' object has no "
@@ -293,8 +293,8 @@ class record(nt.void):
         obj = nt.void.__getitem__(self, indx)
 
         # copy behavior of record.__getattribute__,
-        if isinstance(obj, nt.void) and obj.dtype.fields:
-            return obj.view((self.__class__, obj.dtype.fields))
+        if isinstance(obj, nt.void) and obj.dtype.names is not None:
+            return obj.view((self.__class__, obj.dtype))
         else:
             # return a single element
             return obj
@@ -444,7 +444,7 @@ class recarray(ndarray):
         return self
 
     def __array_finalize__(self, obj):
-        if self.dtype.type is not record and self.dtype.fields:
+        if self.dtype.type is not record and self.dtype.names is not None:
             # if self.dtype is not np.record, invoke __setattr__ which will
             # convert it to a record if it is a void dtype.
             self.dtype = self.dtype
@@ -472,7 +472,7 @@ class recarray(ndarray):
         # with void type convert it to the same dtype.type (eg to preserve
         # numpy.record type if present), since nested structured fields do not
         # inherit type. Don't do this for non-void structures though.
-        if obj.dtype.fields:
+        if obj.dtype.names is not None:
             if issubclass(obj.dtype.type, nt.void):
                 return obj.view(dtype=(self.dtype.type, obj.dtype))
             return obj
@@ -487,7 +487,7 @@ class recarray(ndarray):
 
         # Automatically convert (void) structured types to records
         # (but not non-void structures, subarrays, or non-structured voids)
-        if attr == 'dtype' and issubclass(val.type, nt.void) and val.fields:
+        if attr == 'dtype' and issubclass(val.type, nt.void) and val.names is not None:
             val = sb.dtype((record, val))
 
         newattr = attr not in self.__dict__
@@ -521,7 +521,7 @@ class recarray(ndarray):
         # copy behavior of getattr, except that here
         # we might also be returning a single element
         if isinstance(obj, ndarray):
-            if obj.dtype.fields:
+            if obj.dtype.names is not None:
                 obj = obj.view(type(self))
                 if issubclass(obj.dtype.type, nt.void):
                     return obj.view(dtype=(self.dtype.type, obj.dtype))
@@ -577,7 +577,7 @@ class recarray(ndarray):
 
         if val is None:
             obj = self.getfield(*res)
-            if obj.dtype.fields:
+            if obj.dtype.names is not None:
                 return obj
             return obj.view(ndarray)
         else:
diff --git a/numpy/core/setup_common.py b/numpy/core/setup_common.py
index 6e3109ab5..a3f7acd6d 100644
--- a/numpy/core/setup_common.py
+++ b/numpy/core/setup_common.py
@@ -178,9 +178,10 @@ OPTIONAL_FUNCTION_ATTRIBUTES = [('__attribute__((optimize("unroll-loops")))',
 # gcc 4.8.4 support attributes but not with intrisics
 # tested via "#include<%s> int %s %s(void *){code; return 0;};" % (header, attribute, name, code)
 # function name will be converted to HAVE_<upper-case-name> preprocessor macro
-OPTIONAL_FUNCTION_ATTRIBUTES_WITH_INTRINSICS = [('__attribute__((target("avx2")))',
+OPTIONAL_FUNCTION_ATTRIBUTES_WITH_INTRINSICS = [('__attribute__((target("avx2,fma")))',
                                 'attribute_target_avx2_with_intrinsics',
-                                '__m256 temp = _mm256_set1_ps(1.0)',
+                                '__m256 temp = _mm256_set1_ps(1.0); temp = \
+                                _mm256_fmadd_ps(temp, temp, temp)',
                                 'immintrin.h'),
                                 ('__attribute__((target("avx512f")))',
                                 'attribute_target_avx512f_with_intrinsics',
diff --git a/numpy/core/src/multiarray/alloc.c b/numpy/core/src/multiarray/alloc.c
index addc9f006..a7f34cbe5 100644
--- a/numpy/core/src/multiarray/alloc.c
+++ b/numpy/core/src/multiarray/alloc.c
@@ -25,10 +25,14 @@
 
 #include <assert.h>
 
-#ifdef HAVE_SYS_MMAN_H
+#ifdef NPY_OS_LINUX
 #include <sys/mman.h>
-#if defined MADV_HUGEPAGE && defined HAVE_MADVISE
-#define HAVE_MADV_HUGEPAGE
+#ifndef MADV_HUGEPAGE
+/*
+ * Use code 14 (MADV_HUGEPAGE) if it isn't defined. This gives a chance of
+ * enabling huge pages even if built with linux kernel < 2.6.38
+ */
+#define MADV_HUGEPAGE 14
 #endif
 #endif
 
@@ -74,11 +78,15 @@ _npy_alloc_cache(npy_uintp nelem, npy_uintp esz, npy_uint msz,
 #ifdef _PyPyGC_AddMemoryPressure
         _PyPyPyGC_AddMemoryPressure(nelem * esz);
 #endif
-#ifdef HAVE_MADV_HUGEPAGE
+#ifdef NPY_OS_LINUX
         /* allow kernel allocating huge pages for large arrays */
         if (NPY_UNLIKELY(nelem * esz >= ((1u<<22u)))) {
             npy_uintp offset = 4096u - (npy_uintp)p % (4096u);
             npy_uintp length = nelem * esz - offset;
+            /**
+             * Intentionally not checking for errors that may be returned by
+             * older kernel versions; optimistically tries enabling huge pages.
+             */
             madvise((void*)((npy_uintp)p + offset), length, MADV_HUGEPAGE);
         }
 #endif
diff --git a/numpy/core/src/multiarray/arrayobject.c b/numpy/core/src/multiarray/arrayobject.c
index bbb736fd0..eb939f47c 100644
--- a/numpy/core/src/multiarray/arrayobject.c
+++ b/numpy/core/src/multiarray/arrayobject.c
@@ -483,10 +483,11 @@ array_dealloc(PyArrayObject *self)
             char const * msg = "WRITEBACKIFCOPY detected in array_dealloc. "
                 " Required call to PyArray_ResolveWritebackIfCopy or "
                 "PyArray_DiscardWritebackIfCopy is missing.";
-            Py_INCREF(self); /* hold on to self in next call  since if
-                              * refcount == 0 it will recurse back into
-                              *array_dealloc
-                              */
+            /*
+             * prevent reaching 0 twice and thus recursing into dealloc.
+             * Increasing sys.gettotalrefcount, but path should not be taken.
+             */
+            Py_INCREF(self);
             WARN_IN_DEALLOC(PyExc_RuntimeWarning, msg);
             retval = PyArray_ResolveWritebackIfCopy(self);
             if (retval < 0)
@@ -500,10 +501,11 @@ array_dealloc(PyArrayObject *self)
             char const * msg = "UPDATEIFCOPY detected in array_dealloc. "
                 " Required call to PyArray_ResolveWritebackIfCopy or "
                 "PyArray_DiscardWritebackIfCopy is missing";
-            Py_INCREF(self); /* hold on to self in next call  since if
-                              * refcount == 0 it will recurse back into
-                              *array_dealloc
-                              */
+            /*
+             * prevent reaching 0 twice and thus recursing into dealloc.
+             * Increasing sys.gettotalrefcount, but path should not be taken.
+             */
+            Py_INCREF(self);
             /* 2017-Nov-10 1.14 */
             WARN_IN_DEALLOC(PyExc_DeprecationWarning, msg);
             retval = PyArray_ResolveWritebackIfCopy(self);
@@ -523,12 +525,7 @@ array_dealloc(PyArrayObject *self)
     if ((fa->flags & NPY_ARRAY_OWNDATA) && fa->data) {
         /* Free internal references if an Object array */
         if (PyDataType_FLAGCHK(fa->descr, NPY_ITEM_REFCOUNT)) {
-            Py_INCREF(self); /*hold on to self */
             PyArray_XDECREF(self);
-            /*
-             * Don't need to DECREF -- because we are deleting
-             * self already...
-             */
         }
         npy_free_cache(fa->data, PyArray_NBYTES(self));
     }
diff --git a/numpy/core/src/multiarray/compiled_base.c b/numpy/core/src/multiarray/compiled_base.c
index dc79bfa09..c38067681 100644
--- a/numpy/core/src/multiarray/compiled_base.c
+++ b/numpy/core/src/multiarray/compiled_base.c
@@ -367,6 +367,18 @@ arr_insert(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
 
 #define LIKELY_IN_CACHE_SIZE 8
 
+#ifdef __INTEL_COMPILER
+#pragma intel optimization_level 0
+#endif
+static NPY_INLINE npy_intp
+_linear_search(const npy_double key, const npy_double *arr, const npy_intp len, const npy_intp i0)
+{
+    npy_intp i;
+
+    for (i = i0; i < len && key >= arr[i]; i++);
+    return i - 1;
+}
+
 /** @brief find index of a sorted array such that arr[i] <= key < arr[i + 1].
  *
  * If an starting index guess is in-range, the array values around this
@@ -406,10 +418,7 @@ binary_search_with_guess(const npy_double key, const npy_double *arr,
      * From above we know key >= arr[0] when we start.
      */
     if (len <= 4) {
-        npy_intp i;
-
-        for (i = 1; i < len && key >= arr[i]; ++i);
-        return i - 1;
+        return _linear_search(key, arr, len, 1);
     }
 
     if (guess > len - 3) {
diff --git a/numpy/core/src/multiarray/ctors.c b/numpy/core/src/multiarray/ctors.c
index 53efb1cea..f0a15505d 100644
--- a/numpy/core/src/multiarray/ctors.c
+++ b/numpy/core/src/multiarray/ctors.c
@@ -1849,13 +1849,6 @@ PyArray_GetArrayParamsFromObject(PyObject *op,
             *out_arr = NULL;
             return 0;
         }
-        if (is_object && (requested_dtype != NULL) && 
-                (requested_dtype->type_num != NPY_OBJECT)) {
-            PyErr_SetString(PyExc_ValueError,
-               "cannot create an array from unequal-length (ragged) sequences");
-            Py_DECREF(*out_dtype);
-            return -1;
-        }
         /* If object arrays are forced */
         if (is_object) {
             Py_DECREF(*out_dtype);
@@ -2772,61 +2765,30 @@ PyArray_DescrFromObject(PyObject *op, PyArray_Descr *mintype)
 /* They all zero-out the memory as previously done */
 
 /* steals reference to descr -- and enforces native byteorder on it.*/
+
 /*NUMPY_API
-  Like FromDimsAndData but uses the Descr structure instead of typecode
-  as input.
+  Deprecated, use PyArray_NewFromDescr instead.
 */
 NPY_NO_EXPORT PyObject *
 PyArray_FromDimsAndDataAndDescr(int nd, int *d,
                                 PyArray_Descr *descr,
                                 char *data)
 {
-    PyObject *ret;
-    int i;
-    npy_intp newd[NPY_MAXDIMS];
-    char msg[] = "PyArray_FromDimsAndDataAndDescr: use PyArray_NewFromDescr.";
-
-    if (DEPRECATE(msg) < 0) {
-        /* 2009-04-30, 1.5 */
-        return NULL;
-    }
-    if (!PyArray_ISNBO(descr->byteorder))
-        descr->byteorder = '=';
-    for (i = 0; i < nd; i++) {
-        newd[i] = (npy_intp) d[i];
-    }
-    ret = PyArray_NewFromDescr(&PyArray_Type, descr,
-                               nd, newd,
-                               NULL, data,
-                               (data ? NPY_ARRAY_CARRAY : 0), NULL);
-    return ret;
+    PyErr_SetString(PyExc_NotImplementedError,
+                "PyArray_FromDimsAndDataAndDescr: use PyArray_NewFromDescr.");
+    Py_DECREF(descr);
+    return NULL;
 }
 
 /*NUMPY_API
-  Construct an empty array from dimensions and typenum
+  Deprecated, use PyArray_SimpleNew instead.
 */
 NPY_NO_EXPORT PyObject *
 PyArray_FromDims(int nd, int *d, int type)
 {
-    PyArrayObject *ret;
-    char msg[] = "PyArray_FromDims: use PyArray_SimpleNew.";
-
-    if (DEPRECATE(msg) < 0) {
-        /* 2009-04-30, 1.5 */
-        return NULL;
-    }
-    ret = (PyArrayObject *)PyArray_FromDimsAndDataAndDescr(nd, d,
-                                          PyArray_DescrFromType(type),
-                                          NULL);
-    /*
-     * Old FromDims set memory to zero --- some algorithms
-     * relied on that.  Better keep it the same. If
-     * Object type, then it's already been set to zero, though.
-     */
-    if (ret && (PyArray_DESCR(ret)->type_num != NPY_OBJECT)) {
-        memset(PyArray_DATA(ret), 0, PyArray_NBYTES(ret));
-    }
-    return (PyObject *)ret;
+    PyErr_SetString(PyExc_NotImplementedError,
+                "PyArray_FromDims: use PyArray_SimpleNew.");
+    return NULL;
 }
 
 /* end old calls */
@@ -4000,7 +3962,7 @@ PyArray_FromString(char *data, npy_intp slen, PyArray_Descr *dtype,
         size_t nread = 0;
         char *end;
 
-        if (dtype->f->scanfunc == NULL) {
+        if (dtype->f->fromstr == NULL) {
             PyErr_SetString(PyExc_ValueError,
                             "don't know how to read "       \
                             "character strings with that "  \
diff --git a/numpy/core/src/multiarray/descriptor.c b/numpy/core/src/multiarray/descriptor.c
index 4d22c9ee7..c7db092e6 100644
--- a/numpy/core/src/multiarray/descriptor.c
+++ b/numpy/core/src/multiarray/descriptor.c
@@ -102,6 +102,7 @@ _arraydescr_from_dtype_attr(PyObject *obj, PyArray_Descr **newdescr)
     if (Py_EnterRecursiveCall(
             " while trying to convert the given data type from its "
             "`.dtype` attribute.") != 0) {
+        Py_DECREF(dtypedescr);
         return 1;
     }
 
diff --git a/numpy/core/src/multiarray/item_selection.c b/numpy/core/src/multiarray/item_selection.c
index 762563eb5..19de04a93 100644
--- a/numpy/core/src/multiarray/item_selection.c
+++ b/numpy/core/src/multiarray/item_selection.c
@@ -1456,8 +1456,8 @@ PyArray_LexSort(PyObject *sort_keys, int axis)
 
     /* Now we can check the axis */
     nd = PyArray_NDIM(mps[0]);
-    if ((nd == 0) || (PyArray_SIZE(mps[0]) == 1)) {
-        /* single element case */
+    if ((nd == 0) || (PyArray_SIZE(mps[0]) <= 1)) {
+        /* empty/single element case */
         ret = (PyArrayObject *)PyArray_NewFromDescr(
             &PyArray_Type, PyArray_DescrFromType(NPY_INTP),
             PyArray_NDIM(mps[0]), PyArray_DIMS(mps[0]), NULL, NULL,
@@ -1466,7 +1466,9 @@ PyArray_LexSort(PyObject *sort_keys, int axis)
         if (ret == NULL) {
             goto fail;
         }
-        *((npy_intp *)(PyArray_DATA(ret))) = 0;
+        if (PyArray_SIZE(mps[0]) > 0) {
+            *((npy_intp *)(PyArray_DATA(ret))) = 0;
+        }
         goto finish;
     }
     if (check_and_adjust_axis(&axis, nd) < 0) {
@@ -1516,19 +1518,28 @@ PyArray_LexSort(PyObject *sort_keys, int axis)
         char *valbuffer, *indbuffer;
         int *swaps;
 
-        if (N == 0 || maxelsize == 0 || sizeof(npy_intp) == 0) {
-            goto fail;
+        assert(N > 0);  /* Guaranteed and assumed by indbuffer */
+        int valbufsize = N * maxelsize;
+        if (NPY_UNLIKELY(valbufsize) == 0) {
+            valbufsize = 1;  /* Ensure allocation is not empty */
         }
-        valbuffer = PyDataMem_NEW(N * maxelsize);
+
+        valbuffer = PyDataMem_NEW(valbufsize);
         if (valbuffer == NULL) {
             goto fail;
         }
         indbuffer = PyDataMem_NEW(N * sizeof(npy_intp));
         if (indbuffer == NULL) {
+            PyDataMem_FREE(valbuffer);
+            goto fail;
+        }
+        swaps = malloc(NPY_LIKELY(n > 0) ? n * sizeof(int) : 1);
+        if (swaps == NULL) {
+            PyDataMem_FREE(valbuffer);
             PyDataMem_FREE(indbuffer);
             goto fail;
         }
-        swaps = malloc(n*sizeof(int));
+
         for (j = 0; j < n; j++) {
             swaps[j] = PyArray_ISBYTESWAPPED(mps[j]);
         }
@@ -1557,8 +1568,8 @@ PyArray_LexSort(PyObject *sort_keys, int axis)
 #else
                 if (rcode < 0) {
 #endif
-                    npy_free_cache(valbuffer, N * maxelsize);
-                    npy_free_cache(indbuffer, N * sizeof(npy_intp));
+                    PyDataMem_FREE(valbuffer);
+                    PyDataMem_FREE(indbuffer);
                     free(swaps);
                     goto fail;
                 }
diff --git a/numpy/core/src/multiarray/iterators.c b/numpy/core/src/multiarray/iterators.c
index 83eafaf74..0d7679fe7 100644
--- a/numpy/core/src/multiarray/iterators.c
+++ b/numpy/core/src/multiarray/iterators.c
@@ -116,10 +116,12 @@ get_ptr_simple(PyArrayIterObject* iter, npy_intp *coordinates)
  * This is common initialization code between PyArrayIterObject and
  * PyArrayNeighborhoodIterObject
  *
- * Increase ao refcount
+ * Steals a reference to the array object which gets removed at deallocation,
+ * if the iterator is allocated statically and its dealloc not called, it
+ * can be thought of as borrowing the reference.
  */
-static PyObject *
-array_iter_base_init(PyArrayIterObject *it, PyArrayObject *ao)
+NPY_NO_EXPORT void
+PyArray_RawIterBaseInit(PyArrayIterObject *it, PyArrayObject *ao)
 {
     int nd, i;
 
@@ -131,7 +133,6 @@ array_iter_base_init(PyArrayIterObject *it, PyArrayObject *ao)
     else {
         it->contiguous = 0;
     }
-    Py_INCREF(ao);
     it->ao = ao;
     it->size = PyArray_SIZE(ao);
     it->nd_m1 = nd - 1;
@@ -155,7 +156,7 @@ array_iter_base_init(PyArrayIterObject *it, PyArrayObject *ao)
     it->translate = &get_ptr_simple;
     PyArray_ITER_RESET(it);
 
-    return (PyObject *)it;
+    return;
 }
 
 static void
@@ -170,6 +171,10 @@ array_iter_base_dealloc(PyArrayIterObject *it)
 NPY_NO_EXPORT PyObject *
 PyArray_IterNew(PyObject *obj)
 {
+    /*
+     * Note that internall PyArray_RawIterBaseInit may be called directly on a
+     * statically allocated PyArrayIterObject.
+     */
     PyArrayIterObject *it;
     PyArrayObject *ao;
 
@@ -186,7 +191,8 @@ PyArray_IterNew(PyObject *obj)
         return NULL;
     }
 
-    array_iter_base_init(it, ao);
+    Py_INCREF(ao);  /* PyArray_RawIterBaseInit steals a reference */
+    PyArray_RawIterBaseInit(it, ao);
     return (PyObject *)it;
 }
 
@@ -390,6 +396,10 @@ arrayiter_next(PyArrayIterObject *it)
 static void
 arrayiter_dealloc(PyArrayIterObject *it)
 {
+    /*
+     * Note that it is possible to statically allocate a PyArrayIterObject,
+     * which does not call this function.
+     */
     array_iter_base_dealloc(it);
     PyArray_free(it);
 }
@@ -1779,7 +1789,8 @@ PyArray_NeighborhoodIterNew(PyArrayIterObject *x, npy_intp *bounds,
     }
     PyObject_Init((PyObject *)ret, &PyArrayNeighborhoodIter_Type);
 
-    array_iter_base_init((PyArrayIterObject*)ret, x->ao);
+    Py_INCREF(x->ao);  /* PyArray_RawIterBaseInit steals a reference */
+    PyArray_RawIterBaseInit((PyArrayIterObject*)ret, x->ao);
     Py_INCREF(x);
     ret->_internal_iter = x;
 
diff --git a/numpy/core/src/multiarray/iterators.h b/numpy/core/src/multiarray/iterators.h
index 376dc154a..d942f45b8 100644
--- a/numpy/core/src/multiarray/iterators.h
+++ b/numpy/core/src/multiarray/iterators.h
@@ -7,4 +7,7 @@ NPY_NO_EXPORT PyObject
 NPY_NO_EXPORT int
 iter_ass_subscript(PyArrayIterObject *, PyObject *, PyObject *);
 
+NPY_NO_EXPORT void
+PyArray_RawIterBaseInit(PyArrayIterObject *it, PyArrayObject *ao);
+
 #endif
diff --git a/numpy/core/src/multiarray/mapping.c b/numpy/core/src/multiarray/mapping.c
index add1143b2..9bb85e320 100644
--- a/numpy/core/src/multiarray/mapping.c
+++ b/numpy/core/src/multiarray/mapping.c
@@ -1699,7 +1699,7 @@ array_subscript(PyArrayObject *self, PyObject *op)
                 PyArray_SHAPE(tmp_arr),
                 PyArray_STRIDES(tmp_arr),
                 PyArray_BYTES(tmp_arr),
-                PyArray_FLAGS(self),
+                PyArray_FLAGS(tmp_arr),
                 (PyObject *)self, (PyObject *)tmp_arr);
         Py_DECREF(tmp_arr);
         if (result == NULL) {
diff --git a/numpy/core/src/multiarray/refcount.c b/numpy/core/src/multiarray/refcount.c
index b8230c81a..6033929d9 100644
--- a/numpy/core/src/multiarray/refcount.c
+++ b/numpy/core/src/multiarray/refcount.c
@@ -11,6 +11,7 @@
 #define _MULTIARRAYMODULE
 #include "numpy/arrayobject.h"
 #include "numpy/arrayscalars.h"
+#include "iterators.h"
 
 #include "npy_config.h"
 
@@ -210,21 +211,22 @@ PyArray_XDECREF(PyArrayObject *mp)
     npy_intp i, n;
     PyObject **data;
     PyObject *temp;
-    PyArrayIterObject *it;
+    /*
+     * statically allocating it allows this function to not modify the
+     * reference count of the array for use during dealloc.
+     * (statically is not necessary as such)
+     */
+    PyArrayIterObject it;
 
     if (!PyDataType_REFCHK(PyArray_DESCR(mp))) {
         return 0;
     }
     if (PyArray_DESCR(mp)->type_num != NPY_OBJECT) {
-        it = (PyArrayIterObject *)PyArray_IterNew((PyObject *)mp);
-        if (it == NULL) {
-            return -1;
+        PyArray_RawIterBaseInit(&it, mp);
+        while(it.index < it.size) {
+            PyArray_Item_XDECREF(it.dataptr, PyArray_DESCR(mp));
+            PyArray_ITER_NEXT(&it);
         }
-        while(it->index < it->size) {
-            PyArray_Item_XDECREF(it->dataptr, PyArray_DESCR(mp));
-            PyArray_ITER_NEXT(it);
-        }
-        Py_DECREF(it);
         return 0;
     }
 
@@ -242,16 +244,12 @@ PyArray_XDECREF(PyArrayObject *mp)
         }
     }
     else { /* handles misaligned data too */
-        it = (PyArrayIterObject *)PyArray_IterNew((PyObject *)mp);
-        if (it == NULL) {
-            return -1;
-        }
-        while(it->index < it->size) {
-            NPY_COPY_PYOBJECT_PTR(&temp, it->dataptr);
+        PyArray_RawIterBaseInit(&it, mp);
+        while(it.index < it.size) {
+            NPY_COPY_PYOBJECT_PTR(&temp, it.dataptr);
             Py_XDECREF(temp);
-            PyArray_ITER_NEXT(it);
+            PyArray_ITER_NEXT(&it);
         }
-        Py_DECREF(it);
     }
     return 0;
 }
diff --git a/numpy/core/src/umath/cpuid.c b/numpy/core/src/umath/cpuid.c
index 8673f1736..72c6493e8 100644
--- a/numpy/core/src/umath/cpuid.c
+++ b/numpy/core/src/umath/cpuid.c
@@ -48,6 +48,25 @@ int os_avx512_support(void)
 #endif
 }
 
+static NPY_INLINE
+int cpu_supports_fma(void)
+{
+#ifdef __x86_64__
+    unsigned int feature = 0x01;
+    unsigned int a, b, c, d;
+    __asm__ volatile (
+        "cpuid"				"\n\t"
+	: "=a" (a), "=b" (b), "=c" (c), "=d" (d)
+	: "a" (feature));
+    /*
+     * FMA is the 12th bit of ECX
+     */
+    return (c >> 12) & 1;
+#else
+    return 0;
+#endif
+}
+
 /*
  * Primitive cpu feature detect function
  * Currently only supports checking for avx on gcc compatible compilers.
@@ -63,6 +82,9 @@ npy_cpu_supports(const char * feature)
         return 0;
 #endif
     }
+    else if (strcmp(feature, "fma") == 0) {
+        return cpu_supports_fma() && __builtin_cpu_supports("avx2") && os_avx_support();
+    }
     else if (strcmp(feature, "avx2") == 0) {
         return __builtin_cpu_supports("avx2") && os_avx_support();
     }
diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 1a4885133..2028a5712 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -1594,8 +1594,8 @@ NPY_NO_EXPORT void
 /**end repeat**/
 
 /**begin repeat
- *  #func = exp, log#
- *  #scalarf = npy_expf, npy_logf#
+ *  #func = sin, cos, exp, log#
+ *  #scalarf = npy_sinf, npy_cosf, npy_expf, npy_logf#
  */
 
 NPY_NO_EXPORT NPY_GCC_OPT_3 void
@@ -1610,8 +1610,8 @@ FLOAT_@func@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSE
 /**end repeat**/
 
 /**begin repeat
- * #isa = avx512f, avx2#
- * #ISA = AVX512F, AVX2#
+ * #isa = avx512f, fma#
+ * #ISA = AVX512F, FMA#
  * #CHK = HAVE_ATTRIBUTE_TARGET_AVX512F_WITH_INTRINSICS, HAVE_ATTRIBUTE_TARGET_AVX2_WITH_INTRINSICS#
  */
 
@@ -1642,6 +1642,31 @@ FLOAT_@func@_@isa@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY
 }
 
 /**end repeat1**/
+
+/**begin repeat1
+ *  #func = cos, sin#
+ *  #enum = npy_compute_cos, npy_compute_sin#
+ *  #scalarf = npy_cosf, npy_sinf#
+ */
+
+NPY_NO_EXPORT NPY_GCC_OPT_3 void
+FLOAT_@func@_@isa@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(data))
+{
+    if (!run_unary_@isa@_sincos_FLOAT(args, dimensions, steps, @enum@)) {
+        UNARY_LOOP {
+#if defined @CHK@ && defined NPY_HAVE_SSE2_INTRINSICS
+            @ISA@_sincos_FLOAT((npy_float *)op1, (npy_float *)ip1, 1, steps[0], @enum@);
+#else
+	        const npy_float in1 = *(npy_float *)ip1;
+	        *(npy_float *)op1 = @scalarf@(in1);
+#endif
+        }
+    }
+}
+
+/**end repeat1**/
+
+
 /**end repeat**/
 
 /**begin repeat
diff --git a/numpy/core/src/umath/loops.h.src b/numpy/core/src/umath/loops.h.src
index 7f05a693a..5070ab38b 100644
--- a/numpy/core/src/umath/loops.h.src
+++ b/numpy/core/src/umath/loops.h.src
@@ -178,13 +178,13 @@ NPY_NO_EXPORT void
 /**end repeat**/
 
 /**begin repeat
- *  #func = exp, log#
+ *  #func = sin, cos, exp, log#
  */
 NPY_NO_EXPORT void
 FLOAT_@func@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func));
 
 /**begin repeat1
- * #isa = avx512f, avx2#
+ * #isa = avx512f, fma#
  */
 
 NPY_NO_EXPORT void
diff --git a/numpy/core/src/umath/simd.inc.src b/numpy/core/src/umath/simd.inc.src
index 9816a1da4..7aec1ff49 100644
--- a/numpy/core/src/umath/simd.inc.src
+++ b/numpy/core/src/umath/simd.inc.src
@@ -130,8 +130,9 @@ abs_ptrdiff(char *a, char *b)
  */
 
 /**begin repeat
- * #ISA = AVX2, AVX512F#
- * #isa = avx2, avx512f#
+ * #ISA = FMA, AVX512F#
+ * #isa = fma, avx512f#
+ * #CHK = HAVE_ATTRIBUTE_TARGET_AVX2_WITH_INTRINSICS, HAVE_ATTRIBUTE_TARGET_AVX512F_WITH_INTRINSICS#
  * #REGISTER_SIZE = 32, 64#
  */
 
@@ -141,7 +142,7 @@ abs_ptrdiff(char *a, char *b)
  * #func = exp, log#
  */
 
-#if defined HAVE_ATTRIBUTE_TARGET_@ISA@_WITH_INTRINSICS && defined NPY_HAVE_SSE2_INTRINSICS
+#if defined @CHK@ && defined NPY_HAVE_SSE2_INTRINSICS
 static NPY_INLINE void
 @ISA@_@func@_FLOAT(npy_float *, npy_float *, const npy_intp n, const npy_intp stride);
 #endif
@@ -149,7 +150,7 @@ static NPY_INLINE void
 static NPY_INLINE int
 run_unary_@isa@_@func@_FLOAT(char **args, npy_intp *dimensions, npy_intp *steps)
 {
-#if defined HAVE_ATTRIBUTE_TARGET_@ISA@_WITH_INTRINSICS && defined NPY_HAVE_SSE2_INTRINSICS
+#if defined @CHK@ && defined NPY_HAVE_SSE2_INTRINSICS
     if (IS_OUTPUT_BLOCKABLE_UNARY(sizeof(npy_float), @REGISTER_SIZE@)) {
         @ISA@_@func@_FLOAT((npy_float*)args[1], (npy_float*)args[0], dimensions[0], steps[0]);
         return 1;
@@ -162,6 +163,25 @@ run_unary_@isa@_@func@_FLOAT(char **args, npy_intp *dimensions, npy_intp *steps)
 
 /**end repeat1**/
 
+#if defined @CHK@ && defined NPY_HAVE_SSE2_INTRINSICS
+static NPY_INLINE void
+@ISA@_sincos_FLOAT(npy_float *, npy_float *, const npy_intp n, const npy_intp steps, NPY_TRIG_OP);
+#endif
+
+static NPY_INLINE int
+run_unary_@isa@_sincos_FLOAT(char **args, npy_intp *dimensions, npy_intp *steps, NPY_TRIG_OP my_trig_op)
+{
+#if defined @CHK@ && defined NPY_HAVE_SSE2_INTRINSICS
+    if (IS_OUTPUT_BLOCKABLE_UNARY(sizeof(npy_float), @REGISTER_SIZE@)) {
+        @ISA@_sincos_FLOAT((npy_float*)args[1], (npy_float*)args[0], dimensions[0], steps[0], my_trig_op);
+        return 1;
+    }
+    else
+        return 0;
+#endif
+    return 0;
+}
+
 /**end repeat**/
 
 
@@ -1123,20 +1143,14 @@ sse2_@kind@_@TYPE@(@type@ * ip, @type@ * op, const npy_intp n)
 /* bunch of helper functions used in ISA_exp/log_FLOAT*/
 
 #if defined HAVE_ATTRIBUTE_TARGET_AVX2_WITH_INTRINSICS
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_fmadd(__m256 a, __m256 b, __m256 c)
-{
-    return _mm256_add_ps(_mm256_mul_ps(a, b), c);
-}
-
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_get_full_load_mask(void)
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_get_full_load_mask(void)
 {
     return _mm256_set1_ps(-1.0);
 }
 
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_get_partial_load_mask(const npy_int num_lanes, const npy_int total_elem)
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_get_partial_load_mask(const npy_int num_lanes, const npy_int total_elem)
 {
     float maskint[16] = {-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,
                             1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
@@ -1144,8 +1158,8 @@ avx2_get_partial_load_mask(const npy_int num_lanes, const npy_int total_elem)
     return _mm256_loadu_ps(addr);
 }
 
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_masked_gather(__m256 src,
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_masked_gather(__m256 src,
                    npy_float* addr,
                    __m256i vindex,
                    __m256 mask)
@@ -1153,26 +1167,39 @@ avx2_masked_gather(__m256 src,
     return _mm256_mask_i32gather_ps(src, addr, vindex, mask, 4);
 }
 
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_masked_load(__m256 mask, npy_float* addr)
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_masked_load(__m256 mask, npy_float* addr)
 {
     return _mm256_maskload_ps(addr, _mm256_cvtps_epi32(mask));
 }
 
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_set_masked_lanes(__m256 x, __m256 val, __m256 mask)
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_set_masked_lanes(__m256 x, __m256 val, __m256 mask)
 {
     return _mm256_blendv_ps(x, val, mask);
 }
 
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_blend(__m256 x, __m256 y, __m256 ymask)
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_blend(__m256 x, __m256 y, __m256 ymask)
 {
     return _mm256_blendv_ps(x, y, ymask);
 }
 
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_get_exponent(__m256 x)
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_should_calculate_sine(__m256i k, __m256i andop, __m256i cmp)
+{
+   return _mm256_cvtepi32_ps(
+                _mm256_cmpeq_epi32(_mm256_and_si256(k, andop), cmp));
+}
+
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_should_negate(__m256i k, __m256i andop, __m256i cmp)
+{
+    return fma_should_calculate_sine(k, andop, cmp);
+}
+
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_get_exponent(__m256 x)
 {
     /*
      * Special handling of denormals:
@@ -1198,8 +1225,8 @@ avx2_get_exponent(__m256 x)
     return _mm256_blendv_ps(exp, denorm_exp, denormal_mask);
 }
 
-static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_get_mantissa(__m256 x)
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_FMA __m256
+fma_get_mantissa(__m256 x)
 {
     /*
      * Special handling of denormals:
@@ -1225,7 +1252,7 @@ avx2_get_mantissa(__m256 x)
 }
 
 static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX2 __m256
-avx2_scalef_ps(__m256 poly, __m256 quadrant)
+fma_scalef_ps(__m256 poly, __m256 quadrant)
 {
     /*
      * Handle denormals (which occur when quadrant <= -125):
@@ -1305,6 +1332,18 @@ avx512_blend(__m512 x, __m512 y, __mmask16 ymask)
     return _mm512_mask_mov_ps(x, ymask, y);
 }
 
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX512F __mmask16
+avx512_should_calculate_sine(__m512i k, __m512i andop, __m512i cmp)
+{
+    return _mm512_cmpeq_epi32_mask(_mm512_and_epi32(k, andop), cmp);
+}
+
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX512F __mmask16
+avx512_should_negate(__m512i k, __m512i andop, __m512i cmp)
+{
+    return avx512_should_calculate_sine(k, andop, cmp);
+}
+
 static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_AVX512F __m512
 avx512_get_exponent(__m512 x)
 {
@@ -1325,17 +1364,28 @@ avx512_scalef_ps(__m512 poly, __m512 quadrant)
 #endif
 
 /**begin repeat
- * #ISA = AVX2, AVX512F#
- * #isa = avx2, avx512#
+ * #ISA = FMA, AVX512F#
+ * #isa = fma, avx512#
  * #vtype = __m256, __m512#
  * #vsize = 256, 512#
  * #or = or_ps, kor#
  * #vsub = , _mask#
  * #mask = __m256, __mmask16#
- * #fmadd = avx2_fmadd,_mm512_fmadd_ps#
+ * #fmadd = _mm256_fmadd_ps, _mm512_fmadd_ps#
+ * #CHK = HAVE_ATTRIBUTE_TARGET_AVX2_WITH_INTRINSICS, HAVE_ATTRIBUTE_TARGET_AVX512F_WITH_INTRINSICS#
  **/
 
-#if defined HAVE_ATTRIBUTE_TARGET_@ISA@_WITH_INTRINSICS
+#if defined @CHK@
+
+/*
+ * Vectorized Cody-Waite range reduction technique
+ * Performs the reduction step x* = x - y*C in three steps:
+ * 1) x* = x - y*c1
+ * 2) x* = x - y*c2
+ * 3) x* = x - y*c3
+ * c1, c2 are exact floating points, c3 = C - c1 - c2 simulates higher precision
+ */
+
 static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ @vtype@
 @isa@_range_reduction(@vtype@ x, @vtype@ y, @vtype@ c1, @vtype@ c2, @vtype@ c3)
 {
@@ -1344,12 +1394,56 @@ static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ @vtype@
     reduced_x = @fmadd@(y, c3, reduced_x);
     return reduced_x;
 }
+
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ @mask@
+@isa@_in_range_mask(@vtype@ x, npy_float fmax, npy_float fmin)
+{
+    @mask@ m1 = _mm@vsize@_cmp_ps@vsub@(
+                                x, _mm@vsize@_set1_ps(fmax), _CMP_GT_OQ);
+    @mask@ m2 = _mm@vsize@_cmp_ps@vsub@(
+                                x, _mm@vsize@_set1_ps(fmin), _CMP_LT_OQ);
+    return _mm@vsize@_@or@(m1,m2);
+}
+
+/*
+ * Approximate cosine algorithm for x \in [-PI/4, PI/4]
+ * Maximum ULP across all 32-bit floats = 0.875
+ */
+
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ @vtype@
+@isa@_cosine(@vtype@ x2, @vtype@ invf8, @vtype@ invf6, @vtype@ invf4,
+                                                @vtype@ invf2, @vtype@ invf0)
+{
+    @vtype@ cos = @fmadd@(invf8, x2, invf6);
+    cos = @fmadd@(cos, x2, invf4);
+    cos = @fmadd@(cos, x2, invf2);
+    cos = @fmadd@(cos, x2, invf0);
+    return cos;
+}
+
+/*
+ * Approximate sine algorithm for x \in [-PI/4, PI/4]
+ * Maximum ULP across all 32-bit floats = 0.647
+ */
+
+static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ @vtype@
+@isa@_sine(@vtype@ x, @vtype@ x2, @vtype@ invf9, @vtype@ invf7,
+                                          @vtype@ invf5, @vtype@ invf3,
+                                          @vtype@ zero)
+{
+    @vtype@ sin = @fmadd@(invf9, x2, invf7);
+    sin = @fmadd@(sin, x2, invf5);
+    sin = @fmadd@(sin, x2, invf3);
+    sin = @fmadd@(sin, x2, zero);
+    sin = @fmadd@(sin, x, x);
+    return sin;
+}
 #endif
 /**end repeat**/
 
 /**begin repeat
- * #ISA = AVX2, AVX512F#
- * #isa = avx2, avx512#
+ * #ISA = FMA, AVX512F#
+ * #isa = fma, avx512#
  * #vtype = __m256, __m512#
  * #vsize = 256, 512#
  * #BYTES = 32, 64#
@@ -1358,13 +1452,165 @@ static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ @vtype@
  * #or_masks =_mm256_or_ps, _mm512_kor#
  * #and_masks =_mm256_and_ps, _mm512_kand#
  * #xor_masks =_mm256_xor_ps, _mm512_kxor#
- * #fmadd = avx2_fmadd,_mm512_fmadd_ps#
+ * #fmadd = _mm256_fmadd_ps, _mm512_fmadd_ps#
  * #mask_to_int = _mm256_movemask_ps, #
  * #full_mask= 0xFF, 0xFFFF#
  * #masked_store = _mm256_maskstore_ps, _mm512_mask_storeu_ps#
  * #cvtps_epi32 = _mm256_cvtps_epi32, #
+ * #CHK = HAVE_ATTRIBUTE_TARGET_AVX2_WITH_INTRINSICS, HAVE_ATTRIBUTE_TARGET_AVX512F_WITH_INTRINSICS#
+ */
+
+
+/*
+ * Vectorized approximate sine/cosine algorithms: The following code is a
+ * vectorized version of the algorithm presented here:
+ * https://stackoverflow.com/questions/30463616/payne-hanek-algorithm-implementation-in-c/30465751#30465751
+ * (1) Load data in ZMM/YMM registers and generate mask for elements that are
+ * within range [-71476.0625f, 71476.0625f] for cosine and [-117435.992f,
+ * 117435.992f] for sine.
+ * (2) For elements within range, perform range reduction using Cody-Waite's
+ * method: x* = x - y*PI/2, where y = rint(x*2/PI). x* \in [-PI/4, PI/4].
+ * (3) Map cos(x) to (+/-)sine or (+/-)cosine of x* based on the quadrant k =
+ * int(y).
+ * (4) For elements outside that range, Cody-Waite reduction peforms poorly
+ * leading to catastrophic cancellation. We compute cosine by calling glibc in
+ * a scalar fashion.
+ * (5) Vectorized implementation has a max ULP of 1.49 and performs at least
+ * 5-7x faster than scalar implementations when magnitude of all elements in
+ * the array < 71476.0625f (117435.992f for sine). Worst case performance is
+ * when all the elements are large leading to about 1-2% reduction in
+ * performance.
  */
 
+#if defined @CHK@
+static NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ void
+@ISA@_sincos_FLOAT(npy_float * op,
+                   npy_float * ip,
+                   const npy_intp array_size,
+                   const npy_intp steps,
+                   NPY_TRIG_OP my_trig_op)
+{
+    const npy_intp stride = steps/sizeof(npy_float);
+    const npy_int num_lanes = @BYTES@/sizeof(npy_float);
+    npy_float large_number = 71476.0625f;
+    if (my_trig_op == npy_compute_sin) {
+        large_number = 117435.992f;
+    }
+
+    /* Load up frequently used constants */
+    @vtype@i zeros = _mm@vsize@_set1_epi32(0);
+    @vtype@i ones = _mm@vsize@_set1_epi32(1);
+    @vtype@i twos = _mm@vsize@_set1_epi32(2);
+    @vtype@ two_over_pi = _mm@vsize@_set1_ps(NPY_TWO_O_PIf);
+    @vtype@ codyw_c1 = _mm@vsize@_set1_ps(NPY_CODY_WAITE_PI_O_2_HIGHf);
+    @vtype@ codyw_c2 = _mm@vsize@_set1_ps(NPY_CODY_WAITE_PI_O_2_MEDf);
+    @vtype@ codyw_c3 = _mm@vsize@_set1_ps(NPY_CODY_WAITE_PI_O_2_LOWf);
+    @vtype@ cos_invf0 = _mm@vsize@_set1_ps(NPY_COEFF_INVF0_COSINEf);
+    @vtype@ cos_invf2 = _mm@vsize@_set1_ps(NPY_COEFF_INVF2_COSINEf);
+    @vtype@ cos_invf4 = _mm@vsize@_set1_ps(NPY_COEFF_INVF4_COSINEf);
+    @vtype@ cos_invf6 = _mm@vsize@_set1_ps(NPY_COEFF_INVF6_COSINEf);
+    @vtype@ cos_invf8 = _mm@vsize@_set1_ps(NPY_COEFF_INVF8_COSINEf);
+    @vtype@ sin_invf3 = _mm@vsize@_set1_ps(NPY_COEFF_INVF3_SINEf);
+    @vtype@ sin_invf5 = _mm@vsize@_set1_ps(NPY_COEFF_INVF5_SINEf);
+    @vtype@ sin_invf7 = _mm@vsize@_set1_ps(NPY_COEFF_INVF7_SINEf);
+    @vtype@ sin_invf9 = _mm@vsize@_set1_ps(NPY_COEFF_INVF9_SINEf);
+    @vtype@ cvt_magic = _mm@vsize@_set1_ps(NPY_RINT_CVT_MAGICf);
+    @vtype@ zero_f = _mm@vsize@_set1_ps(0.0f);
+    @vtype@ quadrant, reduced_x, reduced_x2, cos, sin;
+    @vtype@i iquadrant;
+    @mask@ nan_mask, glibc_mask, sine_mask, negate_mask;
+    @mask@ load_mask = @isa@_get_full_load_mask();
+    npy_intp num_remaining_elements = array_size;
+    npy_int indexarr[16];
+    for (npy_int ii = 0; ii < 16; ii++) {
+        indexarr[ii] = ii*stride;
+    }
+    @vtype@i vindex = _mm@vsize@_loadu_si@vsize@((@vtype@i*)&indexarr[0]);
+
+    while (num_remaining_elements > 0) {
+
+        if (num_remaining_elements < num_lanes) {
+            load_mask = @isa@_get_partial_load_mask(num_remaining_elements,
+                                                         num_lanes);
+        }
+
+        @vtype@ x;
+        if (stride == 1) {
+            x = @isa@_masked_load(load_mask, ip);
+        }
+        else {
+            x = @isa@_masked_gather(zero_f, ip, vindex, load_mask);
+        }
+
+        /*
+         * For elements outside of this range, Cody-Waite's range reduction
+         * becomes inaccurate and we will call glibc to compute cosine for
+         * these numbers
+         */
+
+        glibc_mask = @isa@_in_range_mask(x, large_number,-large_number);
+        glibc_mask = @and_masks@(load_mask, glibc_mask);
+        nan_mask = _mm@vsize@_cmp_ps@vsub@(x, x, _CMP_NEQ_UQ);
+        x = @isa@_set_masked_lanes(x, zero_f, @or_masks@(nan_mask, glibc_mask));
+        npy_int iglibc_mask = @mask_to_int@(glibc_mask);
+
+        if (iglibc_mask != @full_mask@) {
+            quadrant = _mm@vsize@_mul_ps(x, two_over_pi);
+
+            /* round to nearest */
+            quadrant = _mm@vsize@_add_ps(quadrant, cvt_magic);
+            quadrant = _mm@vsize@_sub_ps(quadrant, cvt_magic);
+
+            /* Cody-Waite's range reduction algorithm */
+            reduced_x = @isa@_range_reduction(x, quadrant,
+                                                   codyw_c1, codyw_c2, codyw_c3);
+            reduced_x2 = _mm@vsize@_mul_ps(reduced_x, reduced_x);
+
+            /* compute cosine and sine */
+            cos = @isa@_cosine(reduced_x2, cos_invf8, cos_invf6, cos_invf4,
+                                                           cos_invf2, cos_invf0);
+            sin = @isa@_sine(reduced_x, reduced_x2, sin_invf9, sin_invf7,
+                                             sin_invf5, sin_invf3, zero_f);
+
+            iquadrant = _mm@vsize@_cvtps_epi32(quadrant);
+            if (my_trig_op == npy_compute_cos) {
+                iquadrant = _mm@vsize@_add_epi32(iquadrant, ones);
+            }
+
+            /* blend sin and cos based on the quadrant */
+            sine_mask = @isa@_should_calculate_sine(iquadrant, ones, zeros);
+            cos = @isa@_blend(cos, sin, sine_mask);
+
+            /* multiply by -1 for appropriate elements */
+            negate_mask = @isa@_should_negate(iquadrant, twos, twos);
+            cos = @isa@_blend(cos, _mm@vsize@_sub_ps(zero_f, cos), negate_mask);
+            cos = @isa@_set_masked_lanes(cos, _mm@vsize@_set1_ps(NPY_NANF), nan_mask);
+
+            @masked_store@(op, @cvtps_epi32@(load_mask), cos);
+        }
+
+        /* process elements using glibc for large elements */
+        if (my_trig_op == npy_compute_cos) {
+            for (int ii = 0; iglibc_mask != 0; ii++) {
+                if (iglibc_mask & 0x01) {
+                    op[ii] = npy_cosf(ip[ii]);
+                }
+                iglibc_mask  = iglibc_mask >> 1;
+            }
+        }
+        else {
+            for (int ii = 0; iglibc_mask != 0; ii++) {
+                if (iglibc_mask & 0x01) {
+                    op[ii] = npy_sinf(ip[ii]);
+                }
+                iglibc_mask  = iglibc_mask >> 1;
+            }
+        }
+        ip += num_lanes*stride;
+        op += num_lanes;
+        num_remaining_elements -= num_lanes;
+    }
+}
 
 /*
  * Vectorized implementation of exp using AVX2 and AVX512:
@@ -1381,7 +1627,6 @@ static NPY_INLINE NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ @vtype@
  * same x = 0xc2781e37)
  */
 
-#if defined HAVE_ATTRIBUTE_TARGET_@ISA@_WITH_INTRINSICS
 static NPY_GCC_OPT_3 NPY_GCC_TARGET_@ISA@ void
 @ISA@_exp_FLOAT(npy_float * op,
                 npy_float * ip,
diff --git a/numpy/core/tests/test_arrayprint.py b/numpy/core/tests/test_arrayprint.py
index 75a794369..702e68e76 100644
--- a/numpy/core/tests/test_arrayprint.py
+++ b/numpy/core/tests/test_arrayprint.py
@@ -262,11 +262,6 @@ class TestArray2String(object):
         assert_(np.array2string(s, formatter={'numpystr':lambda s: s*2}) ==
                 '[abcabc defdef]')
 
-        # check for backcompat that using FloatFormat works and emits warning
-        with assert_warns(DeprecationWarning):
-            fmt = np.core.arrayprint.FloatFormat(x, 9, 'maxprec', False)
-        assert_equal(np.array2string(x, formatter={'float_kind': fmt}),
-                     '[0. 1. 2.]')
 
     def test_structure_format(self):
         dt = np.dtype([('name', np.str_, 16), ('grades', np.float64, (2,))])
diff --git a/numpy/core/tests/test_deprecations.py b/numpy/core/tests/test_deprecations.py
index 6d71fcbd6..e8aa0c70b 100644
--- a/numpy/core/tests/test_deprecations.py
+++ b/numpy/core/tests/test_deprecations.py
@@ -101,7 +101,7 @@ class _DeprecationTestCase(object):
                         (self.warning_cls.__name__, warning.category))
         if num is not None and num_found != num:
             msg = "%i warnings found but %i expected." % (len(self.log), num)
-            lst = [str(w.category) for w in self.log]
+            lst = [str(w) for w in self.log]
             raise AssertionError("\n".join([msg] + lst))
 
         with warnings.catch_warnings():
@@ -149,16 +149,6 @@ class TestNonTupleNDIndexDeprecation(object):
             a[[0, 1]]
 
 
-class TestRankDeprecation(_DeprecationTestCase):
-    """Test that np.rank is deprecated. The function should simply be
-    removed. The VisibleDeprecationWarning may become unnecessary.
-    """
-
-    def test(self):
-        a = np.arange(10)
-        assert_warns(np.VisibleDeprecationWarning, np.rank, a)
-
-
 class TestComparisonDeprecations(_DeprecationTestCase):
     """This tests the deprecation, for non-element-wise comparison logic.
     This used to mean that when an error occurred during element-wise comparison
@@ -499,6 +489,12 @@ class TestBincount(_DeprecationTestCase):
         self.assert_deprecated(lambda: np.bincount([1, 2, 3], minlength=None))
 
 
+class TestAlen(_DeprecationTestCase):
+    # 2019-08-02, 1.18.0
+    def test_alen(self):
+        self.assert_deprecated(lambda: np.alen(np.array([1, 2, 3])))
+
+
 class TestGeneratorSum(_DeprecationTestCase):
     # 2018-02-25, 1.15.0
     def test_generator_sum(self):
diff --git a/numpy/core/tests/test_indexing.py b/numpy/core/tests/test_indexing.py
index f7485c3f7..70a5a246f 100644
--- a/numpy/core/tests/test_indexing.py
+++ b/numpy/core/tests/test_indexing.py
@@ -617,6 +617,19 @@ class TestSubclasses(object):
         assert_array_equal(s_bool, a[a > 0])
         assert_array_equal(s_bool.base, a[a > 0])
 
+    def test_fancy_on_read_only(self):
+        # Test that fancy indexing on read-only SubClass does not make a
+        # read-only copy (gh-14132)
+        class SubClass(np.ndarray):
+            pass
+
+        a = np.arange(5)
+        s = a.view(SubClass)
+        s.flags.writeable = False
+        s_fancy = s[[0, 1, 2]]
+        assert_(s_fancy.flags.writeable)
+
+
     def test_finalize_gets_full_info(self):
         # Array finalize should be called on the filled array.
         class SubClass(np.ndarray):
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index 53e538f7d..2593045ed 100644
--- a/numpy/core/tests/test_multiarray.py
+++ b/numpy/core/tests/test_multiarray.py
@@ -44,7 +44,7 @@ from numpy.testing import (
     assert_, assert_raises, assert_warns, assert_equal, assert_almost_equal,
     assert_array_equal, assert_raises_regex, assert_array_almost_equal,
     assert_allclose, IS_PYPY, HAS_REFCOUNT, assert_array_less, runstring,
-    temppath, suppress_warnings, break_cycles, assert_raises_regex,
+    temppath, suppress_warnings, break_cycles,
     )
 from numpy.core.tests._locales import CommaDecimalPointLocale
 
@@ -497,9 +497,6 @@ class TestArrayConstruction(object):
         assert_(np.ascontiguousarray(d).flags.c_contiguous)
         assert_(np.asfortranarray(d).flags.f_contiguous)
 
-    def test_ragged(self):
-        assert_raises_regex(ValueError, 'ragged',
-                             np.array, [[1], [2, 3]], dtype=int)
 
 class TestAssignment(object):
     def test_assignment_broadcasting(self):
@@ -6405,20 +6402,22 @@ class TestInner(object):
 
 class TestAlen(object):
     def test_basic(self):
-        m = np.array([1, 2, 3])
-        assert_equal(np.alen(m), 3)
+        with pytest.warns(DeprecationWarning):
+            m = np.array([1, 2, 3])
+            assert_equal(np.alen(m), 3)
 
-        m = np.array([[1, 2, 3], [4, 5, 7]])
-        assert_equal(np.alen(m), 2)
+            m = np.array([[1, 2, 3], [4, 5, 7]])
+            assert_equal(np.alen(m), 2)
 
-        m = [1, 2, 3]
-        assert_equal(np.alen(m), 3)
+            m = [1, 2, 3]
+            assert_equal(np.alen(m), 3)
 
-        m = [[1, 2, 3], [4, 5, 7]]
-        assert_equal(np.alen(m), 2)
+            m = [[1, 2, 3], [4, 5, 7]]
+            assert_equal(np.alen(m), 2)
 
     def test_singleton(self):
-        assert_equal(np.alen(5), 1)
+        with pytest.warns(DeprecationWarning):
+            assert_equal(np.alen(5), 1)
 
 
 class TestChoose(object):
@@ -8080,6 +8079,8 @@ class TestWritebackIfCopy(object):
         arr_wb[...] = 100
         assert_equal(arr, -100)
 
+    @pytest.mark.leaks_references(
+            reason="increments self in dealloc; ignore since deprecated path.")
     def test_dealloc_warning(self):
         with suppress_warnings() as sup:
             sup.record(RuntimeWarning)
diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py
index 3e85054b7..c479a0f6d 100644
--- a/numpy/core/tests/test_numeric.py
+++ b/numpy/core/tests/test_numeric.py
@@ -1341,6 +1341,11 @@ class TestBinaryRepr(object):
             exp = '1' + (width - 1) * '0'
             assert_equal(np.binary_repr(num, width=width), exp)
 
+    def test_large_neg_int64(self):
+        # See gh-14289.
+        assert_equal(np.binary_repr(np.int64(-2**62), width=64),
+                     '11' + '0'*62)
+
 
 class TestBaseRepr(object):
     def test_base3(self):
diff --git a/numpy/core/tests/test_records.py b/numpy/core/tests/test_records.py
index 14413224e..c1b794145 100644
--- a/numpy/core/tests/test_records.py
+++ b/numpy/core/tests/test_records.py
@@ -444,6 +444,48 @@ class TestRecord(object):
         ]
         arr = np.rec.fromarrays(arrays)  # ValueError?
 
+    @pytest.mark.parametrize('nfields', [0, 1, 2])
+    def test_assign_dtype_attribute(self, nfields):
+        dt = np.dtype([('a', np.uint8), ('b', np.uint8), ('c', np.uint8)][:nfields])
+        data = np.zeros(3, dt).view(np.recarray)
+
+        # the original and resulting dtypes differ on whether they are records
+        assert data.dtype.type == np.record
+        assert dt.type != np.record
+
+        # ensure that the dtype remains a record even when assigned
+        data.dtype = dt
+        assert data.dtype.type == np.record
+
+    @pytest.mark.parametrize('nfields', [0, 1, 2])
+    def test_nested_fields_are_records(self, nfields):
+        """ Test that nested structured types are treated as records too """
+        dt = np.dtype([('a', np.uint8), ('b', np.uint8), ('c', np.uint8)][:nfields])
+        dt_outer = np.dtype([('inner', dt)])
+
+        data = np.zeros(3, dt_outer).view(np.recarray)
+        assert isinstance(data, np.recarray)
+        assert isinstance(data['inner'], np.recarray)
+
+        data0 = data[0]
+        assert isinstance(data0, np.record)
+        assert isinstance(data0['inner'], np.record)
+
+    def test_nested_dtype_padding(self):
+        """ test that trailing padding is preserved """
+        # construct a dtype with padding at the end
+        dt = np.dtype([('a', np.uint8), ('b', np.uint8), ('c', np.uint8)])
+        dt_padded_end = dt[['a', 'b']]
+        assert dt_padded_end.itemsize == dt.itemsize
+
+        dt_outer = np.dtype([('inner', dt_padded_end)])
+
+        data = np.zeros(3, dt_outer).view(np.recarray)
+        assert_equal(data['inner'].dtype, dt_padded_end)
+
+        data0 = data[0]
+        assert_equal(data0['inner'].dtype, dt_padded_end)
+
 
 def test_find_duplicate():
     l1 = [1, 2, 3, 4, 5, 6]
diff --git a/numpy/core/tests/test_regression.py b/numpy/core/tests/test_regression.py
index e564ae300..c5289f6ac 100644
--- a/numpy/core/tests/test_regression.py
+++ b/numpy/core/tests/test_regression.py
@@ -436,6 +436,32 @@ class TestRegression(object):
 
         assert_raises(KeyError, np.lexsort, BuggySequence())
 
+    def test_lexsort_zerolen_custom_strides(self):
+        # Ticket #14228
+        xs = np.array([], dtype='i8')
+        assert xs.strides == (8,)
+        assert np.lexsort((xs,)).shape[0] == 0 # Works
+
+        xs.strides = (16,)
+        assert np.lexsort((xs,)).shape[0] == 0 # Was: MemoryError
+
+    def test_lexsort_zerolen_custom_strides_2d(self):
+        xs = np.array([], dtype='i8')
+
+        xs.shape = (0, 2)
+        xs.strides = (16, 16)
+        assert np.lexsort((xs,), axis=0).shape[0] == 0
+
+        xs.shape = (2, 0)
+        xs.strides = (16, 16)
+        assert np.lexsort((xs,), axis=0).shape[0] == 2
+
+    def test_lexsort_zerolen_element(self):
+        dt = np.dtype([])  # a void dtype with no fields
+        xs = np.empty(4, dt)
+
+        assert np.lexsort((xs,)).shape[0] == xs.shape[0]
+
     def test_pickle_py2_bytes_encoding(self):
         # Check that arrays and scalars pickled on Py2 are
         # unpickleable on Py3 using encoding='bytes'
@@ -468,7 +494,7 @@ class TestRegression(object):
                 result = pickle.loads(data, encoding='bytes')
                 assert_equal(result, original)
 
-                if isinstance(result, np.ndarray) and result.dtype.names:
+                if isinstance(result, np.ndarray) and result.dtype.names is not None:
                     for name in result.dtype.names:
                         assert_(isinstance(name, str))
 
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 91d403d98..ef48fed05 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -678,20 +678,49 @@ class TestSpecialFloats(object):
             assert_raises(FloatingPointError, np.log, np.float32(-np.inf))
             assert_raises(FloatingPointError, np.log, np.float32(-1.0))
 
-class TestExpLogFloat32(object):
+    def test_sincos_values(self):
+        with np.errstate(all='ignore'):
+            x = [np.nan,  np.nan, np.nan, np.nan]
+            y = [np.nan, -np.nan, np.inf, -np.inf]
+            for dt in ['f', 'd', 'g']:
+                xf = np.array(x, dtype=dt)
+                yf = np.array(y, dtype=dt)
+                assert_equal(np.sin(yf), xf)
+                assert_equal(np.cos(yf), xf)
+
+        with np.errstate(invalid='raise'):
+            assert_raises(FloatingPointError, np.sin, np.float32(-np.inf))
+            assert_raises(FloatingPointError, np.sin, np.float32(np.inf))
+            assert_raises(FloatingPointError, np.cos, np.float32(-np.inf))
+            assert_raises(FloatingPointError, np.cos, np.float32(np.inf))
+
+
+class TestSIMDFloat32(object):
     def test_exp_float32(self):
         np.random.seed(42)
         x_f32 = np.float32(np.random.uniform(low=0.0,high=88.1,size=1000000))
         x_f64 = np.float64(x_f32)
-        assert_array_max_ulp(np.exp(x_f32), np.float32(np.exp(x_f64)), maxulp=2.6)
+        assert_array_max_ulp(np.exp(x_f32), np.float32(np.exp(x_f64)), maxulp=3)
 
     def test_log_float32(self):
         np.random.seed(42)
         x_f32 = np.float32(np.random.uniform(low=0.0,high=1000,size=1000000))
         x_f64 = np.float64(x_f32)
-        assert_array_max_ulp(np.log(x_f32), np.float32(np.log(x_f64)), maxulp=3.9)
+        assert_array_max_ulp(np.log(x_f32), np.float32(np.log(x_f64)), maxulp=4)
+
+    def test_sincos_float32(self):
+        np.random.seed(42)
+        N = 1000000
+        M = np.int(N/20)
+        index = np.random.randint(low=0, high=N, size=M)
+        x_f32 = np.float32(np.random.uniform(low=-100.,high=100.,size=N))
+        # test coverage for elements > 117435.992f for which glibc is used
+        x_f32[index] = np.float32(10E+10*np.random.rand(M))
+        x_f64 = np.float64(x_f32)
+        assert_array_max_ulp(np.sin(x_f32), np.float32(np.sin(x_f64)), maxulp=2)
+        assert_array_max_ulp(np.cos(x_f32), np.float32(np.cos(x_f64)), maxulp=2)
 
-    def test_strided_exp_log_float32(self):
+    def test_strided_float32(self):
         np.random.seed(42)
         strides = np.random.randint(low=-100, high=100, size=100)
         sizes = np.random.randint(low=1, high=2000, size=100)
@@ -699,9 +728,13 @@ class TestExpLogFloat32(object):
             x_f32 = np.float32(np.random.uniform(low=0.01,high=88.1,size=ii))
             exp_true = np.exp(x_f32)
             log_true = np.log(x_f32)
+            sin_true = np.sin(x_f32)
+            cos_true = np.cos(x_f32)
             for jj in strides:
                 assert_array_almost_equal_nulp(np.exp(x_f32[::jj]), exp_true[::jj], nulp=2)
                 assert_array_almost_equal_nulp(np.log(x_f32[::jj]), log_true[::jj], nulp=2)
+                assert_array_almost_equal_nulp(np.sin(x_f32[::jj]), sin_true[::jj], nulp=2)
+                assert_array_almost_equal_nulp(np.cos(x_f32[::jj]), cos_true[::jj], nulp=2)
 
 class TestLogAddExp(_FilterInvalids):
     def test_logaddexp_values(self):
diff --git a/numpy/ctypeslib.py b/numpy/ctypeslib.py
index 1f842d003..e49f8394f 100644
--- a/numpy/ctypeslib.py
+++ b/numpy/ctypeslib.py
@@ -321,7 +321,7 @@ def ndpointer(dtype=None, ndim=None, shape=None, flags=None):
     # produce a name for the new type
     if dtype is None:
         name = 'any'
-    elif dtype.names:
+    elif dtype.names is not None:
         name = str(id(dtype))
     else:
         name = dtype.str
diff --git a/numpy/distutils/fcompiler/environment.py b/numpy/distutils/fcompiler/environment.py
index 73a5e98e1..bb362d483 100644
--- a/numpy/distutils/fcompiler/environment.py
+++ b/numpy/distutils/fcompiler/environment.py
@@ -59,17 +59,13 @@ class EnvironmentConfig(object):
             if envvar_contents is not None:
                 envvar_contents = convert(envvar_contents)
                 if var and append:
-                    if os.environ.get('NPY_DISTUTILS_APPEND_FLAGS', '0') == '1':
+                    if os.environ.get('NPY_DISTUTILS_APPEND_FLAGS', '1') == '1':
                         var.extend(envvar_contents)
                     else:
+                        # NPY_DISTUTILS_APPEND_FLAGS was explicitly set to 0
+                        # to keep old (overwrite flags rather than append to
+                        # them) behavior
                         var = envvar_contents
-                        if 'NPY_DISTUTILS_APPEND_FLAGS' not in os.environ.keys():
-                            msg = "{} is used as is, not appended ".format(envvar) + \
-                                  "to flags already defined " + \
-                                  "by numpy.distutils! Use NPY_DISTUTILS_APPEND_FLAGS=1 " + \
-                                  "to obtain appending behavior instead (this " + \
-                                  "behavior will become default in a future release)."
-                            warnings.warn(msg, UserWarning, stacklevel=3)
                 else:
                     var = envvar_contents
         if confvar is not None and self._conf:
diff --git a/numpy/distutils/system_info.py b/numpy/distutils/system_info.py
index 2ff0ba7b3..6cfce3b1c 100644
--- a/numpy/distutils/system_info.py
+++ b/numpy/distutils/system_info.py
@@ -456,7 +456,7 @@ class AliasedOptionError(DistutilsError):
 
 class AtlasNotFoundError(NotFoundError):
     """
-    Atlas (http://math-atlas.sourceforge.net/) libraries not found.
+    Atlas (http://github.com/math-atlas/math-atlas) libraries not found.
     Directories to search for the libraries can be specified in the
     numpy/distutils/site.cfg file (section [atlas]) or by setting
     the ATLAS environment variable."""
diff --git a/numpy/distutils/tests/test_fcompiler.py b/numpy/distutils/tests/test_fcompiler.py
index ba19a97ea..6d245fbd4 100644
--- a/numpy/distutils/tests/test_fcompiler.py
+++ b/numpy/distutils/tests/test_fcompiler.py
@@ -45,37 +45,3 @@ def test_fcompiler_flags(monkeypatch):
         else:
             assert_(new_flags == prev_flags + [new_flag])
 
-
-def test_fcompiler_flags_append_warning(monkeypatch):
-    # Test to check that the warning for append behavior changing in future
-    # is triggered.  Need to use a real compiler instance so that we have
-    # non-empty flags to start with (otherwise the "if var and append" check
-    # will always be false).
-    try:
-        with suppress_warnings() as sup:
-            sup.record()
-            fc = numpy.distutils.fcompiler.new_fcompiler(compiler='gnu95')
-            fc.customize()
-    except numpy.distutils.fcompiler.CompilerNotFound:
-        pytest.skip("gfortran not found, so can't execute this test")
-
-    # Ensure NPY_DISTUTILS_APPEND_FLAGS not defined
-    monkeypatch.delenv('NPY_DISTUTILS_APPEND_FLAGS', raising=False)
-
-    for opt, envvar in customizable_flags:
-        new_flag = '-dummy-{}-flag'.format(opt)
-        with suppress_warnings() as sup:
-            sup.record()
-            prev_flags = getattr(fc.flag_vars, opt)
-
-        monkeypatch.setenv(envvar, new_flag)
-        with suppress_warnings() as sup:
-            sup.record()
-            new_flags = getattr(fc.flag_vars, opt)
-            if prev_flags:
-                # Check that warning was issued
-                assert len(sup.log) == 1
-
-        monkeypatch.delenv(envvar)
-        assert_(new_flags == [new_flag])
-
diff --git a/numpy/doc/broadcasting.py b/numpy/doc/broadcasting.py
index f7bd2515b..cb548a0d0 100644
--- a/numpy/doc/broadcasting.py
+++ b/numpy/doc/broadcasting.py
@@ -61,8 +61,7 @@ dimensions are compatible when
 If these conditions are not met, a
 ``ValueError: operands could not be broadcast together`` exception is 
 thrown, indicating that the arrays have incompatible shapes. The size of 
-the resulting array is the maximum size along each dimension of the input 
-arrays.
+the resulting array is the size that is not 1 along each axis of the inputs.
 
 Arrays do not need to have the same *number* of dimensions.  For example,
 if you have a ``256x256x3`` array of RGB values, and you want to scale
diff --git a/numpy/doc/dispatch.py b/numpy/doc/dispatch.py
index 09a3e5134..c9029941b 100644
--- a/numpy/doc/dispatch.py
+++ b/numpy/doc/dispatch.py
@@ -72,7 +72,7 @@ The ``__array_ufunc__`` receives:
 - ``inputs``, which could be a mixture of different types
 - ``kwargs``, keyword arguments passed to the function
 
-For this example we will only handle the method ``'__call__``.
+For this example we will only handle the method ``__call__``.
 
 >>> from numbers import Number
 >>> class DiagonalArray:
@@ -218,12 +218,12 @@ For completeness, to support the usage ``arr.sum()`` add a method ``sum`` that
 calls ``numpy.sum(self)``, and the same for ``mean``.
 
 >>> @implements(np.sum)
-... def sum(a):
+... def sum(arr):
 ...     "Implementation of np.sum for DiagonalArray objects"
 ...     return arr._i * arr._N
 ...
 >>> @implements(np.mean)
-... def sum(a):
+... def mean(arr):
 ...     "Implementation of np.mean for DiagonalArray objects"
 ...     return arr._i / arr._N
 ...
diff --git a/numpy/f2py/cfuncs.py b/numpy/f2py/cfuncs.py
index d59b6301c..17f3861ca 100644
--- a/numpy/f2py/cfuncs.py
+++ b/numpy/f2py/cfuncs.py
@@ -542,7 +542,7 @@ cppmacros[
     'ARRSIZE'] = '#define ARRSIZE(dims,rank) (_PyArray_multiply_list(dims,rank))'
 cppmacros['OLDPYNUM'] = """\
 #ifdef OLDPYNUM
-#error You need to install Numeric Python version 13 or higher. Get it from http:/sourceforge.net/project/?group_id=1369
+#error You need to install NumPy version 13 or higher. See https://scipy.org/install.html
 #endif
 """
 ################# C functions ###############
diff --git a/numpy/f2py/rules.py b/numpy/f2py/rules.py
index 6769f1b1f..1b41498ea 100644..100755
--- a/numpy/f2py/rules.py
+++ b/numpy/f2py/rules.py
@@ -202,7 +202,7 @@ PyMODINIT_FUNC PyInit_#modulename#(void) {
 PyMODINIT_FUNC init#modulename#(void) {
 #endif
 \tint i;
-\tPyObject *m,*d, *s;
+\tPyObject *m,*d, *s, *tmp;
 #if PY_VERSION_HEX >= 0x03000000
 \tm = #modulename#_module = PyModule_Create(&moduledef);
 #else
@@ -224,8 +224,11 @@ PyMODINIT_FUNC init#modulename#(void) {
 \tPyDict_SetItemString(d, \"__doc__\", s);
 \t#modulename#_error = PyErr_NewException (\"#modulename#.error\", NULL, NULL);
 \tPy_DECREF(s);
-\tfor(i=0;f2py_routine_defs[i].name!=NULL;i++)
-\t\tPyDict_SetItemString(d, f2py_routine_defs[i].name,PyFortranObject_NewAsAttr(&f2py_routine_defs[i]));
+\tfor(i=0;f2py_routine_defs[i].name!=NULL;i++) {
+\t\ttmp = PyFortranObject_NewAsAttr(&f2py_routine_defs[i]);
+\t\tPyDict_SetItemString(d, f2py_routine_defs[i].name, tmp);
+\t\tPy_DECREF(tmp);
+\t}
 #initf2pywraphooks#
 #initf90modhooks#
 #initcommonhooks#
diff --git a/numpy/f2py/src/fortranobject.c b/numpy/f2py/src/fortranobject.c
index 4a981bf55..b55385b50 100644
--- a/numpy/f2py/src/fortranobject.c
+++ b/numpy/f2py/src/fortranobject.c
@@ -80,7 +80,10 @@ PyFortranObject_NewAsAttr(FortranDataDef* defs) { /* used for calling F90 module
     PyFortranObject *fp = NULL;
     fp = PyObject_New(PyFortranObject, &PyFortran_Type);
     if (fp == NULL) return NULL;
-    if ((fp->dict = PyDict_New())==NULL) return NULL;
+    if ((fp->dict = PyDict_New())==NULL) {
+        PyObject_Del(fp);
+        return NULL;
+    }
     fp->len = 1;
     fp->defs = defs;
     return (PyObject *)fp;
@@ -91,7 +94,7 @@ PyFortranObject_NewAsAttr(FortranDataDef* defs) { /* used for calling F90 module
 static void
 fortran_dealloc(PyFortranObject *fp) {
     Py_XDECREF(fp->dict);
-    PyMem_Del(fp);
+    PyObject_Del(fp);
 }
 
 
diff --git a/numpy/fft/pocketfft.py b/numpy/fft/pocketfft.py
index b7f6f1434..77ea6e3ba 100644
--- a/numpy/fft/pocketfft.py
+++ b/numpy/fft/pocketfft.py
@@ -271,9 +271,10 @@ def ifft(a, n=None, axis=-1, norm=None):
     a = asarray(a)
     if n is None:
         n = a.shape[axis]
-    fct = 1/n
     if norm is not None and _unitary(norm):
-        fct = 1/sqrt(n)
+        fct = 1/sqrt(max(n, 1))
+    else:
+        fct = 1/max(n, 1)
     output = _raw_fft(a, n, axis, False, False, fct)
     return output
 
diff --git a/numpy/fft/tests/test_pocketfft.py b/numpy/fft/tests/test_pocketfft.py
index db185cb21..453e964fa 100644
--- a/numpy/fft/tests/test_pocketfft.py
+++ b/numpy/fft/tests/test_pocketfft.py
@@ -4,7 +4,7 @@ import numpy as np
 import pytest
 from numpy.random import random
 from numpy.testing import (
-        assert_array_almost_equal, assert_array_equal, assert_raises,
+        assert_array_equal, assert_raises, assert_allclose
         )
 import threading
 import sys
@@ -34,109 +34,115 @@ class TestFFT1D(object):
         x = random(maxlen) + 1j*random(maxlen)
         xr = random(maxlen)
         for i in range(1,maxlen):
-            assert_array_almost_equal(np.fft.ifft(np.fft.fft(x[0:i])), x[0:i],
-                                      decimal=12)
-            assert_array_almost_equal(np.fft.irfft(np.fft.rfft(xr[0:i]),i),
-                                      xr[0:i], decimal=12)
+            assert_allclose(np.fft.ifft(np.fft.fft(x[0:i])), x[0:i],
+                            atol=1e-12)
+            assert_allclose(np.fft.irfft(np.fft.rfft(xr[0:i]),i),
+                            xr[0:i], atol=1e-12)
 
     def test_fft(self):
         x = random(30) + 1j*random(30)
-        assert_array_almost_equal(fft1(x), np.fft.fft(x))
-        assert_array_almost_equal(fft1(x) / np.sqrt(30),
-                                  np.fft.fft(x, norm="ortho"))
+        assert_allclose(fft1(x), np.fft.fft(x), atol=1e-6)
+        assert_allclose(fft1(x) / np.sqrt(30),
+                        np.fft.fft(x, norm="ortho"), atol=1e-6)
 
-    def test_ifft(self):
+    @pytest.mark.parametrize('norm', (None, 'ortho'))
+    def test_ifft(self, norm):
         x = random(30) + 1j*random(30)
-        assert_array_almost_equal(x, np.fft.ifft(np.fft.fft(x)))
-        assert_array_almost_equal(
-            x, np.fft.ifft(np.fft.fft(x, norm="ortho"), norm="ortho"))
+        assert_allclose(
+            x, np.fft.ifft(np.fft.fft(x, norm=norm), norm=norm),
+            atol=1e-6)
+        # Ensure we get the correct error message
+        with pytest.raises(ValueError,
+                           match='Invalid number of FFT data points'):
+            np.fft.ifft([], norm=norm)
 
     def test_fft2(self):
         x = random((30, 20)) + 1j*random((30, 20))
-        assert_array_almost_equal(np.fft.fft(np.fft.fft(x, axis=1), axis=0),
-                                  np.fft.fft2(x))
-        assert_array_almost_equal(np.fft.fft2(x) / np.sqrt(30 * 20),
-                                  np.fft.fft2(x, norm="ortho"))
+        assert_allclose(np.fft.fft(np.fft.fft(x, axis=1), axis=0),
+                        np.fft.fft2(x), atol=1e-6)
+        assert_allclose(np.fft.fft2(x) / np.sqrt(30 * 20),
+                        np.fft.fft2(x, norm="ortho"), atol=1e-6)
 
     def test_ifft2(self):
         x = random((30, 20)) + 1j*random((30, 20))
-        assert_array_almost_equal(np.fft.ifft(np.fft.ifft(x, axis=1), axis=0),
-                                  np.fft.ifft2(x))
-        assert_array_almost_equal(np.fft.ifft2(x) * np.sqrt(30 * 20),
-                                  np.fft.ifft2(x, norm="ortho"))
+        assert_allclose(np.fft.ifft(np.fft.ifft(x, axis=1), axis=0),
+                        np.fft.ifft2(x), atol=1e-6)
+        assert_allclose(np.fft.ifft2(x) * np.sqrt(30 * 20),
+                        np.fft.ifft2(x, norm="ortho"), atol=1e-6)
 
     def test_fftn(self):
         x = random((30, 20, 10)) + 1j*random((30, 20, 10))
-        assert_array_almost_equal(
+        assert_allclose(
             np.fft.fft(np.fft.fft(np.fft.fft(x, axis=2), axis=1), axis=0),
-            np.fft.fftn(x))
-        assert_array_almost_equal(np.fft.fftn(x) / np.sqrt(30 * 20 * 10),
-                                  np.fft.fftn(x, norm="ortho"))
+            np.fft.fftn(x), atol=1e-6)
+        assert_allclose(np.fft.fftn(x) / np.sqrt(30 * 20 * 10),
+                        np.fft.fftn(x, norm="ortho"), atol=1e-6)
 
     def test_ifftn(self):
         x = random((30, 20, 10)) + 1j*random((30, 20, 10))
-        assert_array_almost_equal(
+        assert_allclose(
             np.fft.ifft(np.fft.ifft(np.fft.ifft(x, axis=2), axis=1), axis=0),
-            np.fft.ifftn(x))
-        assert_array_almost_equal(np.fft.ifftn(x) * np.sqrt(30 * 20 * 10),
-                                  np.fft.ifftn(x, norm="ortho"))
+            np.fft.ifftn(x), atol=1e-6)
+        assert_allclose(np.fft.ifftn(x) * np.sqrt(30 * 20 * 10),
+                        np.fft.ifftn(x, norm="ortho"), atol=1e-6)
 
     def test_rfft(self):
         x = random(30)
         for n in [x.size, 2*x.size]:
             for norm in [None, 'ortho']:
-                assert_array_almost_equal(
+                assert_allclose(
                     np.fft.fft(x, n=n, norm=norm)[:(n//2 + 1)],
-                    np.fft.rfft(x, n=n, norm=norm))
-            assert_array_almost_equal(np.fft.rfft(x, n=n) / np.sqrt(n),
-                                      np.fft.rfft(x, n=n, norm="ortho"))
+                    np.fft.rfft(x, n=n, norm=norm), atol=1e-6)
+            assert_allclose(
+                np.fft.rfft(x, n=n) / np.sqrt(n),
+                np.fft.rfft(x, n=n, norm="ortho"), atol=1e-6)
 
     def test_irfft(self):
         x = random(30)
-        assert_array_almost_equal(x, np.fft.irfft(np.fft.rfft(x)))
-        assert_array_almost_equal(
-            x, np.fft.irfft(np.fft.rfft(x, norm="ortho"), norm="ortho"))
+        assert_allclose(x, np.fft.irfft(np.fft.rfft(x)), atol=1e-6)
+        assert_allclose(
+            x, np.fft.irfft(np.fft.rfft(x, norm="ortho"), norm="ortho"), atol=1e-6)
 
     def test_rfft2(self):
         x = random((30, 20))
-        assert_array_almost_equal(np.fft.fft2(x)[:, :11], np.fft.rfft2(x))
-        assert_array_almost_equal(np.fft.rfft2(x) / np.sqrt(30 * 20),
-                                  np.fft.rfft2(x, norm="ortho"))
+        assert_allclose(np.fft.fft2(x)[:, :11], np.fft.rfft2(x), atol=1e-6)
+        assert_allclose(np.fft.rfft2(x) / np.sqrt(30 * 20),
+                        np.fft.rfft2(x, norm="ortho"), atol=1e-6)
 
     def test_irfft2(self):
         x = random((30, 20))
-        assert_array_almost_equal(x, np.fft.irfft2(np.fft.rfft2(x)))
-        assert_array_almost_equal(
-            x, np.fft.irfft2(np.fft.rfft2(x, norm="ortho"), norm="ortho"))
+        assert_allclose(x, np.fft.irfft2(np.fft.rfft2(x)), atol=1e-6)
+        assert_allclose(
+            x, np.fft.irfft2(np.fft.rfft2(x, norm="ortho"), norm="ortho"), atol=1e-6)
 
     def test_rfftn(self):
         x = random((30, 20, 10))
-        assert_array_almost_equal(np.fft.fftn(x)[:, :, :6], np.fft.rfftn(x))
-        assert_array_almost_equal(np.fft.rfftn(x) / np.sqrt(30 * 20 * 10),
-                                  np.fft.rfftn(x, norm="ortho"))
+        assert_allclose(np.fft.fftn(x)[:, :, :6], np.fft.rfftn(x), atol=1e-6)
+        assert_allclose(np.fft.rfftn(x) / np.sqrt(30 * 20 * 10),
+                        np.fft.rfftn(x, norm="ortho"), atol=1e-6)
 
     def test_irfftn(self):
         x = random((30, 20, 10))
-        assert_array_almost_equal(x, np.fft.irfftn(np.fft.rfftn(x)))
-        assert_array_almost_equal(
-            x, np.fft.irfftn(np.fft.rfftn(x, norm="ortho"), norm="ortho"))
+        assert_allclose(x, np.fft.irfftn(np.fft.rfftn(x)), atol=1e-6)
+        assert_allclose(
+            x, np.fft.irfftn(np.fft.rfftn(x, norm="ortho"), norm="ortho"), atol=1e-6)
 
     def test_hfft(self):
         x = random(14) + 1j*random(14)
         x_herm = np.concatenate((random(1), x, random(1)))
         x = np.concatenate((x_herm, x[::-1].conj()))
-        assert_array_almost_equal(np.fft.fft(x), np.fft.hfft(x_herm))
-        assert_array_almost_equal(np.fft.hfft(x_herm) / np.sqrt(30),
-                                  np.fft.hfft(x_herm, norm="ortho"))
+        assert_allclose(np.fft.fft(x), np.fft.hfft(x_herm), atol=1e-6)
+        assert_allclose(np.fft.hfft(x_herm) / np.sqrt(30),
+                        np.fft.hfft(x_herm, norm="ortho"), atol=1e-6)
 
     def test_ihttf(self):
         x = random(14) + 1j*random(14)
         x_herm = np.concatenate((random(1), x, random(1)))
         x = np.concatenate((x_herm, x[::-1].conj()))
-        assert_array_almost_equal(x_herm, np.fft.ihfft(np.fft.hfft(x_herm)))
-        assert_array_almost_equal(
+        assert_allclose(x_herm, np.fft.ihfft(np.fft.hfft(x_herm)), atol=1e-6)
+        assert_allclose(
             x_herm, np.fft.ihfft(np.fft.hfft(x_herm, norm="ortho"),
-                                 norm="ortho"))
+                                 norm="ortho"), atol=1e-6)
 
     @pytest.mark.parametrize("op", [np.fft.fftn, np.fft.ifftn,
                                     np.fft.rfftn, np.fft.irfftn])
@@ -146,7 +152,7 @@ class TestFFT1D(object):
         for a in axes:
             op_tr = op(np.transpose(x, a))
             tr_op = np.transpose(op(x, axes=a), a)
-            assert_array_almost_equal(op_tr, tr_op)
+            assert_allclose(op_tr, tr_op, atol=1e-6)
 
     def test_all_1d_norm_preserving(self):
         # verify that round-trip transforms are norm-preserving
@@ -164,8 +170,8 @@ class TestFFT1D(object):
                 for norm in [None, 'ortho']:
                     tmp = forw(x, n=n, norm=norm)
                     tmp = back(tmp, n=n, norm=norm)
-                    assert_array_almost_equal(x_norm,
-                                              np.linalg.norm(tmp))
+                    assert_allclose(x_norm,
+                                    np.linalg.norm(tmp), atol=1e-6)
 
     @pytest.mark.parametrize("dtype", [np.half, np.single, np.double,
                                        np.longdouble])
@@ -173,8 +179,8 @@ class TestFFT1D(object):
         # make sure that all input precisions are accepted and internally
         # converted to 64bit
         x = random(30).astype(dtype)
-        assert_array_almost_equal(np.fft.ifft(np.fft.fft(x)), x)
-        assert_array_almost_equal(np.fft.irfft(np.fft.rfft(x)), x)
+        assert_allclose(np.fft.ifft(np.fft.fft(x)), x, atol=1e-6)
+        assert_allclose(np.fft.irfft(np.fft.rfft(x)), x, atol=1e-6)
 
 
 @pytest.mark.parametrize(
@@ -190,6 +196,8 @@ def test_fft_with_order(dtype, order, fft):
     # non contiguous arrays
     rng = np.random.RandomState(42)
     X = rng.rand(8, 7, 13).astype(dtype, copy=False)
+    # See discussion in pull/14178
+    _tol = 8.0 * np.sqrt(np.log2(X.size)) * np.finfo(X.dtype).eps
     if order == 'F':
         Y = np.asfortranarray(X)
     else:
@@ -201,7 +209,7 @@ def test_fft_with_order(dtype, order, fft):
         for axis in range(3):
             X_res = fft(X, axis=axis)
             Y_res = fft(Y, axis=axis)
-            assert_array_almost_equal(X_res, Y_res)
+            assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
     elif fft.__name__.endswith(('fft2', 'fftn')):
         axes = [(0, 1), (1, 2), (0, 2)]
         if fft.__name__.endswith('fftn'):
@@ -209,9 +217,9 @@ def test_fft_with_order(dtype, order, fft):
         for ax in axes:
             X_res = fft(X, axes=ax)
             Y_res = fft(Y, axes=ax)
-            assert_array_almost_equal(X_res, Y_res)
+            assert_allclose(X_res, Y_res, atol=_tol, rtol=_tol)
     else:
-        raise ValueError
+        raise ValueError()
 
 
 class TestFFTThreadSafe(object):
diff --git a/numpy/lib/_iotools.py b/numpy/lib/_iotools.py
index 0ebd39b8c..c392929fd 100644
--- a/numpy/lib/_iotools.py
+++ b/numpy/lib/_iotools.py
@@ -121,7 +121,7 @@ def has_nested_fields(ndtype):
 
     """
     for name in ndtype.names or ():
-        if ndtype[name].names:
+        if ndtype[name].names is not None:
             return True
     return False
 
@@ -931,28 +931,27 @@ def easy_dtype(ndtype, names=None, defaultfmt="f%i", **validationargs):
         names = validate(names, nbfields=nbfields, defaultfmt=defaultfmt)
         ndtype = np.dtype(dict(formats=ndtype, names=names))
     else:
-        nbtypes = len(ndtype)
         # Explicit names
         if names is not None:
             validate = NameValidator(**validationargs)
             if isinstance(names, basestring):
                 names = names.split(",")
             # Simple dtype: repeat to match the nb of names
-            if nbtypes == 0:
+            if ndtype.names is None:
                 formats = tuple([ndtype.type] * len(names))
                 names = validate(names, defaultfmt=defaultfmt)
                 ndtype = np.dtype(list(zip(names, formats)))
             # Structured dtype: just validate the names as needed
             else:
-                ndtype.names = validate(names, nbfields=nbtypes,
+                ndtype.names = validate(names, nbfields=len(ndtype.names),
                                         defaultfmt=defaultfmt)
         # No implicit names
-        elif (nbtypes > 0):
+        elif ndtype.names is not None:
             validate = NameValidator(**validationargs)
             # Default initial names : should we change the format ?
-            if ((ndtype.names == tuple("f%i" % i for i in range(nbtypes))) and
+            if ((ndtype.names == tuple("f%i" % i for i in range(len(ndtype.names)))) and
                     (defaultfmt != "f%i")):
-                ndtype.names = validate([''] * nbtypes, defaultfmt=defaultfmt)
+                ndtype.names = validate([''] * len(ndtype.names), defaultfmt=defaultfmt)
             # Explicit initial names : just validate
             else:
                 ndtype.names = validate(ndtype.names, defaultfmt=defaultfmt)
diff --git a/numpy/lib/arraypad.py b/numpy/lib/arraypad.py
index f08d425d6..62330e692 100644
--- a/numpy/lib/arraypad.py
+++ b/numpy/lib/arraypad.py
@@ -323,6 +323,12 @@ def _get_stats(padded, axis, width_pair, length_pair, stat_func):
     if right_length is None or max_length < right_length:
         right_length = max_length
 
+    if (left_length == 0 or right_length == 0) \
+            and stat_func in {np.amax, np.amin}:
+        # amax and amin can't operate on an emtpy array,
+        # raise a more descriptive warning here instead of the default one
+        raise ValueError("stat_length of 0 yields no value for padding")
+
     # Calculate statistic for the left side
     left_slice = _slice_at_axis(
         slice(left_index, left_index + left_length), axis)
@@ -340,6 +346,7 @@ def _get_stats(padded, axis, width_pair, length_pair, stat_func):
     right_chunk = padded[right_slice]
     right_stat = stat_func(right_chunk, axis=axis, keepdims=True)
     _round_if_needed(right_stat, padded.dtype)
+
     return left_stat, right_stat
 
 
@@ -835,7 +842,7 @@ def pad(array, pad_width, mode='constant', **kwargs):
         raise ValueError("unsupported keyword arguments for mode '{}': {}"
                          .format(mode, unsupported_kwargs))
 
-    stat_functions = {"maximum": np.max, "minimum": np.min,
+    stat_functions = {"maximum": np.amax, "minimum": np.amin,
                       "mean": np.mean, "median": np.median}
 
     # Create array with final shape and original values
diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index 9d380e67d..21532838b 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -316,14 +316,17 @@ def average(a, axis=None, weights=None, returned=False):
         The weights array can either be 1-D (in which case its length must be
         the size of `a` along the given axis) or of the same shape as `a`.
         If `weights=None`, then all data in `a` are assumed to have a
-        weight equal to one.
+        weight equal to one.  The 1-D calculation is::
+
+            avg = sum(a * weights) / sum(weights)
+
+        The only constraint on `weights` is that `sum(weights)` must not be 0.
     returned : bool, optional
         Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
         is returned, otherwise only the average is returned.
         If `weights=None`, `sum_of_weights` is equivalent to the number of
         elements over which the average is taken.
 
-
     Returns
     -------
     retval, [sum_of_weights] : array_type or double
diff --git a/numpy/lib/nanfunctions.py b/numpy/lib/nanfunctions.py
index 9a03d0b39..6cffab6ac 100644
--- a/numpy/lib/nanfunctions.py
+++ b/numpy/lib/nanfunctions.py
@@ -1443,7 +1443,7 @@ def nanvar(a, axis=None, dtype=None, out=None, ddof=0, keepdims=np._NoValue):
         the variance of the flattened array.
     dtype : data-type, optional
         Type to use in computing the variance.  For arrays of integer type
-        the default is `float32`; for arrays of float types it is the same as
+        the default is `float64`; for arrays of float types it is the same as
         the array type.
     out : ndarray, optional
         Alternate output array in which to place the result.  It must have
diff --git a/numpy/lib/npyio.py b/numpy/lib/npyio.py
index b9dc444f8..e57a6dd47 100644
--- a/numpy/lib/npyio.py
+++ b/numpy/lib/npyio.py
@@ -506,7 +506,9 @@ def save(file, arr, allow_pickle=True, fix_imports=True):
     Notes
     -----
     For a description of the ``.npy`` format, see :py:mod:`numpy.lib.format`.
-
+        
+    Any data saved to the file is appended to the end of the file. 
+    
     Examples
     --------
     >>> from tempfile import TemporaryFile
@@ -519,6 +521,15 @@ def save(file, arr, allow_pickle=True, fix_imports=True):
     >>> np.load(outfile)
     array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
 
+
+    >>> with open('test.npy', 'wb') as f:
+    ...     np.save(f, np.array([1, 2]))
+    ...     np.save(f, np.array([1, 3]))    
+    >>> with open('test.npy', 'rb') as f:
+    ...     a = np.load(f)
+    ...     b = np.load(f)
+    >>> print(a, b)
+    # [1 2] [1 3]
     """
     own_fid = False
     if hasattr(file, 'write'):
@@ -1776,12 +1787,13 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None,
                                        replace_space=replace_space)
 
         # Skip the first `skip_header` rows
-        for i in range(skip_header):
-            next(fhd)
-
-        # Keep on until we find the first valid values
-        first_values = None
         try:
+            for i in range(skip_header):
+                next(fhd)
+
+            # Keep on until we find the first valid values
+            first_values = None
+
             while not first_values:
                 first_line = _decode_line(next(fhd), encoding)
                 if (names is True) and (comments is not None):
@@ -2168,7 +2180,7 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None,
             outputmask = np.array(masks, dtype=mdtype)
     else:
         # Overwrite the initial dtype names if needed
-        if names and dtype.names:
+        if names and dtype.names is not None:
             dtype.names = names
         # Case 1. We have a structured type
         if len(dtype_flat) > 1:
@@ -2218,7 +2230,7 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None,
             #
             output = np.array(data, dtype)
             if usemask:
-                if dtype.names:
+                if dtype.names is not None:
                     mdtype = [(_, bool) for _ in dtype.names]
                 else:
                     mdtype = bool
diff --git a/numpy/lib/recfunctions.py b/numpy/lib/recfunctions.py
index 6e257bb3f..40060b41a 100644
--- a/numpy/lib/recfunctions.py
+++ b/numpy/lib/recfunctions.py
@@ -72,7 +72,7 @@ def recursive_fill_fields(input, output):
             current = input[field]
         except ValueError:
             continue
-        if current.dtype.names:
+        if current.dtype.names is not None:
             recursive_fill_fields(current, output[field])
         else:
             output[field][:len(current)] = current
@@ -139,11 +139,11 @@ def get_names(adtype):
     names = adtype.names
     for name in names:
         current = adtype[name]
-        if current.names:
+        if current.names is not None:
             listnames.append((name, tuple(get_names(current))))
         else:
             listnames.append(name)
-    return tuple(listnames) or None
+    return tuple(listnames)
 
 
 def get_names_flat(adtype):
@@ -176,9 +176,9 @@ def get_names_flat(adtype):
     for name in names:
         listnames.append(name)
         current = adtype[name]
-        if current.names:
+        if current.names is not None:
             listnames.extend(get_names_flat(current))
-    return tuple(listnames) or None
+    return tuple(listnames)
 
 
 def flatten_descr(ndtype):
@@ -215,8 +215,8 @@ def _zip_dtype(seqarrays, flatten=False):
     else:
         for a in seqarrays:
             current = a.dtype
-            if current.names and len(current.names) <= 1:
-                # special case - dtypes of 0 or 1 field are flattened
+            if current.names is not None and len(current.names) == 1:
+                # special case - dtypes of 1 field are flattened
                 newdtype.extend(_get_fieldspec(current))
             else:
                 newdtype.append(('', current))
@@ -268,7 +268,7 @@ def get_fieldstructure(adtype, lastname=None, parents=None,):
     names = adtype.names
     for name in names:
         current = adtype[name]
-        if current.names:
+        if current.names is not None:
             if lastname:
                 parents[name] = [lastname, ]
             else:
@@ -281,7 +281,7 @@ def get_fieldstructure(adtype, lastname=None, parents=None,):
             elif lastname:
                 lastparent = [lastname, ]
             parents[name] = lastparent or []
-    return parents or None
+    return parents
 
 
 def _izip_fields_flat(iterable):
@@ -435,7 +435,7 @@ def merge_arrays(seqarrays, fill_value=-1, flatten=False,
     if isinstance(seqarrays, (ndarray, np.void)):
         seqdtype = seqarrays.dtype
         # Make sure we have named fields
-        if not seqdtype.names:
+        if seqdtype.names is None:
             seqdtype = np.dtype([('', seqdtype)])
         if not flatten or _zip_dtype((seqarrays,), flatten=True) == seqdtype:
             # Minimal processing needed: just make sure everythng's a-ok
@@ -653,7 +653,7 @@ def rename_fields(base, namemapper):
         for name in ndtype.names:
             newname = namemapper.get(name, name)
             current = ndtype[name]
-            if current.names:
+            if current.names is not None:
                 newdtype.append(
                     (newname, _recursive_rename_fields(current, namemapper))
                     )
@@ -874,16 +874,35 @@ def _get_fields_and_offsets(dt, offset=0):
     scalar fields in the dtype "dt", including nested fields, in left
     to right order.
     """
+
+    # counts up elements in subarrays, including nested subarrays, and returns
+    # base dtype and count
+    def count_elem(dt):
+        count = 1
+        while dt.shape != ():
+            for size in dt.shape:
+                count *= size
+            dt = dt.base
+        return dt, count
+
     fields = []
     for name in dt.names:
         field = dt.fields[name]
-        if field[0].names is None:
-            count = 1
-            for size in field[0].shape:
-                count *= size
-            fields.append((field[0], count, field[1] + offset))
+        f_dt, f_offset = field[0], field[1]
+        f_dt, n = count_elem(f_dt)
+
+        if f_dt.names is None:
+            fields.append((np.dtype((f_dt, (n,))), n, f_offset + offset))
         else:
-            fields.extend(_get_fields_and_offsets(field[0], field[1] + offset))
+            subfields = _get_fields_and_offsets(f_dt, f_offset + offset)
+            size = f_dt.itemsize
+
+            for i in range(n):
+                if i == 0:
+                    # optimization: avoid list comprehension if no subarray
+                    fields.extend(subfields)
+                else:
+                    fields.extend([(d, c, o + i*size) for d, c, o in subfields])
     return fields
 
 
@@ -948,6 +967,12 @@ def structured_to_unstructured(arr, dtype=None, copy=False, casting='unsafe'):
 
     fields = _get_fields_and_offsets(arr.dtype)
     n_fields = len(fields)
+    if n_fields == 0 and dtype is None:
+        raise ValueError("arr has no fields. Unable to guess dtype")
+    elif n_fields == 0:
+        # too many bugs elsewhere for this to work now
+        raise NotImplementedError("arr with no fields is not supported")
+
     dts, counts, offsets = zip(*fields)
     names = ['f{}'.format(n) for n in range(n_fields)]
 
@@ -1039,6 +1064,9 @@ def unstructured_to_structured(arr, dtype=None, names=None, align=False,
     if arr.shape == ():
         raise ValueError('arr must have at least one dimension')
     n_elem = arr.shape[-1]
+    if n_elem == 0:
+        # too many bugs elsewhere for this to work now
+        raise NotImplementedError("last axis with size 0 is not supported")
 
     if dtype is None:
         if names is None:
@@ -1051,7 +1079,11 @@ def unstructured_to_structured(arr, dtype=None, names=None, align=False,
             raise ValueError("don't supply both dtype and names")
         # sanity check of the input dtype
         fields = _get_fields_and_offsets(dtype)
-        dts, counts, offsets = zip(*fields)
+        if len(fields) == 0:
+            dts, counts, offsets = [], [], []
+        else:
+            dts, counts, offsets = zip(*fields)
+
         if n_elem != sum(counts):
             raise ValueError('The length of the last dimension of arr must '
                              'be equal to the number of fields in dtype')
diff --git a/numpy/lib/tests/test_arraypad.py b/numpy/lib/tests/test_arraypad.py
index b7630cdcd..b6dd3b31c 100644
--- a/numpy/lib/tests/test_arraypad.py
+++ b/numpy/lib/tests/test_arraypad.py
@@ -469,6 +469,29 @@ class TestStatistic(object):
             )
         assert_array_equal(a, b)
 
+    @pytest.mark.filterwarnings("ignore:Mean of empty slice:RuntimeWarning")
+    @pytest.mark.filterwarnings(
+        "ignore:invalid value encountered in (true_divide|double_scalars):"
+        "RuntimeWarning"
+    )
+    @pytest.mark.parametrize("mode", ["mean", "median"])
+    def test_zero_stat_length_valid(self, mode):
+        arr = np.pad([1., 2.], (1, 2), mode, stat_length=0)
+        expected = np.array([np.nan, 1., 2., np.nan, np.nan])
+        assert_equal(arr, expected)
+
+    @pytest.mark.parametrize("mode", ["minimum", "maximum"])
+    def test_zero_stat_length_invalid(self, mode):
+        match = "stat_length of 0 yields no value for padding"
+        with pytest.raises(ValueError, match=match):
+            np.pad([1., 2.], 0, mode, stat_length=0)
+        with pytest.raises(ValueError, match=match):
+            np.pad([1., 2.], 0, mode, stat_length=(1, 0))
+        with pytest.raises(ValueError, match=match):
+            np.pad([1., 2.], 1, mode, stat_length=0)
+        with pytest.raises(ValueError, match=match):
+            np.pad([1., 2.], 1, mode, stat_length=(1, 0))
+
 
 class TestConstant(object):
     def test_check_constant(self):
diff --git a/numpy/lib/tests/test_io.py b/numpy/lib/tests/test_io.py
index 78f9f85f3..6ee17c830 100644
--- a/numpy/lib/tests/test_io.py
+++ b/numpy/lib/tests/test_io.py
@@ -1565,6 +1565,13 @@ M   33  21.99
             test = np.genfromtxt(TextIO(data), delimiter=";",
                                  dtype=ndtype, converters=converters)
 
+        # nested but empty fields also aren't supported
+        ndtype = [('idx', int), ('code', object), ('nest', [])]
+        with assert_raises_regex(NotImplementedError,
+                                 'Nested fields.* not supported.*'):
+            test = np.genfromtxt(TextIO(data), delimiter=";",
+                                 dtype=ndtype, converters=converters)
+
     def test_userconverters_with_explicit_dtype(self):
         # Test user_converters w/ explicit (standard) dtype
         data = TextIO('skip,skip,2001-01-01,1.0,skip')
@@ -1681,6 +1688,10 @@ M   33  21.99
             test = np.genfromtxt(data)
             assert_equal(test, np.array([]))
 
+            # when skip_header > 0
+            test = np.genfromtxt(data, skip_header=1)
+            assert_equal(test, np.array([]))
+
     def test_fancy_dtype_alt(self):
         # Check that a nested dtype isn't MIA
         data = TextIO('1,2,3.0\n4,5,6.0\n')
diff --git a/numpy/lib/tests/test_recfunctions.py b/numpy/lib/tests/test_recfunctions.py
index 0126ccaf8..0c839d486 100644
--- a/numpy/lib/tests/test_recfunctions.py
+++ b/numpy/lib/tests/test_recfunctions.py
@@ -115,6 +115,14 @@ class TestRecFunctions(object):
         test = get_names(ndtype)
         assert_equal(test, ('a', ('b', ('ba', 'bb'))))
 
+        ndtype = np.dtype([('a', int), ('b', [])])
+        test = get_names(ndtype)
+        assert_equal(test, ('a', ('b', ())))
+
+        ndtype = np.dtype([])
+        test = get_names(ndtype)
+        assert_equal(test, ())
+
     def test_get_names_flat(self):
         # Test get_names_flat
         ndtype = np.dtype([('A', '|S3'), ('B', float)])
@@ -125,6 +133,14 @@ class TestRecFunctions(object):
         test = get_names_flat(ndtype)
         assert_equal(test, ('a', 'b', 'ba', 'bb'))
 
+        ndtype = np.dtype([('a', int), ('b', [])])
+        test = get_names_flat(ndtype)
+        assert_equal(test, ('a', 'b'))
+
+        ndtype = np.dtype([])
+        test = get_names_flat(ndtype)
+        assert_equal(test, ())
+
     def test_get_fieldstructure(self):
         # Test get_fieldstructure
 
@@ -147,6 +163,11 @@ class TestRecFunctions(object):
                    'BBA': ['B', 'BB'], 'BBB': ['B', 'BB']}
         assert_equal(test, control)
 
+        # 0 fields
+        ndtype = np.dtype([])
+        test = get_fieldstructure(ndtype)
+        assert_equal(test, {})
+
     def test_find_duplicates(self):
         # Test find_duplicates
         a = ma.array([(2, (2., 'B')), (1, (2., 'B')), (2, (2., 'B')),
@@ -248,7 +269,8 @@ class TestRecFunctions(object):
         # including uniform fields with subarrays unpacked
         d = np.array([(1, [2,  3], [[ 4,  5], [ 6,  7]]),
                       (8, [9, 10], [[11, 12], [13, 14]])],
-                     dtype=[('x0', 'i4'), ('x1', ('i4', 2)), ('x2', ('i4', (2, 2)))])
+                     dtype=[('x0', 'i4'), ('x1', ('i4', 2)),
+                            ('x2', ('i4', (2, 2)))])
         dd = structured_to_unstructured(d)
         ddd = unstructured_to_structured(dd, d.dtype)
         assert_(dd.base is d)
@@ -262,6 +284,40 @@ class TestRecFunctions(object):
         assert_equal(res, np.zeros((10, 6), dtype=int))
 
 
+        # test nested combinations of subarrays and structured arrays, gh-13333
+        def subarray(dt, shape):
+            return np.dtype((dt, shape))
+
+        def structured(*dts):
+            return np.dtype([('x{}'.format(i), dt) for i, dt in enumerate(dts)])
+
+        def inspect(dt, dtype=None):
+            arr = np.zeros((), dt)
+            ret = structured_to_unstructured(arr, dtype=dtype)
+            backarr = unstructured_to_structured(ret, dt)
+            return ret.shape, ret.dtype, backarr.dtype
+
+        dt = structured(subarray(structured(np.int32, np.int32), 3))
+        assert_equal(inspect(dt), ((6,), np.int32, dt))
+
+        dt = structured(subarray(subarray(np.int32, 2), 2))
+        assert_equal(inspect(dt), ((4,), np.int32, dt))
+
+        dt = structured(np.int32)
+        assert_equal(inspect(dt), ((1,), np.int32, dt))
+
+        dt = structured(np.int32, subarray(subarray(np.int32, 2), 2))
+        assert_equal(inspect(dt), ((5,), np.int32, dt))
+
+        dt = structured()
+        assert_raises(ValueError, structured_to_unstructured, np.zeros(3, dt))
+
+        # these currently don't work, but we may make it work in the future
+        assert_raises(NotImplementedError, structured_to_unstructured,
+                                           np.zeros(3, dt), dtype=np.int32)
+        assert_raises(NotImplementedError, unstructured_to_structured,
+                                           np.zeros((3,0), dtype=np.int32))
+
     def test_field_assignment_by_name(self):
         a = np.ones(2, dtype=[('a', 'i4'), ('b', 'f8'), ('c', 'u1')])
         newdt = [('b', 'f4'), ('c', 'u1')]
diff --git a/numpy/lib/type_check.py b/numpy/lib/type_check.py
index ac4b03a6c..586824743 100644
--- a/numpy/lib/type_check.py
+++ b/numpy/lib/type_check.py
@@ -395,19 +395,27 @@ def nan_to_num(x, copy=True, nan=0.0, posinf=None, neginf=None):
         in-place (False). The in-place operation only occurs if
         casting to an array does not require a copy.
         Default is True.
+        
+        .. versionadded:: 1.13
     nan : int, float, optional
         Value to be used to fill NaN values. If no value is passed 
         then NaN values will be replaced with 0.0.
+        
+        .. versionadded:: 1.17
     posinf : int, float, optional
         Value to be used to fill positive infinity values. If no value is 
         passed then positive infinity values will be replaced with a very
         large number.
+        
+        .. versionadded:: 1.17
     neginf : int, float, optional
         Value to be used to fill negative infinity values. If no value is 
         passed then negative infinity values will be replaced with a very
         small (or negative) number.
+        
+        .. versionadded:: 1.17
 
-        .. versionadded:: 1.13
+        
 
     Returns
     -------
diff --git a/numpy/lib/utils.py b/numpy/lib/utils.py
index c7dbcc5f9..8bcbd8e86 100644
--- a/numpy/lib/utils.py
+++ b/numpy/lib/utils.py
@@ -1003,93 +1003,6 @@ def _getmembers(item):
                    if hasattr(item, x)]
     return members
 
-#-----------------------------------------------------------------------------
-
-# The following SafeEval class and company are adapted from Michael Spencer's
-# ASPN Python Cookbook recipe: https://code.activestate.com/recipes/364469/
-#
-# Accordingly it is mostly Copyright 2006 by Michael Spencer.
-# The recipe, like most of the other ASPN Python Cookbook recipes was made
-# available under the Python license.
-#   https://en.wikipedia.org/wiki/Python_License
-
-# It has been modified to:
-#   * handle unary -/+
-#   * support True/False/None
-#   * raise SyntaxError instead of a custom exception.
-
-class SafeEval(object):
-    """
-    Object to evaluate constant string expressions.
-
-    This includes strings with lists, dicts and tuples using the abstract
-    syntax tree created by ``compiler.parse``.
-
-    .. deprecated:: 1.10.0
-
-    See Also
-    --------
-    safe_eval
-
-    """
-    def __init__(self):
-        # 2014-10-15, 1.10
-        warnings.warn("SafeEval is deprecated in 1.10 and will be removed.",
-                      DeprecationWarning, stacklevel=2)
-
-    def visit(self, node):
-        cls = node.__class__
-        meth = getattr(self, 'visit' + cls.__name__, self.default)
-        return meth(node)
-
-    def default(self, node):
-        raise SyntaxError("Unsupported source construct: %s"
-                          % node.__class__)
-
-    def visitExpression(self, node):
-        return self.visit(node.body)
-
-    def visitNum(self, node):
-        return node.n
-
-    def visitStr(self, node):
-        return node.s
-
-    def visitBytes(self, node):
-        return node.s
-
-    def visitDict(self, node,**kw):
-        return dict([(self.visit(k), self.visit(v))
-                     for k, v in zip(node.keys, node.values)])
-
-    def visitTuple(self, node):
-        return tuple([self.visit(i) for i in node.elts])
-
-    def visitList(self, node):
-        return [self.visit(i) for i in node.elts]
-
-    def visitUnaryOp(self, node):
-        import ast
-        if isinstance(node.op, ast.UAdd):
-            return +self.visit(node.operand)
-        elif isinstance(node.op, ast.USub):
-            return -self.visit(node.operand)
-        else:
-            raise SyntaxError("Unknown unary op: %r" % node.op)
-
-    def visitName(self, node):
-        if node.id == 'False':
-            return False
-        elif node.id == 'True':
-            return True
-        elif node.id == 'None':
-            return None
-        else:
-            raise SyntaxError("Unknown name: %s" % node.id)
-
-    def visitNameConstant(self, node):
-        return node.value
-
 
 def safe_eval(source):
     """
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index 2a3ff0728..816a200eb 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -1481,6 +1481,12 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
     compute_uv : bool, optional
         Whether or not to compute `u` and `vh` in addition to `s`.  True
         by default.
+    hermitian : bool, optional
+        If True, `a` is assumed to be Hermitian (symmetric if real-valued),
+        enabling a more efficient method for finding singular values.
+        Defaults to False.
+
+        .. versionadded:: 1.17.0
 
     Returns
     -------
@@ -1498,12 +1504,6 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
         size as those of the input `a`. The size of the last two dimensions
         depends on the value of `full_matrices`. Only returned when
         `compute_uv` is True.
-    hermitian : bool, optional
-        If True, `a` is assumed to be Hermitian (symmetric if real-valued),
-        enabling a more efficient method for finding singular values.
-        Defaults to False.
-
-        .. versionadded:: 1.17.0
 
     Raises
     ------
diff --git a/numpy/ma/core.py b/numpy/ma/core.py
index 796677139..bb3788c9a 100644
--- a/numpy/ma/core.py
+++ b/numpy/ma/core.py
@@ -59,14 +59,14 @@ __all__ = [
     'choose', 'clip', 'common_fill_value', 'compress', 'compressed',
     'concatenate', 'conjugate', 'convolve', 'copy', 'correlate', 'cos', 'cosh',
     'count', 'cumprod', 'cumsum', 'default_fill_value', 'diag', 'diagonal',
-    'diff', 'divide', 'dump', 'dumps', 'empty', 'empty_like', 'equal', 'exp',
+    'diff', 'divide', 'empty', 'empty_like', 'equal', 'exp',
     'expand_dims', 'fabs', 'filled', 'fix_invalid', 'flatten_mask',
     'flatten_structured_array', 'floor', 'floor_divide', 'fmod',
     'frombuffer', 'fromflex', 'fromfunction', 'getdata', 'getmask',
     'getmaskarray', 'greater', 'greater_equal', 'harden_mask', 'hypot',
     'identity', 'ids', 'indices', 'inner', 'innerproduct', 'isMA',
     'isMaskedArray', 'is_mask', 'is_masked', 'isarray', 'left_shift',
-    'less', 'less_equal', 'load', 'loads', 'log', 'log10', 'log2',
+    'less', 'less_equal', 'log', 'log10', 'log2',
     'logical_and', 'logical_not', 'logical_or', 'logical_xor', 'make_mask',
     'make_mask_descr', 'make_mask_none', 'mask_or', 'masked',
     'masked_array', 'masked_equal', 'masked_greater',
@@ -77,7 +77,7 @@ __all__ = [
     'maximum_fill_value', 'mean', 'min', 'minimum', 'minimum_fill_value',
     'mod', 'multiply', 'mvoid', 'ndim', 'negative', 'nomask', 'nonzero',
     'not_equal', 'ones', 'outer', 'outerproduct', 'power', 'prod',
-    'product', 'ptp', 'put', 'putmask', 'rank', 'ravel', 'remainder',
+    'product', 'ptp', 'put', 'putmask', 'ravel', 'remainder',
     'repeat', 'reshape', 'resize', 'right_shift', 'round', 'round_',
     'set_fill_value', 'shape', 'sin', 'sinh', 'size', 'soften_mask',
     'sometrue', 'sort', 'sqrt', 'squeeze', 'std', 'subtract', 'sum',
@@ -5887,7 +5887,6 @@ class MaskedArray(ndarray):
         return out[()]
 
     # Array methods
-    clip = _arraymethod('clip', onmask=False)
     copy = _arraymethod('copy')
     diagonal = _arraymethod('diagonal')
     flatten = _arraymethod('flatten')
@@ -7099,23 +7098,6 @@ def resize(x, new_shape):
     return result
 
 
-def rank(obj):
-    """
-    maskedarray version of the numpy function.
-
-    .. note::
-        Deprecated since 1.10.0
-
-    """
-    # 2015-04-12, 1.10.0
-    warnings.warn(
-        "`rank` is deprecated; use the `ndim` function instead. ",
-        np.VisibleDeprecationWarning, stacklevel=2)
-    return np.ndim(getdata(obj))
-
-rank.__doc__ = np.rank.__doc__
-
-
 def ndim(obj):
     """
     maskedarray version of the numpy function.
@@ -7904,93 +7886,6 @@ def _pickle_warn(method):
         stacklevel=3)
 
 
-def dump(a, F):
-    """
-    Pickle a masked array to a file.
-
-    This is a wrapper around ``cPickle.dump``.
-
-    Parameters
-    ----------
-    a : MaskedArray
-        The array to be pickled.
-    F : str or file-like object
-        The file to pickle `a` to. If a string, the full path to the file.
-
-    """
-    _pickle_warn('dump')
-    if not hasattr(F, 'readline'):
-        with open(F, 'w') as F:
-            pickle.dump(a, F)
-    else:
-        pickle.dump(a, F)
-
-
-def dumps(a):
-    """
-    Return a string corresponding to the pickling of a masked array.
-
-    This is a wrapper around ``cPickle.dumps``.
-
-    Parameters
-    ----------
-    a : MaskedArray
-        The array for which the string representation of the pickle is
-        returned.
-
-    """
-    _pickle_warn('dumps')
-    return pickle.dumps(a)
-
-
-def load(F):
-    """
-    Wrapper around ``cPickle.load`` which accepts either a file-like object
-    or a filename.
-
-    Parameters
-    ----------
-    F : str or file
-        The file or file name to load.
-
-    See Also
-    --------
-    dump : Pickle an array
-
-    Notes
-    -----
-    This is different from `numpy.load`, which does not use cPickle but loads
-    the NumPy binary .npy format.
-
-    """
-    _pickle_warn('load')
-    if not hasattr(F, 'readline'):
-        with open(F, 'r') as F:
-            return pickle.load(F)
-    else:
-        return pickle.load(F)
-
-
-def loads(strg):
-    """
-    Load a pickle from the current string.
-
-    The result of ``cPickle.loads(strg)`` is returned.
-
-    Parameters
-    ----------
-    strg : str
-        The string to load.
-
-    See Also
-    --------
-    dumps : Return a string corresponding to the pickling of a masked array.
-
-    """
-    _pickle_warn('loads')
-    return pickle.loads(strg)
-
-
 def fromfile(file, dtype=float, count=-1, sep=''):
     raise NotImplementedError(
         "fromfile() not yet implemented for a MaskedArray.")
diff --git a/numpy/ma/extras.py b/numpy/ma/extras.py
index 639b3dd1f..de1aa3af8 100644
--- a/numpy/ma/extras.py
+++ b/numpy/ma/extras.py
@@ -549,8 +549,11 @@ def average(a, axis=None, weights=None, returned=False):
         The weights array can either be 1-D (in which case its length must be
         the size of `a` along the given axis) or of the same shape as `a`.
         If ``weights=None``, then all data in `a` are assumed to have a
-        weight equal to one.   If `weights` is complex, the imaginary parts
-        are ignored.
+        weight equal to one.  The 1-D calculation is::
+
+            avg = sum(a * weights) / sum(weights)
+
+        The only constraint on `weights` is that `sum(weights)` must not be 0.
     returned : bool, optional
         Flag indicating whether a tuple ``(result, sum of weights)``
         should be returned as output (True), or just the result (False).
diff --git a/numpy/ma/mrecords.py b/numpy/ma/mrecords.py
index 931a7e8b9..826fb0f64 100644
--- a/numpy/ma/mrecords.py
+++ b/numpy/ma/mrecords.py
@@ -208,7 +208,7 @@ class MaskedRecords(MaskedArray, object):
         _localdict = ndarray.__getattribute__(self, '__dict__')
         _data = ndarray.view(self, _localdict['_baseclass'])
         obj = _data.getfield(*res)
-        if obj.dtype.fields:
+        if obj.dtype.names is not None:
             raise NotImplementedError("MaskedRecords is currently limited to"
                                       "simple records.")
         # Get some special attributes
diff --git a/numpy/ma/tests/test_core.py b/numpy/ma/tests/test_core.py
index 9fe550ef8..b72ce56aa 100644
--- a/numpy/ma/tests/test_core.py
+++ b/numpy/ma/tests/test_core.py
@@ -3035,6 +3035,13 @@ class TestMaskedArrayMethods(object):
         assert_equal(clipped._data, x.clip(2, 8))
         assert_equal(clipped._data, mx._data.clip(2, 8))
 
+    def test_clip_out(self):
+        # gh-14140
+        a = np.arange(10)
+        m = np.ma.MaskedArray(a, mask=[0, 1] * 5)
+        m.clip(0, 5, out=m)
+        assert_equal(m.mask, [0, 1] * 5)
+
     def test_compress(self):
         # test compress
         a = masked_array([1., 2., 3., 4., 5.], fill_value=9999)
diff --git a/numpy/random/__init__.py b/numpy/random/__init__.py
index 2d495d67e..e7eecc5cd 100644
--- a/numpy/random/__init__.py
+++ b/numpy/random/__init__.py
@@ -177,7 +177,12 @@ __all__ = [
     'zipf',
 ]
 
-from . import mtrand
+# add these for module-freeze analysis (like PyInstaller)
+from . import _pickle
+from . import common
+from . import bounded_integers
+from . import entropy
+
 from .mtrand import *
 from .generator import Generator, default_rng
 from .bit_generator import SeedSequence
diff --git a/numpy/random/_pickle.py b/numpy/random/_pickle.py
index d20a91ced..3b58f21e8 100644
--- a/numpy/random/_pickle.py
+++ b/numpy/random/_pickle.py
@@ -13,7 +13,7 @@ BitGenerators = {'MT19937': MT19937,
                  }
 
 
-def __generator_ctor(bit_generator_name='mt19937'):
+def __generator_ctor(bit_generator_name='MT19937'):
     """
     Pickling helper function that returns a Generator object
 
@@ -36,7 +36,7 @@ def __generator_ctor(bit_generator_name='mt19937'):
     return Generator(bit_generator())
 
 
-def __bit_generator_ctor(bit_generator_name='mt19937'):
+def __bit_generator_ctor(bit_generator_name='MT19937'):
     """
     Pickling helper function that returns a bit generator object
 
@@ -59,7 +59,7 @@ def __bit_generator_ctor(bit_generator_name='mt19937'):
     return bit_generator()
 
 
-def __randomstate_ctor(bit_generator_name='mt19937'):
+def __randomstate_ctor(bit_generator_name='MT19937'):
     """
     Pickling helper function that returns a legacy RandomState-like object
 
diff --git a/numpy/random/bit_generator.pxd b/numpy/random/bit_generator.pxd
index 79fe69275..984033f17 100644
--- a/numpy/random/bit_generator.pxd
+++ b/numpy/random/bit_generator.pxd
@@ -1,5 +1,5 @@
 
-from .common cimport bitgen_t
+from .common cimport bitgen_t, uint32_t
 cimport numpy as np
 
 cdef class BitGenerator():
@@ -14,9 +14,9 @@ cdef class BitGenerator():
 cdef class SeedSequence():
     cdef readonly object entropy
     cdef readonly tuple spawn_key
-    cdef readonly int pool_size
+    cdef readonly uint32_t pool_size
     cdef readonly object pool
-    cdef readonly int n_children_spawned
+    cdef readonly uint32_t n_children_spawned
 
     cdef mix_entropy(self, np.ndarray[np.npy_uint32, ndim=1] mixer,
                      np.ndarray[np.npy_uint32, ndim=1] entropy_array)
diff --git a/numpy/random/bit_generator.pyx b/numpy/random/bit_generator.pyx
index 6694e5e4d..eb608af6c 100644
--- a/numpy/random/bit_generator.pyx
+++ b/numpy/random/bit_generator.pyx
@@ -116,7 +116,7 @@ def _coerce_to_uint32_array(x):
     Examples
     --------
     >>> import numpy as np
-    >>> from np.random.bit_generator import _coerce_to_uint32_array
+    >>> from numpy.random.bit_generator import _coerce_to_uint32_array
     >>> _coerce_to_uint32_array(12345)
     array([12345], dtype=uint32)
     >>> _coerce_to_uint32_array('12345')
@@ -458,6 +458,8 @@ cdef class SeedSequence():
         -------
         seqs : list of `SeedSequence` s
         """
+        cdef uint32_t i
+
         seqs = []
         for i in range(self.n_children_spawned,
                        self.n_children_spawned + n_children):
diff --git a/numpy/random/common.pyx b/numpy/random/common.pyx
index 6ad5f5b21..74cd5f033 100644
--- a/numpy/random/common.pyx
+++ b/numpy/random/common.pyx
@@ -227,7 +227,7 @@ cdef check_output(object out, object dtype, object size):
         raise ValueError('Supplied output array is not contiguous, writable or aligned.')
     if out_array.dtype != dtype:
         raise TypeError('Supplied output array has the wrong type. '
-                        'Expected {0}, got {0}'.format(dtype, out_array.dtype))
+                        'Expected {0}, got {1}'.format(np.dtype(dtype), out_array.dtype))
     if size is not None:
         try:
             tup_size = tuple(size)
diff --git a/numpy/random/generator.pyx b/numpy/random/generator.pyx
index c7432d8c1..26fd95129 100644
--- a/numpy/random/generator.pyx
+++ b/numpy/random/generator.pyx
@@ -3919,9 +3919,8 @@ cdef class Generator:
         permutation(x)
 
         Randomly permute a sequence, or return a permuted range.
-
         If `x` is a multi-dimensional array, it is only shuffled along its
-        first index.
+        first index. 
 
         Parameters
         ----------
@@ -3950,13 +3949,20 @@ cdef class Generator:
                [0, 1, 2],
                [3, 4, 5]])
 
+        >>> rng.permutation("abc")
+        Traceback (most recent call last):
+            ...
+        numpy.AxisError: x must be an integer or at least 1-dimensional
         """
+
         if isinstance(x, (int, np.integer)):
             arr = np.arange(x)
             self.shuffle(arr)
             return arr
 
         arr = np.asarray(x)
+        if arr.ndim < 1:
+            raise np.AxisError("x must be an integer or at least 1-dimensional")
 
         # shuffle has fast-path for 1-d
         if arr.ndim == 1:
diff --git a/numpy/random/mtrand.pyx b/numpy/random/mtrand.pyx
index 46b6b3388..eb263cd2d 100644
--- a/numpy/random/mtrand.pyx
+++ b/numpy/random/mtrand.pyx
@@ -4134,6 +4134,7 @@ cdef class RandomState:
         out : ndarray
             Permuted sequence or array range.
 
+
         Examples
         --------
         >>> np.random.permutation(10)
@@ -4149,12 +4150,15 @@ cdef class RandomState:
                [3, 4, 5]])
 
         """
+
         if isinstance(x, (int, np.integer)):
             arr = np.arange(x)
             self.shuffle(arr)
             return arr
 
         arr = np.asarray(x)
+        if arr.ndim < 1:
+            raise IndexError("x must be an integer or at least 1-dimensional")
 
         # shuffle has fast-path for 1-d
         if arr.ndim == 1:
diff --git a/numpy/random/setup.py b/numpy/random/setup.py
index a1bf3b83c..a820d326e 100644
--- a/numpy/random/setup.py
+++ b/numpy/random/setup.py
@@ -49,8 +49,8 @@ def configuration(parent_package='', top_path=None):
     elif not is_msvc:
         # Some bit generators require c99
         EXTRA_COMPILE_ARGS += ['-std=c99']
-        INTEL_LIKE = any([val in k.lower() for k in platform.uname()
-                          for val in ('x86', 'i686', 'i386', 'amd64')])
+        INTEL_LIKE = any(arch in platform.machine() 
+                         for arch in ('x86', 'i686', 'i386', 'amd64'))
         if INTEL_LIKE:
             # Assumes GCC or GCC-like compiler
             EXTRA_COMPILE_ARGS += ['-msse2']
diff --git a/numpy/random/tests/test_generator_mt19937.py b/numpy/random/tests/test_generator_mt19937.py
index a962fe84e..853d86fba 100644
--- a/numpy/random/tests/test_generator_mt19937.py
+++ b/numpy/random/tests/test_generator_mt19937.py
@@ -757,6 +757,19 @@ class TestRandomDist(object):
         arr_2d = np.atleast_2d([1, 2, 3, 4, 5, 6, 7, 8, 9, 0]).T
         actual = random.permutation(arr_2d)
         assert_array_equal(actual, np.atleast_2d(desired).T)
+        
+        bad_x_str = "abcd"
+        assert_raises(np.AxisError, random.permutation, bad_x_str)
+
+        bad_x_float = 1.2
+        assert_raises(np.AxisError, random.permutation, bad_x_float)
+
+        random = Generator(MT19937(self.seed))
+        integer_val = 10
+        desired = [3, 0, 8, 7, 9, 4, 2, 5, 1, 6]
+
+        actual = random.permutation(integer_val)
+        assert_array_equal(actual, desired)
 
     def test_beta(self):
         random = Generator(MT19937(self.seed))
diff --git a/numpy/random/tests/test_randomstate.py b/numpy/random/tests/test_randomstate.py
index 3b5a279a3..a0edc5c23 100644
--- a/numpy/random/tests/test_randomstate.py
+++ b/numpy/random/tests/test_randomstate.py
@@ -686,6 +686,21 @@ class TestRandomDist(object):
         actual = random.permutation(arr_2d)
         assert_array_equal(actual, np.atleast_2d(desired).T)
 
+        random.seed(self.seed)
+        bad_x_str = "abcd"
+        assert_raises(IndexError, random.permutation, bad_x_str)
+
+        random.seed(self.seed)
+        bad_x_float = 1.2
+        assert_raises(IndexError, random.permutation, bad_x_float)
+
+        integer_val = 10
+        desired = [9, 0, 8, 5, 1, 3, 4, 7, 6, 2]
+
+        random.seed(self.seed)
+        actual = random.permutation(integer_val)
+        assert_array_equal(actual, desired)
+
     def test_beta(self):
         random.seed(self.seed)
         actual = random.beta(.1, .9, size=(3, 2))
diff --git a/numpy/testing/_private/utils.py b/numpy/testing/_private/utils.py
index 87e66e06f..97a5eac17 100644
--- a/numpy/testing/_private/utils.py
+++ b/numpy/testing/_private/utils.py
@@ -21,7 +21,6 @@ import pprint
 
 from numpy.core import(
      intp, float32, empty, arange, array_repr, ndarray, isnat, array)
-from numpy.lib.utils import deprecate
 
 if sys.version_info[0] >= 3:
     from io import StringIO
@@ -33,7 +32,7 @@ __all__ = [
         'assert_array_equal', 'assert_array_less', 'assert_string_equal',
         'assert_array_almost_equal', 'assert_raises', 'build_err_msg',
         'decorate_methods', 'jiffies', 'memusage', 'print_assert_equal',
-        'raises', 'rand', 'rundocs', 'runstring', 'verbose', 'measure',
+        'raises', 'rundocs', 'runstring', 'verbose', 'measure',
         'assert_', 'assert_array_almost_equal_nulp', 'assert_raises_regex',
         'assert_array_max_ulp', 'assert_warns', 'assert_no_warnings',
         'assert_allclose', 'IgnoreException', 'clear_and_catch_warnings',
@@ -154,22 +153,6 @@ def gisinf(x):
     return st
 
 
-@deprecate(message="numpy.testing.rand is deprecated in numpy 1.11. "
-                   "Use numpy.random.rand instead.")
-def rand(*args):
-    """Returns an array of random numbers with the given shape.
-
-    This only uses the standard library, so it is useful for testing purposes.
-    """
-    import random
-    from numpy.core import zeros, float64
-    results = zeros(args, float64)
-    f = results.flat
-    for i in range(len(f)):
-        f[i] = random.random()
-    return results
-
-
 if os.name == 'nt':
     # Code "stolen" from enthought/debug/memusage.py
     def GetPerformanceAttributes(object, counter, instance=None,
@@ -804,8 +787,11 @@ def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
         # do not trigger a failure (np.ma.masked != True evaluates as
         # np.ma.masked, which is falsy).
         if cond != True:
-            mismatch = 100. * (reduced.size - reduced.sum(dtype=intp)) / ox.size
-            remarks = ['Mismatch: {:.3g}%'.format(mismatch)]
+            n_mismatch = reduced.size - reduced.sum(dtype=intp)
+            percent_mismatch = 100 * n_mismatch / ox.size
+            remarks = [
+                'Mismatched elements: {} / {} ({:.3g}%)'.format(
+                    n_mismatch, ox.size, percent_mismatch)]
 
             with errstate(invalid='ignore', divide='ignore'):
                 # ignore errors for non-numeric types
diff --git a/numpy/testing/tests/test_utils.py b/numpy/testing/tests/test_utils.py
index bf60772d3..4f1b46d4f 100644
--- a/numpy/testing/tests/test_utils.py
+++ b/numpy/testing/tests/test_utils.py
@@ -520,7 +520,7 @@ class TestAlmostEqual(_GenericTest):
         with pytest.raises(AssertionError) as exc_info:
             self._assert_func(x, y, decimal=12)
         msgs = str(exc_info.value).split('\n')
-        assert_equal(msgs[3], 'Mismatch: 100%')
+        assert_equal(msgs[3], 'Mismatched elements: 3 / 3 (100%)')
         assert_equal(msgs[4], 'Max absolute difference: 1.e-05')
         assert_equal(msgs[5], 'Max relative difference: 3.33328889e-06')
         assert_equal(
@@ -536,7 +536,7 @@ class TestAlmostEqual(_GenericTest):
         with pytest.raises(AssertionError) as exc_info:
             self._assert_func(x, y)
         msgs = str(exc_info.value).split('\n')
-        assert_equal(msgs[3], 'Mismatch: 33.3%')
+        assert_equal(msgs[3], 'Mismatched elements: 1 / 3 (33.3%)')
         assert_equal(msgs[4], 'Max absolute difference: 1.e-05')
         assert_equal(msgs[5], 'Max relative difference: 3.33328889e-06')
         assert_equal(msgs[6], ' x: array([1.     , 2.     , 3.00003])')
@@ -548,7 +548,7 @@ class TestAlmostEqual(_GenericTest):
         with pytest.raises(AssertionError) as exc_info:
             self._assert_func(x, y)
         msgs = str(exc_info.value).split('\n')
-        assert_equal(msgs[3], 'Mismatch: 50%')
+        assert_equal(msgs[3], 'Mismatched elements: 1 / 2 (50%)')
         assert_equal(msgs[4], 'Max absolute difference: 1.')
         assert_equal(msgs[5], 'Max relative difference: 1.')
         assert_equal(msgs[6], ' x: array([inf,  0.])')
@@ -560,7 +560,7 @@ class TestAlmostEqual(_GenericTest):
         with pytest.raises(AssertionError) as exc_info:
             self._assert_func(x, y)
         msgs = str(exc_info.value).split('\n')
-        assert_equal(msgs[3], 'Mismatch: 100%')
+        assert_equal(msgs[3], 'Mismatched elements: 2 / 2 (100%)')
         assert_equal(msgs[4], 'Max absolute difference: 2')
         assert_equal(msgs[5], 'Max relative difference: inf')
 
@@ -855,7 +855,8 @@ class TestAssertAllclose(object):
         with pytest.raises(AssertionError) as exc_info:
             assert_allclose(a, b)
         msg = str(exc_info.value)
-        assert_('Mismatch: 25%\nMax absolute difference: 1\n'
+        assert_('Mismatched elements: 1 / 4 (25%)\n'
+                'Max absolute difference: 1\n'
                 'Max relative difference: 0.5' in msg)
 
     def test_equal_nan(self):
diff --git a/numpy/testing/utils.py b/numpy/testing/utils.py
index 98f19e348..1e7d65b89 100644
--- a/numpy/testing/utils.py
+++ b/numpy/testing/utils.py
@@ -19,7 +19,7 @@ __all__ = [
         'assert_array_equal', 'assert_array_less', 'assert_string_equal',
         'assert_array_almost_equal', 'assert_raises', 'build_err_msg',
         'decorate_methods', 'jiffies', 'memusage', 'print_assert_equal',
-        'raises', 'rand', 'rundocs', 'runstring', 'verbose', 'measure',
+        'raises', 'rundocs', 'runstring', 'verbose', 'measure',
         'assert_', 'assert_array_almost_equal_nulp', 'assert_raises_regex',
         'assert_array_max_ulp', 'assert_warns', 'assert_no_warnings',
         'assert_allclose', 'IgnoreException', 'clear_and_catch_warnings',
diff --git a/numpy/tests/test_public_api.py b/numpy/tests/test_public_api.py
index 807c98652..df2fc4802 100644
--- a/numpy/tests/test_public_api.py
+++ b/numpy/tests/test_public_api.py
@@ -1,6 +1,7 @@
 from __future__ import division, absolute_import, print_function
 
 import sys
+import subprocess
 
 import numpy as np
 import pytest
@@ -69,6 +70,28 @@ def test_numpy_namespace():
     assert bad_results == whitelist
 
 
+@pytest.mark.parametrize('name', ['testing', 'Tester'])
+def test_import_lazy_import(name):
+    """Make sure we can actually the the modules we lazy load.
+
+    While not exported as part of the public API, it was accessible.  With the
+    use of __getattr__ and __dir__, this isn't always true It can happen that
+    an infinite recursion may happen.
+
+    This is the only way I found that would force the failure to appear on the
+    badly implemented code.
+
+    We also test for the presence of the lazily imported modules in dir
+
+    """
+    exe = (sys.executable, '-c', "import numpy; numpy." + name)
+    result = subprocess.check_output(exe)
+    assert not result
+
+    # Make sure they are still in the __dir__
+    assert name in dir(np)
+
+
 def test_numpy_linalg():
     bad_results = check_dir(np.linalg)
     assert bad_results == {}