Avoid hangs and debug asserts on invalid parameters for Zipf

[mstoeckl]

• Added a CHANGELOG.md entry

Summary

The Zipf distribution sampler could panic (at debug_assert!(t > F::zero());) when given invalid parameters (like s = inf); or hang in an infinite loop when sampling. This change rejects parameter combinations for which the generalized Zipf distribution is not well defined, and slightly modifies Zipf::new() so that sampling avoids the debug assert and produces the correct value when s = inf.

Motivation

This another simple panic / hang issue that I found by fuzzing input parameters.

Details

The normalization constant (sum from i=1 to n of 1/i^s) of the generalized Zipf distribution needs to be positive and finite for the distribution to be well defined. if s <= 1 and n = inf, then the sum in the constant diverges.

If s = inf, then each term of the constant is 0, so the specific formula in Wikipedia will not work; however, the Zipf distribution approaches the distribution concentrated on i=1 as s -> inf. The previous implementation computed a nan in Zipf::new which triggered the debug assertion (or an infinite loop), but a small change avoids both.

(Edit: originally this PR rejected the case s = inf, but after posting it I changed my mind; in general, if a limiting distribution exists and is close enough to the distribution at the largest finite parameters, then producing that limiting distribution is safe and avoids having the library user write a special case to produce effectively that limiting distribution.)

[dhardy]

I think inv_cdf also needs a special case for s = inf? Perhaps not if it's safe to assume that p <= 1. Maybe that should be a debug_assert.

[mstoeckl]

I think inv_cdf also needs a special case for s = inf? Perhaps not if it's safe to assume that p <= 1. Maybe that should be a debug_assert.

It should be OK. The range of possible inputs for p in inv_cdf should be [0,1]; anything outside this is not a probability value. Zipf::sample is the only caller and provides it a sample from StandardUniform (in [0,1)).

[dhardy]

The case where n.is_finite() && s = inf does, as you say, converge to a result of 1. Is this reasonable, given that infinities often occur due to loss of precision? That depends, but this is certainly not the first case where we have decided to adopt this approach, so it seems acceptable. (If for certain applications this is deemed unacceptable we could look at adding a feature flag for stricter bounds checking.)

The case where n = s = inf is more dubious; I think it would be more reasonable to consider this IllDefined.

Otherwise this PR is acceptable in my opinion.

[dhardy]

Suggested change

	if n.is_infinite() && s <= F::one() {
	if n.is_infinite() && (s.is_infinite() || s <= F::one()) {