[pypy-commit] pypy numpypy-complex2: fix dos line endings
mattip
noreply at buildbot.pypy.org
Mon Oct 1 23:44:33 CEST 2012
Author: mattip <matti.picus at gmail.com>
Branch: numpypy-complex2
Changeset: r57710:e40bc49712a1
Date: 2012-10-01 23:22 +0200
http://bitbucket.org/pypy/pypy/changeset/e40bc49712a1/
Log: fix dos line endings
diff --git a/pypy/rlib/rcomplex.py b/pypy/rlib/rcomplex.py
--- a/pypy/rlib/rcomplex.py
+++ b/pypy/rlib/rcomplex.py
@@ -1,561 +1,561 @@
-import math
-from math import copysign, fabs, pi, e
-from pypy.rlib.rfloat import copysign, asinh, log1p, isinf, isnan
-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
-
-#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 i1 == 0 and i2 == 0 and r1>=0:
- rr = pow(r1, r2)
- ir = 0.
- elif r2 == 0.0 and i2 == 0.0:
- rr, ir = 1, 0
- elif r1 == 1.0 and i1 == 0.0:
- rr, ir = (1.0, 0.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(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))
-
-
-
-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)
-
-
-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)
-
+import math
+from math import copysign, fabs, pi, e
+from pypy.rlib.rfloat import copysign, asinh, log1p, isinf, isnan
+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
+
+#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 i1 == 0 and i2 == 0 and r1>=0:
+ rr = pow(r1, r2)
+ ir = 0.
+ elif r2 == 0.0 and i2 == 0.0:
+ rr, ir = 1, 0
+ elif r1 == 1.0 and i1 == 0.0:
+ rr, ir = (1.0, 0.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(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))
+
+
+
+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)
+
+
+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
--- a/pypy/rlib/test/test_rcomplex.py
+++ b/pypy/rlib/test/test_rcomplex.py
@@ -1,259 +1,259 @@
-from __future__ import with_statement
-
-import pypy.rlib.rcomplex as c
-from pypy.rlib.rfloat import copysign, isnan, isinf
-import os, sys, math, struct
-
-
-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
-
-def parse_testfile2(fname):
- """Parse a file with test values
-
- Empty lines or lines starting with -- are ignored
- yields id, fn, arg1_real, arg1_imag, arg2_real, arg2_imag,
- exp_real, exp_imag where numbers in file may be expressed as floating point or hex
- """
- fname = os.path.join(os.path.dirname(__file__), fname)
- with open(fname) as fp:
- for line in fp:
- # skip comment lines and blank lines
- if line.startswith('--') or not line.strip():
- continue
-
- lhs, rhs = line.split('->')
- lhs_pieces = lhs.split()
- rhs_pieces = rhs.split()
- for i in range(2, len(lhs_pieces)):
- if lhs_pieces[i].lower().startswith('0x'):
- lhs_pieces[i] = struct.unpack('d',
- struct.pack('q',int(lhs_pieces[i])))
- else:
- lhs_pieces[i] = float(lhs_pieces[i])
- for i in range(2):
- if rhs_pieces[i].lower().startswith('0x'):
- rhs_pieces[i] = struct.unpack('d',
- struct.pack('l',int(rhs_pieces[i])))
- else:
- rhs_pieces[i] = float(rhs_pieces[i])
- #id, fn, arg1_real, arg1_imag arg2_real, arg2_imag =
- #exp_real, exp_imag = rhs_pieces[0], rhs_pieces[1]
- flags = rhs_pieces[2:]
- id_f, fn = lhs_pieces[:2]
- if len(lhs_pieces)>4:
- args = (lhs_pieces[2:4], lhs_pieces[4:])
- else:
- args = lhs_pieces[2:]
- yield id_f, fn, args, rhs_pieces[:2], flags
-
-
-
-def parse_testfile(fname):
- """Parse a file with test values
-
- 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
- if line.startswith('--') or not line.strip():
- continue
-
- lhs, rhs = line.split('->')
- id, fn, arg_real, arg_imag = lhs.split()
- rhs_pieces = rhs.split()
- exp_real, exp_imag = rhs_pieces[0], rhs_pieces[1]
- flags = rhs_pieces[2:]
-
- yield (id, fn,
- float(arg_real), float(arg_imag),
- float(exp_real), float(exp_imag),
- flags
- )
-
-def args_to_str(args):
- if isinstance(args[0],(list, tuple)):
- return '(complex(%r, %r), complex(%r, %r))' % \
- (args[0][0], args[0][1], args[1][0], args[1][1])
- else:
- return '(complex(%r, %r))' % (args[0], args[1])
-
-def rAssertAlmostEqual(a, b, rel_err = 2e-15, abs_err = 5e-323, msg=''):
- """Fail if the two floating-point numbers are not almost equal.
-
- Determine whether floating-point values a and b are equal to within
- a (small) rounding error. The default values for rel_err and
- abs_err are chosen to be suitable for platforms where a float is
- represented by an IEEE 754 double. They allow an error of between
- 9 and 19 ulps.
- """
-
- # special values testing
- if isnan(a):
- if isnan(b):
- return
- raise AssertionError(msg + '%r should be nan' % (b,))
-
- if isinf(a):
- if a == b:
- return
- raise AssertionError(msg + 'finite result where infinity expected: '
- 'expected %r, got %r' % (a, b))
-
- # if both a and b are zero, check whether they have the same sign
- # (in theory there are examples where it would be legitimate for a
- # and b to have opposite signs; in practice these hardly ever
- # occur).
- if not a and not b:
- # only check it if we are running on top of CPython >= 2.6
- if sys.version_info >= (2, 6) and copysign(1., a) != copysign(1., b):
- raise AssertionError(msg + 'zero has wrong sign: expected %r, '
- 'got %r' % (a, b))
-
- # if a-b overflows, or b is infinite, return False. Again, in
- # theory there are examples where a is within a few ulps of the
- # max representable float, and then b could legitimately be
- # infinite. In practice these examples are rare.
- try:
- absolute_error = abs(b-a)
- except OverflowError:
- pass
- else:
- # test passes if either the absolute error or the relative
- # error is sufficiently small. The defaults amount to an
- # error of between 9 ulps and 19 ulps on an IEEE-754 compliant
- # machine.
- if absolute_error <= max(abs_err, rel_err * abs(a)):
- return
- raise AssertionError(msg + '%r and %r are not sufficiently close' % (a, b))
-
-def test_specific_values():
- #if not float.__getformat__("double").startswith("IEEE"):
- # return
-
- for id, fn, arg, expected, flags in parse_testfile2('rcomplex_testcases.txt'):
- function = getattr(c, 'c_' + fn)
- #
- if 'divide-by-zero' in flags or 'invalid' in flags:
- try:
- actual = function(*arg)
- except ValueError:
- continue
- else:
- raise AssertionError('ValueError not raised in test '
- '%s: %s%s' % (id, fn, args_to_str(arg)))
- if 'overflow' in flags:
- try:
- actual = function(*arg)
- except OverflowError:
- continue
- else:
- raise AssertionError('OverflowError not raised in test '
- '%s: %s%s' % (id, fn, args_to_str(arg)))
- actual = function(*arg)
-
- if 'ignore-real-sign' in flags:
- actual = (abs(actual[0]), actual[1])
- expected = (abs(expected[0]), expected[1])
- if 'ignore-imag-sign' in flags:
- actual = (actual[0], abs(actual[1]))
- expected = (expected[0], abs(expected[1]))
-
- # for the real part of the log function, we allow an
- # absolute error of up to 2e-15.
- if fn in ('log', 'log10'):
- real_abs_err = 2e-15
- else:
- real_abs_err = 5e-323
-
- error_message = (
- '%s: %s%s\n'
- 'Expected: complex(%r, %r)\n'
- 'Received: complex(%r, %r)\n'
- ) % (id, fn, args_to_str(arg),
- expected[0], expected[1],
- actual[0], actual[1])
-
- rAssertAlmostEqual(expected[0], actual[0],
- abs_err=real_abs_err,
- msg=error_message)
- rAssertAlmostEqual(expected[1], actual[1],
- abs_err=real_abs_err,
- msg=error_message)
-
- for id, fn, a, expected, flags in parse_testfile2('rcomplex_testcases2.txt'):
- function = getattr(c, 'c_' + fn)
- #
- if 'divide-by-zero' in flags or 'invalid' in flags:
- try:
- actual = function(*a)
- except ValueError:
- continue
- else:
- raise AssertionError('ValueError not raised in test '
- '%s: %s%s' % (id, fn, args_to_str(a)))
- if 'overflow' in flags:
- try:
- actual = function(*a)
- except OverflowError:
- continue
- else:
- raise AssertionError('OverflowError not raised in test '
- '%s: %s%s' % (id, fn, args_to_str(a)))
- actual = function(*a)
-
- if 'ignore-real-sign' in flags:
- actual = (abs(actual[0]), actual[1])
- expected = (abs(expected[0]), expected[1])
- if 'ignore-imag-sign' in flags:
- actual = (actual[0], abs(actual[1]))
- expected = (expected[0], abs(expected[1]))
-
- # for the real part of the log function, we allow an
- # absolute error of up to 2e-15.
- if fn in ('log', 'log10'):
- real_abs_err = 2e-15
- else:
- real_abs_err = 5e-323
- error_message = (
- '%s: %s%s\n'
- 'Expected: complex(%r, %r)\n'
- 'Received: complex(%r, %r)\n'
- ) % (id, fn, args_to_str(a),
- expected[0], expected[1],
- actual[0], actual[1])
-
- rAssertAlmostEqual(expected[0], actual[0],
- abs_err=real_abs_err,
- msg=error_message)
- rAssertAlmostEqual(expected[1], actual[1],
- abs_err=real_abs_err,
- msg=error_message)
+from __future__ import with_statement
+
+import pypy.rlib.rcomplex as c
+from pypy.rlib.rfloat import copysign, isnan, isinf
+import os, sys, math, struct
+
+
+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
+
+def parse_testfile2(fname):
+ """Parse a file with test values
+
+ Empty lines or lines starting with -- are ignored
+ yields id, fn, arg1_real, arg1_imag, arg2_real, arg2_imag,
+ exp_real, exp_imag where numbers in file may be expressed as floating point or hex
+ """
+ fname = os.path.join(os.path.dirname(__file__), fname)
+ with open(fname) as fp:
+ for line in fp:
+ # skip comment lines and blank lines
+ if line.startswith('--') or not line.strip():
+ continue
+
+ lhs, rhs = line.split('->')
+ lhs_pieces = lhs.split()
+ rhs_pieces = rhs.split()
+ for i in range(2, len(lhs_pieces)):
+ if lhs_pieces[i].lower().startswith('0x'):
+ lhs_pieces[i] = struct.unpack('d',
+ struct.pack('q',int(lhs_pieces[i])))
+ else:
+ lhs_pieces[i] = float(lhs_pieces[i])
+ for i in range(2):
+ if rhs_pieces[i].lower().startswith('0x'):
+ rhs_pieces[i] = struct.unpack('d',
+ struct.pack('l',int(rhs_pieces[i])))
+ else:
+ rhs_pieces[i] = float(rhs_pieces[i])
+ #id, fn, arg1_real, arg1_imag arg2_real, arg2_imag =
+ #exp_real, exp_imag = rhs_pieces[0], rhs_pieces[1]
+ flags = rhs_pieces[2:]
+ id_f, fn = lhs_pieces[:2]
+ if len(lhs_pieces)>4:
+ args = (lhs_pieces[2:4], lhs_pieces[4:])
+ else:
+ args = lhs_pieces[2:]
+ yield id_f, fn, args, rhs_pieces[:2], flags
+
+
+
+def parse_testfile(fname):
+ """Parse a file with test values
+
+ 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
+ if line.startswith('--') or not line.strip():
+ continue
+
+ lhs, rhs = line.split('->')
+ id, fn, arg_real, arg_imag = lhs.split()
+ rhs_pieces = rhs.split()
+ exp_real, exp_imag = rhs_pieces[0], rhs_pieces[1]
+ flags = rhs_pieces[2:]
+
+ yield (id, fn,
+ float(arg_real), float(arg_imag),
+ float(exp_real), float(exp_imag),
+ flags
+ )
+
+def args_to_str(args):
+ if isinstance(args[0],(list, tuple)):
+ return '(complex(%r, %r), complex(%r, %r))' % \
+ (args[0][0], args[0][1], args[1][0], args[1][1])
+ else:
+ return '(complex(%r, %r))' % (args[0], args[1])
+
+def rAssertAlmostEqual(a, b, rel_err = 2e-15, abs_err = 5e-323, msg=''):
+ """Fail if the two floating-point numbers are not almost equal.
+
+ Determine whether floating-point values a and b are equal to within
+ a (small) rounding error. The default values for rel_err and
+ abs_err are chosen to be suitable for platforms where a float is
+ represented by an IEEE 754 double. They allow an error of between
+ 9 and 19 ulps.
+ """
+
+ # special values testing
+ if isnan(a):
+ if isnan(b):
+ return
+ raise AssertionError(msg + '%r should be nan' % (b,))
+
+ if isinf(a):
+ if a == b:
+ return
+ raise AssertionError(msg + 'finite result where infinity expected: '
+ 'expected %r, got %r' % (a, b))
+
+ # if both a and b are zero, check whether they have the same sign
+ # (in theory there are examples where it would be legitimate for a
+ # and b to have opposite signs; in practice these hardly ever
+ # occur).
+ if not a and not b:
+ # only check it if we are running on top of CPython >= 2.6
+ if sys.version_info >= (2, 6) and copysign(1., a) != copysign(1., b):
+ raise AssertionError(msg + 'zero has wrong sign: expected %r, '
+ 'got %r' % (a, b))
+
+ # if a-b overflows, or b is infinite, return False. Again, in
+ # theory there are examples where a is within a few ulps of the
+ # max representable float, and then b could legitimately be
+ # infinite. In practice these examples are rare.
+ try:
+ absolute_error = abs(b-a)
+ except OverflowError:
+ pass
+ else:
+ # test passes if either the absolute error or the relative
+ # error is sufficiently small. The defaults amount to an
+ # error of between 9 ulps and 19 ulps on an IEEE-754 compliant
+ # machine.
+ if absolute_error <= max(abs_err, rel_err * abs(a)):
+ return
+ raise AssertionError(msg + '%r and %r are not sufficiently close' % (a, b))
+
+def test_specific_values():
+ #if not float.__getformat__("double").startswith("IEEE"):
+ # return
+
+ for id, fn, arg, expected, flags in parse_testfile2('rcomplex_testcases.txt'):
+ function = getattr(c, 'c_' + fn)
+ #
+ if 'divide-by-zero' in flags or 'invalid' in flags:
+ try:
+ actual = function(*arg)
+ except ValueError:
+ continue
+ else:
+ raise AssertionError('ValueError not raised in test '
+ '%s: %s%s' % (id, fn, args_to_str(arg)))
+ if 'overflow' in flags:
+ try:
+ actual = function(*arg)
+ except OverflowError:
+ continue
+ else:
+ raise AssertionError('OverflowError not raised in test '
+ '%s: %s%s' % (id, fn, args_to_str(arg)))
+ actual = function(*arg)
+
+ if 'ignore-real-sign' in flags:
+ actual = (abs(actual[0]), actual[1])
+ expected = (abs(expected[0]), expected[1])
+ if 'ignore-imag-sign' in flags:
+ actual = (actual[0], abs(actual[1]))
+ expected = (expected[0], abs(expected[1]))
+
+ # for the real part of the log function, we allow an
+ # absolute error of up to 2e-15.
+ if fn in ('log', 'log10'):
+ real_abs_err = 2e-15
+ else:
+ real_abs_err = 5e-323
+
+ error_message = (
+ '%s: %s%s\n'
+ 'Expected: complex(%r, %r)\n'
+ 'Received: complex(%r, %r)\n'
+ ) % (id, fn, args_to_str(arg),
+ expected[0], expected[1],
+ actual[0], actual[1])
+
+ rAssertAlmostEqual(expected[0], actual[0],
+ abs_err=real_abs_err,
+ msg=error_message)
+ rAssertAlmostEqual(expected[1], actual[1],
+ abs_err=real_abs_err,
+ msg=error_message)
+
+ for id, fn, a, expected, flags in parse_testfile2('rcomplex_testcases2.txt'):
+ function = getattr(c, 'c_' + fn)
+ #
+ if 'divide-by-zero' in flags or 'invalid' in flags:
+ try:
+ actual = function(*a)
+ except ValueError:
+ continue
+ else:
+ raise AssertionError('ValueError not raised in test '
+ '%s: %s%s' % (id, fn, args_to_str(a)))
+ if 'overflow' in flags:
+ try:
+ actual = function(*a)
+ except OverflowError:
+ continue
+ else:
+ raise AssertionError('OverflowError not raised in test '
+ '%s: %s%s' % (id, fn, args_to_str(a)))
+ actual = function(*a)
+
+ if 'ignore-real-sign' in flags:
+ actual = (abs(actual[0]), actual[1])
+ expected = (abs(expected[0]), expected[1])
+ if 'ignore-imag-sign' in flags:
+ actual = (actual[0], abs(actual[1]))
+ expected = (expected[0], abs(expected[1]))
+
+ # for the real part of the log function, we allow an
+ # absolute error of up to 2e-15.
+ if fn in ('log', 'log10'):
+ real_abs_err = 2e-15
+ else:
+ real_abs_err = 5e-323
+ error_message = (
+ '%s: %s%s\n'
+ 'Expected: complex(%r, %r)\n'
+ 'Received: complex(%r, %r)\n'
+ ) % (id, fn, args_to_str(a),
+ expected[0], expected[1],
+ actual[0], actual[1])
+
+ rAssertAlmostEqual(expected[0], actual[0],
+ abs_err=real_abs_err,
+ msg=error_message)
+ rAssertAlmostEqual(expected[1], actual[1],
+ abs_err=real_abs_err,
+ msg=error_message)
More information about the pypy-commit
mailing list