fix: stop prec::ulps_eq! degenerating into a 1e-9 absolute comparison

[agene0001]

approx's ulps_eq short-circuits on abs_diff_eq(epsilon) before it consults the ULPs distance. The crate's prec::ulps_eq! wrapper paired it with DEFAULT_EPS = 1e-9, which makes the ULPs bound unreachable — every internal ulps_eq!(x, y) was really a 1e-9 absolute comparison.

That matters because the macro is used at ~25 call sites to recognise exact parameter values (p == 1.0, x == x.floor(), shape == 1.0), so anything within 1e-9 of them silently took a degenerate branch:

Binomial(p = 1 − 2⁻³³, n = 100).pmf(99)    0     (true 1.16e-8)
Binomial(p = 1 − 2⁻³³, n = 100).ln_pmf(99) −inf
Geometric(p = 1 − 2⁻³³).max()              1     (true u64::MAX)
Geometric(p = 1 − 2⁻³³).skewness()         inf   (true 92681.9)
Beta(1 + 2⁻³³, 1 + 2⁻³³).pdf(0.0)          1     (true 0)
digamma(−1 + 2⁻³³)                         −inf  (true −8.59e9)

The fix adds DEFAULT_ULPS_EPS = f64::EPSILON and uses it for ulps_eq!/assert_ulps_eq!. The exact values still take the degenerate branches (Binomial(1.0, 5).pmf(5) == 1, Geometric(1.0).max() == 1, digamma at the poles is still −inf) — all covered by existing tests, which pass unchanged.

New tests use 1 − 2⁻³³ rather than 1 − 1e-10 so 1 − p is exact and the mpmath reference values aren't limited by the representation of p.

Also fixes Binomial::entropy, which summed −p·ln p unguarded and returned NaN once any mass underflowed (0 × −inf) — for p = 0.5 that's every n ≳ 1100. Categorical::entropy already filters zero terms; this brings Binomial in line.

Note: PR for #342 is stacked on this one — the ignored test_inverse_cdf_small_p turns out to need both fixes.

🤖 Generated with Claude Code