Merge branch 'master' into isfinite-datetime

22 files changed, 263 insertions, 177 deletions

diff --git a/numpy/core/_methods.py b/numpy/core/_methods.py
index 51362c761..953e7e1b8 100644
--- a/numpy/core/_methods.py
+++ b/numpy/core/_methods.py
@@ -115,10 +115,11 @@ def _var(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):
     # Note that x may not be inexact and that we need it to be an array,
     # not a scalar.
     x = asanyarray(arr - arrmean)
-    if issubclass(arr.dtype.type, nt.complexfloating):
-        x = um.multiply(x, um.conjugate(x), out=x).real
-    else:
+    if issubclass(arr.dtype.type, (nt.floating, nt.integer)):
         x = um.multiply(x, x, out=x)
+    else:
+        x = um.multiply(x, um.conjugate(x), out=x).real
+
     ret = umr_sum(x, axis, dtype, out, keepdims)
 
     # Compute degrees of freedom and make sure it is not negative.
diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py
index 6dd6982df..7591a7952 100644
--- a/numpy/core/code_generators/ufunc_docstrings.py
+++ b/numpy/core/code_generators/ufunc_docstrings.py
@@ -1313,7 +1313,7 @@ add_newdoc('numpy.core.umath', 'floor_divide',
     """
     Return the largest integer smaller or equal to the division of the inputs.
     It is equivalent to the Python ``//`` operator and pairs with the
-    Python ``%`` (`remainder`), function so that ``b = a % b + b * (a // b)``
+    Python ``%`` (`remainder`), function so that ``a = a % b + b * (a // b)``
     up to roundoff.
 
     Parameters
@@ -2607,7 +2607,7 @@ add_newdoc('numpy.core.umath', 'matmul',
 
     >>> a = np.array([[1, 0],
     ...               [0, 1]])
-    >>> b = np.array([[4, 1], 
+    >>> b = np.array([[4, 1],
     ...               [2, 2]])
     >>> np.matmul(a, b)
     array([[4, 1],
diff --git a/numpy/core/src/multiarray/common.c b/numpy/core/src/multiarray/common.c
index addb67732..52694d491 100644
--- a/numpy/core/src/multiarray/common.c
+++ b/numpy/core/src/multiarray/common.c
@@ -343,7 +343,7 @@ PyArray_DTypeFromObjectHelper(PyObject *obj, int maxdims,
             typestr = PyDict_GetItemString(ip, "typestr");
 #if defined(NPY_PY3K)
             /* Allow unicode type strings */
-            if (PyUnicode_Check(typestr)) {
+            if (typestr && PyUnicode_Check(typestr)) {
                 tmp = PyUnicode_AsASCIIString(typestr);
                 typestr = tmp;
             }
diff --git a/numpy/core/src/multiarray/item_selection.c b/numpy/core/src/multiarray/item_selection.c
index ff7c1130a..6065c1df4 100644
--- a/numpy/core/src/multiarray/item_selection.c
+++ b/numpy/core/src/multiarray/item_selection.c
@@ -607,7 +607,8 @@ PyArray_Repeat(PyArrayObject *aop, PyObject *op, int axis)
     else {
         for (j = 0; j < n; j++) {
             if (counts[j] < 0) {
-                PyErr_SetString(PyExc_ValueError, "count < 0");
+                PyErr_SetString(PyExc_ValueError,
+                                "repeats may not contain negative values.");
                 goto fail;
             }
             total += counts[j];
diff --git a/numpy/core/src/multiarray/iterators.c b/numpy/core/src/multiarray/iterators.c
index a3bc8e742..9fcdc91b2 100644
--- a/numpy/core/src/multiarray/iterators.c
+++ b/numpy/core/src/multiarray/iterators.c
@@ -539,6 +539,7 @@ iter_subscript(PyArrayIterObject *self, PyObject *ind)
     char *dptr;
     int size;
     PyObject *obj = NULL;
+    PyObject *new;
     PyArray_CopySwapFunc *copyswap;
 
     if (ind == Py_Ellipsis) {
@@ -640,36 +641,36 @@ iter_subscript(PyArrayIterObject *self, PyObject *ind)
         obj = ind;
     }
 
-    if (PyArray_Check(obj)) {
-        /* Check for Boolean object */
-        if (PyArray_TYPE((PyArrayObject *)obj) == NPY_BOOL) {
-            ret = iter_subscript_Bool(self, (PyArrayObject *)obj);
-            Py_DECREF(indtype);
-        }
-        /* Check for integer array */
-        else if (PyArray_ISINTEGER((PyArrayObject *)obj)) {
-            PyObject *new;
-            new = PyArray_FromAny(obj, indtype, 0, 0,
-                              NPY_ARRAY_FORCECAST | NPY_ARRAY_ALIGNED, NULL);
-            if (new == NULL) {
-                goto fail;
-            }
-            Py_DECREF(obj);
-            obj = new;
-            new = iter_subscript_int(self, (PyArrayObject *)obj);
-            Py_DECREF(obj);
-            return new;
-        }
-        else {
-            goto fail;
-        }
+    /* Any remaining valid input is an array or has been turned into one */
+    if (!PyArray_Check(obj)) {
+        goto fail;
+    }
+
+    /* Check for Boolean array */
+    if (PyArray_TYPE((PyArrayObject *)obj) == NPY_BOOL) {
+        ret = iter_subscript_Bool(self, (PyArrayObject *)obj);
+        Py_DECREF(indtype);
         Py_DECREF(obj);
         return (PyObject *)ret;
     }
-    else {
-        Py_DECREF(indtype);
+
+    /* Only integer arrays left */
+    if (!PyArray_ISINTEGER((PyArrayObject *)obj)) {
+        goto fail;
     }
 
+    Py_INCREF(indtype);
+    new = PyArray_FromAny(obj, indtype, 0, 0,
+                      NPY_ARRAY_FORCECAST | NPY_ARRAY_ALIGNED, NULL);
+    if (new == NULL) {
+        goto fail;
+    }
+    Py_DECREF(indtype);
+    Py_DECREF(obj);
+    ret = (PyArrayObject *)iter_subscript_int(self, (PyArrayObject *)new);
+    Py_DECREF(new);
+    return (PyObject *)ret;
+
 
  fail:
     if (!PyErr_Occurred()) {
diff --git a/numpy/core/src/umath/fast_loop_macros.h b/numpy/core/src/umath/fast_loop_macros.h
index e3cfa1f72..7a1ed66bc 100644
--- a/numpy/core/src/umath/fast_loop_macros.h
+++ b/numpy/core/src/umath/fast_loop_macros.h
@@ -74,73 +74,52 @@
 #define IS_BINARY_CONT(tin, tout) (steps[0] == sizeof(tin) && \
                                    steps[1] == sizeof(tin) && \
                                    steps[2] == sizeof(tout))
+
 /* binary loop input and output contiguous with first scalar */
 #define IS_BINARY_CONT_S1(tin, tout) (steps[0] == 0 && \
                                    steps[1] == sizeof(tin) && \
                                    steps[2] == sizeof(tout))
+
 /* binary loop input and output contiguous with second scalar */
 #define IS_BINARY_CONT_S2(tin, tout) (steps[0] == sizeof(tin) && \
                                    steps[1] == 0 && \
                                    steps[2] == sizeof(tout))
 
-
-/*
- * loop with contiguous specialization
- * op should be the code storing the result in `tout * out`
- * combine with NPY_GCC_OPT_3 to allow autovectorization
- * should only be used where its worthwhile to avoid code bloat
- */
-#define BASE_OUTPUT_LOOP(tout, op) \
-    OUTPUT_LOOP { \
-        tout * out = (tout *)op1; \
-        op; \
-    }
-#define OUTPUT_LOOP_FAST(tout, op) \
-    do { \
-    /* condition allows compiler to optimize the generic macro */ \
-    if (IS_OUTPUT_CONT(tout)) { \
-        BASE_OUTPUT_LOOP(tout, op) \
-    } \
-    else { \
-        BASE_OUTPUT_LOOP(tout, op) \
-    } \
-    } \
-    while (0)
-
 /*
  * loop with contiguous specialization
  * op should be the code working on `tin in` and
- * storing the result in `tout * out`
+ * storing the result in `tout *out`
  * combine with NPY_GCC_OPT_3 to allow autovectorization
  * should only be used where its worthwhile to avoid code bloat
  */
 #define BASE_UNARY_LOOP(tin, tout, op) \
     UNARY_LOOP { \
         const tin in = *(tin *)ip1; \
-        tout * out = (tout *)op1; \
+        tout *out = (tout *)op1; \
         op; \
     }
-#define UNARY_LOOP_FAST(tin, tout, op) \
+
+#define UNARY_LOOP_FAST(tin, tout, op)          \
     do { \
-    /* condition allows compiler to optimize the generic macro */ \
-    if (IS_UNARY_CONT(tin, tout)) { \
-        if (args[0] == args[1]) { \
-            BASE_UNARY_LOOP(tin, tout, op) \
+        /* condition allows compiler to optimize the generic macro */ \
+        if (IS_UNARY_CONT(tin, tout)) { \
+            if (args[0] == args[1]) { \
+                BASE_UNARY_LOOP(tin, tout, op) \
+            } \
+            else { \
+                BASE_UNARY_LOOP(tin, tout, op) \
+            } \
         } \
         else { \
             BASE_UNARY_LOOP(tin, tout, op) \
         } \
     } \
-    else { \
-        BASE_UNARY_LOOP(tin, tout, op) \
-    } \
-    } \
     while (0)
 
 /*
  * loop with contiguous specialization
  * op should be the code working on `tin in1`, `tin in2` and
- * storing the result in `tout * out`
+ * storing the result in `tout *out`
  * combine with NPY_GCC_OPT_3 to allow autovectorization
  * should only be used where its worthwhile to avoid code bloat
  */
@@ -148,9 +127,10 @@
     BINARY_LOOP { \
         const tin in1 = *(tin *)ip1; \
         const tin in2 = *(tin *)ip2; \
-        tout * out = (tout *)op1; \
+        tout *out = (tout *)op1; \
         op; \
     }
+
 /*
  * unfortunately gcc 6/7 regressed and we need to give it additional hints to
  * vectorize inplace operations (PR80198)
@@ -171,59 +151,62 @@
     for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1) { \
         const tin in1 = *(tin *)ip1; \
         const tin in2 = *(tin *)ip2; \
-        tout * out = (tout *)op1; \
+        tout *out = (tout *)op1; \
         op; \
     }
+
 #define BASE_BINARY_LOOP_S(tin, tout, cin, cinp, vin, vinp, op) \
     const tin cin = *(tin *)cinp; \
     BINARY_LOOP { \
         const tin vin = *(tin *)vinp; \
-        tout * out = (tout *)op1; \
+        tout *out = (tout *)op1; \
         op; \
     }
+
 /* PR80198 again, scalar works without the pragma */
 #define BASE_BINARY_LOOP_S_INP(tin, tout, cin, cinp, vin, vinp, op) \
     const tin cin = *(tin *)cinp; \
     BINARY_LOOP { \
         const tin vin = *(tin *)vinp; \
-        tout * out = (tout *)vinp; \
+        tout *out = (tout *)vinp; \
         op; \
     }
-#define BINARY_LOOP_FAST(tin, tout, op) \
+
+#define BINARY_LOOP_FAST(tin, tout, op)         \
     do { \
-    /* condition allows compiler to optimize the generic macro */ \
-    if (IS_BINARY_CONT(tin, tout)) { \
-        if (abs_ptrdiff(args[2], args[0]) == 0 && \
-                abs_ptrdiff(args[2], args[1]) >= NPY_MAX_SIMD_SIZE) { \
-            BASE_BINARY_LOOP_INP(tin, tout, op) \
+        /* condition allows compiler to optimize the generic macro */ \
+        if (IS_BINARY_CONT(tin, tout)) { \
+            if (abs_ptrdiff(args[2], args[0]) == 0 && \
+                    abs_ptrdiff(args[2], args[1]) >= NPY_MAX_SIMD_SIZE) { \
+                BASE_BINARY_LOOP_INP(tin, tout, op) \
+            } \
+            else if (abs_ptrdiff(args[2], args[1]) == 0 && \
+                         abs_ptrdiff(args[2], args[0]) >= NPY_MAX_SIMD_SIZE) { \
+                BASE_BINARY_LOOP_INP(tin, tout, op) \
+            } \
+            else { \
+                BASE_BINARY_LOOP(tin, tout, op) \
+            } \
         } \
-        else if (abs_ptrdiff(args[2], args[1]) == 0 && \
-                     abs_ptrdiff(args[2], args[0]) >= NPY_MAX_SIMD_SIZE) { \
-            BASE_BINARY_LOOP_INP(tin, tout, op) \
+        else if (IS_BINARY_CONT_S1(tin, tout)) { \
+            if (abs_ptrdiff(args[2], args[1]) == 0) { \
+                BASE_BINARY_LOOP_S_INP(tin, tout, in1, args[0], in2, ip2, op) \
+            } \
+            else { \
+                BASE_BINARY_LOOP_S(tin, tout, in1, args[0], in2, ip2, op) \
+            } \
         } \
-        else { \
-            BASE_BINARY_LOOP(tin, tout, op) \
-        } \
-    } \
-    else if (IS_BINARY_CONT_S1(tin, tout)) { \
-        if (abs_ptrdiff(args[2], args[1]) == 0) { \
-            BASE_BINARY_LOOP_S_INP(tin, tout, in1, args[0], in2, ip2, op) \
+        else if (IS_BINARY_CONT_S2(tin, tout)) { \
+            if (abs_ptrdiff(args[2], args[0]) == 0) { \
+                BASE_BINARY_LOOP_S_INP(tin, tout, in2, args[1], in1, ip1, op) \
+            } \
+            else { \
+                BASE_BINARY_LOOP_S(tin, tout, in2, args[1], in1, ip1, op) \
+            }\
         } \
         else { \
-            BASE_BINARY_LOOP_S(tin, tout, in1, args[0], in2, ip2, op) \
-        } \
-    } \
-    else if (IS_BINARY_CONT_S2(tin, tout)) { \
-        if (abs_ptrdiff(args[2], args[0]) == 0) { \
-            BASE_BINARY_LOOP_S_INP(tin, tout, in2, args[1], in1, ip1, op) \
+            BASE_BINARY_LOOP(tin, tout, op) \
         } \
-        else { \
-            BASE_BINARY_LOOP_S(tin, tout, in2, args[1], in1, ip1, op) \
-        }\
-    } \
-    else { \
-        BASE_BINARY_LOOP(tin, tout, op) \
-    } \
     } \
     while (0)
 
diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 1e65acd3b..290a87a33 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -652,7 +652,11 @@ BOOL__ones_like(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UN
 NPY_NO_EXPORT void
 BOOL_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
-    OUTPUT_LOOP_FAST(npy_bool, *out = @val@);
+    /*
+     * The (void)in; suppresses an unused variable warning raised by gcc and allows
+     * us to re-use this macro even though we do not depend on in
+     */
+    UNARY_LOOP_FAST(npy_bool, npy_bool, (void)in; *out = @val@);
 }
 
 /**end repeat**/
@@ -896,7 +900,11 @@ NPY_NO_EXPORT void
 NPY_NO_EXPORT void
 @TYPE@_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
-    OUTPUT_LOOP_FAST(npy_bool, *out = @val@);
+    /*
+     * The (void)in; suppresses an unused variable warning raised by gcc and allows
+     * us to re-use this macro even though we do not depend on in
+     */
+    UNARY_LOOP_FAST(@type@, npy_bool, (void)in; *out = @val@);
 }
 /**end repeat1**/
 
@@ -1019,13 +1027,10 @@ NPY_NO_EXPORT void
  * #c    = u,u,u,ul,ull#
  */
 
-NPY_NO_EXPORT void
+NPY_NO_EXPORT NPY_GCC_OPT_3 void
 @TYPE@_absolute(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
-    UNARY_LOOP {
-        const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = in1;
-    }
+    UNARY_LOOP_FAST(@type@, @type@, *out = in);
 }
 
 NPY_NO_EXPORT NPY_GCC_OPT_3 void
@@ -2232,13 +2237,10 @@ HALF_conjugate(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNU
     }
 }
 
-NPY_NO_EXPORT void
+NPY_NO_EXPORT NPY_GCC_OPT_3 void
 HALF_absolute(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
-    UNARY_LOOP {
-        const npy_half in1 = *(npy_half *)ip1;
-        *((npy_half *)op1) = in1&0x7fffu;
-    }
+    UNARY_LOOP_FAST(npy_half, npy_half, *out = in&0x7fffu);
 }
 
 NPY_NO_EXPORT void
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index ae2fd3cf4..c45029599 100644
--- a/numpy/core/tests/test_multiarray.py
+++ b/numpy/core/tests/test_multiarray.py
@@ -4893,6 +4893,22 @@ class TestFlat(object):
         assert_(e.flags.writebackifcopy is False)
         assert_(f.flags.writebackifcopy is False)
 
+    @pytest.mark.skipif(not HAS_REFCOUNT, reason="Python lacks refcounts")
+    def test_refcount(self):
+        # includes regression test for reference count error gh-13165
+        inds = [np.intp(0), np.array([True]*self.a.size), np.array([0]), None]
+        indtype = np.dtype(np.intp)
+        rc_indtype = sys.getrefcount(indtype)
+        for ind in inds:
+            rc_ind = sys.getrefcount(ind)
+            for _ in range(100):
+                try:
+                    self.a.flat[ind]
+                except IndexError:
+                    pass
+            assert_(abs(sys.getrefcount(ind) - rc_ind) < 50)
+            assert_(abs(sys.getrefcount(indtype) - rc_indtype) < 50)
+
 
 class TestResize(object):
     def test_basic(self):
diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py
index 90ac43a56..1822a7adf 100644
--- a/numpy/core/tests/test_numeric.py
+++ b/numpy/core/tests/test_numeric.py
@@ -216,6 +216,9 @@ class TestNonarrayArgs(object):
             assert_(np.isnan(np.var([])))
             assert_(w[0].category is RuntimeWarning)
 
+        B = np.array([None, 0])
+        B[0] = 1j
+        assert_almost_equal(np.var(B), 0.25)
 
 class TestIsscalar(object):
     def test_isscalar(self):
diff --git a/numpy/core/tests/test_regression.py b/numpy/core/tests/test_regression.py
index 1ba0cda21..4b551d8aa 100644
--- a/numpy/core/tests/test_regression.py
+++ b/numpy/core/tests/test_regression.py
@@ -2448,3 +2448,9 @@ class TestRegression(object):
         for proto in range(2, pickle.HIGHEST_PROTOCOL + 1):
             dumped = pickle.dumps(arr, protocol=proto)
             assert_equal(pickle.loads(dumped), arr)
+
+    def test_bad_array_interface(self):
+        class T(object):
+            __array_interface__ = {}
+
+        np.array([T()])
diff --git a/numpy/doc/byteswapping.py b/numpy/doc/byteswapping.py
index f9491ed43..7a749c8d5 100644
--- a/numpy/doc/byteswapping.py
+++ b/numpy/doc/byteswapping.py
@@ -31,16 +31,16 @@ Let's say the two integers were in fact 1 and 770.  Because 770 = 256 *
 3 + 2, the 4 bytes in memory would contain respectively: 0, 1, 3, 2.
 The bytes I have loaded from the file would have these contents:
 
->>> big_end_str = chr(0) + chr(1) + chr(3) + chr(2)
->>> big_end_str
-'\\x00\\x01\\x03\\x02'
+>>> big_end_buffer = bytearray([0,1,3,2])
+>>> big_end_buffer
+bytearray(b'\\x00\\x01\\x03\\x02')
 
 We might want to use an ``ndarray`` to access these integers.  In that
 case, we can create an array around this memory, and tell numpy that
 there are two integers, and that they are 16 bit and big-endian:
 
 >>> import numpy as np
->>> big_end_arr = np.ndarray(shape=(2,),dtype='>i2', buffer=big_end_str)
+>>> big_end_arr = np.ndarray(shape=(2,),dtype='>i2', buffer=big_end_buffer)
 >>> big_end_arr[0]
 1
 >>> big_end_arr[1]
@@ -53,7 +53,7 @@ integer, the dtype string would be ``<u4``.
 
 In fact, why don't we try that?
 
->>> little_end_u4 = np.ndarray(shape=(1,),dtype='<u4', buffer=big_end_str)
+>>> little_end_u4 = np.ndarray(shape=(1,),dtype='<u4', buffer=big_end_buffer)
 >>> little_end_u4[0] == 1 * 256**1 + 3 * 256**2 + 2 * 256**3
 True
 
@@ -97,7 +97,7 @@ Data and dtype endianness don't match, change dtype to match data
 
 We make something where they don't match:
 
->>> wrong_end_dtype_arr = np.ndarray(shape=(2,),dtype='<i2', buffer=big_end_str)
+>>> wrong_end_dtype_arr = np.ndarray(shape=(2,),dtype='<i2', buffer=big_end_buffer)
 >>> wrong_end_dtype_arr[0]
 256
 
@@ -110,7 +110,7 @@ the correct endianness:
 
 Note the array has not changed in memory:
 
->>> fixed_end_dtype_arr.tobytes() == big_end_str
+>>> fixed_end_dtype_arr.tobytes() == big_end_buffer
 True
 
 Data and type endianness don't match, change data to match dtype
@@ -126,7 +126,7 @@ that needs a certain byte ordering.
 
 Now the array *has* changed in memory:
 
->>> fixed_end_mem_arr.tobytes() == big_end_str
+>>> fixed_end_mem_arr.tobytes() == big_end_buffer
 False
 
 Data and dtype endianness match, swap data and dtype
@@ -140,7 +140,7 @@ the previous operations:
 >>> swapped_end_arr = big_end_arr.byteswap().newbyteorder()
 >>> swapped_end_arr[0]
 1
->>> swapped_end_arr.tobytes() == big_end_str
+>>> swapped_end_arr.tobytes() == big_end_buffer
 False
 
 An easier way of casting the data to a specific dtype and byte ordering
@@ -149,7 +149,7 @@ can be achieved with the ndarray astype method:
 >>> swapped_end_arr = big_end_arr.astype('<i2')
 >>> swapped_end_arr[0]
 1
->>> swapped_end_arr.tobytes() == big_end_str
+>>> swapped_end_arr.tobytes() == big_end_buffer
 False
 
 """
diff --git a/numpy/f2py/src/fortranobject.c b/numpy/f2py/src/fortranobject.c
index 78b06f066..4a981bf55 100644
--- a/numpy/f2py/src/fortranobject.c
+++ b/numpy/f2py/src/fortranobject.c
@@ -135,7 +135,7 @@ format_def(char *buf, Py_ssize_t size, FortranDataDef def)
 
     if (def.data == NULL) {
         static const char notalloc[] = ", not allocated";
-        if (size < sizeof(notalloc)) {
+        if ((size_t) size < sizeof(notalloc)) {
             return -1;
         }
         memcpy(p, notalloc, sizeof(notalloc));
diff --git a/numpy/lib/tests/test_type_check.py b/numpy/lib/tests/test_type_check.py
index 2982ca31a..b3f114b92 100644
--- a/numpy/lib/tests/test_type_check.py
+++ b/numpy/lib/tests/test_type_check.py
@@ -360,6 +360,14 @@ class TestNanToNum(object):
         assert_(vals[1] == 0)
         assert_all(vals[2] > 1e10) and assert_all(np.isfinite(vals[2]))
         assert_equal(type(vals), np.ndarray)
+        
+        # perform the same tests but with nan, posinf and neginf keywords
+        with np.errstate(divide='ignore', invalid='ignore'):
+            vals = nan_to_num(np.array((-1., 0, 1))/0., 
+                              nan=10, posinf=20, neginf=30)
+        assert_equal(vals, [30, 10, 20])
+        assert_all(np.isfinite(vals[[0, 2]]))
+        assert_equal(type(vals), np.ndarray)
 
         # perform the same test but in-place
         with np.errstate(divide='ignore', invalid='ignore'):
@@ -371,26 +379,48 @@ class TestNanToNum(object):
         assert_(vals[1] == 0)
         assert_all(vals[2] > 1e10) and assert_all(np.isfinite(vals[2]))
         assert_equal(type(vals), np.ndarray)
+        
+        # perform the same test but in-place
+        with np.errstate(divide='ignore', invalid='ignore'):
+            vals = np.array((-1., 0, 1))/0.
+        result = nan_to_num(vals, copy=False, nan=10, posinf=20, neginf=30)
+
+        assert_(result is vals)
+        assert_equal(vals, [30, 10, 20])
+        assert_all(np.isfinite(vals[[0, 2]]))
+        assert_equal(type(vals), np.ndarray)
 
     def test_array(self):
         vals = nan_to_num([1])
         assert_array_equal(vals, np.array([1], int))
         assert_equal(type(vals), np.ndarray)
+        vals = nan_to_num([1], nan=10, posinf=20, neginf=30)
+        assert_array_equal(vals, np.array([1], int))
+        assert_equal(type(vals), np.ndarray)
 
     def test_integer(self):
         vals = nan_to_num(1)
         assert_all(vals == 1)
         assert_equal(type(vals), np.int_)
+        vals = nan_to_num(1, nan=10, posinf=20, neginf=30)
+        assert_all(vals == 1)
+        assert_equal(type(vals), np.int_)
 
     def test_float(self):
         vals = nan_to_num(1.0)
         assert_all(vals == 1.0)
         assert_equal(type(vals), np.float_)
+        vals = nan_to_num(1.1, nan=10, posinf=20, neginf=30)
+        assert_all(vals == 1.1)
+        assert_equal(type(vals), np.float_)
 
     def test_complex_good(self):
         vals = nan_to_num(1+1j)
         assert_all(vals == 1+1j)
         assert_equal(type(vals), np.complex_)
+        vals = nan_to_num(1+1j, nan=10, posinf=20, neginf=30)
+        assert_all(vals == 1+1j)
+        assert_equal(type(vals), np.complex_)
 
     def test_complex_bad(self):
         with np.errstate(divide='ignore', invalid='ignore'):
@@ -414,6 +444,16 @@ class TestNanToNum(object):
         # !! inf.  Comment out for now, and see if it
         # !! changes
         #assert_all(vals.real < -1e10) and assert_all(np.isfinite(vals))
+    
+    def test_do_not_rewrite_previous_keyword(self):
+        # This is done to test that when, for instance, nan=np.inf then these 
+        # values are not rewritten by posinf keyword to the posinf value.
+        with np.errstate(divide='ignore', invalid='ignore'):
+            vals = nan_to_num(np.array((-1., 0, 1))/0., nan=np.inf, posinf=999)
+        assert_all(np.isfinite(vals[[0, 2]]))
+        assert_all(vals[0] < -1e10)
+        assert_equal(vals[[1, 2]], [np.inf, 999])
+        assert_equal(type(vals), np.ndarray)
 
 
 class TestRealIfClose(object):
diff --git a/numpy/lib/type_check.py b/numpy/lib/type_check.py
index f55517732..2b254b6c0 100644
--- a/numpy/lib/type_check.py
+++ b/numpy/lib/type_check.py
@@ -363,18 +363,23 @@ def _getmaxmin(t):
     return f.max, f.min
 
 
-def _nan_to_num_dispatcher(x, copy=None):
+def _nan_to_num_dispatcher(x, copy=None, nan=None, posinf=None, neginf=None):
     return (x,)
 
 
 @array_function_dispatch(_nan_to_num_dispatcher)
-def nan_to_num(x, copy=True):
+def nan_to_num(x, copy=True, nan=0.0, posinf=None, neginf=None):
     """
-    Replace NaN with zero and infinity with large finite numbers.
+    Replace NaN with zero and infinity with large finite numbers (default
+    behaviour) or with the numbers defined by the user using the `nan`, 
+    `posinf` and/or `neginf` keywords.
 
-    If `x` is inexact, NaN is replaced by zero, and infinity and -infinity
-    replaced by the respectively largest and most negative finite floating
-    point values representable by ``x.dtype``.
+    If `x` is inexact, NaN is replaced by zero or by the user defined value in
+    `nan` keyword, infinity is replaced by the largest finite floating point 
+    values representable by ``x.dtype`` or by the user defined value in 
+    `posinf` keyword and -infinity is replaced by the most negative finite 
+    floating point values representable by ``x.dtype`` or by the user defined 
+    value in `neginf` keyword.
 
     For complex dtypes, the above is applied to each of the real and
     imaginary components of `x` separately.
@@ -390,6 +395,17 @@ def nan_to_num(x, copy=True):
         in-place (False). The in-place operation only occurs if
         casting to an array does not require a copy.
         Default is True.
+    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.
+    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.
+    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.13
 
@@ -424,6 +440,9 @@ def nan_to_num(x, copy=True):
     >>> np.nan_to_num(x)
     array([ 1.79769313e+308, -1.79769313e+308,  0.00000000e+000, # may vary
            -1.28000000e+002,  1.28000000e+002])
+    >>> np.nan_to_num(x, nan=-9999, posinf=33333333, neginf=33333333)
+    array([ 3.3333333e+07,  3.3333333e+07, -9.9990000e+03, 
+           -1.2800000e+02,  1.2800000e+02])
     >>> y = np.array([complex(np.inf, np.nan), np.nan, complex(np.nan, np.inf)])
     array([  1.79769313e+308,  -1.79769313e+308,   0.00000000e+000, # may vary
          -1.28000000e+002,   1.28000000e+002])
@@ -431,6 +450,8 @@ def nan_to_num(x, copy=True):
     array([  1.79769313e+308 +0.00000000e+000j, # may vary
              0.00000000e+000 +0.00000000e+000j,
              0.00000000e+000 +1.79769313e+308j])
+    >>> np.nan_to_num(y, nan=111111, posinf=222222)
+    array([222222.+111111.j, 111111.     +0.j, 111111.+222222.j])
     """
     x = _nx.array(x, subok=True, copy=copy)
     xtype = x.dtype.type
@@ -444,10 +465,17 @@ def nan_to_num(x, copy=True):
 
     dest = (x.real, x.imag) if iscomplex else (x,)
     maxf, minf = _getmaxmin(x.real.dtype)
+    if posinf is not None:
+        maxf = posinf
+    if neginf is not None:
+        minf = neginf
     for d in dest:
-        _nx.copyto(d, 0.0, where=isnan(d))
-        _nx.copyto(d, maxf, where=isposinf(d))
-        _nx.copyto(d, minf, where=isneginf(d))
+        idx_nan = isnan(d)
+        idx_posinf = isposinf(d)
+        idx_neginf = isneginf(d)
+        _nx.copyto(d, nan, where=idx_nan)
+        _nx.copyto(d, maxf, where=idx_posinf)
+        _nx.copyto(d, minf, where=idx_neginf)
     return x[()] if isscalar else x
 
 #-----------------------------------------------------------------------------
diff --git a/numpy/ma/core.py b/numpy/ma/core.py
index 1f5bc2a51..90aff6ec8 100644
--- a/numpy/ma/core.py
+++ b/numpy/ma/core.py
@@ -3450,7 +3450,7 @@ class MaskedArray(ndarray):
         # We could try to force a reshape, but that wouldn't work in some
         # cases.
         return self._mask
-    
+
     @mask.setter
     def mask(self, value):
         self.__setmask__(value)
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 836f47363..c37b2c6d6 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -1743,7 +1743,8 @@ def chebroots(c):
     if len(c) == 2:
         return np.array([-c[0]/c[1]])
 
-    m = chebcompanion(c)
+    # rotated companion matrix reduces error
+    m = chebcompanion(c)[::-1,::-1]
     r = la.eigvals(m)
     r.sort()
     return r
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index dcd0a2b4d..3870d1054 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -1476,7 +1476,8 @@ def hermroots(c):
     if len(c) == 2:
         return np.array([-.5*c[0]/c[1]])
 
-    m = hermcompanion(c)
+    # rotated companion matrix reduces error
+    m = hermcompanion(c)[::-1,::-1]
     r = la.eigvals(m)
     r.sort()
     return r
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index 48e5eab20..b12c0d792 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -1471,7 +1471,8 @@ def hermeroots(c):
     if len(c) == 2:
         return np.array([-c[0]/c[1]])
 
-    m = hermecompanion(c)
+    # rotated companion matrix reduces error
+    m = hermecompanion(c)[::-1,::-1]
     r = la.eigvals(m)
     r.sort()
     return r
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 41d15142a..e103dde92 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -1475,7 +1475,8 @@ def lagroots(c):
     if len(c) == 2:
         return np.array([1 + c[0]/c[1]])
 
-    m = lagcompanion(c)
+    # rotated companion matrix reduces error
+    m = lagcompanion(c)[::-1,::-1]
     r = la.eigvals(m)
     r.sort()
     return r
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index d56e2ae2c..9eec9740d 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -1505,7 +1505,8 @@ def legroots(c):
     if len(c) == 2:
         return np.array([-c[0]/c[1]])
 
-    m = legcompanion(c)
+    # rotated companion matrix reduces error
+    m = legcompanion(c)[::-1,::-1]
     r = la.eigvals(m)
     r.sort()
     return r
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index e3eb7ec5a..afa330502 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -1428,7 +1428,8 @@ def polyroots(c):
     if len(c) == 2:
         return np.array([-c[0]/c[1]])
 
-    m = polycompanion(c)
+    # rotated companion matrix reduces error
+    m = polycompanion(c)[::-1,::-1]
     r = la.eigvals(m)
     r.sort()
     return r
diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx
index 2426dbaa4..c97c166aa 100644
--- a/numpy/random/mtrand/mtrand.pyx
+++ b/numpy/random/mtrand/mtrand.pyx
@@ -866,10 +866,9 @@ cdef class RandomState:
         """
         tomaxint(size=None)
 
-        Random integers between 0 and ``sys.maxint``, inclusive.
-
         Return a sample of uniformly distributed random integers in the interval
-        [0, ``sys.maxint``].
+        [0, ``np.iinfo(np.int).max``]. The np.int type translates to the C long
+        integer type and its precision is platform dependent.
 
         Parameters
         ----------
@@ -891,16 +890,15 @@ cdef class RandomState:
 
         Examples
         --------
-        >>> RS = np.random.mtrand.RandomState() # need a RandomState object
-        >>> RS.tomaxint((2,2,2))
-        array([[[1170048599, 1600360186],
+        >>> rs = np.random.RandomState() # need a RandomState object
+        >>> rs.tomaxint((2,2,2))
+        array([[[1170048599, 1600360186], # random
                 [ 739731006, 1947757578]],
                [[1871712945,  752307660],
                 [1601631370, 1479324245]]])
-        >>> import sys
-        >>> sys.maxint
+        >>> np.iinfo(np.int).max
         2147483647
-        >>> RS.tomaxint((2,2,2)) < sys.maxint
+        >>> rs.tomaxint((2,2,2)) < np.iinfo(np.int).max
         array([[[ True,  True],
                 [ True,  True]],
                [[ True,  True],
@@ -992,7 +990,7 @@ cdef class RandomState:
             raise ValueError("high is out of bounds for %s" % dtype)
         if ilow >= ihigh and np.prod(size) != 0:
             raise ValueError("Range cannot be empty (low >= high) unless no samples are taken")
- 
+
         with self.lock:
             ret = randfunc(ilow, ihigh - 1, size, self.state_address)
 
@@ -1040,7 +1038,7 @@ cdef class RandomState:
                 .. versionadded:: 1.7.0
 
         Parameters
-        -----------
+        ----------
         a : 1-D array-like or int
             If an ndarray, a random sample is generated from its elements.
             If an int, the random sample is generated as if a were np.arange(a)
@@ -1056,12 +1054,12 @@ cdef class RandomState:
             entries in a.
 
         Returns
-        --------
+        -------
         samples : single item or ndarray
             The generated random samples
 
         Raises
-        -------
+        ------
         ValueError
             If a is an int and less than zero, if a or p are not 1-dimensional,
             if a is an array-like of size 0, if p is not a vector of
@@ -1070,11 +1068,11 @@ cdef class RandomState:
             size
 
         See Also
-        ---------
+        --------
         randint, shuffle, permutation
 
         Examples
-        ---------
+        --------
         Generate a uniform random sample from np.arange(5) of size 3:
 
         >>> np.random.choice(5, 3)
@@ -1421,6 +1419,7 @@ cdef class RandomState:
         >>> 3 + 2.5 * np.random.randn(2, 4)
         array([[-4.49401501,  4.00950034, -1.81814867,  7.29718677],   # random
                [ 0.39924804,  4.68456316,  4.99394529,  4.84057254]])  # random
+
         """
         if len(args) == 0:
             return self.standard_normal()
@@ -1436,8 +1435,8 @@ cdef class RandomState:
         Return random integers of type np.int from the "discrete uniform"
         distribution in the closed interval [`low`, `high`].  If `high` is
         None (the default), then results are from [1, `low`]. The np.int
-        type translates to the C long type used by Python 2 for "short"
-        integers and its precision is platform dependent.
+        type translates to the C long integer type and its precision
+        is platform dependent.
 
         This function has been deprecated. Use randint instead.
 
@@ -1708,7 +1707,7 @@ cdef class RandomState:
         .. math:: f(x; a,b) = \\frac{1}{B(\\alpha, \\beta)} x^{\\alpha - 1}
                                                          (1 - x)^{\\beta - 1},
 
-        where the normalisation, B, is the beta function,
+        where the normalization, B, is the beta function,
 
         .. math:: B(\\alpha, \\beta) = \\int_0^1 t^{\\alpha - 1}
                                      (1 - t)^{\\beta - 1} dt.
@@ -2329,7 +2328,7 @@ cdef class RandomState:
 
         Draw samples from a noncentral chi-square distribution.
 
-        The noncentral :math:`\\chi^2` distribution is a generalisation of
+        The noncentral :math:`\\chi^2` distribution is a generalization of
         the :math:`\\chi^2` distribution.
 
         Parameters
@@ -2359,7 +2358,7 @@ cdef class RandomState:
 
         .. math:: P(x;df,nonc) = \\sum^{\\infty}_{i=0}
                                \\frac{e^{-nonc/2}(nonc/2)^{i}}{i!}
-                               \\P_{Y_{df+2i}}(x),
+                               P_{Y_{df+2i}}(x),
 
         where :math:`Y_{q}` is the Chi-square with q degrees of freedom.
 
@@ -2395,6 +2394,7 @@ cdef class RandomState:
         >>> values = plt.hist(np.random.noncentral_chisquare(3, 20, 100000),
         ...                   bins=200, density=True)
         >>> plt.show()
+
         """
         cdef ndarray odf, ononc
         cdef double fdf, fnonc
@@ -2921,7 +2921,7 @@ cdef class RandomState:
         Parameters
         ----------
         a : float or array_like of floats
-            Parameter of the distribution. Should be greater than zero.
+            Parameter of the distribution. Must be non-negative.
         size : int or tuple of ints, optional
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  If size is ``None`` (default),
@@ -3434,7 +3434,7 @@ cdef class RandomState:
         >>> # values, drawn from a normal distribution.
         >>> b = []
         >>> for i in range(1000):
-        ...    a = 10. + np.random.random(100)
+        ...    a = 10. + np.random.standard_normal(100)
         ...    b.append(np.product(a))
 
         >>> b = np.array(b) / np.min(b) # scale values to be positive
@@ -3779,8 +3779,8 @@ cdef class RandomState:
         product p*n <=5, where p = population proportion estimate, and n =
         number of samples, in which case the binomial distribution is used
         instead. For example, a sample of 15 people shows 4 who are left
-        handed, and 11 who are right handed. Then p = 4/15 = 27%. 0.27*15 = 4,
-        so the binomial distribution should be used in this case.
+        handed, and 11 who are right handed. Then p = 4/15 = 27%. Since 
+        0.27*15 = 4, the binomial distribution should be used in this case.
 
         References
         ----------
@@ -4052,7 +4052,7 @@ cdef class RandomState:
         Parameters
         ----------
         a : float or array_like of floats
-            Distribution parameter. Should be greater than 1.
+            Distribution parameter. Must be greater than 1.
         size : int or tuple of ints, optional
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  If size is ``None`` (default),
@@ -4590,7 +4590,7 @@ cdef class RandomState:
 
         Draw samples from a multinomial distribution.
 
-        The multinomial distribution is a multivariate generalisation of the
+        The multinomial distribution is a multivariate generalization of the
         binomial distribution.  Take an experiment with one of ``p``
         possible outcomes.  An example of such an experiment is throwing a dice,
         where the outcome can be 1 through 6.  Each sample drawn from the
@@ -4706,7 +4706,7 @@ cdef class RandomState:
         Draw `size` samples of dimension k from a Dirichlet distribution. A
         Dirichlet-distributed random variable can be seen as a multivariate
         generalization of a Beta distribution. The Dirichlet distribution
-        is a conjugate prior of a multinomial distribution in Bayesian 
+        is a conjugate prior of a multinomial distribution in Bayesian
         inference.
 
         Parameters
@@ -4731,23 +4731,22 @@ cdef class RandomState:
 
         Notes
         -----
-
-        The Dirichlet distribution is a distribution over vectors 
-        :math:`x` that fulfil the conditions :math:`x_i>0` and 
+        The Dirichlet distribution is a distribution over vectors
+        :math:`x` that fulfil the conditions :math:`x_i>0` and
         :math:`\\sum_{i=1}^k x_i = 1`.
 
-        The probability density function :math:`p` of a 
-        Dirichlet-distributed random vector :math:`X` is 
+        The probability density function :math:`p` of a
+        Dirichlet-distributed random vector :math:`X` is
         proportional to
 
         .. math:: p(x) \\propto \\prod_{i=1}^{k}{x^{\\alpha_i-1}_i},
 
-        where :math:`\\alpha` is a vector containing the positive 
+        where :math:`\\alpha` is a vector containing the positive
         concentration parameters.
 
         The method uses the following property for computation: let :math:`Y`
-        be a random vector which has components that follow a standard gamma 
-        distribution, then :math:`X = \\frac{1}{\\sum_{i=1}^k{Y_i}} Y` 
+        be a random vector which has components that follow a standard gamma
+        distribution, then :math:`X = \\frac{1}{\\sum_{i=1}^k{Y_i}} Y`
         is Dirichlet-distributed
 
         References
@@ -4962,7 +4961,7 @@ cdef class RandomState:
             return arr
 
         arr = np.asarray(x)
-    
+
         # shuffle has fast-path for 1-d
         if arr.ndim == 1:
             # Return a copy if same memory
@@ -4975,7 +4974,7 @@ cdef class RandomState:
         idx = np.arange(arr.shape[0], dtype=np.intp)
         self.shuffle(idx)
         return arr[idx]
-        
+
 
 _rand = RandomState()
 seed = _rand.seed