Question: Weirdness with the inverse Erf function

Dear Maple Experts,

  I'm fairly new to Maple.  I've been trying to compute Gaussian tail probabilities accurately.  As I understand it, Maple has an "erf" function, but not an inverse erf.  Other posts have suggested doing something like:

  inv_erf := x->solve(erf(y/sqrt(2.0))=x,y);

At first ; for instance, about 68% of a Gaussian's probability is within 1 standard deviation, so:

> erf(1/sqrt(2))

    0.6826894920

> evalf(inv_erf(0.6826894920));

    0.9999999994

So far, so good.  Unfortunately, my interest lies in fairly large deviations (up to about 9 standard deviations).  All of a sudden, I start to get complex (!) answers:

> evalf(inv_erf(1-(10^(-16))));
                                 8.304785423
> evalf(inv_erf(1-(10^(-17))));
                         -3.937718724 - 4.575448791 I
 

I've tried using "erfc" instead of "erf", under the logic that it might be better at estimating the tails, but it doesn't seem to resolve the problem.  What's going on and how do I fix it?

  Thanks for any help.