Post Reply 
(Free42) roundoff for complex SQRT
02-12-2021, 02:45 PM
Post: #21
RE: (Free42) roundoff for complex SQRT
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
Find all posts by this user
Quote this message in a reply
02-12-2021, 07:44 PM
Post: #22
RE: (Free42) roundoff for complex SQRT
(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.
Visit this user's website Find all posts by this user
Quote this message in a reply
02-12-2021, 08:12 PM
Post: #23
RE: (Free42) roundoff for complex SQRT
according to the thread, you implemented a specific Re()=0 test?
Find all posts by this user
Quote this message in a reply
02-12-2021, 09:35 PM
Post: #24
RE: (Free42) roundoff for complex SQRT
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;
}
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)