Issue 1:
From https://stackoverflow.com/questions/68927829/how-to-use-use-numpy-random-choice-to-get-progressively-longer-sequences-with-th

When generating random numbers with np.random.Generator.choice and having probabilities, the chosen elements in the first positions will depend on the number of elements to generate, leading to inconsistent results if one tries to replicate results but with varying input lengths:

import numpy as np
n = 100
prob = np.random.choice(100, n)
prob = prob / np.sum(prob)
for i in np.arange(10, 40, 10):
rng = np.random.default_rng(seed=123)
selected = rng.choice(n, size=i, replace=False, p=prob)
print(str(i) + ': ' + str(selected))

10: [71 5 20 16 81 92 26 82 88 52]
20: [71 5 20 16 81 92 26 82 88 52 24 19 75 65 22 80 53 25 13 48]
30: [71 5 20 16 81 92 26 82 88 52 24 19 75 65 22 80 53 15 50 60 1 48 74 91 94 87 18 73 31 79]

Which is quite unexpected and does not happen in e.g. R's sample.

I would assume the optimal algorithm for drawing these random numbers amounts to having some binary tree saving weights of its left and right subtrees, and then drawing random numbers in the tree's range and passing down the tree until reaching the bottom and removing that weight, which should not make any choice at position n depend on what's going to happen at position n+1.

Issue 2:
In the FFTPACK code, we have some "clever" code that uses the pattern NPY_SIGINT_ON, then run some inner loop, then NPY_SIGINT_OFF. And what this does is:

NPY_SIGINT_ON uses setjmp to make a jump buffer target, saves the current signal handler for SIGINT (control-C), and then registers a new SIGINT handler that attempts to handle control-C by longjmping out to the setjmp buffer. The idea is that this allows us to cancel long-running calculations.

The actual implementations of these macros is:

# define NPY_SIGINT_ON { \
PyOS_sighandler_t _npy_sig_save; \
_npy_sig_save = PyOS_setsig(SIGINT, _PyArray_SigintHandler);\
\
if (NPY_SIGSETJMP(*((NPY_SIGJMP_BUF *)_PyArray_GetSigintBuf(\
)), \
1) == 0) {

# define NPY_SIGINT_OFF } \
PyOS_setsig(SIGINT, _npy_sig_save); \
}

So this has two race conditions:

Minor problem: If we receive a control-C in between the call to PyOS_setsig and the call to NPY_SIGSETJMP, then we'll longjmp out of the signal handler into a uninitialized buffer, and will not go to space today.

Major problem: suppose that the following sequence occurs:

thread 1 enters NPY_SIGINT_ON, stashes Python's default sigint handler in its _npy_sig_save

thread 2 enters NPY_SIGINT_ON, stashes our sigint handler in its _npy_sig_save

thread 1 leaves via NPY_SIGINT_OFF, restores Python's default sigint handler from its _npy_sig_save

thread 2 leaves via NPY_SIGINT_OFF, then restores our sigint handler from its _npy_sig_save

Now our signal handler gets left installed indefinitely, and eventually when control-C gets hit it will attempt to jump out to a stack frame that was reused ages ago. This is a bad problem and you will not go to space today.

By adding some printf checks to NPY_SIGINT_OFF, I've confirmed that this is the cause of the segfault that Oscar reported on numpy-discussion. It's very easy to replicate even from Python code:

import numpy
numpy.test()
import os
import signal
os.kill(os.getpid(), signal.SIGINT)
# -> segfault

Also, the whole architecture here is somewhat screwed up, because the Python model is that all signals are directed to the main thread (by setting them to masked on all other threads). So it doesn't even make sense to be installing a signal handler to interrupt operations on a thread that can't receive signals.

(If you look at the definition of _PyArray_SigintHandler and _PyArray_GetSigintBuf in multiarraymodule.c, then you'll see that there's some basic defenses against this: we keep our jump buffer in TLS, and also a flag sigint_buf_init, so that idea is that if we start running the signal handler in a thread that hasn't actually allocated a jump buffer, then we throw away the signal. This is also pretty buggy: (1) nothing ever sets the flag back to 0, so if a thread has ever had a valid jump buffer on its stack then we will happily jump to it, even if it's no longer valid. (2) we shouldn't throw away the signal, we should pass it on to Python.)

So the simple solution is (1) check whether we are actually the one thread that can receive signals, (2) if so, then register our handler etc.; (3) if not, then make these macros into no-ops.

The one wrinkle here is that in theory it's possible for someone to call sigprocmask themselves and enable signal receiving on a non-default thread. We might want to be robust to this. A nice extra constraint is that there's no C API to check whether we are in the main thread; the only way I can think of to check that is actually to call sigprocmask and check whether SIGINT delivery is enabled for this thread.