"""This module provides exponential, logarithmic, trigonometric, inverse trigonometric, hyperbolic and inverse hyperbolic functions for complex numbers. It is intended as a drop-in replacement for the cmath module in the standard library. For the most part, it uses formulas and methods described by W. Kahan in his `Branch cuts for elementary functions' paper. Design goals ------------ - make all functions numerically sound; both the real part and the imaginary part should be within a few ulps of the true result, where this is reasonable. (But see the note on accuracy below.) - avoid unnecessary overflows in intermediate expressions, when the final result is representable. - fix buggy branch cuts in asinh - do the 'right thing' with respect to signed zeros on platforms that follow IEEE754 and C99/IEC 60559. In particular, all branch cuts should be continuous from both sides in this case. - don't do anything unreasonable on platforms that don't support signed zeros, or have signed zero support that doesn't comply fully with the above standards. With no signed zeros, continuity at branch cuts should match the documentation. Behaviour here is untested. Non-design goals ---------------- - do the right thing with NaNs and infinities. It's hard to do this portably. I believe that many of the routines below do actually do the right thing in the presence of NaNs, but this is mostly accidental. - (related to the above): give sensible and consistent exceptions. Again, it seems difficult to do this across platforms. - produce results that are provably accurate to within 1ulp. Note on accuracy ---------------- In an ideal world, the complex-valued function f(z) would return the closest representable complex number to the true mathematical value of f(z): that is, both the real and imaginary part of the result would be accurate to within <= 0.5ulp. Achieving this level of accuracy is very hard---most C math libraries don't manage it. (But see the crlibm and MPFR libraries.) A slightly more realistic goal is <1ulp. In practice, one might hope that the returned real and imaginary parts are always within a few ulps (say 10 or 20) of those of the closest representable value. But even this is an unrealistic goal in some situtations. For example let z be the complex number: 0.79999999999999993338661852249060757458209991455078125 + 0.600000000000000088817841970012523233890533447265625j which is exactly representable, assuming IEEE doubles are being used for the real and imaginary parts. The nearest representable complex number to the natural logarithm of z is: 6.16297582203915472977912941627176741932192527428924222476780414581298828125e-33 + 2.4980915447965088560522417537868022918701171875j It would take a lot of effort to get the real part anywhere near correct here. Other problems occur in computing trigonometric functions for complex numbers z with large real part, or hyperbolic trig functions for z with large imaginary part. Notes on signed zeros --------------------- There are many subtle difficulties here: for example, the expression "1 + z" should add 1 to the real part of z and leave the imaginary part untouched. But in Python, if the imaginary part of z is -0. then the imaginary part of 1+z will be +0: 1 gets coerced to the complex number 1 + 0j, and then the imaginary parts get added to give -0j + 0j = 0j. But z - 1 always does the right thing: subtracting +0. won't change the sign of zero. Similarly, z + (-1.) is fine. So we can either work with the underlying reals directly, or rewrite the erroneous 1+z as -(-z-1), which works. Similarly, 1-z should be rewritten as -(z-1). An alternative fix is to use the complex number 1 - 0j (note negative sign) in place of 1 in the expressions 1+z and 1-z. Similarly, the expression i*z is special-cased so that i*(x+i*y) = -y + i*x; see the function mul_by_j. The code below should `do the right thing' regardless of whether signed zeros are present. In particular: - on a platform (hardware + C math library) that supports signed zeros, so that for example: atan2(-0., positive) gives -0. atan2(0., positive) gives 0. atan2(-0., negative) gives -pi atan2(0., negative) gives pi sin(-0.) gives -0. 0. * -0. gives -0. and so on, the functions below should all be continuous on *both sides* of every branch cut, and `do the right thing' with respect to signed zeros. (For example sqrt(4 - 0j) = 2 - 0j.) - if atan2(-0., negative) gives a positive result then continuity at branch cuts is as specified in the docstrings, and there's no attempt to do the right thing with -0., if it even makes sense. the above two possibilities should cover the vast majority of cases. But if neither of the above hold---for example atan2 returns values in [-pi, pi), or does something else nonstandard, then behaviour on branch cuts may be unpredictable. Bugs ---- - atanh (and atan) has not been fully checked, and may be buggy - many routines don't accept floats. e.g. >>> from cmath_py import tan >>> tan(0.) Traceback (most recent call last): File "", line 1, in File "cmath_py.py", line 394, in tan return -mul_by_j(tanh(mul_by_j(z))) File "cmath_py.py", line 261, in mul_by_j return complex_with_signed_zero(-z.imag, z.real) AttributeError: 'float' object has no attribute 'imag' """ from __future__ import division from math import exp as fexp, log as flog, sqrt as fsqrt, atan as fatan from math import sin as fsin, cos as fcos, tan as ftan from math import sinh as fsinh, cosh as fcosh, tanh as ftanh from math import atan2, hypot, pi, e __all__ = ['acos', 'acosh', 'asin', 'asinh', 'atan', 'atanh', 'cos', 'cosh', 'e', 'exp', 'log', 'log10', 'pi', 'sin', 'sinh', 'sqrt', 'tan', 'tanh'] # the constant half_pi should be less than the true mathematical value # of pi/2; in other words, it should be the biggest floating-point # value in the interval (0, pi/2). Otherwise we'll get strange # behaviours on boundaries when tan(half_pi) evaluates to something # negative. half_pi = pi/2 assert ftan(half_pi) > 0. def complexify(z): """Interpret the Python object z as a complex number.""" # we look for, in order, # (1) existence of real and imag attributes: in this case z is # probably already complex, or usable as such # (2) existence of a __complex__ method # (3) existence of a __float__ method. # Note that we don't want to just do # z = complex(z) # since this would allow things like sqrt("2") to work. if not (hasattr(z, "real") and hasattr(z, "imag")): if hasattr(z, "__complex__"): z = z.__complex__() elif hasattr(z, "__float__"): z = complex(z.__float__(), 0.) else: raise TypeError("A complex number is expected") return z # Here are three floating-point functions; (these are present in the # C99 standard, and so should be present in most but not all C math # libraries; in would be better to use the library version when it # exists). def copysign(x, y): """float with the magnitude of x but the sign bit of y""" # evil atan2 hack to determine sign of y if y == 0. if y > 0. or y == 0. and atan2(y, -1.) > 0: return abs(x) else: return -abs(x) def fasinh(x): """Inverse hyperbolic sine of the float x.""" # this function isn't part of the standard C89 library, so we have # to roll our own. The usual formula is: # # asinh(x) = copysign(log(abs(x) + sqrt(1+x*x)), x) # # For abs(x) >= 1. we rewrite as # # asinh(x) = log(abs(x)) + log(1 + sqrt(1 + 1/(x*x))) # # which avoids overflow errors for large x. For abs(x) < 1., the # argument of the log in the usual formula is close to 1, giving a # loss of precision. In this case we rewrite as: # # asinh(x) = log1p(abs(x) + x*x/(sqrt(1+x*x))) ax = abs(x) if ax < 1.: result = flog1p(ax*(1.+ax/(hypot(1., ax)+1.))) else: result = flog(ax)+flog(1.+hypot(1., 1./ax)) return copysign(result, x) def flog1p(x): """Compute log(1+x) accurately for a floating-point argument x, even when x is close to 0.""" # If x <= -0.5 or x >= 1 then either x is huge or 1+x will be # exactly representable. Either way, we can just compute log(1+x) # directly. # For -0.5 <= x <= 1, we do the following: (1) let y = 1 + x, # which will be a floating-point approximation to the true # mathematical value of 1+x. (2) compute the error in this # approximation: e = (y-1.)-x. This should give an exact # result---that is, the true mathematical value of (y-1.)-x will # be exactly representable as a float. Now 1+x = e+y, so log1p(x) # = log(e+y) = log(y(1+e/y)) = log(y) + log(1+e/y). Since e/y # will be tiny, in general smaller than the machine epsilon, # log(1+e/y) can be well approximated by e/y. So we get the # approximation log(1+x) ~ log(y) + e/y. # Warning: this approach can be catastrophically bad in the # absence of a sensible hardware rounding scheme: suppose that we # try to compute log(1+x) for some negative tiny x (-2**-150, # say), and that the value y=1+x, instead of being rounded to the # nearest representable float (1.), gets rounded *down* to the # largest representable float that's strictly less than 1. Then # flog(y) and (y-1.-x)/y will both have size around the machine # epsilon, and the result of the subtraction will be either 0. or # have size around the square of the machine epsilon; the # information from the original x will have disappeared entirely. # A solution in this case is to explicitly test for small x # (smaller than the machine epsilon), and simply return x directly # in this case. But this requires knowing, or computing, the # machine epsilon, which is going to be difficult to do reliably # and portably. if x == 0.: return x if x <= -0.5 or x >= 1.: return flog(1.+x) else: y = 1.+x return flog(y)-(y-1.-x)/y # Python's complex constructor doesn't preserve the signs of zeros. For example: # # >>> complex(-0., -0.).real # okay # -0.0 # >>> complex(-0., -0.).imag # negative sign lost # 0.0 # # Here's a replacement for the complex constructor which does respect the sign. def complex_with_signed_zero(x, y): """Workaround for the fact that complex(x, y) loses the sign on y if y is zero.""" if y == 0. and atan2(y, -1) < 1.: return -complex(-x, -y) else: return complex(x, y) def mul_by_j(z): """Return the result of multiplying the complex number z by the square root 1j of -1. The only difference between mul_by_j(z) and 1j*z is in dealing with signed zeros. """ return complex_with_signed_zero(-z.imag, z.real) def log(z, base = None): """Natural logarithm of the complex argument z. There is a single branch cut, running along the real axis from 0 to -infinity, and the function is continuous from above on this branch cut. The log is undefined at 0. The optional argument base gives the base of the logarithm, which can be any complex number not equal to 0 or 1. It defaults to e. """ z = complexify(z) if z == 0.: raise ValueError("Math domain error: logarithm of 0 is undefined") # It's not enough to return the real part as log(x*x+y*y)/2, for # two reasons: first, the x*x+y*y might overflow or underflow, and # second, this loses valuable bits of y when for example x is # close to 1 and y is small, or vice versa, so that the argument # of log is close to 1. Similar objections apply to log(abs(z)) # or log(hypot(x, y)), though the danger of overflow here is far # smaller. Instead, we write sqrt(x*x+y*y) as # |x|*sqrt(1+(y/x)*(y/x)) if |y| < |x| and vice versa if |x| <= # |y|. Then the log can be expressed as log(|x|) + # log1p((y/x)*(y/x))/2, or log(|y|) + log1p((x/y)*(x/y))/2. x, y = z.real, z.imag ax = abs(x) ay = abs(y) am = max(ax, ay) anoam = min(ax, ay)/am real_part = flog(am) + flog1p(anoam*anoam)/2. result = complex_with_signed_zero(real_part, atan2(y, x)) if base is None: return result else: return result/log(base) def log10(z): """Logarithm to base 10 of the complex argument z. The branch cut is exactly the same as for log. """ return log(z)/flog(10) def pow(z, w): """complex number z raised to the complex power w. For any fixed w, the branch cut for z runs along the real axis from 0 to -infinity. For a fixed z, the function is continuous in w over the whole complex plane.""" return exp(w*log(z)) def sqrt(z): """Square root of the complex argument z. There is a single branch cut, running along the real axis from 0 to -infinity, and the function is continuous from above on this branch cut. """ # Formula: if u+iv = sqrt(x+iy) # then # # u = sqrt((x + sqrt(x^2 + y^2))/2), v = (y/2)/u. # # there's a possibility of cancellation error in x + sqrt(x^2 + # y^2) if x is negative, so for x < 0 use the alternative formula # # v = sqrt((-x + sqrt(x^2 + y^2))/2), u = (y/2)/v # # In most cases hypot(x, y) can be used to compute sqrt(x^2 + y^2) # without overflow or underflow. Even so, we still have problems # at the extremes of the floating point range: for x and y very # large either hypot(x, y) or the sum abs(x) + hypot(x, y) can # overflow. At the other end of the scale, if both x and y are # denormalized then hypot(x, y) will also be denormalized, so will # lack full precision. Both problems can be resolved by scaling # by the inverse of max(abs(x), abs(y)). z = complexify(z) x, y = z.real, z.imag if x == y == 0.: # preserve sign of y in result return complex_with_signed_zero(0., y) ax, ay = abs(x), abs(y) if ax >= ay: s = fsqrt(ax) * fsqrt((1. + hypot(1., ay/ax))/2.) else: xony = ax/ay s = fsqrt(ay) * fsqrt((xony + hypot(1., xony))/2.) d = ay/(2.*s) if x >= 0: return complex_with_signed_zero(s, copysign(d, y)) else: return complex_with_signed_zero(d, copysign(s, y)) def exp(z): """Exponential function for the complex argument z.""" # signs: check that if y == 0. then result has imaginary part 0. # while if y == -0. then result has imaginary part -0. # This preserves the relation exp(z*) = (exp(z))* # there's a very small possibility that if x is large then exp(x) # could overflow even though exp(z) is representable. In that case, # the simplest thing to do seems to be to compute exp(x-1) in place of # x; this won't incur significant extra error, since x-1 will be # exactly representable. z = complexify(z) x, y = z.real, z.imag if x > 700: exp_xm1 = fexp(x-1) return complex_with_signed_zero(exp_xm1*fcos(y)*e, exp_xm1*fsin(y)*e) else: exp_x = fexp(x) return complex_with_signed_zero(exp_x*fcos(y), exp_x*fsin(y)) def sin(z): """Sine of the complex argument z.""" return -mul_by_j(sinh(mul_by_j(z))) def cos(z): """Cosine of the complex argument z.""" return cosh(mul_by_j(z)) def tan(z): """Tangent of the complex argument z.""" return -mul_by_j(tanh(mul_by_j(z))) def cosh(z): """Hyperbolic cosine of the complex argument z.""" z = complexify(z) x, y = z.real, z.imag return complex_with_signed_zero(fcosh(x)*fcos(y), fsinh(x)*fsin(y)) def sinh(z): """Hyperbolic sine of the complex argument z.""" z = complexify(z) x, y = z.real, z.imag return complex_with_signed_zero(fsinh(x)*fcos(y), fcosh(x)*fsin(y)) def tanh(z): """Hyperbolic tangent of the complex argument z.""" # If abs(x) is large then computation of cosh(x) can cause # overflow; so in this case we use an approximation. The cutoff # value of 25 is chosen so that exp(-2*25) is less than the # machine epsilon (and hence exp(x) = cosh(x) and tanh(x) = +/- 1 # to within machine precision). z = complexify(z) x, y = z.real, z.imag if abs(x) > 25: return complex_with_signed_zero(copysign(1., x), fsin(2*y)*fexp(-2*abs(x))) else: t = ftan(y) b = 1+t*t s = fsinh(x) s2 = s*s r = fsqrt(1+s2) d = 1+b*s2 return complex_with_signed_zero(b*r*s/d, t/d) def asin(z): """Inverse sine of the complex argument z, as defined by the formula asin(z) = -1j*log(1j*z + sqrt(1-z*z)) There are two branch cuts: one running from 1 along the positive real axis to infinity, continuous from below, and one running from -1 along the negative real axis to -infinity, continuous from above. """ z = complexify(z) # Formula suggested by Kahan x = z.real # the conjugation in the next line is important: it makes sure # that the right thing happens when the imaginary part of z is # -0., conj1 = complex(1).conjugate() s1 = sqrt(conj1-z) s2 = sqrt(conj1+z) # no possibility of cancellation error here: # both s1.real and s2.real are positive, and # s1.imag and s2.imag have opposite signs real_part = atan2(x, s1.real*s2.real-s1.imag*s2.imag) imag_part = fasinh(s1.real*s2.imag-s2.real*s1.imag) return complex_with_signed_zero(real_part, imag_part) def acos(z): """Inverse cosine of the complex argument z, as defined by the formula acos(z) = -1j*log(z + 1j*sqrt(1-z*z)) There are two branch cuts: one running from 1 along the positive real axis to infinity, continuous from below, and one running from -1 along the negative real axis to -infinity, continuous from above. """ z = complexify(z) # Formula suggested by Kahan # conj1 used as in asin conj1 = complex(1).conjugate() s1 = sqrt(conj1-z) s2 = sqrt(conj1+z) real_part = 2. * atan2(s1.real, s2.real) imag_part = fasinh(s2.real*s1.imag - s2.imag*s1.real) return complex_with_signed_zero(real_part, imag_part) def atan(z): """Inverse tangent of the complex argument z, as defined by the formula atan(z) = (log(1+1j*z) - log(1-1j*z))/(2j) Equivalently, atan(z) = atanh(1j*z)/1j. There are two branch cuts: one along the positive imaginary axis from j to j*infinity, continuous from the right, and one along the negative imaginary axis from -j to -j*infinity, continuous from the left. The function atan is undefined at 1j and -1j. """ return -mul_by_j(atanh(mul_by_j(z))) def asinh(z): """Inverse hyperbolic sine of the complex argument z, as defined by the formula asinh(z) = log(z + sqrt(1+z*z)) There are two branch cuts: one running from 1j up the imaginary axis to +infinity j, continuous from the right, and one running from -1j down to -infinity j, continuous from the left. """ return -mul_by_j(asin(mul_by_j(z))) def acosh(z): """Inverse hyperbolic cosine of the complex argument z, as defined by the formula acosh(z) = 2*log(sqrt((z+1)/2) + sqrt((z-1)/2)) There is a single branch cut, running from 1 to -infinity along the real axis. The function is continuous from above on this branch cut. """ z = complexify(z) # Formula suggested by Kahan s1 = sqrt(z-1) s2 = sqrt(z+complex(1).conjugate()) real_part = fasinh(s1.real*s2.real+s1.imag*s2.imag) imag_part = 2*atan2(s1.imag, s2.real) return complex_with_signed_zero(real_part, imag_part) def atanh(z): """Inverse hyperbolic tangent of the complex argument z, as defined by the formula atanh(z) = (log(1+z) - log(1-z))/2 There are two branch cuts: one along the positive real axis from 1 to infinity, continuous from below, and one along the negative real axis from -1 to -infinity, continuous from above. The function atanh is undefined at 1 and -1. """ # This is the trickiest of the inverse hyperbolic functions to get # right. Danger areas are: z near 0, where the precision of z # gets lost in 1+z and 1-z, z near 1 or -1, where one or other of # the logs might blow up, and z near infinity, where we get # subtractive cancellation error in the formula above. # Formula suggested by Kahan: # Kahan suggests that the real part should be computed as # log1p(4*x/((1-x)*(1-x) + y*y))/4 # this works fine for most x >= 0. and y. When x and y are large # attention must be paid to computing the argument 4*x/((1-x)*(1-x)+y*y) # whilst avoiding overflow. # When x == 1 and y is tiny, so that y*y=0, we need another way. # In this case the real part is log(1+4/y^2)/4 = # log((y^2+4)/y^2)/4 = (log(hypot(y,2)) - log(abs(y)))/2, and we # can use this formula directly. # In general, the real part should be # fatanh(2x/(1+x^2+y^2))/2 # Assuming that x >= 0, this is okay except when the argument is close to 1. # If the argument is not close to 1, then our problems are reduced to: # (1) evaluating 2x/(1+x^2+y^2) without unnecessary overflow or underflow, and # (2) writing a real version of atanh # For (1), if |x| >= max(|y|, 2): write as 2/(1+(1/x)^2 + (y/x)^2)/x. # if |y| >= max(|x|, 2): write as 2(x/y)/(1+(x/y)^2 + (1/y)^2)/y; # It's still subject to overflow for very large x and y, though. # It fails when x == 1 and y*y underflows. It's also subject to # overflow when x and y are large. # The obvious formula # # log(hypot(1+x, y)/hypot(1-x, y))/2 # # works except when hypot(1+x, y) and hypot(1-x, y) are almost equal; # that is, when x and y are large. Similarly for # # (hypot(1+x, y) - hypot(1-x, y))/2 z = complexify(z) x, y = z.real, z.imag if x < 0.: return -atanh(-z) # avoid overflow when x or y is large # The formula given by Kahan for x == 1 loses significant accuracy # for y large. However, the usual formula works just fine in this # case. if x == 1. and abs(y) < 1.1547: if abs(y) < 1.1547: real_part = log(hypot(y, 2.)/abs(y))/2. imag_part = copysign(atan2(2., -abs(y))/2, y) return complex_with_signed_zero(real_part, imag_part) else: real_part = flog1p(4/(y*y))/4 imag_part = copysign(atan2(2., -abs(y))/2, y) else: real_part = flog1p(4*x/((1-x)*(1-x) + y*y))/4 imag_part = -atan2(-2*y, (1-x)*(1+x) - y*y)/2 return complex_with_signed_zero(real_part, imag_part) # Problems with existing cmath module # # (1) branch cuts for asinh # (2) continuity of branch cuts for atan, atanh # (3) asinh catastrophically inaccurate for small values (and possibly some large ones too) # (4) asin overflows unnecessarily for 1e160, 1e-160 # (5) tanh overflow/underflow problems (cmath.tanh(1000) gives the amusing response (nannanj); the correct value is close to 1.) # For everything to work correctly with signed zeros, we need: # multiplication does the correct thing with respect to signs # atan2 does the correct thing when x = +/- 0. For example: # atan2(-5., 0.) = pi, atan2(-5., -0.) = -pi # atan2(5., 0.) = 0., atan2(5., -0.) = -0. # abs does the correct thing