+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not quite HP Calculators - but related (/forum-8.html)
+--- Thread: (Free42) roundoff for complex SQRT (/thread-10440.html)
Pages: 1 2
RE: (Free42) roundoff for complex SQRT - Werner - 02-12-2021 02:45 PM
Resurrecting this old thread ...
When X = 0+iB, SQRT(X) should have equal real and imaginary parts (apart from the sign).
But that is (still) not the case.
The 48G algorithm I posted gives the correct result - in the 48, with its 3 extra digits, but not in Free42.
Example:
-1E-12
SQRT
SQRT
COMPLEX
-
equals -1E-37, exactly what you get when you follow the 48G algorithm.
Cheers, Werner
RE: (Free42) roundoff for complex SQRT - Thomas Okken - 02-12-2021 07:44 PM
(02-12-2021 02:45 PM)Werner Wrote: When X = 0+iB, SQRT(X) should have equal real and imaginary parts (apart from the sign).
But that is (still) not the case.
I don't think that is something I ever tested. I implemented your algorithm in 2.0.21 and verified that it gave the desired result for 0+iB when B/2 is a perfect square... I suppose this means an explicit check for Re(X)=0 will be needed after all.
RE: (Free42) roundoff for complex SQRT - Werner - 02-12-2021 08:12 PM
according to the thread, you implemented a specific Re()=0 test?
RE: (Free42) roundoff for complex SQRT - Thomas Okken - 02-12-2021 09:35 PM
That's what I was thinking of doing at first, but then you suggested the 48G algorithm, and I implemented that instead. There is no check for pure imaginary.
free42/common/core_commands6.cc
Code:
static int mappable_sqrt_c(phloat xre, phloat xim, phloat *yre, phloat *yim) {
if (xre == 0 && xim == 0) {
*yre = 0;
*yim = 0;
return ERR_NONE;
}
phloat r = hypot(xre, xim);
phloat a = sqrt((r + fabs(xre)) / 2);
phloat b = xim / (a * 2);
if (p_isinf(a)) {
xre /= 100;
xim /= 100;
r = hypot(xre, xim);
a = sqrt((r + fabs(xre)) / 2);
b = xim / (a * 2);
a *= 10;
b *= 10;
}
if (xre >= 0) {
*yre = a;
*yim = b;
} else if (b >= 0) {
*yre = b;
*yim = a;
} else {
*yre = -b;
*yim = -a;
}
return ERR_NONE;
}