Better inverse Erfc function for the Normal distribution

The Normal distribution uses InverseErfc for the inverse CDF. This is a wrapper for InverseErf and as such it is limited in precision as p -> 0 since the method computes using (2p - 1). When p is < 1e-16 then the inverse CDF returns -infinity. The smallest z is around -8.2
as computed by matlab:

norminv(1e-16)
-8.222082216130437

scipy, R, matlab, octave, Mathematica all have methods to invert small p. Here is output from matlab for norminv(5e-324): -38.472346342768788. It is accurate down to the minimum value for a 64-bit float.

There are two free libraries that seem to be used for this: Boost and Cephes (function ntdri). The Boost licence is compatible with Apache and there are some Boost derived works already in Commons Numbers. Incorporating the Boost Erf functions into numbers would be useful.

The inability to invert the Erfc for all p affects the Normal, Truncated normal and Levy distributions.

Boost error function inverses:
https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_erf/error_inv.html