perf: batch array-valued distribution parameters into a single fill

[vchamarthi]

What

When a distribution parameter is an array rather than a scalar,
vec_cont1_array / vec_cont2_array invoke the fill once per element:

for i from 0 <= i < cnp_PyArray_MultiIter_SIZE(multi):
    func(state, 1, &array_data[i], oa_data[0], ob_data[0])

That is 100,000 VSL calls for a 100,000-element array.

For location and scale families the parameters can be applied after the draw
instead, so one call covers the whole request:

normal(loc, scale)     == loc + scale * N(0, 1)
exponential(scale)     == scale * Exp(1)
lognormal(mean, sigma) == exp(mean + sigma * N(0, 1))

These are exact identities, not approximations.

Routed: normal, uniform, exponential, laplace, gumbel, logistic,
rayleigh, lognormal.

Not routed: gamma, beta, f, chisquare, vonmises, wald,
triangular, pareto, weibull, power, and the discrete distributions.
Their parameters change the shape, so no affine step recovers them. gamma is
kept in the benchmark below as a control.

Numbers

Xeon Gold 6338, single thread, N=100k, 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.

distribution	before	after
gumbel	242.6	3.44	70.6x
rayleigh	182.2	2.63	69.4x
exponential	134.3	2.16	62.2x
normal	141.4	3.14	45.1x
lognormal	176.2	4.11	42.8x
laplace	145.4	3.86	37.6x
logistic	71.9	9.85	7.3x
uniform	40.9	8.37	4.9x
gamma (control)	186.4	185.4	1.0x

Scalar-parameter paths are untouched and not measured here.

Streams change

Array-parameter streams change: the draw is now standardised and transformed
rather than parameterised up front. Scalar-parameter paths are byte-identical.
Noted in the changelog.

Testing

• Suite: 190 -> 211 passed, 3 skipped. The 21 new tests also pass on unpatched
  master, so they are invariant guards rather than fitted to this change.
• Constant-valued arrays match the scalar path (assert_allclose).
• KS test against each analytic CDF with per-element varying parameters, D vs
  crit 0.00364: all pass.
• Broadcast shapes, size handling and error cases compared against stock
  numpy: agree, including incompatible size raising ValueError.
• Verified on upstream master 6b7962f.

Notes for reviewers

• _param_out_shape reimplements the shape and compatibility logic rather than
  reusing the multi-iterator, since the batched path has no iterator. It raises
  ValueError("size is not compatible with inputs") where the old path
  surfaced numpy's broadcast error text for some inputs. Same exception type,
  different message.
• lognormal is built on the normal fill, not the lognormal one: its parameters
  sit inside the exponential, so an affine step cannot apply them to a
  standardised lognormal.
• Existing test_uniform_array_bounds_return_ndarray covers this path and still
  passes.

[antonwolfy]

Tests relating note is not something the user might be interesting in. I'd propose to avoid adding such things in the changelog.

[antonwolfy]

We are testing only with default method=ICDF, and so no BoxMuller paths coverage.

[antonwolfy]

The per-element test is statistically weak.
If per-element scale were ignored, that comparison is a ~50/50 coin flip, so the test can pass even when the bug it targets is present. It also never sweeps a per-element loc. Strengthen to a magnitude check.

[antonwolfy]

Suggested change

	* The random streams for the array-valued-parameter paths of the distributions above have changed: with a fixed seed these now produce different (but equally valid) samples. Scalar-parameter paths are unaffected. [gh-171](https://github.com/IntelPython/mkl_random/pull/171)
	* Changed the random streams for the array-valued-parameter paths of the distributions above: with a fixed seed these now produce different (but equally valid) samples. Scalar-parameter paths are unaffected. `lognormal`'s Box-Muller path now draws through MKL's Gaussian `BOXMULLER2` generator. [gh-171](https://github.com/IntelPython/mkl_random/pull/171)

[antonwolfy]

The behavior changed now on array path:

• Old path (vec_cont2_array → irk_uniform_vec(low, high)): MKL generated the values directly on the interval via VSL_RNG_METHOD_UNIFORM_STD_ACCURATE with the actual low/high. MKL's uniform is generated on [low, high), so results stayed strictly below high.
• New path (vec_uniform_array): fills a standard U[0, 1) and reconstructs with low + u·(high − low) via np.subtract + np.multiply + np.add. Each of those three steps rounds to nearest double.

The new path is aligned how NumPy implements that, but need to align documentation also:

Suggested change

	Upper boundary of the output interval. All values generated will be
	less than high. The high limit may be included in the returned array of
	floats due to floating-point rounding in the equation
	``low + (high-low) * random_sample()``. high - low must be
	non-negative. The default value is 1.0.

[antonwolfy]

Also needs to add a note to the CHANGELOG that array-valued bounds may include high due to floating-point rounding

[antonwolfy]

It is stricter than the legacy check. This matches NumPy's own semantics and is one-directional (no previously-invalid call now passes). It's arguably a fix, but it's an acceptance-behavior change worth a changelog line.