impl/python.py

# Derived from fdlibm e_exp.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: k * LN2_HI is exact for every k this range needs.
LN2 = 0.6931471805599453
LN2_HI = 6.93147180369123816490e-01
LN2_LO = 1.90821492927058770002e-10


def exp(x: float) -> float:
    """e^x using only +, -, * and /, in the same order as the TypeScript and
    Rust versions, so all three return the same double. See the README."""
    if isinstance(x, bool) or not isinstance(x, (int, float)) or not math.isfinite(x) or x < -708 or x > 709:
        raise ValueError("exp needs x between -708 and 709, received %r" % (x,))
    x = float(x)
    k = math.floor(x / LN2 + 0.5)
    r = x - k * LN2_HI - k * LN2_LO
    p = 1 / 20922789888000
    p = 1 / 1307674368000 + r * p
    p = 1 / 87178291200 + r * p
    p = 1 / 6227020800 + r * p
    p = 1 / 479001600 + r * p
    p = 1 / 39916800 + r * p
    p = 1 / 3628800 + r * p
    p = 1 / 362880 + r * p
    p = 1 / 40320 + r * p
    p = 1 / 5040 + r * p
    p = 1 / 720 + r * p
    p = 1 / 120 + r * p
    p = 1 / 24 + r * p
    p = 1 / 6 + r * p
    p = 1 / 2 + r * p
    p = 1 + r * p
    result = 1 + r * p
    for _ in range(k):
        result *= 2
    for _ in range(-k):
        result /= 2
    return result