#5743 Bigfloat erf and erfc are slow and fail at high fpprec

[David Scherfgen]

By Claude:

Bigfloat erf and erfc fail at high fpprec for many arguments, for example erf(0.55b0 - 1.25b0*%i) at fpprec : 100 and erfc(1.23b0) at fpprec : 220:

(%i1) fpprec : 100$
(%i2) erf(0.55b0 - 1.25b0*%i);
gamma_incomplete: continued fractions failed for gamma_incomplete(5.0b-1, -1.375b0*%i-1.26b0).
 -- an error. To debug this try: debugmode(true);
(%i3) fpprec : 220$
(%i4) erfc(1.23b0);
gamma_incomplete: continued fractions failed for gamma_incomplete(5.0b-1, 1.5129b0).
 -- an error. To debug this try: debugmode(true);

At fpprec : 200 they fail for about one in seven random complex arguments in [-2, 2]^2. Well below that, they are slow. At fpprec : 200, erf(1.5b0 + 1.5b0*%i) takes over a third of a second, and at the default fpprec, erf(0.5b0 - 1.2b0*%i) takes 11 ms:

(%i5) fpprec : 200$
(%i6) showtime : true$
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%i7) erf(1.5b0 + 1.5b0*%i)$
Evaluation took 0.3679 seconds (0.3760 elapsed) using 373.283 MB.
(%i8) showtime : false$
(%i9) fpprec : 16$
(%i10) block([t0 : elapsed_real_time()], for k thru 100 do erf(0.5b0 - 1.2b0*%i),
         elapsed_real_time() - t0);
(%o10)                        1.1440000000000001

Bigfloat erf and erfc go through gamma_incomplete(1/2, z^2), and its continued fraction needs a number of iterations that grows like fpprec^2, and like 1/realpart(z)^2, and gives up after 10000 of them. erfi, inverse_erf and inverse_erfc call erf or erfc and are slow for the same reason.

[Raymond Toy]

Fixed in commit [71a379] merging branch bug-5743-bfloat-erf-erfc-slow

Basically stopped using gamma_incomplete for bigfloats because that's slow. Instead use the power series or asymptotic series to compute erf/erfc. Use increased precision for these calculations.