impl/python.py

# Derived from fdlibm e_log.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 split of ln 2: the high part has trailing zero bits, so e * LN2_HI
# is exact for any exponent a double can have.
LN2_HI = 6.93147180369123816490e-01
LN2_LO = 1.90821492927058770002e-10
SQRT2 = 1.4142135623730951


def ln(x: float) -> float:
    """Natural logarithm using only +, -, * and /.

    math.log comes from the platform's C library, whose last bit may differ
    from JavaScript's and Rust's. This follows the same steps in the same
    order as the TypeScript and Rust versions, so all three return the same
    double. See the README for the method.
    """
    if isinstance(x, bool) or not isinstance(x, (int, float)) or not math.isfinite(x) or x <= 0:
        raise ValueError("ln is only defined for finite numbers greater than zero, received %r" % (x,))
    m = float(x)
    e = 0
    while m >= 2:
        m /= 2
        e += 1
    while m < 1:
        m *= 2
        e -= 1
    if m > SQRT2:
        m /= 2
        e += 1
    f = m - 1
    s = f / (2 + f)
    z = s * s
    p = 1 / 25
    p = 1 / 23 + z * p
    p = 1 / 21 + z * p
    p = 1 / 19 + z * p
    p = 1 / 17 + z * p
    p = 1 / 15 + z * p
    p = 1 / 13 + z * p
    p = 1 / 11 + z * p
    p = 1 / 9 + z * p
    p = 1 / 7 + z * p
    p = 1 / 5 + z * p
    p = 1 / 3 + z * p
    lnm = 2 * s + 2 * s * (z * p)
    return e * LN2_HI + (e * LN2_LO + lnm)