Describe the issue:

np.unravel_index silently returns incorrect coordinate values when indices is a sufficiently large two-dimensional integer array with a trailing singleton dimension, such as shape (N, 1).

According to the documentation, indices may be an integer array, and every returned coordinate array should have the same shape as indices. Therefore, (N, 1) is a valid input shape and should be processed element by element.

For example, the following arrays contain the same flat-index values:

indices_1d = np.arange(10_000, dtype=np.intp)
indices_2d = indices_1d[:, None]

The results should satisfy:

coords_2d[i].ravel() == coords_1d[i]

for every coordinate dimension. Equivalently, applying np.ravel_multi_index to the result of np.unravel_index should recover the original indices.

Instead, with the (N, 1) input, coordinates become incorrect after an apparent internal buffer boundary. In the example below, values up to index 8192 are handled correctly, but later indices are mapped to repeated coordinates.

For example, with shape (1000, 1000):

flat index 8192 -> expected (8, 192), returned (8, 192)
flat index 8193 -> expected (8, 193), returned (8, 192)
flat index 8194 -> expected (8, 194), returned (8, 192)

The one-dimensional version of the same indices works correctly.

This is silent data corruption: the output arrays have the expected shape and no exception is raised, but some coordinate values are wrong.

I reproduced this behavior with NumPy 2.3.0, 2.3.5, and 2.5.1. The same example works correctly with NumPy 2.2.6.

This may be related to issue #29690 and PR #29780, which involved unravel_index, buffered NpyIter operation, and an axis of size one. However, that issue concerned a debug-build assertion. The problem reported here concerns silently incorrect returned values when the input crosses a buffer boundary.

Reproduce the code example:

import numpy as np

shape = (1000, 1000)
first_bad = 8193

# Same flat-index values, represented with two different input shapes.
indices_1d = np.arange(first_bad + 1, dtype=np.intp)
indices_2d = indices_1d[:, None].copy()

# Prove that the values are identical.
np.testing.assert_array_equal(indices_2d.ravel(), indices_1d)

expected_row_1d = indices_1d // shape[1]
expected_col_1d = indices_1d % shape[1]

expected_row_2d = expected_row_1d[:, None]
expected_col_2d = expected_col_1d[:, None]

# Control 1: one-dimensional input works correctly.
row_1d, col_1d = np.unravel_index(indices_1d, shape)

np.testing.assert_array_equal(row_1d, expected_row_1d)
np.testing.assert_array_equal(col_1d, expected_col_1d)

restored_1d = np.ravel_multi_index((row_1d, col_1d), shape)
np.testing.assert_array_equal(restored_1d, indices_1d)

# Affected case: the same values stored in an (N, 1) array.
row_2d, col_2d = np.unravel_index(indices_2d, shape)

# Control 2: independently expected coordinates are valid.
np.testing.assert_array_equal(
    np.ravel_multi_index((expected_row_2d, expected_col_2d), shape),
    indices_2d,
)

# Actual round trip using coordinates returned from np.unravel_index.
restored_2d = np.ravel_multi_index((row_2d, col_2d), shape)

mismatch = np.flatnonzero(restored_2d.ravel() != indices_1d)

print("NumPy version:", np.__version__)
print("1-D input shape:", indices_1d.shape)
print("2-D input shape:", indices_2d.shape)
print("Input dtype:", indices_1d.dtype)
print("2-D input strides:", indices_2d.strides)
print("2-D input C-contiguous:", indices_2d.flags.c_contiguous)
print("Number of round-trip mismatches for 2-D input:", mismatch.size)

if mismatch.size:
    pos = int(mismatch[0])

    print("First mismatch position:", pos)
    print("Flat index at that position:", int(indices_1d[pos]))

    print(
        "1-D returned coordinate:",
        int(row_1d[pos]),
        int(col_1d[pos]),
    )
    print(
        "2-D returned coordinate:",
        int(row_2d[pos, 0]),
        int(col_2d[pos, 0]),
    )
    print(
        "Expected coordinate:",
        int(expected_row_1d[pos]),
        int(expected_col_1d[pos]),
    )
    print("2-D round-tripped value:", int(restored_2d[pos, 0]))

# Directly compare the same values under different input shapes.
np.testing.assert_array_equal(row_2d.ravel(), row_1d)
np.testing.assert_array_equal(col_2d.ravel(), col_1d)

# Equivalent form of the same correctness requirement.
np.testing.assert_array_equal(row_2d, expected_row_2d)
np.testing.assert_array_equal(col_2d, expected_col_2d)
np.testing.assert_array_equal(restored_2d, indices_2d)

Error message:

NumPy version: 2.3.4
1-D input shape: (8194,)
2-D input shape: (8194, 1)
Input dtype: int64
2-D input strides: (8, 8)
2-D input C-contiguous: True
Number of round-trip mismatches for 2-D input: 1
First mismatch position: 8193
Flat index at that position: 8193
1-D returned coordinate: 8 193
2-D returned coordinate: 8 192
Expected coordinate: 8 193
2-D round-tripped value: 8192

AssertionError: Arrays are not equal

Mismatched elements: 1 / 8194 (0.0122%)
ACTUAL:   ... 191, 192, 192
DESIRED:  ... 191, 192, 193

Python and NumPy Versions:

'numpy_version': '2.3.4', 'python': '3.12.12 (main, Apr 4 2026, 15:58:07) [GCC 8.5.0 20210514 (Red ' 'Hat 8.5.0-28)

P.S. Additional version testing: I reproduced the issue across multiple Python versions. With NumPy 2.5.1, the reproducer reports the same 1,807 incorrect indices on CPython 3.12.12, 3.13.0, and 3.14.0. The first mismatch is consistently at flat index 8193: the expected coordinate is (8, 193), but np.unravel_index returns (8, 192), which maps back to 8192. I also reproduced the same behavior on CPython 3.11.15 with NumPy 2.4.6, the latest NumPy release compatible with Python 3.11. This therefore does not appear to be specific to a single Python version. NumPy 2.2.6 remains the last known good release, and NumPy 2.3.0 remains the first known bad release.

Runtime Environment:

[
{
'numpy_version': '2.3.4',
'python': '3.12.12 (main, Apr 4 2026, 15:58:07) [GCC 8.5.0 20210514 (Red '
'Hat 8.5.0-28)]',
'uname': uname_result(system='Linux',
node='',
release='4.18.0-553.117.1.el8_10.x86_64',
version='#1 SMP Wed Apr 8 10:49:10 UTC 2026',
machine='x86_64')
},
{
'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'],
'found': ['SSSE3',
'SSE41',
'POPCNT',
'SSE42',
'AVX',
'F16C',
'FMA3',
'AVX2'],
'not_found': ['AVX512F',
'AVX512CD',
'AVX512_KNL',
'AVX512_KNM',
'AVX512_SKX',
'AVX512_CLX',
'AVX512_CNL',
'AVX512_ICL',
'AVX512_SPR']
}
},
{
'architecture': 'Haswell',
'filepath': '/numpy.libs/libscipy_openblas64_-8fb3d286.so',
'internal_api': 'openblas',
'num_threads': 64,
'prefix': 'libscipy_openblas',
'threading_layer': 'pthreads',
'user_api': 'blas',
'version': '0.3.30'
}
]

How does this issue affect you or how did you find it:

This issue has a concrete downstream impact on MRI reconstruction software.

I encountered it while investigating why pygrappa, a Python implementation of the classical GRAPPA parallel MRI reconstruction algorithm, was leaving many missing k-space locations unfilled.

In affected cases, the reconstruction completes without raising an exception, and the acquired k-space samples remain unchanged. However, a substantial fraction of the missing k-space locations are never filled and remain exactly zero.

This produces incomplete GRAPPA reconstructions and can lead to residual aliasing or other image artifacts in the reconstructed MRI data. The problem is difficult to detect because the function returns arrays with the expected shape and emits no warning or error. Flattening the input indices before calling np.unravel_index restores complete target coverage in pygrappa.

I first observed the behavior on real multicoil MRI data, then reproduced it using fully synthetic pygrappa inputs, and finally reduced it to the self-contained NumPy-only reproducer included above.