Opened 9 months ago

Last modified 6 months ago

#14933 new bug

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 (4)

comment:1 by waddlesplash, 9 months ago

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

comment:2 by pulkomandy, 9 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, 8 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, 6 months ago

Component: - GeneralSystem/
Note: See TracTickets for help on using tickets.