fix: 64-bit randint could return high; correct the integer fills

[vchamarthi]

What

Four fixes in the integer fills, all in mkl_distributions.cpp.

1. randint could return high. In the masked branch of the 64-bit fill,
rng was incremented before the mask was derived from it:

rng = hi - lo;
rng++;                      /* now a count, not a maximum */
mask = smear_bits(rng);     /* one bit too wide */
if (value <= rng) accept;   /* admits lo + rng, which is hi + 1 */

Fixed by leaving rng as the inclusive maximum and using rng + 1 only where
viRngUniform wants an exclusive bound.

The same surplus mask bit also forced roughly half of all draws to be rejected
and redrawn whenever rng + 1 was a power of two. With the mask correct,
mask == rng means nothing can be rejected, so those requests now fill the
result buffer directly.

2. The integer fills handled at most two MKL_INT_MAX chunks. They used
if where the other 43 chunking sites in the same file use while, so a larger
request recursed once, fell through, and passed a count that overflows MKL's
32-bit int. NDEBUG compiles out the assert, so the error was swallowed and
the tail was never written.

3. multinomial decremented its chunk counter by elements, not draws
(len -= k * MKL_INT_MAX), skipping work and leaving the tail of large outputs
unwritten.

4. The fills narrower than int staged the whole request. viRngUniform
only emits 32-bit integers, so a uint8 fill allocated four bytes per element
to produce one, then read that buffer back to narrow it: 40 MB allocated to
write 10 MB of output. Now staged one cache-resident tile at a time.

Numbers

Xeon Gold 6338, single thread, ns/element, min of 2 runs. Intel-channel
mkl / mkl-devel 2026.1.0-intel_236, meson build. Baseline and patched
installed in separate conda envs.

case	before	after
uint64, power-of-two range 2**32	25.4	1.39	18.3x
uint64, power-of-two range 2**48	25.4	1.39	18.3x
uint64, iinfo.max upper bound	9.69	1.93	5.0x
bool, 10M	3.85	1.11	3.4x
uint8, 10M	3.79	1.11	3.4x
uint16, 10M	3.93	1.12	3.5x
uint32, 10M (control)	1.90	1.89	unchanged

Streams change

Only power-of-two ranges at or above INT_MAX change, and only because the old
mask there was wrong. Of 12 seeded outputs compared, exactly one hash moved:
randint(0, 2**32, dtype='uint64'). The other 11 are bit-identical, including
every narrow dtype, uint32, standard_normal and multinomial. The
narrow-int tiling is bit-identical by construction: a VSL stream is sequential,
so tiling draws the same values in the same order.

Testing

• Suite: 190 -> 209 passed, 3 skipped. The 19 new tests also pass on unpatched
  master, so they are invariant guards rather than fitted to this change.
• Out-of-range scan on unpatched master, randint(0, 2**31), 8e9 draws per
  dtype: 3 hits on default int, 3 on int64, 2 on uint64. After the fix, 0
  on all three.
• Verified on upstream master 6b7962f.

Notes for reviewers

• test_randint_in_bounds_fuzz only uses high in {4, 8, 16}, all of which
  take the 32-bit path, so the masked 64-bit branch had no coverage. The new
  test covers ranges at and above INT_MAX.
• The out-of-range failure fires about once per 2**32 draws, so the added test
  asserts the bound invariant rather than trying to reproduce it.
• Fixes 2 and 3 are reasoned from the source, not reproduced: they need requests
  above two MKL_INT_MAX chunks (~4.4 GB for bool, ~26 GB for
  multinomial). Worth a second opinion on whether the while conversion is
  correct in all eight fills.
• irk_discrete_uniform_long_vec already derived its mask from the inclusive
  maximum, which is the pattern this restores in the 64-bit fill.

[ndgrigorian]

@vlad-perevezentsev
Is this also fixed in #168 ? Or separate?

[vlad-perevezentsev]

@vlad-perevezentsev Is this also fixed in #168 ? Or separate?

Looks separate
irk_rand_uint32_vec/irk_rand_uint64_vec still call viRngUniformBits32/64 (and this PR adds one more viRngUniformBits64 call in mask == rng branch)
For WH/MCG31/R250/MRG32K3A those return -1014 so the full-range scalar draw is still broken for these BRNGs

[antonwolfy]

There is near-twin irk_discrete_uniform_long_vec (used by randint_untyped) keeps the old single-path rejection loop and the pre-fix VSL_RNG_METHOD_UNIFORM_STD constant.
It'd be nice to research on performance improvement there also in the follow-up PR.

[ndgrigorian]

@vlad-perevezentsev Is this also fixed in #168 ? Or separate?

Looks separate irk_rand_uint32_vec/irk_rand_uint64_vec still call viRngUniformBits32/64 (and this PR adds one more viRngUniformBits64 call in mask == rng branch) For WH/MCG31/R250/MRG32K3A those return -1014 so the full-range scalar draw is still broken for these BRNGs

Can you @vchamarthi work with @vlad-perevezentsev to fix the full-range scalar draw also for the other BRNGs?