/rpython/rtyper/lltypesystem/module/ll_math.py

https://bitbucket.org/pypy/pypy/ · Python · 429 lines · 302 code · 54 blank · 73 comment · 101 complexity · 74f27ffe1a3c50640b53545ef532c199 MD5 · raw file

  1. import errno
  2. import math
  3. import py
  4. import sys
  5. from rpython.translator import cdir
  6. from rpython.rlib import jit, rposix
  7. from rpython.rlib.rfloat import INFINITY, NAN, isfinite, isinf, isnan
  8. from rpython.rlib.rposix import UNDERSCORE_ON_WIN32
  9. from rpython.rtyper.lltypesystem import lltype, rffi
  10. from rpython.tool.sourcetools import func_with_new_name
  11. from rpython.translator.tool.cbuild import ExternalCompilationInfo
  12. from rpython.translator.platform import platform
  13. use_library_isinf_isnan = False
  14. if sys.platform == "win32":
  15. if platform.name == "msvc":
  16. # When compiled with /O2 or /Oi (enable intrinsic functions)
  17. # It's no more possible to take the address of some math functions.
  18. # Ensure that the compiler chooses real functions instead.
  19. eci = ExternalCompilationInfo(
  20. includes=['math.h', 'float.h'],
  21. post_include_bits=['#pragma function(floor)'],
  22. )
  23. use_library_isinf_isnan = True
  24. else:
  25. eci = ExternalCompilationInfo()
  26. # Some math functions are C99 and not defined by the Microsoft compiler
  27. cdir = py.path.local(cdir)
  28. math_eci = ExternalCompilationInfo(
  29. include_dirs=[cdir],
  30. includes=['src/ll_math.h'],
  31. separate_module_files=[cdir.join('src', 'll_math.c')],
  32. )
  33. math_prefix = '_pypy_math_'
  34. else:
  35. eci = ExternalCompilationInfo(
  36. libraries=['m'])
  37. math_eci = eci
  38. math_prefix = ''
  39. def llexternal(name, ARGS, RESULT, **kwargs):
  40. return rffi.llexternal(name, ARGS, RESULT, compilation_info=eci,
  41. sandboxsafe=True, **kwargs)
  42. def math_llexternal(name, ARGS, RESULT, **kwargs):
  43. return rffi.llexternal(math_prefix + name, ARGS, RESULT,
  44. compilation_info=math_eci,
  45. sandboxsafe=True, **kwargs)
  46. math_fabs = llexternal('fabs', [rffi.DOUBLE], rffi.DOUBLE)
  47. math_log = llexternal('log', [rffi.DOUBLE], rffi.DOUBLE)
  48. math_log10 = llexternal('log10', [rffi.DOUBLE], rffi.DOUBLE)
  49. math_log1p = math_llexternal('log1p', [rffi.DOUBLE], rffi.DOUBLE)
  50. math_copysign = llexternal(UNDERSCORE_ON_WIN32 + 'copysign',
  51. [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE,
  52. elidable_function=True)
  53. math_atan2 = llexternal('atan2', [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE)
  54. math_frexp = llexternal('frexp', [rffi.DOUBLE, rffi.INTP], rffi.DOUBLE)
  55. math_modf = llexternal('modf', [rffi.DOUBLE, rffi.DOUBLEP], rffi.DOUBLE)
  56. math_ldexp = llexternal('ldexp', [rffi.DOUBLE, rffi.INT], rffi.DOUBLE,
  57. save_err=rffi.RFFI_FULL_ERRNO_ZERO)
  58. math_pow = llexternal('pow', [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE,
  59. save_err=rffi.RFFI_FULL_ERRNO_ZERO)
  60. math_fmod = llexternal('fmod', [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE,
  61. save_err=rffi.RFFI_FULL_ERRNO_ZERO)
  62. math_hypot = llexternal(UNDERSCORE_ON_WIN32 + 'hypot',
  63. [rffi.DOUBLE, rffi.DOUBLE], rffi.DOUBLE,
  64. save_err=rffi.RFFI_FULL_ERRNO_ZERO)
  65. math_floor = llexternal('floor', [rffi.DOUBLE], rffi.DOUBLE, elidable_function=True)
  66. math_sqrt = llexternal('sqrt', [rffi.DOUBLE], rffi.DOUBLE)
  67. math_sin = llexternal('sin', [rffi.DOUBLE], rffi.DOUBLE, elidable_function=True)
  68. math_cos = llexternal('cos', [rffi.DOUBLE], rffi.DOUBLE, elidable_function=True)
  69. @jit.elidable
  70. def sqrt_nonneg(x):
  71. return math_sqrt(x)
  72. sqrt_nonneg.oopspec = "math.sqrt_nonneg(x)"
  73. # ____________________________________________________________
  74. #
  75. # Error handling functions
  76. ERANGE = errno.ERANGE
  77. EDOM = errno.EDOM
  78. def _likely_raise(errno, x):
  79. """Call this with errno != 0. It usually raises the proper RPython
  80. exception, but may also just ignore it and return in case of underflow.
  81. """
  82. assert errno
  83. if errno == ERANGE:
  84. # We consider underflow to not be an error, like CPython.
  85. # On some platforms (Ubuntu/ia64) it seems that errno can be
  86. # set to ERANGE for subnormal results that do *not* underflow
  87. # to zero. So to be safe, we'll ignore ERANGE whenever the
  88. # function result is less than one in absolute value.
  89. if math_fabs(x) < 1.0:
  90. return
  91. raise OverflowError("math range error")
  92. else:
  93. raise ValueError("math domain error")
  94. # ____________________________________________________________
  95. #
  96. # Custom implementations
  97. VERY_LARGE_FLOAT = 1.0
  98. while VERY_LARGE_FLOAT * 100.0 != INFINITY:
  99. VERY_LARGE_FLOAT *= 64.0
  100. _lib_isnan = llexternal('_isnan', [lltype.Float], lltype.Signed)
  101. _lib_finite = llexternal('_finite', [lltype.Float], lltype.Signed)
  102. def ll_math_isnan(y):
  103. # By not calling into the external function the JIT can inline this.
  104. # Floats are awesome.
  105. if use_library_isinf_isnan and not jit.we_are_jitted():
  106. return bool(_lib_isnan(y))
  107. return y != y
  108. def ll_math_isinf(y):
  109. if jit.we_are_jitted():
  110. return (y + VERY_LARGE_FLOAT) == y
  111. elif use_library_isinf_isnan:
  112. return not _lib_finite(y) and not _lib_isnan(y)
  113. else:
  114. return y == INFINITY or y == -INFINITY
  115. def ll_math_isfinite(y):
  116. # Use a custom hack that is reasonably well-suited to the JIT.
  117. # Floats are awesome (bis).
  118. if use_library_isinf_isnan and not jit.we_are_jitted():
  119. return bool(_lib_finite(y))
  120. return (y - y) == 0.0 # (y - y) is NaN if y is an infinite or NaN
  121. ll_math_floor = math_floor
  122. ll_math_copysign = math_copysign
  123. def ll_math_atan2(y, x):
  124. """wrapper for atan2 that deals directly with special cases before
  125. delegating to the platform libm for the remaining cases. This
  126. is necessary to get consistent behaviour across platforms.
  127. Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
  128. always follow C99.
  129. """
  130. if isnan(x):
  131. return NAN
  132. if not isfinite(y):
  133. if isnan(y):
  134. return NAN
  135. if isinf(x):
  136. if math_copysign(1.0, x) == 1.0:
  137. # atan2(+-inf, +inf) == +-pi/4
  138. return math_copysign(0.25 * math.pi, y)
  139. else:
  140. # atan2(+-inf, -inf) == +-pi*3/4
  141. return math_copysign(0.75 * math.pi, y)
  142. # atan2(+-inf, x) == +-pi/2 for finite x
  143. return math_copysign(0.5 * math.pi, y)
  144. if isinf(x) or y == 0.0:
  145. if math_copysign(1.0, x) == 1.0:
  146. # atan2(+-y, +inf) = atan2(+-0, +x) = +-0.
  147. return math_copysign(0.0, y)
  148. else:
  149. # atan2(+-y, -inf) = atan2(+-0., -x) = +-pi.
  150. return math_copysign(math.pi, y)
  151. return math_atan2(y, x)
  152. # XXX Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
  153. # log(-ve), log(NaN). For now I'm ignoring this issue as these are a bit
  154. # more marginal platforms for us.
  155. def ll_math_frexp(x):
  156. # deal with special cases directly, to sidestep platform differences
  157. if not isfinite(x) or not x:
  158. mantissa = x
  159. exponent = 0
  160. else:
  161. exp_p = lltype.malloc(rffi.INTP.TO, 1, flavor='raw')
  162. try:
  163. mantissa = math_frexp(x, exp_p)
  164. exponent = rffi.cast(lltype.Signed, exp_p[0])
  165. finally:
  166. lltype.free(exp_p, flavor='raw')
  167. return (mantissa, exponent)
  168. INT_MAX = int(2**31-1)
  169. INT_MIN = int(-2**31)
  170. def ll_math_ldexp(x, exp):
  171. if x == 0.0 or not isfinite(x):
  172. return x # NaNs, zeros and infinities are returned unchanged
  173. if exp > INT_MAX:
  174. # overflow (64-bit platforms only)
  175. r = math_copysign(INFINITY, x)
  176. errno = ERANGE
  177. elif exp < INT_MIN:
  178. # underflow to +-0 (64-bit platforms only)
  179. r = math_copysign(0.0, x)
  180. errno = 0
  181. else:
  182. r = math_ldexp(x, exp)
  183. errno = rposix.get_saved_errno()
  184. if isinf(r):
  185. errno = ERANGE
  186. if errno:
  187. _likely_raise(errno, r)
  188. return r
  189. def ll_math_modf(x):
  190. # some platforms don't do the right thing for NaNs and
  191. # infinities, so we take care of special cases directly.
  192. if not isfinite(x):
  193. if isnan(x):
  194. return (x, x)
  195. else: # isinf(x)
  196. return (math_copysign(0.0, x), x)
  197. intpart_p = lltype.malloc(rffi.DOUBLEP.TO, 1, flavor='raw')
  198. try:
  199. fracpart = math_modf(x, intpart_p)
  200. intpart = intpart_p[0]
  201. finally:
  202. lltype.free(intpart_p, flavor='raw')
  203. return (fracpart, intpart)
  204. def ll_math_fmod(x, y):
  205. # fmod(x, +/-Inf) returns x for finite x.
  206. if isinf(y) and isfinite(x):
  207. return x
  208. r = math_fmod(x, y)
  209. errno = rposix.get_saved_errno()
  210. if isnan(r):
  211. if isnan(x) or isnan(y):
  212. errno = 0
  213. else:
  214. errno = EDOM
  215. if errno:
  216. _likely_raise(errno, r)
  217. return r
  218. def ll_math_hypot(x, y):
  219. # hypot(x, +/-Inf) returns Inf, even if x is a NaN.
  220. if isinf(x):
  221. return math_fabs(x)
  222. if isinf(y):
  223. return math_fabs(y)
  224. r = math_hypot(x, y)
  225. errno = rposix.get_saved_errno()
  226. if not isfinite(r):
  227. if isnan(r):
  228. if isnan(x) or isnan(y):
  229. errno = 0
  230. else:
  231. errno = EDOM
  232. else: # isinf(r)
  233. if isfinite(x) and isfinite(y):
  234. errno = ERANGE
  235. else:
  236. errno = 0
  237. if errno:
  238. _likely_raise(errno, r)
  239. return r
  240. def ll_math_pow(x, y):
  241. # deal directly with IEEE specials, to cope with problems on various
  242. # platforms whose semantics don't exactly match C99
  243. if isnan(y):
  244. if x == 1.0:
  245. return 1.0 # 1**Nan = 1
  246. return y
  247. if not isfinite(x):
  248. if isnan(x):
  249. if y == 0.0:
  250. return 1.0 # NaN**0 = 1
  251. return x
  252. else: # isinf(x)
  253. odd_y = not isinf(y) and math_fmod(math_fabs(y), 2.0) == 1.0
  254. if y > 0.0:
  255. if odd_y:
  256. return x
  257. return math_fabs(x)
  258. elif y == 0.0:
  259. return 1.0
  260. else: # y < 0.0
  261. if odd_y:
  262. return math_copysign(0.0, x)
  263. return 0.0
  264. if isinf(y):
  265. if math_fabs(x) == 1.0:
  266. return 1.0
  267. elif y > 0.0 and math_fabs(x) > 1.0:
  268. return y
  269. elif y < 0.0 and math_fabs(x) < 1.0:
  270. if x == 0.0:
  271. raise ValueError("0**-inf: divide by zero")
  272. return -y # result is +inf
  273. else:
  274. return 0.0
  275. r = math_pow(x, y)
  276. errno = rposix.get_saved_errno()
  277. if not isfinite(r):
  278. if isnan(r):
  279. # a NaN result should arise only from (-ve)**(finite non-integer)
  280. errno = EDOM
  281. else: # isinf(r)
  282. # an infinite result here arises either from:
  283. # (A) (+/-0.)**negative (-> divide-by-zero)
  284. # (B) overflow of x**y with x and y finite
  285. if x == 0.0:
  286. errno = EDOM
  287. else:
  288. errno = ERANGE
  289. if errno:
  290. _likely_raise(errno, r)
  291. return r
  292. def ll_math_sqrt(x):
  293. if x < 0.0:
  294. raise ValueError("math domain error")
  295. if isfinite(x):
  296. return sqrt_nonneg(x)
  297. return x # +inf or nan
  298. def ll_math_log(x):
  299. if x <= 0:
  300. raise ValueError("math domain error")
  301. return math_log(x)
  302. def ll_math_log10(x):
  303. if x <= 0:
  304. raise ValueError("math domain error")
  305. return math_log10(x)
  306. def ll_math_log1p(x):
  307. if x == 0.0:
  308. return x # returns 0.0 or -0.0
  309. if x <= -1.0:
  310. if x == -1:
  311. raise OverflowError("math range error")
  312. raise ValueError("math domain error")
  313. return math_log1p(x)
  314. def ll_math_sin(x):
  315. if isinf(x):
  316. raise ValueError("math domain error")
  317. return math_sin(x)
  318. def ll_math_cos(x):
  319. if isinf(x):
  320. raise ValueError("math domain error")
  321. return math_cos(x)
  322. # ____________________________________________________________
  323. #
  324. # Default implementations
  325. def new_unary_math_function(name, can_overflow, c99):
  326. if sys.platform == 'win32' and c99:
  327. c_func = math_llexternal(name, [rffi.DOUBLE], rffi.DOUBLE,
  328. save_err=rffi.RFFI_FULL_ERRNO_ZERO)
  329. else:
  330. c_func = llexternal(name, [rffi.DOUBLE], rffi.DOUBLE,
  331. save_err=rffi.RFFI_FULL_ERRNO_ZERO)
  332. def ll_math(x):
  333. r = c_func(x)
  334. # Error checking fun. Copied from CPython 2.6
  335. errno = rposix.get_saved_errno()
  336. if not isfinite(r):
  337. if isnan(r):
  338. if isnan(x):
  339. errno = 0
  340. else:
  341. errno = EDOM
  342. else: # isinf(r)
  343. if not isfinite(x):
  344. errno = 0
  345. elif can_overflow:
  346. errno = ERANGE
  347. else:
  348. errno = EDOM
  349. if errno:
  350. _likely_raise(errno, r)
  351. return r
  352. return func_with_new_name(ll_math, 'll_math_' + name)
  353. # ____________________________________________________________
  354. unary_math_functions = [
  355. 'acos', 'asin', 'atan',
  356. 'ceil', 'cosh', 'exp', 'fabs',
  357. 'sinh', 'tan', 'tanh',
  358. 'acosh', 'asinh', 'atanh', 'expm1',
  359. ]
  360. unary_math_functions_can_overflow = [
  361. 'cosh', 'exp', 'sinh', 'expm1',
  362. ]
  363. unary_math_functions_c99 = [
  364. 'acosh', 'asinh', 'atanh', 'expm1',
  365. ]
  366. for name in unary_math_functions:
  367. can_overflow = name in unary_math_functions_can_overflow
  368. c99 = name in unary_math_functions_c99
  369. globals()['ll_math_' + name] = new_unary_math_function(name, can_overflow,
  370. c99)