blob: 64c56d46538aaf11e14f74635b09cebd01a74e4b [file] [log] [blame]
Guido van Rossume8769491992-08-13 12:14:11 +00001# Complex numbers
Guido van Rossum81a12bc1994-10-08 18:56:41 +00002# ---------------
3
Guido van Rossum72ba6161996-07-30 19:02:01 +00004# [Now that Python has a complex data type built-in, this is not very
5# useful, but it's still a nice example class]
6
Guido van Rossum81a12bc1994-10-08 18:56:41 +00007# This module represents complex numbers as instances of the class Complex.
8# A Complex instance z has two data attribues, z.re (the real part) and z.im
9# (the imaginary part). In fact, z.re and z.im can have any value -- all
10# arithmetic operators work regardless of the type of z.re and z.im (as long
11# as they support numerical operations).
12#
13# The following functions exist (Complex is actually a class):
14# Complex([re [,im]) -> creates a complex number from a real and an imaginary part
15# IsComplex(z) -> true iff z is a complex number (== has .re and .im attributes)
16# ToComplex(z) -> a complex number equal to z; z itself if IsComplex(z) is true
17# if z is a tuple(re, im) it will also be converted
18# PolarToComplex([r [,phi [,fullcircle]]]) ->
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000019# the complex number z for which r == z.radius() and phi == z.angle(fullcircle)
20# (r and phi default to 0)
Guido van Rossum72ba6161996-07-30 19:02:01 +000021# exp(z) -> returns the complex exponential of z. Equivalent to pow(math.e,z).
Guido van Rossum81a12bc1994-10-08 18:56:41 +000022#
23# Complex numbers have the following methods:
24# z.abs() -> absolute value of z
25# z.radius() == z.abs()
26# z.angle([fullcircle]) -> angle from positive X axis; fullcircle gives units
27# z.phi([fullcircle]) == z.angle(fullcircle)
28#
29# These standard functions and unary operators accept complex arguments:
30# abs(z)
31# -z
32# +z
33# not z
34# repr(z) == `z`
35# str(z)
36# hash(z) -> a combination of hash(z.re) and hash(z.im) such that if z.im is zero
37# the result equals hash(z.re)
38# Note that hex(z) and oct(z) are not defined.
39#
40# These conversions accept complex arguments only if their imaginary part is zero:
41# int(z)
Guido van Rossum81a12bc1994-10-08 18:56:41 +000042# float(z)
43#
44# The following operators accept two complex numbers, or one complex number
45# and one real number (int, long or float):
46# z1 + z2
47# z1 - z2
48# z1 * z2
49# z1 / z2
50# pow(z1, z2)
51# cmp(z1, z2)
52# Note that z1 % z2 and divmod(z1, z2) are not defined,
53# nor are shift and mask operations.
54#
55# The standard module math does not support complex numbers.
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000056# The cmath modules should be used instead.
Guido van Rossum81a12bc1994-10-08 18:56:41 +000057#
58# Idea:
59# add a class Polar(r, phi) and mixed-mode arithmetic which
60# chooses the most appropriate type for the result:
61# Complex for +,-,cmp
62# Polar for *,/,pow
Guido van Rossume8769491992-08-13 12:14:11 +000063
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000064import math
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +000065import sys
Guido van Rossume8769491992-08-13 12:14:11 +000066
Guido van Rossum81a12bc1994-10-08 18:56:41 +000067twopi = math.pi*2.0
68halfpi = math.pi/2.0
Guido van Rossume8769491992-08-13 12:14:11 +000069
Guido van Rossum81a12bc1994-10-08 18:56:41 +000070def IsComplex(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000071 return hasattr(obj, 're') and hasattr(obj, 'im')
Guido van Rossume8769491992-08-13 12:14:11 +000072
Guido van Rossum81a12bc1994-10-08 18:56:41 +000073def ToComplex(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000074 if IsComplex(obj):
75 return obj
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000076 elif isinstance(obj, tuple):
77 return Complex(*obj)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000078 else:
79 return Complex(obj)
Guido van Rossum7565b931993-12-17 14:23:52 +000080
Guido van Rossum81a12bc1994-10-08 18:56:41 +000081def PolarToComplex(r = 0, phi = 0, fullcircle = twopi):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000082 phi = phi * (twopi / fullcircle)
83 return Complex(math.cos(phi)*r, math.sin(phi)*r)
Guido van Rossum81a12bc1994-10-08 18:56:41 +000084
85def Re(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000086 if IsComplex(obj):
87 return obj.re
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000088 return obj
Guido van Rossum81a12bc1994-10-08 18:56:41 +000089
90def Im(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000091 if IsComplex(obj):
92 return obj.im
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000093 return 0
Guido van Rossum81a12bc1994-10-08 18:56:41 +000094
95class Complex:
96
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000097 def __init__(self, re=0, im=0):
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +000098 _re = 0
99 _im = 0
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000100 if IsComplex(re):
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000101 _re = re.re
102 _im = re.im
103 else:
104 _re = re
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000105 if IsComplex(im):
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000106 _re = _re - im.im
107 _im = _im + im.re
108 else:
109 _im = _im + im
110 # this class is immutable, so setting self.re directly is
111 # not possible.
112 self.__dict__['re'] = _re
113 self.__dict__['im'] = _im
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000114
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000115 def __setattr__(self, name, value):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000116 raise TypeError('Complex numbers are immutable')
Guido van Rossume8769491992-08-13 12:14:11 +0000117
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000118 def __hash__(self):
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000119 if not self.im:
120 return hash(self.re)
121 return hash((self.re, self.im))
Guido van Rossume8769491992-08-13 12:14:11 +0000122
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000123 def __repr__(self):
124 if not self.im:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000125 return 'Complex(%r)' % (self.re,)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000126 else:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000127 return 'Complex(%r, %r)' % (self.re, self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000128
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000129 def __str__(self):
130 if not self.im:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000131 return repr(self.re)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000132 else:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000133 return 'Complex(%r, %r)' % (self.re, self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000134
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000135 def __neg__(self):
136 return Complex(-self.re, -self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000137
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000138 def __pos__(self):
139 return self
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000140
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000141 def __abs__(self):
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000142 return math.hypot(self.re, self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000143
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000144 def __int__(self):
145 if self.im:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000146 raise ValueError("can't convert Complex with nonzero im to int")
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000147 return int(self.re)
Guido van Rossume8769491992-08-13 12:14:11 +0000148
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000149 def __float__(self):
150 if self.im:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000151 raise ValueError("can't convert Complex with nonzero im to float")
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000152 return float(self.re)
Guido van Rossume8769491992-08-13 12:14:11 +0000153
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000154 def __cmp__(self, other):
155 other = ToComplex(other)
156 return cmp((self.re, self.im), (other.re, other.im))
Guido van Rossume8769491992-08-13 12:14:11 +0000157
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000158 def __rcmp__(self, other):
159 other = ToComplex(other)
160 return cmp(other, self)
Guido van Rossume8769491992-08-13 12:14:11 +0000161
Jack Diederich4dafcc42006-11-28 19:15:13 +0000162 def __bool__(self):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000163 return not (self.re == self.im == 0)
Guido van Rossume8769491992-08-13 12:14:11 +0000164
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000165 abs = radius = __abs__
Guido van Rossume8769491992-08-13 12:14:11 +0000166
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000167 def angle(self, fullcircle = twopi):
168 return (fullcircle/twopi) * ((halfpi - math.atan2(self.re, self.im)) % twopi)
Guido van Rossume8769491992-08-13 12:14:11 +0000169
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000170 phi = angle
Guido van Rossume8769491992-08-13 12:14:11 +0000171
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000172 def __add__(self, other):
173 other = ToComplex(other)
174 return Complex(self.re + other.re, self.im + other.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000175
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000176 __radd__ = __add__
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000177
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000178 def __sub__(self, other):
179 other = ToComplex(other)
180 return Complex(self.re - other.re, self.im - other.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000181
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000182 def __rsub__(self, other):
183 other = ToComplex(other)
184 return other - self
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000185
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000186 def __mul__(self, other):
187 other = ToComplex(other)
188 return Complex(self.re*other.re - self.im*other.im,
189 self.re*other.im + self.im*other.re)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000190
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000191 __rmul__ = __mul__
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000192
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000193 def __div__(self, other):
194 other = ToComplex(other)
195 d = float(other.re*other.re + other.im*other.im)
Collin Winter6f2df4d2007-07-17 20:59:35 +0000196 if not d: raise ZeroDivisionError('Complex division')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000197 return Complex((self.re*other.re + self.im*other.im) / d,
198 (self.im*other.re - self.re*other.im) / d)
199
200 def __rdiv__(self, other):
201 other = ToComplex(other)
202 return other / self
203
204 def __pow__(self, n, z=None):
205 if z is not None:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000206 raise TypeError('Complex does not support ternary pow()')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000207 if IsComplex(n):
208 if n.im:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000209 if self.im: raise TypeError('Complex to the Complex power')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000210 else: return exp(math.log(self.re)*n)
211 n = n.re
212 r = pow(self.abs(), n)
213 phi = n*self.angle()
214 return Complex(math.cos(phi)*r, math.sin(phi)*r)
215
216 def __rpow__(self, base):
217 base = ToComplex(base)
218 return pow(base, self)
219
Guido van Rossum72ba6161996-07-30 19:02:01 +0000220def exp(z):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000221 r = math.exp(z.re)
222 return Complex(math.cos(z.im)*r,math.sin(z.im)*r)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000223
224
225def checkop(expr, a, b, value, fuzz = 1e-6):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000226 print(' ', a, 'and', b, end=' ')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000227 try:
228 result = eval(expr)
229 except:
Neal Norwitzac3625f2006-03-17 05:49:33 +0000230 result = sys.exc_info()[0]
Collin Winter6f2df4d2007-07-17 20:59:35 +0000231 print('->', result)
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000232 if isinstance(result, str) or isinstance(value, str):
233 ok = (result == value)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000234 else:
235 ok = abs(result - value) <= fuzz
236 if not ok:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000237 print('!!\t!!\t!! should be', value, 'diff', abs(result - value))
Guido van Rossume8769491992-08-13 12:14:11 +0000238
Guido van Rossume8769491992-08-13 12:14:11 +0000239def test():
Collin Winter6f2df4d2007-07-17 20:59:35 +0000240 print('test constructors')
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000241 constructor_test = (
242 # "expect" is an array [re,im] "got" the Complex.
243 ( (0,0), Complex() ),
244 ( (0,0), Complex() ),
245 ( (1,0), Complex(1) ),
246 ( (0,1), Complex(0,1) ),
247 ( (1,2), Complex(Complex(1,2)) ),
248 ( (1,3), Complex(Complex(1,2),1) ),
249 ( (0,0), Complex(0,Complex(0,0)) ),
250 ( (3,4), Complex(3,Complex(4)) ),
251 ( (-1,3), Complex(1,Complex(3,2)) ),
252 ( (-7,6), Complex(Complex(1,2),Complex(4,8)) ) )
253 cnt = [0,0]
254 for t in constructor_test:
255 cnt[0] += 1
256 if ((t[0][0]!=t[1].re)or(t[0][1]!=t[1].im)):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000257 print(" expected", t[0], "got", t[1])
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000258 cnt[1] += 1
Collin Winter6f2df4d2007-07-17 20:59:35 +0000259 print(" ", cnt[1], "of", cnt[0], "tests failed")
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000260 # test operators
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000261 testsuite = {
262 'a+b': [
263 (1, 10, 11),
264 (1, Complex(0,10), Complex(1,10)),
265 (Complex(0,10), 1, Complex(1,10)),
266 (Complex(0,10), Complex(1), Complex(1,10)),
267 (Complex(1), Complex(0,10), Complex(1,10)),
268 ],
269 'a-b': [
270 (1, 10, -9),
271 (1, Complex(0,10), Complex(1,-10)),
272 (Complex(0,10), 1, Complex(-1,10)),
273 (Complex(0,10), Complex(1), Complex(-1,10)),
274 (Complex(1), Complex(0,10), Complex(1,-10)),
275 ],
276 'a*b': [
277 (1, 10, 10),
278 (1, Complex(0,10), Complex(0, 10)),
279 (Complex(0,10), 1, Complex(0,10)),
280 (Complex(0,10), Complex(1), Complex(0,10)),
281 (Complex(1), Complex(0,10), Complex(0,10)),
282 ],
283 'a/b': [
284 (1., 10, 0.1),
285 (1, Complex(0,10), Complex(0, -0.1)),
286 (Complex(0, 10), 1, Complex(0, 10)),
287 (Complex(0, 10), Complex(1), Complex(0, 10)),
288 (Complex(1), Complex(0,10), Complex(0, -0.1)),
289 ],
290 'pow(a,b)': [
291 (1, 10, 1),
292 (1, Complex(0,10), 1),
293 (Complex(0,10), 1, Complex(0,10)),
294 (Complex(0,10), Complex(1), Complex(0,10)),
295 (Complex(1), Complex(0,10), 1),
296 (2, Complex(4,0), 16),
297 ],
298 'cmp(a,b)': [
299 (1, 10, -1),
300 (1, Complex(0,10), 1),
301 (Complex(0,10), 1, -1),
302 (Complex(0,10), Complex(1), -1),
303 (Complex(1), Complex(0,10), 1),
304 ],
305 }
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000306 for expr in sorted(testsuite):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000307 print(expr + ':')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000308 t = (expr,)
309 for item in testsuite[expr]:
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000310 checkop(*(t+item))
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000311
Guido van Rossume8769491992-08-13 12:14:11 +0000312
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000313if __name__ == '__main__':
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000314 test()