PageRenderTime 82ms CodeModel.GetById 24ms app.highlight 51ms RepoModel.GetById 1ms app.codeStats 0ms

/pypy/module/math/interp_math.py

https://bitbucket.org/dac_io/pypy
Python | 604 lines | 428 code | 58 blank | 118 comment | 78 complexity | c37546f094a8f5e1e92cd583a22eafa9 MD5 | raw file
  1import math
  2import sys
  3
  4from pypy.rlib import rfloat, unroll
  5from pypy.interpreter.error import OperationError
  6from pypy.interpreter.gateway import NoneNotWrapped
  7
  8class State:
  9    def __init__(self, space):
 10        self.w_e = space.wrap(math.e)
 11        self.w_pi = space.wrap(math.pi)
 12def get(space):
 13    return space.fromcache(State)
 14
 15def _get_double(space, w_x):
 16    if space.is_w(space.type(w_x), space.w_float):
 17        return space.float_w(w_x)
 18    else:
 19        return space.float_w(space.float(w_x))
 20
 21def math1(space, f, w_x):
 22    x = _get_double(space, w_x)
 23    try:
 24        y = f(x)
 25    except OverflowError:
 26        raise OperationError(space.w_OverflowError,
 27                             space.wrap("math range error"))
 28    except ValueError:
 29        raise OperationError(space.w_ValueError,
 30                             space.wrap("math domain error"))
 31    return space.wrap(y)
 32math1._annspecialcase_ = 'specialize:arg(1)'
 33
 34def math1_w(space, f, w_x):
 35    x = _get_double(space, w_x)
 36    try:
 37        r = f(x)
 38    except OverflowError:
 39        raise OperationError(space.w_OverflowError,
 40                             space.wrap("math range error"))
 41    except ValueError:
 42        raise OperationError(space.w_ValueError,
 43                             space.wrap("math domain error"))
 44    return r
 45math1_w._annspecialcase_ = 'specialize:arg(1)'
 46
 47def math2(space, f, w_x, w_snd):
 48    x = _get_double(space, w_x)
 49    snd = _get_double(space, w_snd)
 50    try:
 51        r = f(x, snd)
 52    except OverflowError:
 53        raise OperationError(space.w_OverflowError,
 54                             space.wrap("math range error"))
 55    except ValueError:
 56        raise OperationError(space.w_ValueError,
 57                             space.wrap("math domain error"))
 58    return space.wrap(r)
 59math2._annspecialcase_ = 'specialize:arg(1)'
 60
 61def trunc(space, w_x):
 62    """Truncate x."""
 63    return space.trunc(w_x)
 64
 65def copysign(space, w_x, w_y):
 66    """Return x with the sign of y."""
 67    # No exceptions possible.
 68    x = _get_double(space, w_x)
 69    y = _get_double(space, w_y)
 70    return space.wrap(rfloat.copysign(x, y))
 71
 72def isinf(space, w_x):
 73    """Return True if x is infinity."""
 74    return space.wrap(rfloat.isinf(_get_double(space, w_x)))
 75
 76def isnan(space, w_x):
 77    """Return True if x is not a number."""
 78    return space.wrap(rfloat.isnan(_get_double(space, w_x)))
 79
 80def pow(space, w_x, w_y):
 81    """pow(x,y)
 82
 83       Return x**y (x to the power of y).
 84    """
 85    return math2(space, math.pow, w_x, w_y)
 86
 87def cosh(space, w_x):
 88    """cosh(x)
 89
 90       Return the hyperbolic cosine of x.
 91    """
 92    return math1(space, math.cosh, w_x)
 93
 94def ldexp(space, w_x,  w_i):
 95    """ldexp(x, i) -> x * (2**i)
 96    """
 97    x = _get_double(space, w_x)
 98    if (space.isinstance_w(w_i, space.w_int) or
 99        space.isinstance_w(w_i, space.w_long)):
100        try:
101            exp = space.int_w(w_i)
102        except OperationError, e:
103            if not e.match(space, space.w_OverflowError):
104                raise
105            if space.is_true(space.lt(w_i, space.wrap(0))):
106                exp = -sys.maxint
107            else:
108                exp = sys.maxint
109    else:
110        raise OperationError(space.w_TypeError,
111                             space.wrap("integer required for second argument"))
112    try:
113        r = math.ldexp(x, exp)
114    except OverflowError:
115        raise OperationError(space.w_OverflowError,
116                             space.wrap("math range error"))
117    except ValueError:
118        raise OperationError(space.w_ValueError,
119                             space.wrap("math domain error"))
120    return space.wrap(r)
121
122def hypot(space, w_x, w_y):
123    """hypot(x,y)
124
125       Return the Euclidean distance, sqrt(x*x + y*y).
126    """
127    return math2(space, math.hypot, w_x, w_y)
128
129def tan(space, w_x):
130    """tan(x)
131
132       Return the tangent of x (measured in radians).
133    """
134    return math1(space, math.tan, w_x)
135
136def asin(space, w_x):
137    """asin(x)
138
139       Return the arc sine (measured in radians) of x.
140    """
141    return math1(space, math.asin, w_x)
142
143def fabs(space, w_x):
144    """fabs(x)
145
146       Return the absolute value of the float x.
147    """
148    return math1(space, math.fabs, w_x)
149
150def floor(space, w_x):
151    """floor(x)
152
153       Return the floor of x as a float.
154       This is the largest integral value <= x.
155    """
156    x = _get_double(space, w_x)
157    return space.wrap(math.floor(x))
158
159def sqrt(space, w_x):
160    """sqrt(x)
161
162       Return the square root of x.
163    """
164    return math1(space, math.sqrt, w_x)
165
166def frexp(space, w_x):
167    """frexp(x)
168
169       Return the mantissa and exponent of x, as pair (m, e).
170       m is a float and e is an int, such that x = m * 2.**e.
171       If x is 0, m and e are both 0.  Else 0.5 <= abs(m) < 1.0.
172    """
173    mant, expo = math1_w(space, math.frexp, w_x)
174    return space.newtuple([space.wrap(mant), space.wrap(expo)])
175
176degToRad = math.pi / 180.0
177
178def degrees(space, w_x):
179    """degrees(x) -> converts angle x from radians to degrees
180    """
181    return space.wrap(_get_double(space, w_x) / degToRad)
182
183def _log_any(space, w_x, base):
184    # base is supposed to be positive or 0.0, which means we use e
185    try:
186        if space.is_true(space.isinstance(w_x, space.w_long)):
187            # special case to support log(extremely-large-long)
188            num = space.bigint_w(w_x)
189            result = num.log(base)
190        else:
191            x = _get_double(space, w_x)
192            if base == 10.0:
193                result = math.log10(x)
194            else:
195                result = math.log(x)
196                if base != 0.0:
197                    den = math.log(base)
198                    result /= den
199    except OverflowError:
200        raise OperationError(space.w_OverflowError,
201                             space.wrap('math range error'))
202    except ValueError:
203        raise OperationError(space.w_ValueError,
204                             space.wrap('math domain error'))
205    return space.wrap(result)
206
207def log(space, w_x, w_base=NoneNotWrapped):
208    """log(x[, base]) -> the logarithm of x to the given base.
209       If the base not specified, returns the natural logarithm (base e) of x.
210    """
211    if w_base is None:
212        base = 0.0
213    else:
214        base = _get_double(space, w_base)
215        if base <= 0.0:
216            # just for raising the proper errors
217            return math1(space, math.log, w_base)
218    return _log_any(space, w_x, base)
219
220def log10(space, w_x):
221    """log10(x) -> the base 10 logarithm of x.
222    """
223    return _log_any(space, w_x, 10.0)
224
225def fmod(space, w_x, w_y):
226    """fmod(x,y)
227
228       Return fmod(x, y), according to platform C.  x % y may differ.
229    """
230    return math2(space, math.fmod, w_x, w_y)
231
232def atan(space, w_x):
233    """atan(x)
234
235       Return the arc tangent (measured in radians) of x.
236    """
237    return math1(space, math.atan, w_x)
238
239def ceil(space, w_x):
240    """ceil(x)
241
242       Return the ceiling of x as a float.
243       This is the smallest integral value >= x.
244    """
245    return math1(space, math.ceil, w_x)
246
247def sinh(space, w_x):
248    """sinh(x)
249
250       Return the hyperbolic sine of x.
251    """
252    return math1(space, math.sinh, w_x)
253
254def cos(space, w_x):
255    """cos(x)
256
257       Return the cosine of x (measured in radians).
258    """
259    return math1(space, math.cos, w_x)
260
261def tanh(space, w_x):
262    """tanh(x)
263
264       Return the hyperbolic tangent of x.
265    """
266    return math1(space, math.tanh, w_x)
267
268def radians(space, w_x):
269    """radians(x) -> converts angle x from degrees to radians
270    """
271    return space.wrap(_get_double(space, w_x) * degToRad)
272
273def sin(space, w_x):
274    """sin(x)
275
276       Return the sine of x (measured in radians).
277    """
278    return math1(space, math.sin, w_x)
279
280def atan2(space, w_y, w_x):
281    """atan2(y, x)
282
283       Return the arc tangent (measured in radians) of y/x.
284       Unlike atan(y/x), the signs of both x and y are considered.
285    """
286    return math2(space, math.atan2, w_y,  w_x)
287
288def modf(space, w_x):
289    """modf(x)
290
291       Return the fractional and integer parts of x.  Both results carry the sign
292       of x.  The integer part is returned as a real.
293    """
294    frac, intpart = math1_w(space, math.modf, w_x)
295    return space.newtuple([space.wrap(frac), space.wrap(intpart)])
296
297def exp(space, w_x):
298    """exp(x)
299
300       Return e raised to the power of x.
301    """
302    return math1(space, math.exp, w_x)
303
304def acos(space, w_x):
305    """acos(x)
306
307       Return the arc cosine (measured in radians) of x.
308    """
309    return math1(space, math.acos, w_x)
310
311def fsum(space, w_iterable):
312    """Sum an iterable of floats, trying to keep precision."""
313    w_iter = space.iter(w_iterable)
314    inf_sum = special_sum = 0.0
315    partials = []
316    while True:
317        try:
318            w_value = space.next(w_iter)
319        except OperationError, e:
320            if not e.match(space, space.w_StopIteration):
321                raise
322            break
323        v = _get_double(space, w_value)
324        original = v
325        added = 0
326        for y in partials:
327            if abs(v) < abs(y):
328                v, y = y, v
329            hi = v + y
330            yr = hi - v
331            lo = y - yr
332            if lo != 0.0:
333                partials[added] = lo
334                added += 1
335            v = hi
336        del partials[added:]
337        if v != 0.0:
338            if rfloat.isinf(v) or rfloat.isnan(v):
339                if (not rfloat.isinf(original) and
340                    not rfloat.isnan(original)):
341                    raise OperationError(space.w_OverflowError,
342                                         space.wrap("intermediate overflow"))
343                if rfloat.isinf(original):
344                    inf_sum += original
345                special_sum += original
346                del partials[:]
347            else:
348                partials.append(v)
349    if special_sum != 0.0:
350        if rfloat.isnan(special_sum):
351            raise OperationError(space.w_ValueError, space.wrap("-inf + inf"))
352        return space.wrap(special_sum)
353    hi = 0.0
354    if partials:
355        hi = partials[-1]
356        j = 0
357        lo = 0
358        for j in range(len(partials) - 2, -1, -1):
359            v = hi
360            y = partials[j]
361            assert abs(y) < abs(v)
362            hi = v + y
363            yr = hi - v
364            lo = y - yr
365            if lo != 0.0:
366                break
367        if j > 0 and (lo < 0.0 and partials[j - 1] < 0.0 or
368                      lo > 0.0 and partials[j - 1] > 0.0):
369            y = lo * 2.0
370            v = hi + y
371            yr = v - hi
372            if y == yr:
373                hi = v
374    return space.wrap(hi)
375
376def log1p(space, w_x):
377    """Find log(x + 1)."""
378    return math1(space, rfloat.log1p, w_x)
379
380def acosh(space, w_x):
381    """Inverse hyperbolic cosine"""
382    return math1(space, rfloat.acosh, w_x)
383
384def asinh(space, w_x):
385    """Inverse hyperbolic sine"""
386    return math1(space, rfloat.asinh, w_x)
387
388def atanh(space, w_x):
389    """Inverse hyperbolic tangent"""
390    return math1(space, rfloat.atanh, w_x)
391
392def expm1(space, w_x):
393    """exp(x) - 1"""
394    return math1(space, rfloat.expm1, w_x)
395
396def erf(space, w_x):
397    """The error function"""
398    return math1(space, _erf, w_x)
399
400def erfc(space, w_x):
401    """The complementary error function"""
402    return math1(space, _erfc, w_x)
403
404def gamma(space, w_x):
405    """Compute the gamma function for x."""
406    return math1(space, _gamma, w_x)
407
408def lgamma(space, w_x):
409    """Compute the natural logarithm of the gamma function for x."""
410    return math1(space, _lgamma, w_x)
411
412# Implementation of the error function, the complimentary error function, the
413# gamma function, and the natural log of the gamma function.  These exist in
414# libm, but I hear those implementations are horrible.
415
416ERF_SERIES_CUTOFF = 1.5
417ERF_SERIES_TERMS = 25
418ERFC_CONTFRAC_CUTOFF = 30.
419ERFC_CONTFRAC_TERMS = 50
420_sqrtpi = 1.772453850905516027298167483341145182798
421
422def _erf_series(x):
423    x2 = x * x
424    acc = 0.
425    fk = ERF_SERIES_TERMS + .5
426    for i in range(ERF_SERIES_TERMS):
427        acc = 2.0 + x2 * acc / fk
428        fk -= 1.
429    return acc * x * math.exp(-x2) / _sqrtpi
430
431def _erfc_contfrac(x):
432    if x >= ERFC_CONTFRAC_CUTOFF:
433        return 0.
434    x2 = x * x
435    a = 0.
436    da = .5
437    p = 1.
438    p_last = 0.
439    q = da + x2
440    q_last = 1.
441    for i in range(ERFC_CONTFRAC_TERMS):
442        a += da
443        da += 2.
444        b = da + x2
445        p_last, p = p, b * p - a * p_last
446        q_last, q = q, b * q - a * q_last
447    return p / q * x * math.exp(-x2) / _sqrtpi
448
449def _erf(x):
450    if rfloat.isnan(x):
451        return x
452    absx = abs(x)
453    if absx < ERF_SERIES_CUTOFF:
454        return _erf_series(x)
455    else:
456        cf = _erfc_contfrac(absx)
457        return 1. - cf if x > 0. else cf - 1.
458
459def _erfc(x):
460    if rfloat.isnan(x):
461        return x
462    absx = abs(x)
463    if absx < ERF_SERIES_CUTOFF:
464        return 1. - _erf_series(x)
465    else:
466        cf = _erfc_contfrac(absx)
467        return cf if x > 0. else 2. - cf
468
469def _sinpi(x):
470    y = math.fmod(abs(x), 2.)
471    n = int(rfloat.round_away(2. * y))
472    if n == 0:
473        r = math.sin(math.pi * y)
474    elif n == 1:
475        r = math.cos(math.pi * (y - .5))
476    elif n == 2:
477        r = math.sin(math.pi * (1. - y))
478    elif n == 3:
479        r = -math.cos(math.pi * (y - 1.5))
480    elif n == 4:
481        r = math.sin(math.pi * (y - 2.))
482    else:
483        raise AssertionError("should not reach")
484    return rfloat.copysign(1., x) * r
485
486_lanczos_g = 6.024680040776729583740234375
487_lanczos_g_minus_half = 5.524680040776729583740234375
488_lanczos_num_coeffs = [
489    23531376880.410759688572007674451636754734846804940,
490    42919803642.649098768957899047001988850926355848959,
491    35711959237.355668049440185451547166705960488635843,
492    17921034426.037209699919755754458931112671403265390,
493    6039542586.3520280050642916443072979210699388420708,
494    1439720407.3117216736632230727949123939715485786772,
495    248874557.86205415651146038641322942321632125127801,
496    31426415.585400194380614231628318205362874684987640,
497    2876370.6289353724412254090516208496135991145378768,
498    186056.26539522349504029498971604569928220784236328,
499    8071.6720023658162106380029022722506138218516325024,
500    210.82427775157934587250973392071336271166969580291,
501    2.5066282746310002701649081771338373386264310793408
502]
503_lanczos_den_coeffs = [
504    0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
505    13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0]
506LANCZOS_N = len(_lanczos_den_coeffs)
507_lanczos_n_iter = unroll.unrolling_iterable(range(LANCZOS_N))
508_lanczos_n_iter_back = unroll.unrolling_iterable(range(LANCZOS_N - 1, -1, -1))
509_gamma_integrals = [
510    1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
511    3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
512    1307674368000.0, 20922789888000.0, 355687428096000.0,
513    6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
514    51090942171709440000.0, 1124000727777607680000.0]
515
516def _lanczos_sum(x):
517    num = 0.
518    den = 0.
519    assert x > 0.
520    if x < 5.:
521        for i in _lanczos_n_iter_back:
522            num = num * x + _lanczos_num_coeffs[i]
523            den = den * x + _lanczos_den_coeffs[i]
524    else:
525        for i in _lanczos_n_iter:
526            num = num / x + _lanczos_num_coeffs[i]
527            den = den / x + _lanczos_den_coeffs[i]
528    return num / den
529
530def _gamma(x):
531    if rfloat.isnan(x) or (rfloat.isinf(x) and x > 0.):
532        return x
533    if rfloat.isinf(x):
534        raise ValueError("math domain error")
535    if x == 0.:
536        raise ValueError("math domain error")
537    if x == math.floor(x):
538        if x < 0.:
539            raise ValueError("math domain error")
540        if x < len(_gamma_integrals):
541            return _gamma_integrals[int(x) - 1]
542    absx = abs(x)
543    if absx < 1e-20:
544        r = 1. / x
545        if rfloat.isinf(r):
546            raise OverflowError("math range error")
547        return r
548    if absx > 200.:
549        if x < 0.:
550            return 0. / -_sinpi(x)
551        else:
552            raise OverflowError("math range error")
553    y = absx + _lanczos_g_minus_half
554    if absx > _lanczos_g_minus_half:
555        q = y - absx
556        z = q - _lanczos_g_minus_half
557    else:
558        q = y - _lanczos_g_minus_half
559        z = q - absx
560    z = z * _lanczos_g / y
561    if x < 0.:
562        r = -math.pi / _sinpi(absx) / absx * math.exp(y) / _lanczos_sum(absx)
563        r -= z * r
564        if absx < 140.:
565            r /= math.pow(y, absx - .5)
566        else:
567            sqrtpow = math.pow(y, absx / 2. - .25)
568            r /= sqrtpow
569            r /= sqrtpow
570    else:
571        r = _lanczos_sum(absx) / math.exp(y)
572        r += z * r
573        if absx < 140.:
574            r *= math.pow(y, absx - .5)
575        else:
576            sqrtpow = math.pow(y, absx / 2. - .25)
577            r *= sqrtpow
578            r *= sqrtpow
579    if rfloat.isinf(r):
580        raise OverflowError("math range error")
581    return r
582
583def _lgamma(x):
584    if rfloat.isnan(x):
585        return x
586    if rfloat.isinf(x):
587        return rfloat.INFINITY
588    if x == math.floor(x) and x <= 2.:
589        if x <= 0.:
590            raise ValueError("math range error")
591        return 0.
592    absx = abs(x)
593    if absx < 1e-20:
594        return -math.log(absx)
595    if x > 0.:
596        r = (math.log(_lanczos_sum(x)) - _lanczos_g + (x - .5) *
597             (math.log(x + _lanczos_g - .5) - 1))
598    else:
599        r = (math.log(math.pi) - math.log(abs(_sinpi(absx))) - math.log(absx) -
600             (math.log(_lanczos_sum(absx)) - _lanczos_g +
601              (absx - .5) * (math.log(absx + _lanczos_g - .5) - 1)))
602    if rfloat.isinf(r):
603        raise OverflowError("math domain error")
604    return r