Opened 12 months ago

Closed 4 weeks ago

#14933 closed bug (fixed)

atan2 returns wrong result in some case

Reported by: hanya Owned by: nobody
Priority: normal Milestone: Unscheduled
Component: System/ Version: R1/Development
Keywords: Cc:
Blocked By: Blocking:
Has a Patch: no Platform: x86-64


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)

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.

Change History (9)

comment:1 by waddlesplash, 12 months ago

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

comment:2 by pulkomandy, 12 months ago

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!

comment:3 by hanya, 11 months ago

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

comment:4 by waddlesplash, 9 months ago

Component: - GeneralSystem/

comment:5 by pulkomandy, 7 weeks ago

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:

comment:6 by pulkomandy, 6 weeks ago

Resolution: fixed
Status: newclosed

Fixed in hrev53683.

comment:7 by waddlesplash, 6 weeks ago

Resolution: fixed
Status: closedreopened

That just made everything crash. Reverted in hrev53689.

comment:9 by waddlesplash, 4 weeks ago

Resolution: fixed
Status: reopenedclosed

Merged in hrev53728.

Note: See TracTickets for help on using tickets.