Dear Axel and Acer,

  Thanks to your help, I think I've got it working now.  For the sake of future readers, here's how:

(1) I used fsolve instead of solve,

(2) I used "fulldigits", which forces fsolve to use full precision in its internal computations, and

(3) I gave an initial guess.

(I was lucky enough in my problem that getting a good initial guess didn't require any extra effort.)

To demonstrate the difference (and the syntax), here's my original function and the better function.

> Digits:=75:
> inv_erf_bad := x->fsolve(erf(y/sqrt(2.0))=x,y):
> inv_erf_better := (x,initial_guess)->fsolve(erf(y/sqrt(2.0))=x,y=initial_guess, fulldigits);

> inv_erf_bad(.999999999999999999);
fsolve(erf(0.707106781186547524400844362104849039284835937688474036588339868995366239232 y) =
    0.999999999999999999, y)
> inv_erf_better(.999999999999999999,9.0);
          8.83510978817539578676690418014833167528733167476311521056292913279530244841
 

Dear Axel and Acer,

  Thanks to your help, I think I've got it working now.  For the sake of future readers, here's how:

(1) I used fsolve instead of solve,

(2) I used "fulldigits", which forces fsolve to use full precision in its internal computations, and

(3) I gave an initial guess.

(I was lucky enough in my problem that getting a good initial guess didn't require any extra effort.)

To demonstrate the difference (and the syntax), here's my original function and the better function.

> Digits:=75:
> inv_erf_bad := x->fsolve(erf(y/sqrt(2.0))=x,y):
> inv_erf_better := (x,initial_guess)->fsolve(erf(y/sqrt(2.0))=x,y=initial_guess, fulldigits);

> inv_erf_bad(.999999999999999999);
fsolve(erf(0.707106781186547524400844362104849039284835937688474036588339868995366239232 y) =
    0.999999999999999999, y)
> inv_erf_better(.999999999999999999,9.0);
          8.83510978817539578676690418014833167528733167476311521056292913279530244841
 

Dear Acer,

  Thanks for introducing me to the Statistics package, I'll definitely have to play around with that!  However, I'm not sure that it addresses the problem with the tails.  I think the analog to my original weird example would be:

> restart:with(Statistics):

> X:=RandomVariable('Normal'(0,1)):

> evalf[50](Quantile(X,1-10^(-18)));
-3.9377187255341233799589551305177585550470808666755 -
    4.5754487927888755526460423197296824251197610676463 I

which is still strangely complex valued.  I should have mentioned that I had previously tried setting the "Digits:=100:", and it didn't resolve the issue, so I don't think it's a problem with precision per se.  Incidentally, if I examine the corresponding left tail, I get:

> evalf[50](Quantile(X,10^(-18)));  
1.4142135623730950488016887242096980785696718753769
    RootOf(500000000000000000 erf(_Z) + 499999999999999999)

I don't even know how to parse that answer.  (The correct answer should be somewhere around -8 or -9.)

It's weird, but if I evaluate "erf" or "erfc" directly (instead of trying to invert them), they work fine for these values.  For example:

> evalf[50](erfc(9/sqrt(2)));
                                                                 -18
          0.22571768119076812954710041519374945159600838016441 10

What should I do?  It seems silly to invert the CDF by hand.

Dear Acer,

  Thanks for introducing me to the Statistics package, I'll definitely have to play around with that!  However, I'm not sure that it addresses the problem with the tails.  I think the analog to my original weird example would be:

> restart:with(Statistics):

> X:=RandomVariable('Normal'(0,1)):

> evalf[50](Quantile(X,1-10^(-18)));
-3.9377187255341233799589551305177585550470808666755 -
    4.5754487927888755526460423197296824251197610676463 I

which is still strangely complex valued.  I should have mentioned that I had previously tried setting the "Digits:=100:", and it didn't resolve the issue, so I don't think it's a problem with precision per se.  Incidentally, if I examine the corresponding left tail, I get:

> evalf[50](Quantile(X,10^(-18)));  
1.4142135623730950488016887242096980785696718753769
    RootOf(500000000000000000 erf(_Z) + 499999999999999999)

I don't even know how to parse that answer.  (The correct answer should be somewhere around -8 or -9.)

It's weird, but if I evaluate "erf" or "erfc" directly (instead of trying to invert them), they work fine for these values.  For example:

> evalf[50](erfc(9/sqrt(2)));
                                                                 -18
          0.22571768119076812954710041519374945159600838016441 10

What should I do?  It seems silly to invert the CDF by hand.