Author: Jasper Schulz <[email protected]>
Branch: numpypy-complex2
Changeset: r55862:4a0dbf4567f8
Date: 2012-06-27 16:29 +0200
http://bitbucket.org/pypy/pypy/changeset/4a0dbf4567f8/
Log: started to outsource complex methods
diff --git a/pypy/rlib/rcomplex.py b/pypy/rlib/rcomplex.py
new file mode 100644
--- /dev/null
+++ b/pypy/rlib/rcomplex.py
@@ -0,0 +1,539 @@
+import math
+from math import copysign
+from pypy.module.cmath.special_value import isfinite
+
+#binary
+
+def c_add(x, y):
+ (r1, i1), (r2, i2) = x, y
+ r = r1 + r2
+ i = i1 + i2
+ return (r, i)
+
+def c_sub(x, y):
+ (r1, i1), (r2, i2) = x, y
+ r = r1 - r2
+ i = i1 - i2
+ return (r, i)
+
+def c_mul(x, y):
+ (r1, i1), (r2, i2) = x, y
+ r = r1 * r2 - i1 * i2
+ i = r1 * i2 + i1 * r2
+ return (r, i)
+
+def c_div(x, y): #x/y
+ (r1, i1), (r2, i2) = x, y
+ if r2 < 0:
+ abs_r2 = -r2
+ else:
+ abs_r2 = r2
+ if i2 < 0:
+ abs_i2 = -i2
+ else:
+ abs_i2 = i2
+ if abs_r2 >= abs_i2:
+ if abs_r2 == 0.0:
+ raise ZeroDivisionError
+ else:
+ ratio = i2 / r2
+ denom = r2 + i2 * ratio
+ rr = (r1 + i1 * ratio) / denom
+ ir = (i1 - r1 * ratio) / denom
+ else:
+ ratio = r2 / i2
+ denom = r2 * ratio + i2
+ assert i2 != 0.0
+ rr = (r1 * ratio + i1) / denom
+ ir = (i1 * ratio - r1) / denom
+ return (rr, ir)
+
+def c_pow(x, y):
+ (r1, i1), (r2, i2) = x, y
+ if r2 == 0.0 and i2 == 0.0:
+ rr, ir = 1, 0
+ elif r1 == 0.0 and i1 == 0.0:
+ if i2 != 0.0 or r2 < 0.0:
+ raise ZeroDivisionError
+ rr, ir = (0.0, 0.0)
+ else:
+ vabs = math.hypot(r1,i1)
+ len = math.pow(vabs,r2)
+ at = math.atan2(i1,r1)
+ phase = at * r2
+ if i2 != 0.0:
+ len /= math.exp(at * i2)
+ phase += i2 * math.log(vabs)
+ rr = len * math.cos(phase)
+ ir = len * math.sin(phase)
+ return (rr, ir)
+
+#unary
+
+def c_neg(r, i):
+ return (-r, -i)
+
+
+def c_sqrt(r, i):
+ '''
+ 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(r) or not isfinite(i):
+ return sqrt_special_values[special_type(r)][special_type(i)]
+
+ if r == 0. and i == 0.:
+ return (0., y)
+
+ ar = fabs(r)
+ ai = fabs(i)
+
+ if ar < DBL_MIN and ai < DBL_MIN and (ar > 0. or ai > 0.):
+ # here we catch cases where hypot(ar, ai) is subnormal
+ ar = math.ldexp(ar, CM_SCALE_UP)
+ ai1= math.ldexp(ai, CM_SCALE_UP)
+ s = math.ldexp(math.sqrt(ar + math.hypot(ar, ai1)),
+ CM_SCALE_DOWN)
+ else:
+ ar /= 8.
+ s = 2.*math.sqrt(ar + math.hypot(ar, ai/8.))
+
+ d = ai/(2.*s)
+
+ if x >= 0.:
+ return (s, copysign(d, i))
+ else:
+ return (d, copysign(s, i))
+
+
+def c_acos(r, i):
+ if not isfinite(r) or not isfinite(i):
+ return acos_special_values[special_type(r)][special_type(i)]
+
+ if fabs(r) > CM_LARGE_DOUBLE or fabs(i) > CM_LARGE_DOUBLE:
+ # avoid unnecessary overflow for large arguments
+ real = math.atan2(fabs(i), r)
+ # split into cases to make sure that the branch cut has the
+ # correct continuity on systems with unsigned zeros
+ if r < 0.:
+ imag = -copysign(math.log(math.hypot(r/2., i/2.)) +
+ M_LN2*2., i)
+ else:
+ imag = copysign(math.log(math.hypot(r/2., i/2.)) +
+ M_LN2*2., -i)
+ else:
+ s1r, s1i = c_sqrt(1.-r, -i)
+ s2r, s2i = c_sqrt(1.+r, i)
+ real = 2.*math.atan2(s1r, s2r)
+ imag = asinh(s2r*s1i - s2i*s1r)
+ return (real, imag)
+
+
+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)
+
+
+def c_asin(x, y):
+ # asin(z) = -i asinh(iz)
+ sx, sy = c_asinh(-y, x)
+ return (sy, -sx)
+
+
+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)
+
+
+def c_atan(x, y):
+ # atan(z) = -i atanh(iz)
+ sx, sy = c_atanh(-y, x)
+ return (sy, -sx)
+
+
+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)
+
+
+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)
+
+
+def c_log10(x, y):
+ rx, ry = c_log(x, y)
+ return (rx / M_LN10, ry / M_LN10)
+
+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
+
+
+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
+
+
+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
+
+
+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
+
+
+def c_cos(r, i):
+ # cos(z) = cosh(iz)
+ return c_cosh(-i, r)
+
+def c_sin(r, i):
+ # sin(z) = -i sinh(iz)
+ sr, si = c_sinh(-i, r)
+ return si, -sr
+
+def c_tan(r, i):
+ # tan(z) = -i tanh(iz)
+ sr, si = c_tanh(-i, r)
+ return si, -sr
+
+
+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
+
+
+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)
+
+
+def c_abs(r, i):
+ if not isfinite(r) or not isfinite(i):
+ # 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(r):
+ return INF
+ if isinf(i):
+ return INF
+
+ # either the real or imaginary part is a NaN,
+ # and neither is infinite. Result should be NaN.
+ return NAN
+
+ result = math.hypot(r, i)
+ if not isfinite(result):
+ raise OverflowError("math range error")
+ return result
+
+
+def c_polar(r, i):
+ real = c_abs(r, i)
+ phi = c_phase(r, i)
+ return real, phi
+
+
+def c_isinf(r, i):
+ return isinf(r) or isinf(i)
+
+
+def c_isnan(r, i):
+ return isnan(r) or isnan(i)
+
diff --git a/pypy/rlib/test/test_rcomplex.py b/pypy/rlib/test/test_rcomplex.py
new file mode 100644
--- /dev/null
+++ b/pypy/rlib/test/test_rcomplex.py
@@ -0,0 +1,31 @@
+
+import pypy.rlib.rcomplex as c
+
+
+def test_add():
+ for c1, c2, result in [
+ ((0, 0), (0, 0), (0, 0)),
+ ((1, 0), (2, 0), (3, 0)),
+ ((0, 3), (0, 2), (0, 5)),
+ ((10., -3.), (-5, 7), (5, 4)),
+ ]:
+ assert c.c_add(c1, c2) == result
+
+def test_sub():
+ for c1, c2, result in [
+ ((0, 0), (0, 0), (0, 0)),
+ ((1, 0), (2, 0), (-1, 0)),
+ ((0, 3), (0, 2), (0, 1)),
+ ((10, -3), (-5, 7), (15, -10)),
+ ((42, 0.3), (42, 0.3), (0, 0))
+ ]:
+ assert c.c_sub(c1, c2) == result
+
+def test_mul():
+ for c1, c2, result in [
+ ((0, 0), (0, 0), (0, 0)),
+ ((1, 0), (2, 0), (2, 0)),
+ ((0, 3), (0, 2), (-6, 0)),
+ ((0, -3), (-5, 0), (0, 15)),
+ ]:
+ assert c.c_mul(c1, c2) == result
\ No newline at end of file
_______________________________________________
pypy-commit mailing list
[email protected]
http://mail.python.org/mailman/listinfo/pypy-commit
