Describe the issue:

np.random.default_rng(0).poisson(lam) (as well as the legacy API
np.random.poisson) draws samples whose empirical variance is substantially
larger than the theoretical variance (which should be lam), once lam is large.
The sample mean stays correct, but var(x) / lam grows well above 1.0.

Measured var/lam on current main / 2.4.4 (size = 2_000_000):

Root cause: the transformed-rejection sampler random_poisson_ptrs
(numpy/random/src/distributions/distributions.c,
around L600) evaluates the log pmf: -lam + k*loglam - loggam(k+1), in
float64.
For example: for promising candidates (that is, those that should be accepted at
L600) k ≈ lam at lam = 1e17 the above three terms are ~1e17, ~3.91e18,
and ~3.81e18. They nearly cancel to an exact result of about ~ -20. But a
float64 computation returns 0.0 here due to rounding error, thereby never accepting the promising candidate k.

Reproduce the code example:

import numpy as np
 
lam = 1e17
size = 2_000_000
rng = np.random.default_rng()
x = rng.poisson(lam, size=size)
variance = x.var()
ratio = variance / lam

print(f"sample mean       = {x.mean()}")
print(f"variance / lambda = {ratio}

Error message:

No exception. The sampler silently returns samples from a wrong distribution.

Python and NumPy Versions:

NumPy: 2.5.0
Python: 3.13.2

Runtime Environment:

[{'numpy_version': '2.5.0',
  'python': '3.13.2 (main, Feb  4 2025, 14:51:09) [Clang 16.0.0 '
            '(clang-1600.0.26.6)]',
  'uname': uname_result(system='Darwin', node='infra04-wg313.visitor.eduroam.wireless.umn.edu', release='24.6.0', version='Darwin Kernel Version 24.6.0: Mon Jan 19 21:58:37 PST 2026; root:xnu-11417.140.69.708.3~1/RELEASE_ARM64_T8103', machine='arm64')},
 {'simd_extensions': {'baseline': ['NEON', 'NEON_FP16', 'NEON_VFPV4', 'ASIMD'],
                      'found': ['ASIMDHP', 'ASIMDDP'],
                      'not_found': ['ASIMDFHM']}},
 {'ignore_floating_point_errors_in_matmul': True}]

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

We use np.random.poisson as a backend in our model counting tool:
[pepin](https://github.com/meelgroup/pepin). For large model counts we received
wrong counts, which we traced back to this variance inflation in the Poisson
sampler at large lam.