Merge pull request #3104 from seberg/nditer-allow-0d

Make AdvancedNew iter more 0-d aware

12 files changed, 188 insertions, 137 deletions

diff --git a/doc/release/1.8.0-notes.rst b/doc/release/1.8.0-notes.rst
index f41c8e716..1750b5d14 100644
--- a/doc/release/1.8.0-notes.rst
+++ b/doc/release/1.8.0-notes.rst
@@ -37,6 +37,13 @@ compiler, then it's possible you will encounter problems. If so, please
 file a bug and as a temporary workaround you can re-enable the old build
 system by exporting the shell variable NPY_SEPARATE_COMPILATION=0.
 
+For the AdvancedNew iterator the ``oa_ndim`` flag should now be -1 to indicate
+that no ``op_axes`` and ``itershape`` are passed in. The ``oa_ndim == 0``
+case, now indicates a 0-D iteration and ``op_axes`` being NULL and the old 
+usage is deprecated. This does not effect the ``NpyIter_New`` or
+``NpyIter_MultiNew`` functions.
+
+
 New features
 ============
 
diff --git a/doc/source/reference/c-api.iterator.rst b/doc/source/reference/c-api.iterator.rst
index 7e2900bcc..1e3565bc1 100644
--- a/doc/source/reference/c-api.iterator.rst
+++ b/doc/source/reference/c-api.iterator.rst
@@ -634,12 +634,12 @@ Construction and Destruction
     Extends :cfunc:`NpyIter_MultiNew` with several advanced options providing
     more control over broadcasting and buffering.
 
-    If 0/NULL values are passed to ``oa_ndim``, ``op_axes``, ``itershape``,
+    If -1/NULL values are passed to ``oa_ndim``, ``op_axes``, ``itershape``,
     and ``buffersize``, it is equivalent to :cfunc:`NpyIter_MultiNew`.
 
-    The parameter ``oa_ndim``, when non-zero, specifies the number of
+    The parameter ``oa_ndim``, when not zero or -1, specifies the number of
     dimensions that will be iterated with customized broadcasting.
-    If it is provided, ``op_axes`` and/or ``itershape`` must also be provided.
+    If it is provided, ``op_axes`` must and ``itershape`` can also be provided.
     The ``op_axes`` parameter let you control in detail how the
     axes of the operand arrays get matched together and iterated.
     In ``op_axes``, you must provide an array of ``nop`` pointers
@@ -649,6 +649,11 @@ Construction and Destruction
     -1 which means ``newaxis``.  Within each ``op_axes[j]`` array, axes
     may not be repeated.  The following example is how normal broadcasting
     applies to a 3-D array, a 2-D array, a 1-D array and a scalar.
+    
+    **Note**: Before NumPy 1.8 ``oa_ndim == 0` was used for signalling that
+    that ``op_axes`` and ``itershape`` are unused. This is deprecated and
+    should be replaced with -1. Better backward compatibility may be
+    achieved by using :cfunc:`NpyIter_MultiNew` for this case.
 
     .. code-block:: c
 
diff --git a/numpy/core/src/multiarray/nditer_api.c b/numpy/core/src/multiarray/nditer_api.c
index 09e572f10..007958800 100644
--- a/numpy/core/src/multiarray/nditer_api.c
+++ b/numpy/core/src/multiarray/nditer_api.c
@@ -134,12 +134,10 @@ NpyIter_RemoveAxis(NpyIter *iter, int axis)
     axisdata = NIT_INDEX_AXISDATA(axisdata_del, 1);
     memmove(axisdata_del, axisdata, (ndim-1-xdim)*sizeof_axisdata);
 
-    /* If there is more than one dimension, shrink the iterator */
-    if (ndim > 1) {
-        NIT_NDIM(iter) = ndim-1;
-    }
-    /* Otherwise convert it to a singleton dimension */
-    else {
+    /* Shrink the iterator */
+    NIT_NDIM(iter) = ndim - 1;
+    /* If it is now 0-d fill the singleton dimension */
+    if (ndim == 1) {
         npy_intp *strides = NAD_STRIDES(axisdata_del);
         NAD_SHAPE(axisdata_del) = 1;
         for (iop = 0; iop < nop; ++iop) {
@@ -642,6 +640,9 @@ NpyIter_GetIterIndex(NpyIter *iter)
         npy_intp sizeof_axisdata;
 
         iterindex = 0;
+        if (ndim == 0) {
+            return 0;
+        }
         sizeof_axisdata = NIT_AXISDATA_SIZEOF(itflags, ndim, nop);
         axisdata = NIT_INDEX_AXISDATA(NIT_AXISDATA(iter), ndim-1);
 
@@ -1750,6 +1751,8 @@ npyiter_goto_iterindex(NpyIter *iter, npy_intp iterindex)
 
     NIT_ITERINDEX(iter) = iterindex;
 
+    ndim = ndim ? ndim : 1;
+
     if (iterindex == 0) {
         dataptr = NIT_RESETDATAPTR(iter);
 
diff --git a/numpy/core/src/multiarray/nditer_constr.c b/numpy/core/src/multiarray/nditer_constr.c
index cfbaea321..a40cbc7bc 100644
--- a/numpy/core/src/multiarray/nditer_constr.c
+++ b/numpy/core/src/multiarray/nditer_constr.c
@@ -54,8 +54,7 @@ static int
 npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, npyiter_opitflags *op_itflags,
                     char **op_dataptr,
                     npy_uint32 *op_flags, int **op_axes,
-                    npy_intp *itershape,
-                    int output_scalars);
+                    npy_intp *itershape);
 static void
 npyiter_replace_axisdata(NpyIter *iter, int iop,
                       PyArrayObject *op,
@@ -75,7 +74,7 @@ static PyArray_Descr *
 npyiter_get_common_dtype(int nop, PyArrayObject **op,
                         npyiter_opitflags *op_itflags, PyArray_Descr **op_dtype,
                         PyArray_Descr **op_request_dtypes,
-                        int only_inputs, int output_scalars);
+                        int only_inputs);
 static PyArrayObject *
 npyiter_new_temp_array(NpyIter *iter, PyTypeObject *subtype,
                 npy_uint32 flags, npyiter_opitflags *op_itflags,
@@ -86,7 +85,7 @@ npyiter_allocate_arrays(NpyIter *iter,
                         npy_uint32 flags,
                         PyArray_Descr **op_dtype, PyTypeObject *subtype,
                         npy_uint32 *op_flags, npyiter_opitflags *op_itflags,
-                        int **op_axes, int output_scalars);
+                        int **op_axes);
 static void
 npyiter_get_priority_subtype(int nop, PyArrayObject **op,
                             npyiter_opitflags *op_itflags,
@@ -122,8 +121,7 @@ NpyIter_AdvancedNew(int nop, PyArrayObject **op_in, npy_uint32 flags,
 
     npy_int8 *perm;
     NpyIter_BufferData *bufferdata = NULL;
-    int any_allocate = 0, any_missing_dtypes = 0,
-            output_scalars = 0, need_subtype = 0;
+    int any_allocate = 0, any_missing_dtypes = 0, need_subtype = 0;
 
     /* The subtype for automatically allocated outputs */
     double subtype_priority = NPY_PRIORITY;
@@ -158,6 +156,22 @@ NpyIter_AdvancedNew(int nop, PyArrayObject **op_in, npy_uint32 flags,
         return NULL;
     }
 
+    /*
+     * Before 1.8, if `oa_ndim == 0`, this meant `op_axes != NULL` was an error.
+     * With 1.8, `oa_ndim == -1` takes this role, while op_axes in that case
+     * enforces a 0-d iterator. Using `oa_ndim == 0` with `op_axes == NULL`
+     * is thus deprecated with version 1.8.
+     */
+    if ((oa_ndim == 0) && (op_axes == NULL)) {
+        char* mesg = "using `oa_ndim == 0` when `op_axes` is NULL is "
+                     "deprecated. Use `oa_ndim == -1` or the MultiNew "
+                     "iterator for NumPy <1.8 compatibility";
+        if (DEPRECATE(mesg) < 0) {
+            return NULL;
+        }
+        oa_ndim = -1;
+    }
+
     /* Error check 'oa_ndim' and 'op_axes', which must be used together */
     if (!npyiter_check_op_axes(nop, oa_ndim, op_axes, itershape)) {
         return NULL;
@@ -175,12 +189,6 @@ NpyIter_AdvancedNew(int nop, PyArrayObject **op_in, npy_uint32 flags,
     /* Calculate how many dimensions the iterator should have */
     ndim = npyiter_calculate_ndim(nop, op_in, oa_ndim);
 
-    /* If 'ndim' is zero, any outputs should be scalars */
-    if (ndim == 0) {
-        output_scalars = 1;
-        ndim = 1;
-    }
-
     NPY_IT_TIME_POINT(c_calculate_ndim);
 
     /* Allocate memory for the iterator */
@@ -231,8 +239,7 @@ NpyIter_AdvancedNew(int nop, PyArrayObject **op_in, npy_uint32 flags,
 
     /* Fill in the AXISDATA arrays and set the ITERSIZE field */
     if (!npyiter_fill_axisdata(iter, flags, op_itflags, op_dataptr,
-                                        op_flags, op_axes, itershape,
-                                        output_scalars)) {
+                                        op_flags, op_axes, itershape)) {
         NpyIter_Deallocate(iter);
         return NULL;
     }
@@ -338,8 +345,7 @@ NpyIter_AdvancedNew(int nop, PyArrayObject **op_in, npy_uint32 flags,
         dtype = npyiter_get_common_dtype(nop, op,
                                     op_itflags, op_dtype,
                                     op_request_dtypes,
-                                    only_inputs,
-                                    output_scalars);
+                                    only_inputs);
         if (dtype == NULL) {
             NpyIter_Deallocate(iter);
             return NULL;
@@ -389,7 +395,7 @@ NpyIter_AdvancedNew(int nop, PyArrayObject **op_in, npy_uint32 flags,
      * done now using a memory layout matching the iterator.
      */
     if (!npyiter_allocate_arrays(iter, flags, op_dtype, subtype, op_flags,
-                            op_itflags, op_axes, output_scalars)) {
+                            op_itflags, op_axes)) {
         NpyIter_Deallocate(iter);
         return NULL;
     }
@@ -504,7 +510,7 @@ NpyIter_MultiNew(int nop, PyArrayObject **op_in, npy_uint32 flags,
 {
     return NpyIter_AdvancedNew(nop, op_in, flags, order, casting,
                             op_flags, op_request_dtypes,
-                            0, NULL, NULL, 0);
+                            -1, NULL, NULL, 0);
 }
 
 /*NUMPY_API
@@ -521,7 +527,7 @@ NpyIter_New(PyArrayObject *op, npy_uint32 flags,
 
     return NpyIter_AdvancedNew(1, &op, flags, order, casting,
                             &op_flags, &dtype,
-                            0, NULL, NULL, 0);
+                            -1, NULL, NULL, 0);
 }
 
 /*NUMPY_API
@@ -758,53 +764,60 @@ npyiter_check_op_axes(int nop, int oa_ndim, int **op_axes,
     char axes_dupcheck[NPY_MAXDIMS];
     int iop, idim;
 
-    if (oa_ndim == 0 && (op_axes != NULL || itershape != NULL)) {
-        PyErr_Format(PyExc_ValueError,
-                "If 'op_axes' or 'itershape' is not NULL in the"
-                "iterator constructor, 'oa_ndim' must be greater than zero");
-        return 0;
-    }
-    else if (oa_ndim > 0) {
-        if (oa_ndim > NPY_MAXDIMS) {
+    if (oa_ndim < 0) {
+        /*
+         * If `oa_ndim < 0`, `op_axes` and `itershape` are signalled to
+         * be unused and should be NULL. (Before NumPy 1.8 this was
+         * signalled by `oa_ndim == 0`.)
+         */
+        if (op_axes != NULL || itershape != NULL) {
             PyErr_Format(PyExc_ValueError,
+                    "If 'op_axes' or 'itershape' is not NULL in the iterator "
+                    "constructor, 'oa_ndim' must be zero or greater");
+            return 0;
+        }
+        return 1;
+    }
+    if (oa_ndim > NPY_MAXDIMS) {
+        PyErr_Format(PyExc_ValueError,
                 "Cannot construct an iterator with more than %d dimensions "
                 "(%d were requested for op_axes)",
                 (int)NPY_MAXDIMS, oa_ndim);
-            return 0;
-        }
-        else if (op_axes == NULL) {
-            PyErr_Format(PyExc_ValueError,
-                    "If 'oa_ndim' is greater than zero in the iterator "
-                    "constructor, then op_axes cannot be NULL");
-            return 0;
-        }
+        return 0;
+    }
+    if (op_axes == NULL) {
+        PyErr_Format(PyExc_ValueError,
+                "If 'oa_ndim' is zero or greater in the iterator "
+                "constructor, then op_axes cannot be NULL");
+        return 0;
+    }
 
-        /* Check that there are no duplicates in op_axes */
-        for (iop = 0; iop < nop; ++iop) {
-            int *axes = op_axes[iop];
-            if (axes != NULL) {
-                memset(axes_dupcheck, 0, NPY_MAXDIMS);
-                for (idim = 0; idim < oa_ndim; ++idim) {
-                    npy_intp i = axes[idim];
-                    if (i >= 0) {
-                        if (i >= NPY_MAXDIMS) {
-                            PyErr_Format(PyExc_ValueError,
-                                    "The 'op_axes' provided to the iterator "
-                                    "constructor for operand %d "
-                                    "contained invalid "
-                                    "values %d", (int)iop, (int)i);
-                            return 0;
-                        } else if(axes_dupcheck[i] == 1) {
-                            PyErr_Format(PyExc_ValueError,
-                                    "The 'op_axes' provided to the iterator "
-                                    "constructor for operand %d "
-                                    "contained duplicate "
-                                    "value %d", (int)iop, (int)i);
-                            return 0;
-                        }
-                        else {
-                            axes_dupcheck[i] = 1;
-                        }
+    /* Check that there are no duplicates in op_axes */
+    for (iop = 0; iop < nop; ++iop) {
+        int *axes = op_axes[iop];
+        if (axes != NULL) {
+            memset(axes_dupcheck, 0, NPY_MAXDIMS);
+            for (idim = 0; idim < oa_ndim; ++idim) {
+                npy_intp i = axes[idim];
+                if (i >= 0) {
+                    if (i >= NPY_MAXDIMS) {
+                        PyErr_Format(PyExc_ValueError,
+                                "The 'op_axes' provided to the iterator "
+                                "constructor for operand %d "
+                                "contained invalid "
+                                "values %d", (int)iop, (int)i);
+                        return 0;
+                    }
+                    else if (axes_dupcheck[i] == 1) {
+                        PyErr_Format(PyExc_ValueError,
+                                "The 'op_axes' provided to the iterator "
+                                "constructor for operand %d "
+                                "contained duplicate "
+                                "value %d", (int)iop, (int)i);
+                        return 0;
+                    }
+                    else {
+                        axes_dupcheck[i] = 1;
                     }
                 }
             }
@@ -819,7 +832,7 @@ npyiter_calculate_ndim(int nop, PyArrayObject **op_in,
                        int oa_ndim)
 {
     /* If 'op_axes' is being used, force 'ndim' */
-    if (oa_ndim > 0 ) {
+    if (oa_ndim >= 0 ) {
         return oa_ndim;
     }
     /* Otherwise it's the maximum 'ndim' from the operands */
@@ -1439,8 +1452,7 @@ static int
 npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, npyiter_opitflags *op_itflags,
                     char **op_dataptr,
                     npy_uint32 *op_flags, int **op_axes,
-                    npy_intp *itershape,
-                    int output_scalars)
+                    npy_intp *itershape)
 {
     npy_uint32 itflags = NIT_ITFLAGS(iter);
     int idim, ndim = NIT_NDIM(iter);
@@ -1540,6 +1552,13 @@ npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, npyiter_opitflags *op_itf
     axisdata = NIT_AXISDATA(iter);
     sizeof_axisdata = NIT_AXISDATA_SIZEOF(itflags, ndim, nop);
 
+    if (ndim == 0) {
+        /* Need to fill the first axisdata, even if the iterator is 0-d */
+        NAD_SHAPE(axisdata) = 1;
+        NAD_INDEX(axisdata) = 0;
+        memcpy(NAD_PTRS(axisdata), op_dataptr, NPY_SIZEOF_INTP*nop);
+    }
+
     /* Now process the operands, filling in the axisdata */
     for (idim = 0; idim < ndim; ++idim) {
         npy_intp bshape = broadcast_shape[ndim-idim-1];
@@ -1560,7 +1579,7 @@ npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, npyiter_opitflags *op_itf
                     ondim = PyArray_NDIM(op_cur);
                     if (bshape == 1) {
                         strides[iop] = 0;
-                        if (idim >= ondim && !output_scalars &&
+                        if (idim >= ondim &&
                                     (op_flags[iop] & NPY_ITER_NO_BROADCAST)) {
                             goto operand_different_than_broadcast;
                         }
@@ -1681,8 +1700,8 @@ npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, npyiter_opitflags *op_itf
     }
 
     /* Now fill in the ITERSIZE member */
-    NIT_ITERSIZE(iter) = broadcast_shape[0];
-    for (idim = 1; idim < ndim; ++idim) {
+    NIT_ITERSIZE(iter) = 1;
+    for (idim = 0; idim < ndim; ++idim) {
         NIT_ITERSIZE(iter) *= broadcast_shape[idim];
     }
     /* The range defaults to everything */
@@ -2003,7 +2022,10 @@ npyiter_replace_axisdata(NpyIter *iter, int iop,
     NIT_RESETDATAPTR(iter)[iop] = op_dataptr;
     NIT_BASEOFFSETS(iter)[iop] = baseoffset;
     axisdata = axisdata0;
-    for (idim = 0; idim < ndim; ++idim, NIT_ADVANCE_AXISDATA(axisdata, 1)) {
+    /* Fill at least one axisdata, for the 0-d case */
+    NAD_PTRS(axisdata)[iop] = op_dataptr;
+    NIT_ADVANCE_AXISDATA(axisdata, 1);
+    for (idim = 1; idim < ndim; ++idim, NIT_ADVANCE_AXISDATA(axisdata, 1)) {
         NAD_PTRS(axisdata)[iop] = op_dataptr;
     }
 }
@@ -2029,7 +2051,7 @@ npyiter_compute_index_strides(NpyIter *iter, npy_uint32 flags)
     /*
      * If there is only one element being iterated, we just have
      * to touch the first AXISDATA because nothing will ever be
-     * incremented.
+     * incremented. This also initializes the data for the 0-d case.
      */
     if (NIT_ITERSIZE(iter) == 1) {
         if (itflags & NPY_ITFLAG_HASINDEX) {
@@ -2399,7 +2421,7 @@ static PyArray_Descr *
 npyiter_get_common_dtype(int nop, PyArrayObject **op,
                         npyiter_opitflags *op_itflags, PyArray_Descr **op_dtype,
                         PyArray_Descr **op_request_dtypes,
-                        int only_inputs, int output_scalars)
+                        int only_inputs)
 {
     int iop;
     npy_intp narrs = 0, ndtypes = 0;
@@ -2698,7 +2720,7 @@ npyiter_allocate_arrays(NpyIter *iter,
                         npy_uint32 flags,
                         PyArray_Descr **op_dtype, PyTypeObject *subtype,
                         npy_uint32 *op_flags, npyiter_opitflags *op_itflags,
-                        int **op_axes, int output_scalars)
+                        int **op_axes)
 {
     npy_uint32 itflags = NIT_ITFLAGS(iter);
     int idim, ndim = NIT_NDIM(iter);
@@ -2729,7 +2751,7 @@ npyiter_allocate_arrays(NpyIter *iter,
         if (op[iop] == NULL) {
             PyArrayObject *out;
             PyTypeObject *op_subtype;
-            int ondim = output_scalars ? 0 : ndim;
+            int ondim = ndim;
 
             /* Check whether the subtype was disabled */
             op_subtype = (op_flags[iop] & NPY_ITER_NO_SUBTYPE) ?
@@ -2902,7 +2924,7 @@ npyiter_allocate_arrays(NpyIter *iter,
         if ((itflags & NPY_ITFLAG_BUFFER) &&
                                 !(op_itflags[iop] & NPY_OP_ITFLAG_CAST)) {
             NpyIter_AxisData *axisdata = NIT_AXISDATA(iter);
-            if (ndim == 1) {
+            if (ndim <= 1) {
                 op_itflags[iop] |= NPY_OP_ITFLAG_BUFNEVER;
                 NBF_STRIDES(bufferdata)[iop] = NAD_STRIDES(axisdata)[iop];
             }
diff --git a/numpy/core/src/multiarray/nditer_impl.h b/numpy/core/src/multiarray/nditer_impl.h
index 1251baa6e..ae24f46e6 100644
--- a/numpy/core/src/multiarray/nditer_impl.h
+++ b/numpy/core/src/multiarray/nditer_impl.h
@@ -294,7 +294,7 @@ struct NpyIter_AD {
 #define NIT_SIZEOF_ITERATOR(itflags, ndim, nop) ( \
         sizeof(struct NpyIter_InternalOnly) + \
         NIT_AXISDATA_OFFSET(itflags, ndim, nop) + \
-        NIT_AXISDATA_SIZEOF(itflags, ndim, nop)*(ndim))
+        NIT_AXISDATA_SIZEOF(itflags, ndim, nop)*(ndim ? ndim : 1))
 
 /* Internal helper functions shared between implementation files */
 NPY_NO_EXPORT void
diff --git a/numpy/core/src/multiarray/nditer_pywrap.c b/numpy/core/src/multiarray/nditer_pywrap.c
index 4621491a3..61f0c42b6 100644
--- a/numpy/core/src/multiarray/nditer_pywrap.c
+++ b/numpy/core/src/multiarray/nditer_pywrap.c
@@ -95,7 +95,6 @@ NpyIter_GlobalFlagsConverter(PyObject *flags_in, npy_uint32 *flags)
     npy_uint32 flag;
 
     if (flags_in == NULL || flags_in == Py_None) {
-        *flags = 0;
         return 1;
     }
 
@@ -526,7 +525,7 @@ npyiter_convert_op_axes(PyObject *op_axes_in, npy_intp nop,
         return 0;
     }
 
-    *oa_ndim = 0;
+    *oa_ndim = -1;
 
     /* Copy the tuples into op_axes */
     for (iop = 0; iop < nop; ++iop) {
@@ -545,13 +544,8 @@ npyiter_convert_op_axes(PyObject *op_axes_in, npy_intp nop,
                 Py_DECREF(a);
                 return 0;
             }
-            if (*oa_ndim == 0) {
+            if (*oa_ndim == -1) {
                 *oa_ndim = PySequence_Size(a);
-                if (*oa_ndim == 0) {
-                    PyErr_SetString(PyExc_ValueError,
-                            "op_axes must have at least one dimension");
-                    return 0;
-                }
                 if (*oa_ndim > NPY_MAXDIMS) {
                     PyErr_SetString(PyExc_ValueError,
                             "Too many dimensions in op_axes");
@@ -575,7 +569,7 @@ npyiter_convert_op_axes(PyObject *op_axes_in, npy_intp nop,
                     op_axes[iop][idim] = -1;
                 }
                 else {
-                    op_axes[iop][idim] = PyInt_AsLong(v);
+                    op_axes[iop][idim] = PyArray_PyIntAsInt(v);
                     if (op_axes[iop][idim]==-1 &&
                                                 PyErr_Occurred()) {
                         Py_DECREF(a);
@@ -589,7 +583,7 @@ npyiter_convert_op_axes(PyObject *op_axes_in, npy_intp nop,
         }
     }
 
-    if (*oa_ndim == 0) {
+    if (*oa_ndim == -1) {
         PyErr_SetString(PyExc_ValueError,
                 "If op_axes is provided, at least one list of axes "
                 "must be contained within it");
@@ -726,7 +720,7 @@ npyiter_init(NewNpyArrayIterObject *self, PyObject *args, PyObject *kwds)
     NPY_CASTING casting = NPY_SAFE_CASTING;
     npy_uint32 op_flags[NPY_MAXARGS];
     PyArray_Descr *op_request_dtypes[NPY_MAXARGS];
-    int oa_ndim = 0;
+    int oa_ndim = -1;
     int op_axes_arrays[NPY_MAXARGS][NPY_MAXDIMS];
     int *op_axes[NPY_MAXARGS];
     PyArray_Dims itershape = {NULL, 0};
@@ -784,7 +778,7 @@ npyiter_init(NewNpyArrayIterObject *self, PyObject *args, PyObject *kwds)
     }
 
     if (itershape.len > 0) {
-        if (oa_ndim == 0) {
+        if (oa_ndim == -1) {
             oa_ndim = itershape.len;
             memset(op_axes, 0, sizeof(op_axes[0]) * nop);
         }
@@ -800,10 +794,9 @@ npyiter_init(NewNpyArrayIterObject *self, PyObject *args, PyObject *kwds)
         itershape.ptr = NULL;
     }
 
-
     self->iter = NpyIter_AdvancedNew(nop, op, flags, order, casting, op_flags,
                                   op_request_dtypes,
-                                  oa_ndim, oa_ndim > 0 ? op_axes : NULL,
+                                  oa_ndim, oa_ndim >= 0 ? op_axes : NULL,
                                   itershape.ptr,
                                   buffersize);
 
@@ -860,7 +853,7 @@ NpyIter_NestedIters(PyObject *NPY_UNUSED(self),
 
     int iop, nop = 0, inest, nnest = 0;
     PyArrayObject *op[NPY_MAXARGS];
-    npy_uint32 flags = 0, flags_inner = 0;
+    npy_uint32 flags = 0, flags_inner;
     NPY_ORDER order = NPY_KEEPORDER;
     NPY_CASTING casting = NPY_SAFE_CASTING;
     npy_uint32 op_flags[NPY_MAXARGS], op_flags_inner[NPY_MAXARGS];
diff --git a/numpy/core/src/private/lowlevel_strided_loops.h b/numpy/core/src/private/lowlevel_strided_loops.h
index 94c6a2121..c9fd1248f 100644
--- a/numpy/core/src/private/lowlevel_strided_loops.h
+++ b/numpy/core/src/private/lowlevel_strided_loops.h
@@ -256,6 +256,7 @@ PyArray_CastRawArrays(npy_intp count,
  * 'stransfer' with the provided dst_stride/src_stride and
  * dst_strides[0]/src_strides[0], so the caller can use those values to
  * specialize the function.
+ * Note that even if ndim == 0, everything needs to be set as if ndim == 1.
  *
  * The return value is the number of elements it couldn't copy.  A return value
  * of 0 means all elements were copied, a larger value means the end of
diff --git a/numpy/core/src/umath/reduction.c b/numpy/core/src/umath/reduction.c
index e6ed04e99..f69aea2d0 100644
--- a/numpy/core/src/umath/reduction.c
+++ b/numpy/core/src/umath/reduction.c
@@ -513,7 +513,7 @@ PyUFunc_ReduceWrapper(PyArrayObject *operand, PyArrayObject *out,
                                NPY_KEEPORDER, casting,
                                op_flags,
                                op_dtypes,
-                               0, NULL, NULL, buffersize);
+                               -1, NULL, NULL, buffersize);
     if (iter == NULL) {
         goto fail;
     }
diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c
index 124185bfd..9c499d322 100644
--- a/numpy/core/src/umath/ufunc_object.c
+++ b/numpy/core/src/umath/ufunc_object.c
@@ -1211,7 +1211,7 @@ iterator_loop(PyUFuncObject *ufunc,
                         NPY_ITER_DELAY_BUFALLOC,
                         order, NPY_UNSAFE_CASTING,
                         op_flags, dtype,
-                        0, NULL, NULL, buffersize);
+                        -1, NULL, NULL, buffersize);
     if (iter == NULL) {
         return -1;
     }
@@ -1509,7 +1509,7 @@ execute_fancy_ufunc_loop(PyUFuncObject *ufunc,
                         NPY_ITER_GROWINNER,
                         order, NPY_UNSAFE_CASTING,
                         op_flags, dtypes,
-                        0, NULL, NULL, buffersize);
+                        -1, NULL, NULL, buffersize);
     if (iter == NULL) {
         return -1;
     }
@@ -1976,18 +1976,6 @@ PyUFunc_GeneralizedFunction(PyUFuncObject *ufunc,
                       NPY_ITER_NO_BROADCAST;
     }
 
-    /*
-     * If there are no iteration dimensions, create a fake one
-     * so that the scalar edge case works right.
-     */
-    if (iter_ndim == 0) {
-        iter_ndim = 1;
-        iter_shape[0] = 1;
-        for (i = 0; i < nop; ++i) {
-            op_axes[i][0] = -1;
-        }
-    }
-
     /* Create the iterator */
     iter = NpyIter_AdvancedNew(nop, op, NPY_ITER_MULTI_INDEX|
                                       NPY_ITER_REFS_OK|
diff --git a/numpy/core/tests/test_einsum.py b/numpy/core/tests/test_einsum.py
index 4fba533db..a60d9a415 100644
--- a/numpy/core/tests/test_einsum.py
+++ b/numpy/core/tests/test_einsum.py
@@ -241,6 +241,7 @@ class TestEinSum(TestCase):
             assert_equal(np.einsum(a, [0,0]), np.trace(a).astype(dtype))
 
         # multiply(a, b)
+        assert_equal(np.einsum("..., ...", 3, 4), 12) # scalar case
         for n in range(1,17):
             a = np.arange(3*n, dtype=dtype).reshape(3,n)
             b = np.arange(2*3*n, dtype=dtype).reshape(2,3,n)
diff --git a/numpy/core/tests/test_nditer.py b/numpy/core/tests/test_nditer.py
index 219629191..dd0ed8cfa 100644
--- a/numpy/core/tests/test_nditer.py
+++ b/numpy/core/tests/test_nditer.py
@@ -2472,5 +2472,59 @@ def test_iter_allocated_array_dtypes():
         c[1,1] = a / b
     assert_equal(it.operands[2], [[8, 12], [20, 5]])
 
+
+def test_0d_iter():
+    # Basic test for iteration of 0-d arrays:
+    i = nditer([2, 3], ['multi_index'], [['readonly']]*2)
+    assert_equal(i.ndim, 0)
+    assert_equal(i.next(), (2, 3))
+    assert_equal(i.multi_index, ())
+    assert_equal(i.iterindex, 0)
+    assert_raises(StopIteration, i.next)
+    # test reset:
+    i.reset()
+    assert_equal(i.next(), (2, 3))
+    assert_raises(StopIteration, i.next)
+
+    # test forcing to 0-d
+    i = nditer(np.arange(5), ['multi_index'], [['readonly']], op_axes=[()])
+    assert_equal(i.ndim, 0)
+    assert_equal(len(i), 1)
+    # note that itershape=(), still behaves like None due to the conversions
+
+    # Test a more complex buffered casting case (same as another test above)
+    sdt = [('a', 'f4'), ('b', 'i8'), ('c', 'c8', (2,3)), ('d', 'O')]
+    a = np.array(0.5, dtype='f4')
+    i = nditer(a, ['buffered','refs_ok'], ['readonly'],
+                    casting='unsafe', op_dtypes=sdt)
+    vals = i.next()
+    assert_equal(vals['a'], 0.5)
+    assert_equal(vals['b'], 0)
+    assert_equal(vals['c'], [[(0.5)]*3]*2)
+    assert_equal(vals['d'], 0.5)
+
+
+def test_0d_nested_iter():
+    a = np.arange(12).reshape(2,3,2)
+    i, j = np.nested_iters(a, [[],[1,0,2]])
+    vals = []
+    for x in i:
+        vals.append([y for y in j])
+    assert_equal(vals, [[0,1,2,3,4,5,6,7,8,9,10,11]])
+
+    i, j = np.nested_iters(a, [[1,0,2],[]])
+    vals = []
+    for x in i:
+        vals.append([y for y in j])
+    assert_equal(vals, [[0],[1],[2],[3],[4],[5],[6],[7],[8],[9],[10],[11]])
+
+    i, j, k = np.nested_iters(a, [[2,0], [] ,[1]])
+    vals = []
+    for x in i:
+        for y in j:
+            vals.append([z for z in k])
+    assert_equal(vals, [[0,2,4],[1,3,5],[6,8,10],[7,9,11]])
+
+
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/lib/index_tricks.py b/numpy/lib/index_tricks.py
index 07e7c1e39..dc57c1048 100644
--- a/numpy/lib/index_tricks.py
+++ b/numpy/lib/index_tricks.py
@@ -532,30 +532,7 @@ class ndindex(object):
     (2, 0, 0)
     (2, 1, 0)
 
-    """
-    # This is a hack to handle 0-d arrays correctly.
-    # Fixing nditer would be more work but should be done eventually,
-    #  and then this entire __new__ method can be removed.
-    def __new__(cls, *shape):
-        if len(shape) == 1 and isinstance(shape[0], tuple):
-           shape = shape[0]
-        if len(shape) == 0:
-            class zero_dim_iter(object):
-                def __init__(self):
-                    self._N = 1
-                def __iter__(self):
-                    return self
-                def ndincr(self):
-                    self.next()
-                def next(self):
-                    if self._N > 0:
-                        self._N -= 1
-                        return ()
-                    raise StopIteration
-            return zero_dim_iter()
-        else:
-            return super(ndindex, cls).__new__(cls)
-            
+    """     
     def __init__(self, *shape):
         if len(shape) == 1 and isinstance(shape[0], tuple):
             shape = shape[0]