impl/python.py

# Derived from fdlibm s_atan.c.
# Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
# Developed at SunSoft, a Sun Microsystems, Inc. business.
# Permission to use, copy, modify, and distribute this
# software is freely granted, provided that this notice
# is preserved.

import math

# fdlibm s_atan.c (Sun Microsystems, 1993). atan of 0.5, 1, 1.5 and infinity,
# each split into a high and a low part so the sum carries more than 53 bits.
ATAN_HI = (
    4.63647609000806093515e-01,
    7.85398163397448278999e-01,
    9.82793723247329054082e-01,
    1.57079632679489655800e00,
)
ATAN_LO = (
    2.26987774529616870924e-17,
    3.06161699786838301793e-17,
    1.39033110312309984516e-17,
    6.12323399573676603587e-17,
)
AT = (
    3.33333333333329318027e-01,
    -1.99999999998764832476e-01,
    1.42857142725034663711e-01,
    -1.11111104054623557880e-01,
    9.09088713343650656196e-02,
    -7.69187620504482999495e-02,
    6.66107313738753120669e-02,
    -5.83357013379057348645e-02,
    4.97687799461593236017e-02,
    -3.65315727442169155270e-02,
    1.62858201153657823623e-02,
)
TWO_POW_66 = 73786976294838206464.0
TWO_POW_M29 = 1.862645149230957e-09


def atan(x: float) -> float:
    """Arctangent in radians using only +, -, * and /, in the same order as the
    TypeScript and Rust versions, so all three return the same double.
    math.atan comes from the C library and may differ in the last bit.

    fdlibm's method: reduce |x| against atan(0.5), atan(1), atan(1.5) or
    pi/2, then an odd polynomial on the remainder. fdlibm picks the interval
    from the high word of the double; comparing |x| with the same boundaries
    (7/16, 11/16, 19/16, 39/16, 2^66, 2^-29) is exactly equivalent."""
    if isinstance(x, bool) or not isinstance(x, (int, float)) or not math.isfinite(x):
        raise ValueError("x must be a finite number, received %r" % (x,))
    x = float(x)
    negative = x < 0
    a = -x if negative else x
    if a >= TWO_POW_66:
        z = ATAN_HI[3] + ATAN_LO[3]
        return -z if negative else z
    if a < 0.4375:
        if a < TWO_POW_M29:
            # atan(x) rounds to x here; + 0.0 turns -0 into 0.
            return x + 0.0
        i = -1
    else:
        x = a
        if a < 1.1875:
            if a < 0.6875:
                i = 0
                x = (2.0 * x - 1.0) / (2.0 + x)
            else:
                i = 1
                x = (x - 1.0) / (x + 1.0)
        elif a < 2.4375:
            i = 2
            x = (x - 1.5) / (1.0 + 1.5 * x)
        else:
            i = 3
            x = -1.0 / x
    z = x * x
    w = z * z
    s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10])))))
    s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9]))))
    if i < 0:
        return x - x * (s1 + s2)
    z = ATAN_HI[i] - ((x * (s1 + s2) - ATAN_LO[i]) - x)
    return -z if negative else z