/pypy/rpython/lltypesystem/module/ll_math.py
Python | 430 lines | 321 code | 46 blank | 63 comment | 97 complexity | a2823afb42ed13dd46004d3925d574c3 MD5 | raw file
Possible License(s): AGPL-3.0, BSD-3-Clause, Apache-2.0
- import math
- import errno
- import py
- import sys
- from pypy.rpython.lltypesystem import lltype, rffi
- from pypy.tool.sourcetools import func_with_new_name
- from pypy.tool.autopath import pypydir
- from pypy.rlib import jit, rposix
- from pypy.translator.tool.cbuild import ExternalCompilationInfo
- from pypy.translator.platform import platform
- from pypy.rlib.rfloat import isfinite, isinf, isnan, INFINITY, NAN
- use_library_isinf_isnan = False
- if sys.platform == "win32":
- if platform.name == "msvc":
- # When compiled with /O2 or /Oi (enable intrinsic functions)
- # It's no more possible to take the address of some math functions.
- # Ensure that the compiler chooses real functions instead.
- eci = ExternalCompilationInfo(
- includes = ['math.h', 'float.h'],
- post_include_bits = ['#pragma function(floor)'],
- )
- use_library_isinf_isnan = True
- else:
- eci = ExternalCompilationInfo()
- # Some math functions are C99 and not defined by the Microsoft compiler
- cdir = py.path.local(pypydir).join('translator', 'c')
- math_eci = ExternalCompilationInfo(
- include_dirs = [cdir],
- includes = ['src/ll_math.h'],
- separate_module_files=[cdir.join('src', 'll_math.c')],
- export_symbols=['_pypy_math_acosh', '_pypy_math_asinh',
- '_pypy_math_atanh',
- '_pypy_math_expm1', '_pypy_math_log1p'],
- )
- math_prefix = '_pypy_math_'
- else:
- eci = ExternalCompilationInfo(
- libraries=['m'])
- math_eci = eci
- math_prefix = ''
- def llexternal(name, ARGS, RESULT, **kwargs):
- return rffi.llexternal(name, ARGS, RESULT, compilation_info=eci,
- sandboxsafe=True, **kwargs)
- def math_llexternal(name, ARGS, RESULT):
- return rffi.llexternal(math_prefix + name, ARGS, RESULT,
- compilation_info=math_eci,
- sandboxsafe=True)
- if sys.platform == 'win32':
- underscore = '_'
- else:
- underscore = ''
- math_fabs = llexternal('fabs', [rffi.DOUBLE], rffi.DOUBLE)
- math_log = llexternal('log', [rffi.DOUBLE], rffi.DOUBLE)
- math_log10 = llexternal('log10', [rffi.DOUBLE], rffi.DOUBLE)
- math_copysign = llexternal(underscore + 'copysign',
- [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE,
- elidable_function=True)
- math_atan2 = llexternal('atan2', [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE)
- math_frexp = llexternal('frexp', [rffi.DOUBLE, rffi.INTP], rffi.DOUBLE)
- math_modf = llexternal('modf', [rffi.DOUBLE, rffi.DOUBLEP], rffi.DOUBLE)
- math_ldexp = llexternal('ldexp', [rffi.DOUBLE, rffi.INT], rffi.DOUBLE)
- math_pow = llexternal('pow', [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE)
- math_fmod = llexternal('fmod', [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE)
- math_hypot = llexternal(underscore + 'hypot',
- [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE)
- math_floor = llexternal('floor', [rffi.DOUBLE], rffi.DOUBLE, elidable_function=True)
- math_sqrt = llexternal('sqrt', [rffi.DOUBLE], rffi.DOUBLE)
- math_sin = llexternal('sin', [rffi.DOUBLE], rffi.DOUBLE)
- math_cos = llexternal('cos', [rffi.DOUBLE], rffi.DOUBLE)
- @jit.elidable
- def sqrt_nonneg(x):
- return math_sqrt(x)
- sqrt_nonneg.oopspec = "math.sqrt_nonneg(x)"
- # ____________________________________________________________
- #
- # Error handling functions
- ERANGE = errno.ERANGE
- EDOM = errno.EDOM
- def _error_reset():
- rposix.set_errno(0)
- def _likely_raise(errno, x):
- """Call this with errno != 0. It usually raises the proper RPython
- exception, but may also just ignore it and return in case of underflow.
- """
- assert errno
- if errno == ERANGE:
- # We consider underflow to not be an error, like CPython.
- # On some platforms (Ubuntu/ia64) it seems that errno can be
- # set to ERANGE for subnormal results that do *not* underflow
- # to zero. So to be safe, we'll ignore ERANGE whenever the
- # function result is less than one in absolute value.
- if math_fabs(x) < 1.0:
- return
- raise OverflowError("math range error")
- else:
- raise ValueError("math domain error")
- # ____________________________________________________________
- #
- # Custom implementations
- VERY_LARGE_FLOAT = 1.0
- while VERY_LARGE_FLOAT * 100.0 != INFINITY:
- VERY_LARGE_FLOAT *= 64.0
- _lib_isnan = rffi.llexternal("_isnan", [lltype.Float], lltype.Signed,
- compilation_info=eci)
- _lib_finite = rffi.llexternal("_finite", [lltype.Float], lltype.Signed,
- compilation_info=eci)
- def ll_math_isnan(y):
- # By not calling into the external function the JIT can inline this.
- # Floats are awesome.
- if use_library_isinf_isnan and not jit.we_are_jitted():
- return bool(_lib_isnan(y))
- return y != y
- def ll_math_isinf(y):
- if jit.we_are_jitted():
- return (y + VERY_LARGE_FLOAT) == y
- elif use_library_isinf_isnan:
- return not _lib_finite(y) and not _lib_isnan(y)
- else:
- return y == INFINITY or y == -INFINITY
- def ll_math_isfinite(y):
- # Use a custom hack that is reasonably well-suited to the JIT.
- # Floats are awesome (bis).
- if use_library_isinf_isnan and not jit.we_are_jitted():
- return bool(_lib_finite(y))
- z = 0.0 * y
- return z == z # i.e.: z is not a NaN
- ll_math_floor = math_floor
- ll_math_copysign = math_copysign
- def ll_math_atan2(y, x):
- """wrapper for atan2 that deals directly with special cases before
- delegating to the platform libm for the remaining cases. This
- is necessary to get consistent behaviour across platforms.
- Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
- always follow C99.
- """
- if isnan(x):
- return NAN
- if not isfinite(y):
- if isnan(y):
- return NAN
- if isinf(x):
- if math_copysign(1.0, x) == 1.0:
- # atan2(+-inf, +inf) == +-pi/4
- return math_copysign(0.25 * math.pi, y)
- else:
- # atan2(+-inf, -inf) == +-pi*3/4
- return math_copysign(0.75 * math.pi, y)
- # atan2(+-inf, x) == +-pi/2 for finite x
- return math_copysign(0.5 * math.pi, y)
- if isinf(x) or y == 0.0:
- if math_copysign(1.0, x) == 1.0:
- # atan2(+-y, +inf) = atan2(+-0, +x) = +-0.
- return math_copysign(0.0, y)
- else:
- # atan2(+-y, -inf) = atan2(+-0., -x) = +-pi.
- return math_copysign(math.pi, y)
- return math_atan2(y, x)
- # XXX Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
- # log(-ve), log(NaN). For now I'm ignoring this issue as these are a bit
- # more marginal platforms for us.
- def ll_math_frexp(x):
- # deal with special cases directly, to sidestep platform differences
- if not isfinite(x) or not x:
- mantissa = x
- exponent = 0
- else:
- exp_p = lltype.malloc(rffi.INTP.TO, 1, flavor='raw')
- try:
- mantissa = math_frexp(x, exp_p)
- exponent = rffi.cast(lltype.Signed, exp_p[0])
- finally:
- lltype.free(exp_p, flavor='raw')
- return (mantissa, exponent)
- INT_MAX = int(2**31-1)
- INT_MIN = int(-2**31)
- def ll_math_ldexp(x, exp):
- if x == 0.0 or not isfinite(x):
- return x # NaNs, zeros and infinities are returned unchanged
- if exp > INT_MAX:
- # overflow (64-bit platforms only)
- r = math_copysign(INFINITY, x)
- errno = ERANGE
- elif exp < INT_MIN:
- # underflow to +-0 (64-bit platforms only)
- r = math_copysign(0.0, x)
- errno = 0
- else:
- _error_reset()
- r = math_ldexp(x, exp)
- errno = rposix.get_errno()
- if isinf(r):
- errno = ERANGE
- if errno:
- _likely_raise(errno, r)
- return r
- def ll_math_modf(x):
- # some platforms don't do the right thing for NaNs and
- # infinities, so we take care of special cases directly.
- if not isfinite(x):
- if isnan(x):
- return (x, x)
- else: # isinf(x)
- return (math_copysign(0.0, x), x)
- intpart_p = lltype.malloc(rffi.DOUBLEP.TO, 1, flavor='raw')
- try:
- fracpart = math_modf(x, intpart_p)
- intpart = intpart_p[0]
- finally:
- lltype.free(intpart_p, flavor='raw')
- return (fracpart, intpart)
- def ll_math_fmod(x, y):
- # fmod(x, +/-Inf) returns x for finite x.
- if isinf(y) and isfinite(x):
- return x
- _error_reset()
- r = math_fmod(x, y)
- errno = rposix.get_errno()
- if isnan(r):
- if isnan(x) or isnan(y):
- errno = 0
- else:
- errno = EDOM
- if errno:
- _likely_raise(errno, r)
- return r
- def ll_math_hypot(x, y):
- # hypot(x, +/-Inf) returns Inf, even if x is a NaN.
- if isinf(x):
- return math_fabs(x)
- if isinf(y):
- return math_fabs(y)
- _error_reset()
- r = math_hypot(x, y)
- errno = rposix.get_errno()
- if not isfinite(r):
- if isnan(r):
- if isnan(x) or isnan(y):
- errno = 0
- else:
- errno = EDOM
- else: # isinf(r)
- if isfinite(x) and isfinite(y):
- errno = ERANGE
- else:
- errno = 0
- if errno:
- _likely_raise(errno, r)
- return r
- def ll_math_pow(x, y):
- # deal directly with IEEE specials, to cope with problems on various
- # platforms whose semantics don't exactly match C99
- if isnan(y):
- if x == 1.0:
- return 1.0 # 1**Nan = 1
- return y
- if not isfinite(x):
- if isnan(x):
- if y == 0.0:
- return 1.0 # NaN**0 = 1
- return x
- else: # isinf(x)
- odd_y = not isinf(y) and math_fmod(math_fabs(y), 2.0) == 1.0
- if y > 0.0:
- if odd_y:
- return x
- return math_fabs(x)
- elif y == 0.0:
- return 1.0
- else: # y < 0.0
- if odd_y:
- return math_copysign(0.0, x)
- return 0.0
- if isinf(y):
- if math_fabs(x) == 1.0:
- return 1.0
- elif y > 0.0 and math_fabs(x) > 1.0:
- return y
- elif y < 0.0 and math_fabs(x) < 1.0:
- if x == 0.0:
- raise ValueError("0**-inf: divide by zero")
- return -y # result is +inf
- else:
- return 0.0
- _error_reset()
- r = math_pow(x, y)
- errno = rposix.get_errno()
- if not isfinite(r):
- if isnan(r):
- # a NaN result should arise only from (-ve)**(finite non-integer)
- errno = EDOM
- else: # isinf(r)
- # an infinite result here arises either from:
- # (A) (+/-0.)**negative (-> divide-by-zero)
- # (B) overflow of x**y with x and y finite
- if x == 0.0:
- errno = EDOM
- else:
- errno = ERANGE
- if errno:
- _likely_raise(errno, r)
- return r
- def ll_math_sqrt(x):
- if x < 0.0:
- raise ValueError, "math domain error"
- if isfinite(x):
- return sqrt_nonneg(x)
- return x # +inf or nan
- def ll_math_log(x):
- if x <= 0:
- raise ValueError("math domain error")
- return math_log(x)
- def ll_math_log10(x):
- if x <= 0:
- raise ValueError("math domain error")
- return math_log10(x)
- def ll_math_sin(x):
- if isinf(x):
- raise ValueError("math domain error")
- return math_sin(x)
- def ll_math_cos(x):
- if isinf(x):
- raise ValueError("math domain error")
- return math_cos(x)
- # ____________________________________________________________
- #
- # Default implementations
- def new_unary_math_function(name, can_overflow, c99):
- if sys.platform == 'win32' and c99:
- c_func = math_llexternal(name, [rffi.DOUBLE], rffi.DOUBLE)
- else:
- c_func = llexternal(name, [rffi.DOUBLE], rffi.DOUBLE)
- def ll_math(x):
- _error_reset()
- r = c_func(x)
- # Error checking fun. Copied from CPython 2.6
- errno = rposix.get_errno()
- if not isfinite(r):
- if isnan(r):
- if isnan(x):
- errno = 0
- else:
- errno = EDOM
- else: # isinf(r)
- if not isfinite(x):
- errno = 0
- elif can_overflow:
- errno = ERANGE
- else:
- errno = EDOM
- if errno:
- _likely_raise(errno, r)
- return r
- return func_with_new_name(ll_math, 'll_math_' + name)
- # ____________________________________________________________
- unary_math_functions = [
- 'acos', 'asin', 'atan',
- 'ceil', 'cosh', 'exp', 'fabs',
- 'sinh', 'tan', 'tanh',
- 'acosh', 'asinh', 'atanh', 'log1p', 'expm1',
- ]
- unary_math_functions_can_overflow = [
- 'cosh', 'exp', 'log1p', 'sinh', 'expm1',
- ]
- unary_math_functions_c99 = [
- 'acosh', 'asinh', 'atanh', 'log1p', 'expm1',
- ]
- for name in unary_math_functions:
- can_overflow = name in unary_math_functions_can_overflow
- c99 = name in unary_math_functions_c99
- globals()['ll_math_' + name] = new_unary_math_function(name, can_overflow, c99)