Special functions and `libm` compatibility

[YeungOnion]

statrs should route all float ops through num_traits which means we'll rely on libm or std based.
For a brief period, I thought, we should consider using libm since that's a growing default for the no_std environment, but then how complete does a math function backend need to be for consistency?
libm implements special functions that aren't available in std and there are other crates that maintain special functions, this exposes cases where statrs and other functions aren't the same.

The question I have on this is,

What footguns can we mitigate for users who need some statistics functions that wouldn't already be mitigating if they had a crash course in floating point?

I think the community convergence onto a special functions crate could reduce some of this, since the number of places needed to make consistent choices could be reduced. Right now num_traits has this role for the overlap between std and libm. Perhaps a crate would be able to implement sane trait-level defaults while exposing specific algorithms as functions where the backend selection would only control the trait-level would be useful. Then we can contribute upstream to broader numerics needs as well.

special exposes its functions
as traits on float and double atop libm
puruspe is also in this space, author is that of peroxide (which to me, it feels like Matlab in Rust).

I don't know what direction to take yet, but some input on if this would be an issue someone actually stumbles on would be nice to know. If this is fruitful, would also try to get input from stainless-steel and Axect.

---

Here's a summary on what we actually use and whether its available elsewhere. Note that no_std environment will typically imply libm version of it could be consumed by a user of statrs.

Function statrs needs	Where used	Hand-rolled?	In libm?	In puruspe?	In special?
erf / erfc	function/erf.rs	Boost-derived rational/continued-fraction (erf_impl)	libm::erf, libm::erfc	erf, erfc	Error::error, Error::compl_error
erf_inv / erfc_inv	function/erf.rs	Boost-derived (erf_inv_impl)	-	inverf, inverfc	Error::inv_error only (no erfc_inv)
ln_gamma	function/gamma.rs	Lanczos approx (Pugh, 2004)	libm::lgamma / lgamma_r	ln_gamma	Gamma::ln_gamma
gamma (Γ)	function/gamma.rs	Lanczos approx (Pugh, 2004)	libm::tgamma	gamma	Gamma::gamma
digamma / inv_digamma	function/gamma.rs	Algorithm AS 103	-	-	Gamma::digamma only (no inv_digamma; also adds trigamma)
incomplete gamma (gamma_ui/li/ur/lr)	function/gamma.rs	continued-fraction (Cephes-style, uncited)	-	gammp, gammq, invgammp	Gamma::inc_gamma
ln_beta, beta	function/beta.rs	built from ln_gamma	-	beta	Beta::ln_beta only (no plain beta)
incomplete/regularized beta (beta_inc, beta_reg, inv_beta_reg)	function/beta.rs	Algorithm AS 63/64/109	-	betai, invbetai	Beta::inc_beta, Beta::inv_inc_beta
logistic, logit	function/logistic.rs	trivial, built from exp/ln	-	-	-
factorial, ln_factorial, binomial, multinomial	function/factorial.rs	lookup table + ln_gamma	-	-	-
harmonic, gen_harmonic	function/harmonic.rs	definition	-	-	-