BinomialFunction: wrong samples with Cython for p > 0.5, and non-terminating generated code with approximate=False #1890

When the normal approximation is not used (n*p <= 5 or n*(1-p) <= 5, or approximate=False), BinomialFunction generates code for the sequential inversion algorithm taken from numpy: it draws X ~ B(n, P) with P = min(p, 1-p), starting from qn = (1-P)**n, and returns n - X for p > 0.5. The generated code has three problems.

1. Cython only, wrong distribution for p > 0.5. The Cython template receives the original p instead of the reduced P (binomial.py#L72; the C++ template uses P, L133). For p > 0.5, q = 1 - P = p, so the factor p/q that takes the place of P/q in the recurrence equals 1, the "probabilities" sum to (2p)**n, and the samples pile up near n. This is the default path (approximate=True) whenever p > 0.5 and n*(1-p) <= 5:

from brian2 import *
prefs.codegen.target = "cython"
seed(20261010)
bino = BinomialFunction(10, 0.8)  # n*(1-p) = 2 <= 5: exact (inversion) branch
G = NeuronGroup(100000, "x : 1")
G.run_regularly("x = bino()")
run(defaultclock.dt)
print(G.x[:].mean(), G.x[:].var())  # B(10, 0.8): mean 8, variance 1.6

On Linux (gcc 13.3, Cython 3.1.3), with both 240e2f0 and the 2.10.1 release, this prints 9.10688 0.09545666560000002; with prefs.codegen.target = "numpy" it prints 8.00088 1.5984592256. C++ standalone samples this case correctly (mean 8.0009 over 100 000 samples). The substitution dates back to the first implementation (304621c, #428).

2. Cython and C++, constants rounded to 15 decimal places (L72-L75, L133-L136). With approximate=False and qn < 5e-16, the generated code starts from px = 0.000000000000000 and while U > px never terminates: BinomialFunction(100, 0.3, approximate=False) (qn = 3.2e-16) did not finish within 120 s on either target, while numpy samples it normally. Just above this threshold qn loses most of its digits: for BinomialFunction(50, 0.5, approximate=False), qn = 2**-50 = 8.9e-16 is written as 0.000000000000001, and 100 000 samples have mean 24.24 and variance 8.70 instead of 25 and 12.5, on both targets. With approximate=True the inversion is only used when n*P <= 5, so qn >= 0.5**10 and the rounding is negligible.

3. Even when written exactly, qn underflows in double precision for large n*P: for p = 0.5 it is subnormal, with reduced precision, from n = 1023, and zero from n = 1075; with qn = 0 the loop of problem 2 does not terminate either. For example, BinomialFunction(2000, 0.5, approximate=False) has qn = 0 and did not finish within 120 s on either target. numpy avoids this regime by switching to BTPE above n*min(p, 1-p) = 30; the templates have no such alternative.

A patch with tests is ready and I will open a PR referencing this issue. It uses P in the Cython template, writes the constants of the inversion algorithm with repr(float(x)) (as NodeRenderer.render_Constant does), and, when qn is not a normal double, registers an availability_check that raises a NotImplementedError at code generation with an explanation, so that the numpy target keeps working with the same parameters. The tests compare the sampled distribution with the exact pmf on numpy, Cython and C++ standalone.

For problem 3, is an explicit error acceptable to you, or would you prefer a fallback to the normal approximation, or a port of BTPE to the templates?

I used an AI assistant (Claude Code) to investigate the bug, run the reproductions and measurements, and draft this text.