(Free42) roundoff for complex SQRT - Printable Version +- 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; }```