atan2 returns wrong result in some case

[hanya]

atan2 C function returns wrong result in the following case.

// atan2.c
#include <math.h>
#include <stdio.h>

int main(int argc, char* argv[]) {
    // from Python/Lib/test/cmath_testcases.txt acos0200
    // acos0200 acos 1.6206860518683021e+308 1.0308426226285283e+308 -> 
    // 0.56650826093826223 -710.54206874241561

    double real = 1.6206860518683021e+308;
    double imag = 1.0308426226285283e+308;
    double r = atan2(imag, real); // real part
    printf("expected: 0.56650826093826223\n");
    printf("  result: %1.17lf\n", r);
    //x86_64, result: 0.08332244049906361
}

// gcc -o atan2 atan2.c -g -O0

Python contains some test cases for its complex calculations. It contains near infinity acos0200 test executed by test_cmath. The test can be executed as follows:

python -m test test_cmath

The same calculation in the Python interactive:

>>> import cmath
>>> c = complex(1.6206860518683021e+308, 1.0308426226285283e+308)
>>> cmath.acos(c)
(0.08332244049906361-710.5420687424156j)

The imaginary part is correct.

The calculation of the cmath.acos is done by cmath_acos_impl function in Modules/cmathmodule.c. The real part is calculated as follows:

r.real = atan2(fabs(z.imag), z.real);

fabs call is removed from the C code above since the value is positive.

The result on the x86_gcc2(Python2) and secondary x86(Python3) are correct.

[waddlesplash]

Probably we should just update/swap the glibc math code (#10913), likely with musl's math library.

[pulkomandy]

I would prefer that we replace only parts we know are broken, for now. And usually I grab the replacements from FreeBSD, but musl is fine as well.

Also, it would be great to turn your example into an actual testcase to put in src/tests and using cppunit!

[hanya]

I built Python 3.7.2 with libm.a from musl and the test_cmath all passed on x86_64.

[pulkomandy]

So, I checked what we're doing. The atan2 code we're currently using is in the arch/generic/ directory and was imported for PowerPC. I suspect it uses PowerPC floating point format, which is not compatible with x86.

I checked what glibc does, and there, the same e_atan2.S file is used for both x86 and x86_64. So I think we should do just that. However I have no x86_64 system to test this on.

If someone can confirm this builds and fixes the problem: ​https://review.haiku-os.org/c/haiku/+/2065

[pulkomandy]

Fixed in hrev53683.

[waddlesplash]

That just made everything crash. Reverted in hrev53689.

[waddlesplash]

​https://review.haiku-os.org/c/haiku/+/2077 should fix this.

[waddlesplash]

Merged in hrev53728.

[nielx]

Assign tickets with status=closed and resolution=fixed within the R1/beta2 development window to the R1/beta2 Milestone