[pypy-commit] pypy numpypy-complex2: remove duplicate code
mattip
noreply at buildbot.pypy.org
Wed Aug 29 00:16:30 CEST 2012
Author: mattip <matti.picus at gmail.com>
Branch: numpypy-complex2
Changeset: r56922:0ed79855240e
Date: 2012-08-29 01:13 +0300
http://bitbucket.org/pypy/pypy/changeset/0ed79855240e/
Log: remove duplicate code
diff --git a/pypy/module/cmath/interp_cmath.py b/pypy/module/cmath/interp_cmath.py
--- a/pypy/module/cmath/interp_cmath.py
+++ b/pypy/module/cmath/interp_cmath.py
@@ -1,28 +1,10 @@
import math
-from math import fabs
from pypy.rlib.objectmodel import specialize
-from pypy.rlib.rfloat import copysign, asinh, log1p, isinf, isnan
from pypy.tool.sourcetools import func_with_new_name
from pypy.interpreter.error import OperationError
from pypy.interpreter.gateway import NoneNotWrapped
from pypy.module.cmath import names_and_docstrings
-from pypy.rlib.constant import DBL_MIN, CM_SCALE_UP, CM_SCALE_DOWN
-from pypy.rlib.constant import CM_LARGE_DOUBLE, DBL_MANT_DIG
-from pypy.rlib.constant import M_LN2, M_LN10
-from pypy.rlib.constant import CM_SQRT_LARGE_DOUBLE, CM_SQRT_DBL_MIN
-from pypy.rlib.constant import CM_LOG_LARGE_DOUBLE
-from pypy.rlib.special_value import isfinite, special_type, INF, NAN
-from pypy.rlib.special_value import sqrt_special_values
-from pypy.rlib.special_value import acos_special_values
-from pypy.rlib.special_value import acosh_special_values
-from pypy.rlib.special_value import asinh_special_values
-from pypy.rlib.special_value import atanh_special_values
-from pypy.rlib.special_value import log_special_values
-from pypy.rlib.special_value import exp_special_values
-from pypy.rlib.special_value import cosh_special_values
-from pypy.rlib.special_value import sinh_special_values
-from pypy.rlib.special_value import tanh_special_values
-from pypy.rlib.special_value import rect_special_values
+from pypy.rlib import rcomplex
pi = math.pi
e = math.e
@@ -56,235 +38,40 @@
def c_neg(x, y):
- return (-x, -y)
+ return rcomplex.c_neg(x,y)
@unaryfn
def c_sqrt(x, y):
- # Method: use symmetries to reduce to the case when x = z.real and y
- # = z.imag are nonnegative. Then the real part of the result is
- # given by
- #
- # s = sqrt((x + hypot(x, y))/2)
- #
- # and the imaginary part is
- #
- # d = (y/2)/s
- #
- # If either x or y is very large then there's a risk of overflow in
- # computation of the expression x + hypot(x, y). We can avoid this
- # by rewriting the formula for s as:
- #
- # s = 2*sqrt(x/8 + hypot(x/8, y/8))
- #
- # This costs us two extra multiplications/divisions, but avoids the
- # overhead of checking for x and y large.
- #
- # If both x and y are subnormal then hypot(x, y) may also be
- # subnormal, so will lack full precision. We solve this by rescaling
- # x and y by a sufficiently large power of 2 to ensure that x and y
- # are normal.
-
- if not isfinite(x) or not isfinite(y):
- return sqrt_special_values[special_type(x)][special_type(y)]
-
- if x == 0. and y == 0.:
- return (0., y)
-
- ax = fabs(x)
- ay = fabs(y)
-
- if ax < DBL_MIN and ay < DBL_MIN and (ax > 0. or ay > 0.):
- # here we catch cases where hypot(ax, ay) is subnormal
- ax = math.ldexp(ax, CM_SCALE_UP)
- ay1= math.ldexp(ay, CM_SCALE_UP)
- s = math.ldexp(math.sqrt(ax + math.hypot(ax, ay1)),
- CM_SCALE_DOWN)
- else:
- ax /= 8.
- s = 2.*math.sqrt(ax + math.hypot(ax, ay/8.))
-
- d = ay/(2.*s)
-
- if x >= 0.:
- return (s, copysign(d, y))
- else:
- return (d, copysign(s, y))
-
+ return rcomplex.c_sqrt(x,y)
@unaryfn
def c_acos(x, y):
- if not isfinite(x) or not isfinite(y):
- return acos_special_values[special_type(x)][special_type(y)]
-
- if fabs(x) > CM_LARGE_DOUBLE or fabs(y) > CM_LARGE_DOUBLE:
- # avoid unnecessary overflow for large arguments
- real = math.atan2(fabs(y), x)
- # split into cases to make sure that the branch cut has the
- # correct continuity on systems with unsigned zeros
- if x < 0.:
- imag = -copysign(math.log(math.hypot(x/2., y/2.)) +
- M_LN2*2., y)
- else:
- imag = copysign(math.log(math.hypot(x/2., y/2.)) +
- M_LN2*2., -y)
- else:
- s1x, s1y = c_sqrt(1.-x, -y)
- s2x, s2y = c_sqrt(1.+x, y)
- real = 2.*math.atan2(s1x, s2x)
- imag = asinh(s2x*s1y - s2y*s1x)
- return (real, imag)
-
+ return rcomplex.c_acos(x,y)
@unaryfn
def c_acosh(x, y):
- # XXX the following two lines seem unnecessary at least on Linux;
- # the tests pass fine without them
- if not isfinite(x) or not isfinite(y):
- return acosh_special_values[special_type(x)][special_type(y)]
-
- if fabs(x) > CM_LARGE_DOUBLE or fabs(y) > CM_LARGE_DOUBLE:
- # avoid unnecessary overflow for large arguments
- real = math.log(math.hypot(x/2., y/2.)) + M_LN2*2.
- imag = math.atan2(y, x)
- else:
- s1x, s1y = c_sqrt(x - 1., y)
- s2x, s2y = c_sqrt(x + 1., y)
- real = asinh(s1x*s2x + s1y*s2y)
- imag = 2.*math.atan2(s1y, s2x)
- return (real, imag)
-
+ return rcomplex.c_acosh(x,y)
@unaryfn
def c_asin(x, y):
- # asin(z) = -i asinh(iz)
- sx, sy = c_asinh(-y, x)
- return (sy, -sx)
-
+ return rcomplex.c_asin(x,y)
@unaryfn
def c_asinh(x, y):
- if not isfinite(x) or not isfinite(y):
- return asinh_special_values[special_type(x)][special_type(y)]
-
- if fabs(x) > CM_LARGE_DOUBLE or fabs(y) > CM_LARGE_DOUBLE:
- if y >= 0.:
- real = copysign(math.log(math.hypot(x/2., y/2.)) +
- M_LN2*2., x)
- else:
- real = -copysign(math.log(math.hypot(x/2., y/2.)) +
- M_LN2*2., -x)
- imag = math.atan2(y, fabs(x))
- else:
- s1x, s1y = c_sqrt(1.+y, -x)
- s2x, s2y = c_sqrt(1.-y, x)
- real = asinh(s1x*s2y - s2x*s1y)
- imag = math.atan2(y, s1x*s2x - s1y*s2y)
- return (real, imag)
-
+ return rcomplex.c_asinh(x,y)
@unaryfn
def c_atan(x, y):
- # atan(z) = -i atanh(iz)
- sx, sy = c_atanh(-y, x)
- return (sy, -sx)
-
+ return rcomplex.c_atan(x,y)
@unaryfn
def c_atanh(x, y):
- if not isfinite(x) or not isfinite(y):
- return atanh_special_values[special_type(x)][special_type(y)]
-
- # Reduce to case where x >= 0., using atanh(z) = -atanh(-z).
- if x < 0.:
- return c_neg(*c_atanh(*c_neg(x, y)))
-
- ay = fabs(y)
- if x > CM_SQRT_LARGE_DOUBLE or ay > CM_SQRT_LARGE_DOUBLE:
- # if abs(z) is large then we use the approximation
- # atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign
- # of y
- h = math.hypot(x/2., y/2.) # safe from overflow
- real = x/4./h/h
- # the two negations in the next line cancel each other out
- # except when working with unsigned zeros: they're there to
- # ensure that the branch cut has the correct continuity on
- # systems that don't support signed zeros
- imag = -copysign(math.pi/2., -y)
- elif x == 1. and ay < CM_SQRT_DBL_MIN:
- # C99 standard says: atanh(1+/-0.) should be inf +/- 0i
- if ay == 0.:
- raise ValueError("math domain error")
- #real = INF
- #imag = y
- else:
- real = -math.log(math.sqrt(ay)/math.sqrt(math.hypot(ay, 2.)))
- imag = copysign(math.atan2(2., -ay) / 2, y)
- else:
- real = log1p(4.*x/((1-x)*(1-x) + ay*ay))/4.
- imag = -math.atan2(-2.*y, (1-x)*(1+x) - ay*ay) / 2.
- return (real, imag)
-
+ return rcomplex.c_atanh(x,y)
@unaryfn
def c_log(x, y):
- # The usual formula for the real part is log(hypot(z.real, z.imag)).
- # There are four situations where this formula is potentially
- # problematic:
- #
- # (1) the absolute value of z is subnormal. Then hypot is subnormal,
- # so has fewer than the usual number of bits of accuracy, hence may
- # have large relative error. This then gives a large absolute error
- # in the log. This can be solved by rescaling z by a suitable power
- # of 2.
- #
- # (2) the absolute value of z is greater than DBL_MAX (e.g. when both
- # z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
- # Again, rescaling solves this.
- #
- # (3) the absolute value of z is close to 1. In this case it's
- # difficult to achieve good accuracy, at least in part because a
- # change of 1ulp in the real or imaginary part of z can result in a
- # change of billions of ulps in the correctly rounded answer.
- #
- # (4) z = 0. The simplest thing to do here is to call the
- # floating-point log with an argument of 0, and let its behaviour
- # (returning -infinity, signaling a floating-point exception, setting
- # errno, or whatever) determine that of c_log. So the usual formula
- # is fine here.
-
- # XXX the following two lines seem unnecessary at least on Linux;
- # the tests pass fine without them
- if not isfinite(x) or not isfinite(y):
- return log_special_values[special_type(x)][special_type(y)]
-
- ax = fabs(x)
- ay = fabs(y)
-
- if ax > CM_LARGE_DOUBLE or ay > CM_LARGE_DOUBLE:
- real = math.log(math.hypot(ax/2., ay/2.)) + M_LN2
- elif ax < DBL_MIN and ay < DBL_MIN:
- if ax > 0. or ay > 0.:
- # catch cases where hypot(ax, ay) is subnormal
- real = math.log(math.hypot(math.ldexp(ax, DBL_MANT_DIG),
- math.ldexp(ay, DBL_MANT_DIG)))
- real -= DBL_MANT_DIG*M_LN2
- else:
- # log(+/-0. +/- 0i)
- raise ValueError("math domain error")
- #real = -INF
- #imag = atan2(y, x)
- else:
- h = math.hypot(ax, ay)
- if 0.71 <= h and h <= 1.73:
- am = max(ax, ay)
- an = min(ax, ay)
- real = log1p((am-1)*(am+1) + an*an) / 2.
- else:
- real = math.log(h)
- imag = math.atan2(y, x)
- return (real, imag)
-
+ return rcomplex.c_log(x,y)
_inner_wrapped_log = wrapped_log
@@ -300,196 +87,38 @@
@unaryfn
def c_log10(x, y):
- rx, ry = c_log(x, y)
- return (rx / M_LN10, ry / M_LN10)
-
+ return rcomplex.c_log10(x,y)
@unaryfn
def c_exp(x, y):
- if not isfinite(x) or not isfinite(y):
- if isinf(x) and isfinite(y) and y != 0.:
- if x > 0:
- real = copysign(INF, math.cos(y))
- imag = copysign(INF, math.sin(y))
- else:
- real = copysign(0., math.cos(y))
- imag = copysign(0., math.sin(y))
- r = (real, imag)
- else:
- r = exp_special_values[special_type(x)][special_type(y)]
-
- # need to raise ValueError if y is +/- infinity and x is not
- # a NaN and not -infinity
- if isinf(y) and (isfinite(x) or (isinf(x) and x > 0)):
- raise ValueError("math domain error")
- return r
-
- if x > CM_LOG_LARGE_DOUBLE:
- l = math.exp(x-1.)
- real = l * math.cos(y) * math.e
- imag = l * math.sin(y) * math.e
- else:
- l = math.exp(x)
- real = l * math.cos(y)
- imag = l * math.sin(y)
- if isinf(real) or isinf(imag):
- raise OverflowError("math range error")
- return real, imag
-
+ return rcomplex.c_exp(x,y)
@unaryfn
def c_cosh(x, y):
- if not isfinite(x) or not isfinite(y):
- if isinf(x) and isfinite(y) and y != 0.:
- if x > 0:
- real = copysign(INF, math.cos(y))
- imag = copysign(INF, math.sin(y))
- else:
- real = copysign(INF, math.cos(y))
- imag = -copysign(INF, math.sin(y))
- r = (real, imag)
- else:
- r = cosh_special_values[special_type(x)][special_type(y)]
-
- # need to raise ValueError if y is +/- infinity and x is not
- # a NaN
- if isinf(y) and not isnan(x):
- raise ValueError("math domain error")
- return r
-
- if fabs(x) > CM_LOG_LARGE_DOUBLE:
- # deal correctly with cases where cosh(x) overflows but
- # cosh(z) does not.
- x_minus_one = x - copysign(1., x)
- real = math.cos(y) * math.cosh(x_minus_one) * math.e
- imag = math.sin(y) * math.sinh(x_minus_one) * math.e
- else:
- real = math.cos(y) * math.cosh(x)
- imag = math.sin(y) * math.sinh(x)
- if isinf(real) or isinf(imag):
- raise OverflowError("math range error")
- return real, imag
-
+ return rcomplex.c_cosh(x,y)
@unaryfn
def c_sinh(x, y):
- # special treatment for sinh(+/-inf + iy) if y is finite and nonzero
- if not isfinite(x) or not isfinite(y):
- if isinf(x) and isfinite(y) and y != 0.:
- if x > 0:
- real = copysign(INF, math.cos(y))
- imag = copysign(INF, math.sin(y))
- else:
- real = -copysign(INF, math.cos(y))
- imag = copysign(INF, math.sin(y))
- r = (real, imag)
- else:
- r = sinh_special_values[special_type(x)][special_type(y)]
-
- # need to raise ValueError if y is +/- infinity and x is not
- # a NaN
- if isinf(y) and not isnan(x):
- raise ValueError("math domain error")
- return r
-
- if fabs(x) > CM_LOG_LARGE_DOUBLE:
- x_minus_one = x - copysign(1., x)
- real = math.cos(y) * math.sinh(x_minus_one) * math.e
- imag = math.sin(y) * math.cosh(x_minus_one) * math.e
- else:
- real = math.cos(y) * math.sinh(x)
- imag = math.sin(y) * math.cosh(x)
- if isinf(real) or isinf(imag):
- raise OverflowError("math range error")
- return real, imag
-
+ return rcomplex.c_sinh(x,y)
@unaryfn
def c_tanh(x, y):
- # Formula:
- #
- # tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /
- # (1+tan(y)^2 tanh(x)^2)
- #
- # To avoid excessive roundoff error, 1-tanh(x)^2 is better computed
- # as 1/cosh(x)^2. When abs(x) is large, we approximate 1-tanh(x)^2
- # by 4 exp(-2*x) instead, to avoid possible overflow in the
- # computation of cosh(x).
-
- if not isfinite(x) or not isfinite(y):
- if isinf(x) and isfinite(y) and y != 0.:
- if x > 0:
- real = 1.0 # vv XXX why is the 2. there?
- imag = copysign(0., 2. * math.sin(y) * math.cos(y))
- else:
- real = -1.0
- imag = copysign(0., 2. * math.sin(y) * math.cos(y))
- r = (real, imag)
- else:
- r = tanh_special_values[special_type(x)][special_type(y)]
-
- # need to raise ValueError if y is +/-infinity and x is finite
- if isinf(y) and isfinite(x):
- raise ValueError("math domain error")
- return r
-
- if fabs(x) > CM_LOG_LARGE_DOUBLE:
- real = copysign(1., x)
- imag = 4. * math.sin(y) * math.cos(y) * math.exp(-2.*fabs(x))
- else:
- tx = math.tanh(x)
- ty = math.tan(y)
- cx = 1. / math.cosh(x)
- txty = tx * ty
- denom = 1. + txty * txty
- real = tx * (1. + ty*ty) / denom
- imag = ((ty / denom) * cx) * cx
- return real, imag
-
+ return rcomplex.c_tanh(x,y)
@unaryfn
def c_cos(x, y):
- # cos(z) = cosh(iz)
- return c_cosh(-y, x)
+ return rcomplex.c_cos(x,y)
@unaryfn
def c_sin(x, y):
- # sin(z) = -i sinh(iz)
- sx, sy = c_sinh(-y, x)
- return sy, -sx
+ return rcomplex.c_sin(x,y)
@unaryfn
def c_tan(x, y):
- # tan(z) = -i tanh(iz)
- sx, sy = c_tanh(-y, x)
- return sy, -sx
-
+ return rcomplex.c_tan(x,y)
def c_rect(r, phi):
- if not isfinite(r) or not isfinite(phi):
- # if r is +/-infinity and phi is finite but nonzero then
- # result is (+-INF +-INF i), but we need to compute cos(phi)
- # and sin(phi) to figure out the signs.
- if isinf(r) and isfinite(phi) and phi != 0.:
- if r > 0:
- real = copysign(INF, math.cos(phi))
- imag = copysign(INF, math.sin(phi))
- else:
- real = -copysign(INF, math.cos(phi))
- imag = -copysign(INF, math.sin(phi))
- z = (real, imag)
- else:
- z = rect_special_values[special_type(r)][special_type(phi)]
-
- # need to raise ValueError if r is a nonzero number and phi
- # is infinite
- if r != 0. and not isnan(r) and isinf(phi):
- raise ValueError("math domain error")
- return z
-
- real = r * math.cos(phi)
- imag = r * math.sin(phi)
- return real, imag
+ return rcomplex.c_rect(r,phi)
def wrapped_rect(space, w_x, w_y):
x = space.float_w(w_x)
@@ -500,28 +129,7 @@
def c_phase(x, y):
- # Windows screws up atan2 for inf and nan, and alpha Tru64 5.1 doesn't
- # follow C99 for atan2(0., 0.).
- if isnan(x) or isnan(y):
- return NAN
- if isinf(y):
- if isinf(x):
- if copysign(1., x) == 1.:
- # atan2(+-inf, +inf) == +-pi/4
- return copysign(0.25 * math.pi, y)
- else:
- # atan2(+-inf, -inf) == +-pi*3/4
- return copysign(0.75 * math.pi, y)
- # atan2(+-inf, x) == +-pi/2 for finite x
- return copysign(0.5 * math.pi, y)
- if isinf(x) or y == 0.:
- if copysign(1., x) == 1.:
- # atan2(+-y, +inf) = atan2(+-0, +x) = +-0.
- return copysign(0., y)
- else:
- # atan2(+-y, -inf) = atan2(+-0., -x) = +-pi.
- return copysign(math.pi, y)
- return math.atan2(y, x)
+ return rcomplex.c_phase(x,y)
def wrapped_phase(space, w_z):
x, y = space.unpackcomplex(w_z)
@@ -531,28 +139,10 @@
def c_abs(x, y):
- if not isfinite(x) or not isfinite(y):
- # C99 rules: if either the real or the imaginary part is an
- # infinity, return infinity, even if the other part is a NaN.
- if isinf(x):
- return INF
- if isinf(y):
- return INF
-
- # either the real or imaginary part is a NaN,
- # and neither is infinite. Result should be NaN.
- return NAN
-
- result = math.hypot(x, y)
- if not isfinite(result):
- raise OverflowError("math range error")
- return result
-
+ return rcomplex.c_abs(x,y)
def c_polar(x, y):
- phi = c_phase(x, y)
- r = c_abs(x, y)
- return r, phi
+ return rcomplex.c_polar(x,y)
def wrapped_polar(space, w_z):
x, y = space.unpackcomplex(w_z)
@@ -562,7 +152,7 @@
def c_isinf(x, y):
- return isinf(x) or isinf(y)
+ return rcomplex.c_isinf(x,y)
def wrapped_isinf(space, w_z):
x, y = space.unpackcomplex(w_z)
@@ -572,7 +162,7 @@
def c_isnan(x, y):
- return isnan(x) or isnan(y)
+ return rcomplex.c_isnan(x,y)
def wrapped_isnan(space, w_z):
x, y = space.unpackcomplex(w_z)
diff --git a/pypy/module/cmath/test/test_cmath.py b/pypy/module/cmath/test/test_cmath.py
--- a/pypy/module/cmath/test/test_cmath.py
+++ b/pypy/module/cmath/test/test_cmath.py
@@ -115,7 +115,6 @@
Empty lines or lines starting with -- are ignored
yields id, fn, arg_real, arg_imag, exp_real, exp_imag
"""
- fname = os.path.join(os.path.dirname(__file__), fname)
with open(fname) as fp:
for line in fp:
# skip comment lines and blank lines
@@ -186,8 +185,10 @@
def test_specific_values():
#if not float.__getformat__("double").startswith("IEEE"):
# return
-
- for id, fn, ar, ai, er, ei, flags in parse_testfile('cmath_testcases.txt'):
+
+ # too fragile...
+ fname = os.path.join(os.path.dirname(__file__), '../../../rlib/test', 'rcomplex_testcases.txt')
+ for id, fn, ar, ai, er, ei, flags in parse_testfile(fname):
arg = (ar, ai)
expected = (er, ei)
function = getattr(interp_cmath, 'c_' + fn)
More information about the pypy-commit
mailing list