HP Prime bug for approximate(erfc(z)), z as (0+b*i), or erf(b*i)

[Gil]

Try for instance on the HP PRIME approximate(erfc(5*i))

ie with no real part in the complex number input (argument).

It returns -82982738 79.76*i (a complex number with no real part).
According to Wolfram Alpha, the result should be
1 - 8.2982738 806 7680 ×10^9 i
(with the number 1 as the real part of the resulting complex number).

Fazit
1) As no real part of the HP PRIME result is shown, that part has to be considered being equal to 0, and the PRIME result seems therefore clearly to be wrong here (a bug?).
2) As a basic program on the HP50G allows to find (1.,-82982738 80.61), we see that the imaginary part of the HP PRIME is not so well "calibrated", letting suppose that the internal algorithm of the HP PRIME for erfc function could or should be improved.

Logically, on HP PRIME, approximate(erf(5*i)) returns
1 -82982738 79.76*i (a complex number with the number 1 as its real part, instead of 0).

Strange.

Or did I miss something?

Regards

Gil

[Albert Chan]

(11-15-2023 04:00 PM)Gil Wrote:  Logically, on HP PRIME, approximate(erf(5*i)) returns
1 -82982738 79.76*i (a complex number with the number 1 as its real part, instead of 0).

Strange.

erf(x) = 2/sqrt(pi) * x * hyp1f1(1/2, 3/2, -x*x)

erf(x)/x = 2/sqrt(pi) * hyp1f1(1/2, 3/2, -x*x)

If x is purely real or imaginery, RHS is real --> LHS is real too.

Without computation, erf(5i) is purely imaginery.

[Albert Chan]

(11-15-2023 04:53 PM)Albert Chan Wrote:  Without computation, erf(5i) is purely imaginery.

I would not consider this a bug though.
Difference is very tiny, especially compared in polar form.

p2> from mpmath import *
p2> x = erf(5j)
p2> x
mpc(real='0.0', imag='8298273880.6768036')
p2> polar(x)
(mpf('8298273880.6768036'), mpf('1.5707963267948966'))
p2> polar(x + 1)
(mpf('8298273880.6768036'), mpf('1.5707963266743896'))

[Brian Zilli]

The cutoff appears to be 3; that is, for any a>3, erf(a*i) will return a real part of 1. Near 3, the relative error is orders of magnitude larger.

It also leads to wonky results when computing purely real integrals in the CAS (and numerical Home in the 2020 firmware of the Mac emulator but not on a physical G2 with the latest firmware). e.g., int(e^(x^2),x,0,3.001) returns 1452.55199576+0.886226925453*I, when it should be approximately 1452.67.

[Brian Zilli]

Also interesting: try erf(3.16845945394*i) or int(e^(x^2),x,0,3.16845945394) in the CAS.

[parisse]

I confirm a small bug for erfc if |z| is between 3 and 6.5 and re(z) is small (missed a 1 in real part). In this area, erf is computed with a continued fraction expansion
2*exp(z^2)*int(exp(-t^2),t=z..inf)=1/(z+1/2/(z+1/(z+3/2/(z+...))))
There is a warning of low accuracy if re(z) is small : I must indeed stop the fraction expansion and since we do not have multiprecision floats there are some wrong digits.

[Albert Chan]

Cas> a := 3.16845945394*i
Cas> int(exp(-x*x), x=0..a)       → 0.886226925453-51940029439.9*i
Cas> Ans * 2/sqrt(pi)                → 1.-58608047158.3*i
Cas> erf(a)                               → 1.-58608047158.3*i 

Try numerical integration instead (double exponential method, quad6)

Cas> quad(x -> exp(-x*x), 0, a)       → [3833.13077361*i, 1.16978024776e−12]
Cas> Ans[1] * 2/sqrt(pi)                   → 4325.2249097*i == erf(3.16845945394*i)

[parisse]

The continued fraction expansion "low accuracy" is probably full wrong for some values of z near the imaginary axis...

[parisse]

I have disabled the continued fraction, it should now work correctly (it will just be a little bit slower if |z| is between 4 and 6.5).

[Gil]

I have the Android version for mobile phone.

Should be I desinstall it and install it again?

Or wait for an automatic update?

Right now erf (5i)returns 1+b×i.

[Gil]

An interesting approach for erf(z) is the one I found in
https://math.stackexchange.com/questions...inary-part

for which the resulting calculations (for a number of about 30 summations) seem to give better results (in the last digit) in my HP50G than the built-in function in HP PRIME.

[StephenG1CMZ]

(11-16-2023 01:33 PM)Gil Wrote:  I have the Android version for mobile phone.

Should be I desinstall it and install it again?

Or wait for an automatic update?

Right now erf (5i)returns 1+b×i.

I would not recommend uninstalling and reinstalling because Google Play does not reliably list the HP Prime Pro app.
So at best, you might get back what you had. At worst, there is a risk that you might be unable to reinstall.
https://www.hpmuseum.org/forum/thread-20...t+and+miss
Eventually, a new version of the App might be published, and that might include the updates described.
KlassKuperus said in September that an update was being worked on.

[Albert Chan]

(11-15-2023 08:23 PM)parisse Wrote:  There is a warning of low accuracy if re(z) is small : I must indeed stop the fraction expansion and
since we do not have multiprecision floats there are some wrong digits.

Throwing more precisions may not help.
Problem is more to do with how deep "fraction expansion" default we pick.

Let n = "last" numerator --> erfc(z, n=1) = exp(-z*z)/√(pi) / (z + (1/2)/(z + 1/z))

Let z = 3.16845945394*I, precisions = 1000 digits, results rounded to integer.

erfc(z, n=0.5) = -4292 I
erfc(z, n=1.0) = -4317 I
erfc(z, n=1.5) = -4323 I
erfc(z, n=2.0) = -4324 I
erfc(z, n=2.5) = -4325 I
erfc(z, n=3.0) = -4325 I
erfc(z, n=3.5) = -4326 I
erfc(z, n=4.0) = -4311 I
erfc(z, n=4.5) = -4324 I
erfc(z, n=5.0) = -4325 I
...
erfc(z, n=19.0) = -4324 I
erfc(z, n=19.5) = -4326 I
erfc(z, n=20.0) = 58957691655 I
erfc(z, n=20.5) = -4325 I
erfc(z, n=21.0) = -4326 I
...

For z=3.16845945394*I, n=20 is really bad.
But changing n will not help. It would just shift problematic z somewhere else.

---

More precisions can't help because CF formula denominator is close to 0
But, erfc(z) is continuous ... something is going to break.

1 / ±0 = ±∞

Cas> erfc(3.16845945394*i)      → +58608047158.3*i
Cas> erfc(3.16845945395*i)      → −14620291004.7*i

CF Denomaintor = x * p(x²)/q(x²), rational polynomial with all positive terms
--> erfc CF formula for non-zero real argument will not get divide-by-0 issue.

(04-06-2024 11:04 PM)Albert Chan Wrote:  denominator may get to 0 only if re(z)=0