integral competition HP50g vs. DM42

[John Keith]

(08-23-2020 07:38 PM)Albert Chan Wrote:  >>> k = 2/pi**0.5
>>> erf = lambda x: k*x * hyp1f1(1/2, 3/2, -x*x)
>>> CiS = lambda x, erf=erf: (1+1j)/2 * erf((1-1j)/k * x)
>>> CiS(3.9)
(0.42233270272270274 + 0.47520241152278453j)

CiS(3.9) called erf(z = (1-1j) * 3.4562850092657564), |z| ≈ 5 > 3
As expected, we lose many digits precision due to massive cancellation.

...

>>> CiS(3.9, erf = lambda x: 1 - erfc(x))
(0.42233271026093344 + 0.47520240235068834j)

Bringing up an old thread with a question and a comment. The question: How do you get an accurate erfc(z) from 1-erf(z)? erf(3.9) = 0.999999965208, and that number subtracted from 1 is 0.000000034792 which has only 5 significant digits.

Now the comment: Much of the loss of precision in the first example comes from the large negative x in hyp1f1(1/2, 3/2, -x*x) which involves summing large alternating positive and negative numbers resulting in cancellation losses. The result, 0.22723763473 has only 7 correct digits. J-M Baillard's Special Function package uses a different formula, erf x = (2x/Pi^(1/2)) * exp(-x*x) * 1F1( 1 , 3/2 , x^2 ), which returns 0.999999965204, which is within 4 ULPs of the exact result. The positive x seems to result in more stable behavior of the hypergeometric 1F1 program.

Additional note: Baillard's program for C(x) and S(x) also uses a continued fraction if x > sqrt(pi).