feat: add GeneralizedExtremeValueDistribution (GEV) — unifies Gumbel/Fréchet/Weibull-max families (depends on #54)

[OldCrow]

Summary

Add the Generalized Extreme Value (GEV) distribution, which unifies three extreme value families (Gumbel, Fréchet, Weibull) under a single shape parameter ξ.

Background

GEV(μ, σ, ξ) is the limit distribution of properly normalised maxima of i.i.d. sequences:

• ξ = 0: Gumbel (Type I, double-exp) — see issue feat: add LogisticDistribution and GumbelDistribution — closed-form SIMD via vector_exp #54
• ξ > 0: Fréchet (Type II, heavy-tailed power law)
• ξ < 0: Weibull maximum (Type III, bounded upper tail)

Scipy: genextreme. Boost.Math does not have GEV as a named distribution.

Formulas

Let z = (x - μ) / σ.

For ξ ≠ 0:

t(x)    = (1 + ξ·z)^(-1/ξ)          [requires ξ·z > -1]
CDF(x)  = exp(-t(x))
PDF(x)  = (1/σ) · t(x)^(ξ+1) · exp(-t(x))
log PDF = -log(σ) + (ξ+1)·log(t(x)) - t(x)

For ξ = 0 (Gumbel limit):
Delegate to GumbelDistribution(μ, σ) from issue #54.

The log-space computation (ξ+1)·log(t) uses vector_log and vector_exp — the same SIMD primitives as Pareto and Weibull.

Implementation notes

• Prerequisite: issue feat: add LogisticDistribution and GumbelDistribution — closed-form SIMD via vector_exp #54 (GumbelDistribution) should be implemented first. GEV delegates to it at ξ = 0 (using a small epsilon threshold, e.g. |ξ| < 1e-10).
• MLE: Maximum likelihood for GEV is well-defined and standard; numerical optimisation via Newton–Raphson or profile likelihood.
• Mean = μ + σ(Γ(1-ξ)-1)/ξ for ξ < 1, ξ ≠ 0; μ + σγ for ξ = 0 (γ = Euler–Mascheroni constant).
• Domain: support is (μ - σ/ξ, +∞) for ξ > 0; (-∞, μ - σ/ξ) for ξ < 0; (-∞, +∞) for ξ = 0.
• is_delegation_wrapper = false (the Fréchet and Weibull-max cases are standalone; the Gumbel case delegates).
• Dispatch thresholds: NEVER in all tables until profiled. The log+exp pipeline suggests thresholds similar to LogNormal or Pareto.
• Kurtosis is undefined (infinite) for ξ ≥ 1/4; return NaN for those parameter values.