impl/python.py

# Derived from fdlibm e_rem_pio2.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

from .math_sin_cos_types import SinCos

# pi/2 in two parts (fdlibm): PIO2_1 has 33 significant bits, so k * PIO2_1 is
# exact for every k up to 2^20, which covers the permitted range.
PIO2_1 = 1.57079632673412561417e+00
PIO2_1T = 6.07710050650619224932e-11
INV_PIO2 = 6.36619772367581382433e-01
LIMIT = 1000000


def _sin_small(r: float) -> float:
    s = r * r
    p = 1 / 355687428096000
    p = -1 / 1307674368000 + s * p
    p = 1 / 6227020800 + s * p
    p = -1 / 39916800 + s * p
    p = 1 / 362880 + s * p
    p = -1 / 5040 + s * p
    p = 1 / 120 + s * p
    p = -1 / 6 + s * p
    return r + r * s * p


def _cos_small(r: float) -> float:
    s = r * r
    p = -1 / 6402373705728000
    p = 1 / 20922789888000 + s * p
    p = -1 / 87178291200 + s * p
    p = 1 / 479001600 + s * p
    p = -1 / 3628800 + s * p
    p = 1 / 40320 + s * p
    p = -1 / 720 + s * p
    p = 1 / 24 + s * p
    p = -1 / 2 + s * p
    return 1 + s * p


def sin_cos(radians: float) -> SinCos:
    """Sine and cosine using only +, -, * and /, in the same order as the
    TypeScript and Rust versions, so all three return the same doubles.
    math.sin and math.cos come from the C library and may differ in the last
    bit. See the README for the method."""
    if (
        isinstance(radians, bool)
        or not isinstance(radians, (int, float))
        or not math.isfinite(radians)
        or radians < -LIMIT
        or radians > LIMIT
    ):
        raise ValueError("radians must be a finite number between -1000000 and 1000000, received %r" % (radians,))
    x = float(radians)
    k = math.floor(x * INV_PIO2 + 0.5)
    r = x - k * PIO2_1 - k * PIO2_1T
    s = _sin_small(r)
    c = _cos_small(r)
    quadrant = k % 4
    if quadrant == 0:
        return SinCos(sin=s + 0.0, cos=c + 0.0)
    if quadrant == 1:
        return SinCos(sin=c + 0.0, cos=-s + 0.0)
    if quadrant == 2:
        return SinCos(sin=-s + 0.0, cos=-c + 0.0)
    return SinCos(sin=-c + 0.0, cos=s + 0.0)