blob: b6c4b10e006f012470316898d7b369ad4916091c [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)
42# long(z)
43# float(z)
44#
45# The following operators accept two complex numbers, or one complex number
46# and one real number (int, long or float):
47# z1 + z2
48# z1 - z2
49# z1 * z2
50# z1 / z2
51# pow(z1, z2)
52# cmp(z1, z2)
53# Note that z1 % z2 and divmod(z1, z2) are not defined,
54# nor are shift and mask operations.
55#
56# The standard module math does not support complex numbers.
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000057# The cmath modules should be used instead.
Guido van Rossum81a12bc1994-10-08 18:56:41 +000058#
59# Idea:
60# add a class Polar(r, phi) and mixed-mode arithmetic which
61# chooses the most appropriate type for the result:
62# Complex for +,-,cmp
63# Polar for *,/,pow
Guido van Rossume8769491992-08-13 12:14:11 +000064
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000065import math
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +000066import sys
Guido van Rossume8769491992-08-13 12:14:11 +000067
Guido van Rossum81a12bc1994-10-08 18:56:41 +000068twopi = math.pi*2.0
69halfpi = math.pi/2.0
Guido van Rossume8769491992-08-13 12:14:11 +000070
Guido van Rossum81a12bc1994-10-08 18:56:41 +000071def IsComplex(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000072 return hasattr(obj, 're') and hasattr(obj, 'im')
Guido van Rossume8769491992-08-13 12:14:11 +000073
Guido van Rossum81a12bc1994-10-08 18:56:41 +000074def ToComplex(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000075 if IsComplex(obj):
76 return obj
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000077 elif isinstance(obj, tuple):
78 return Complex(*obj)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000079 else:
80 return Complex(obj)
Guido van Rossum7565b931993-12-17 14:23:52 +000081
Guido van Rossum81a12bc1994-10-08 18:56:41 +000082def PolarToComplex(r = 0, phi = 0, fullcircle = twopi):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000083 phi = phi * (twopi / fullcircle)
84 return Complex(math.cos(phi)*r, math.sin(phi)*r)
Guido van Rossum81a12bc1994-10-08 18:56:41 +000085
86def Re(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000087 if IsComplex(obj):
88 return obj.re
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000089 return obj
Guido van Rossum81a12bc1994-10-08 18:56:41 +000090
91def Im(obj):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000092 if IsComplex(obj):
93 return obj.im
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +000094 return 0
Guido van Rossum81a12bc1994-10-08 18:56:41 +000095
96class Complex:
97
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +000098 def __init__(self, re=0, im=0):
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +000099 _re = 0
100 _im = 0
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000101 if IsComplex(re):
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000102 _re = re.re
103 _im = re.im
104 else:
105 _re = re
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000106 if IsComplex(im):
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000107 _re = _re - im.im
108 _im = _im + im.re
109 else:
110 _im = _im + im
111 # this class is immutable, so setting self.re directly is
112 # not possible.
113 self.__dict__['re'] = _re
114 self.__dict__['im'] = _im
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000115
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000116 def __setattr__(self, name, value):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000117 raise TypeError('Complex numbers are immutable')
Guido van Rossume8769491992-08-13 12:14:11 +0000118
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000119 def __hash__(self):
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000120 if not self.im:
121 return hash(self.re)
122 return hash((self.re, self.im))
Guido van Rossume8769491992-08-13 12:14:11 +0000123
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000124 def __repr__(self):
125 if not self.im:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000126 return 'Complex(%r)' % (self.re,)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000127 else:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000128 return 'Complex(%r, %r)' % (self.re, self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000129
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000130 def __str__(self):
131 if not self.im:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000132 return repr(self.re)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000133 else:
Walter Dörwald70a6b492004-02-12 17:35:32 +0000134 return 'Complex(%r, %r)' % (self.re, self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000135
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000136 def __neg__(self):
137 return Complex(-self.re, -self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000138
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000139 def __pos__(self):
140 return self
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000141
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000142 def __abs__(self):
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000143 return math.hypot(self.re, self.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000144
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000145 def __int__(self):
146 if self.im:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000147 raise ValueError("can't convert Complex with nonzero im to int")
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000148 return int(self.re)
Guido van Rossume8769491992-08-13 12:14:11 +0000149
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000150 def __long__(self):
151 if self.im:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000152 raise ValueError("can't convert Complex with nonzero im to long")
153 return int(self.re)
Guido van Rossume8769491992-08-13 12:14:11 +0000154
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000155 def __float__(self):
156 if self.im:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000157 raise ValueError("can't convert Complex with nonzero im to float")
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000158 return float(self.re)
Guido van Rossume8769491992-08-13 12:14:11 +0000159
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000160 def __cmp__(self, other):
161 other = ToComplex(other)
162 return cmp((self.re, self.im), (other.re, other.im))
Guido van Rossume8769491992-08-13 12:14:11 +0000163
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000164 def __rcmp__(self, other):
165 other = ToComplex(other)
166 return cmp(other, self)
Guido van Rossume8769491992-08-13 12:14:11 +0000167
Jack Diederich4dafcc42006-11-28 19:15:13 +0000168 def __bool__(self):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000169 return not (self.re == self.im == 0)
Guido van Rossume8769491992-08-13 12:14:11 +0000170
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000171 abs = radius = __abs__
Guido van Rossume8769491992-08-13 12:14:11 +0000172
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000173 def angle(self, fullcircle = twopi):
174 return (fullcircle/twopi) * ((halfpi - math.atan2(self.re, self.im)) % twopi)
Guido van Rossume8769491992-08-13 12:14:11 +0000175
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000176 phi = angle
Guido van Rossume8769491992-08-13 12:14:11 +0000177
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000178 def __add__(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 __radd__ = __add__
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000183
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000184 def __sub__(self, other):
185 other = ToComplex(other)
186 return Complex(self.re - other.re, self.im - other.im)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000187
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000188 def __rsub__(self, other):
189 other = ToComplex(other)
190 return other - self
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000191
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000192 def __mul__(self, other):
193 other = ToComplex(other)
194 return Complex(self.re*other.re - self.im*other.im,
195 self.re*other.im + self.im*other.re)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000196
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000197 __rmul__ = __mul__
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000198
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000199 def __div__(self, other):
200 other = ToComplex(other)
201 d = float(other.re*other.re + other.im*other.im)
Collin Winter6f2df4d2007-07-17 20:59:35 +0000202 if not d: raise ZeroDivisionError('Complex division')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000203 return Complex((self.re*other.re + self.im*other.im) / d,
204 (self.im*other.re - self.re*other.im) / d)
205
206 def __rdiv__(self, other):
207 other = ToComplex(other)
208 return other / self
209
210 def __pow__(self, n, z=None):
211 if z is not None:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000212 raise TypeError('Complex does not support ternary pow()')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000213 if IsComplex(n):
214 if n.im:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000215 if self.im: raise TypeError('Complex to the Complex power')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000216 else: return exp(math.log(self.re)*n)
217 n = n.re
218 r = pow(self.abs(), n)
219 phi = n*self.angle()
220 return Complex(math.cos(phi)*r, math.sin(phi)*r)
221
222 def __rpow__(self, base):
223 base = ToComplex(base)
224 return pow(base, self)
225
Guido van Rossum72ba6161996-07-30 19:02:01 +0000226def exp(z):
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000227 r = math.exp(z.re)
228 return Complex(math.cos(z.im)*r,math.sin(z.im)*r)
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000229
230
231def checkop(expr, a, b, value, fuzz = 1e-6):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000232 print(' ', a, 'and', b, end=' ')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000233 try:
234 result = eval(expr)
235 except:
Neal Norwitzac3625f2006-03-17 05:49:33 +0000236 result = sys.exc_info()[0]
Collin Winter6f2df4d2007-07-17 20:59:35 +0000237 print('->', result)
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000238 if isinstance(result, str) or isinstance(value, str):
239 ok = (result == value)
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000240 else:
241 ok = abs(result - value) <= fuzz
242 if not ok:
Collin Winter6f2df4d2007-07-17 20:59:35 +0000243 print('!!\t!!\t!! should be', value, 'diff', abs(result - value))
Guido van Rossume8769491992-08-13 12:14:11 +0000244
Guido van Rossume8769491992-08-13 12:14:11 +0000245def test():
Collin Winter6f2df4d2007-07-17 20:59:35 +0000246 print('test constructors')
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000247 constructor_test = (
248 # "expect" is an array [re,im] "got" the Complex.
249 ( (0,0), Complex() ),
250 ( (0,0), Complex() ),
251 ( (1,0), Complex(1) ),
252 ( (0,1), Complex(0,1) ),
253 ( (1,2), Complex(Complex(1,2)) ),
254 ( (1,3), Complex(Complex(1,2),1) ),
255 ( (0,0), Complex(0,Complex(0,0)) ),
256 ( (3,4), Complex(3,Complex(4)) ),
257 ( (-1,3), Complex(1,Complex(3,2)) ),
258 ( (-7,6), Complex(Complex(1,2),Complex(4,8)) ) )
259 cnt = [0,0]
260 for t in constructor_test:
261 cnt[0] += 1
262 if ((t[0][0]!=t[1].re)or(t[0][1]!=t[1].im)):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000263 print(" expected", t[0], "got", t[1])
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000264 cnt[1] += 1
Collin Winter6f2df4d2007-07-17 20:59:35 +0000265 print(" ", cnt[1], "of", cnt[0], "tests failed")
Martin v. Löwis4a1e48c2005-04-09 10:51:19 +0000266 # test operators
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000267 testsuite = {
268 'a+b': [
269 (1, 10, 11),
270 (1, Complex(0,10), Complex(1,10)),
271 (Complex(0,10), 1, Complex(1,10)),
272 (Complex(0,10), Complex(1), Complex(1,10)),
273 (Complex(1), Complex(0,10), Complex(1,10)),
274 ],
275 'a-b': [
276 (1, 10, -9),
277 (1, Complex(0,10), Complex(1,-10)),
278 (Complex(0,10), 1, Complex(-1,10)),
279 (Complex(0,10), Complex(1), Complex(-1,10)),
280 (Complex(1), Complex(0,10), Complex(1,-10)),
281 ],
282 'a*b': [
283 (1, 10, 10),
284 (1, Complex(0,10), Complex(0, 10)),
285 (Complex(0,10), 1, Complex(0,10)),
286 (Complex(0,10), Complex(1), Complex(0,10)),
287 (Complex(1), Complex(0,10), Complex(0,10)),
288 ],
289 'a/b': [
290 (1., 10, 0.1),
291 (1, Complex(0,10), Complex(0, -0.1)),
292 (Complex(0, 10), 1, Complex(0, 10)),
293 (Complex(0, 10), Complex(1), Complex(0, 10)),
294 (Complex(1), Complex(0,10), Complex(0, -0.1)),
295 ],
296 'pow(a,b)': [
297 (1, 10, 1),
298 (1, Complex(0,10), 1),
299 (Complex(0,10), 1, Complex(0,10)),
300 (Complex(0,10), Complex(1), Complex(0,10)),
301 (Complex(1), Complex(0,10), 1),
302 (2, Complex(4,0), 16),
303 ],
304 'cmp(a,b)': [
305 (1, 10, -1),
306 (1, Complex(0,10), 1),
307 (Complex(0,10), 1, -1),
308 (Complex(0,10), Complex(1), -1),
309 (Complex(1), Complex(0,10), 1),
310 ],
311 }
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000312 for expr in sorted(testsuite):
Collin Winter6f2df4d2007-07-17 20:59:35 +0000313 print(expr + ':')
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000314 t = (expr,)
315 for item in testsuite[expr]:
Raymond Hettingerbdaad8c2005-04-09 14:55:07 +0000316 checkop(*(t+item))
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000317
Guido van Rossume8769491992-08-13 12:14:11 +0000318
Guido van Rossum81a12bc1994-10-08 18:56:41 +0000319if __name__ == '__main__':
Andrew M. Kuchling946c53e2003-04-24 17:13:18 +0000320 test()