Issue 1:
I noticed that in many of my codes, seemingly harmless lines like

numpy.add.at(target, idx, vals)

take a large share of the runtime. I investigated and found that one gets a speed-up of a factor of 40 (!) by simply moving the critical code to C++.

(Not sure what's actually done in numpy.ufunc.) I put this into a very simple module, https://github.com/nschloe/fastfunc, so feel free to take some code out of there.

Issue 2:
numpy 1.18.4 and 1.18.5 test errors

run numpy.test()

_______________ TestRandomDist.test_multivariate_normal[eigh] _________________

self = <numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d19bcc40>
method = 'eigh'

@pytest.mark.parametrize("method", ["svd", "eigh", "cholesky"])
def test_multivariate_normal(self, method):
random = Generator(MT19937(self.seed))
mean = (.123456789, 10)
cov = [[1, 0], [0, 1]]
size = (3, 2)
actual = random.multivariate_normal(mean, cov, size, method=method)
desired = np.array([[[-1.747478062846581, 11.25613495182354 ],
[-0.9967333370066214, 10.342002097029821 ]],
[[ 0.7850019631242964, 11.181113712443013 ],
[ 0.8901349653255224, 8.873825399642492 ]],
[[ 0.7130260107430003, 9.551628690083056 ],
[ 0.7127098726541128, 11.991709234143173 ]]])

> assert_array_almost_equal(actual, desired, decimal=15)
E AssertionError:
E Arrays are not almost equal to 15 decimals
E
E Mismatched elements: 6 / 12 (50%)
E Max absolute difference: 3.7418697
E Max relative difference: 2.2477228
E x: array([[[ 1.994391640846581, 11.25613495182354 ],
E [ 1.243646915006621, 10.342002097029821]],
E ...
E y: array([[[-1.747478062846581, 11.25613495182354 ],
E [-0.996733337006621, 10.342002097029821]],
E ...

actual = array([[[ 1.99439164, 11.25613495],
[ 1.24364692, 10.3420021 ]],

[[-0.53808839, 11.18111371],
[-0.64322139, 8.8738254 ]],

[[-0.46611243, 9.55162869],
[-0.46579629, 11.99170923]]])
cov = [[1, 0], [0, 1]]
desired = array([[[-1.74747806, 11.25613495],
[-0.99673334, 10.3420021 ]],

[[ 0.78500196, 11.18111371],
[ 0.89013497, 8.8738254 ]],

[[ 0.71302601, 9.55162869],
[ 0.71270987, 11.99170923]]])
mean = (0.123456789, 10)
method = 'eigh'
random = Generator(MT19937) at 0x7F77D1CD1D60
self = <numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d19bcc40>
size = (3, 2)

lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py:1264: AssertionError
________ TestRandomDist.test_multivariate_normal_basic_stats[cholesky] _________

self = <numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d18b54f0>
method = 'cholesky'

@pytest.mark.parametrize("method", ["svd", "eigh", "cholesky"])
def test_multivariate_normal_basic_stats(self, method):
random = Generator(MT19937(self.seed))
n_s = 1000
mean = np.array([1, 2])
cov = np.array([[2, 1], [1, 2]])
s = random.multivariate_normal(mean, cov, size=(n_s,), method=method)
s_center = s - mean
cov_emp = (s_center.T @ s_center) / (n_s - 1)
# these are pretty loose and are only designed to detect major errors
assert np.all(np.abs(s_center.mean(-2)) < 0.1)
> assert np.all(np.abs(cov_emp - cov) < 0.2)
E AssertionError: assert False
E + where False = <function all at 0x7f7861f6b3a0>(array([[0.55018097, 0.10366964],\n [0.10366964, 0.51195134]]) < 0.2)
E + where <function all at 0x7f7861f6b3a0> = np.all
E + and array([[0.55018097, 0.10366964],\n [0.10366964, 0.51195134]]) = <ufunc 'absolute'>((array([[2.55018097, 0.89633036],\n [0.89633036, 1.48804866]]) - array([[2, 1],\n [1, 2]])))
E + where <ufunc 'absolute'> = np.abs

cov = array([[2, 1],
[1, 2]])
cov_emp = array([[2.55018097, 0.89633036],
[0.89633036, 1.48804866]])
mean = array([1, 2])
method = 'cholesky'
n_s = 1000
random = Generator(MT19937) at 0x7F77D1CD1A90
s = array([[-0.7576799 , 3.53844484],
[-0.34235607, 2.41886531],
[ 2.77073967, 3.44656296],
...,
[-0.51773839, 0.49648485],
[ 0.57610204, 3.15188553],
[ 1.10139309, 1.47737657]])
s_center = array([[-1.7576799 , 1.53844484],
[-1.34235607, 0.41886531],
[ 1.77073967, 1.44656296],
...,
[-1.51773839, -1.50351515],
[-0.42389796, 1.15188553],
[ 0.10139309, -0.52262343]])
self = <numpy.random.tests.test_generator_mt19937.TestRandomDist object at 0x7f77d18b54f0>

lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py:1335: AssertionError
=========================== short test summary info ============================
FAILED lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py::TestRandomDist::test_multivariate_normal[eigh]
FAILED lib/python3.8/site-packages/numpy/random/tests/test_generator_mt19937.py::TestRandomDist::test_multivariate_normal_basic_stats[cholesky]
2 failed, 10264 passed, 68 skipped, 180 deselected, 17 xfailed, 3 xpassed in 69.53s (0:01:09)
False

OS is ubuntu, same failure on python-3.7.6 and 3.8.2.

numpy installed via "pip3 install --no-binary numy numpy" compiled from source with MKL.

But same test failures persist if link numpy to OpenBLAS (version 0.3.9)

numpy-1.17.5 (installed from source the same way) produces no test error in both.

>>> import numpy as np
>>> np.show_config()
blas_mkl_info:
libraries = ['mkl_rt', 'pthread']
library_dirs = ['/opt/intel/mkl/lib/intel64/']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['/opt/intel/mkl/include']
blas_opt_info:
libraries = ['mkl_rt', 'pthread']
library_dirs = ['/opt/intel/mkl/lib/intel64/']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['/opt/intel/mkl/include']
lapack_mkl_info:
libraries = ['mkl_rt', 'pthread']
library_dirs = ['/opt/intel/mkl/lib/intel64/']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['/opt/intel/mkl/include']
lapack_opt_info:
libraries = ['mkl_rt', 'pthread']
library_dirs = ['/opt/intel/mkl/lib/intel64/']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['/opt/intel/mkl/include']

One of the failures is the same as #15546 for numpy-1.18.1

but this time it fails also in python 3.7.6