PageRenderTime 105ms CodeModel.GetById 40ms app.highlight 39ms RepoModel.GetById 21ms app.codeStats 0ms

/pypy/rpython/lltypesystem/module/ll_math.py

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