| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | # Complex numbers | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | # --------------- | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1996-07-30 19:02:01 +00:00
										 |  |  | # [Now that Python has a complex data type built-in, this is not very | 
					
						
							|  |  |  | # useful, but it's still a nice example class] | 
					
						
							|  |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | # This module represents complex numbers as instances of the class Complex. | 
					
						
							|  |  |  | # A Complex instance z has two data attribues, z.re (the real part) and z.im | 
					
						
							|  |  |  | # (the imaginary part).  In fact, z.re and z.im can have any value -- all | 
					
						
							|  |  |  | # arithmetic operators work regardless of the type of z.re and z.im (as long | 
					
						
							|  |  |  | # as they support numerical operations). | 
					
						
							|  |  |  | # | 
					
						
							|  |  |  | # The following functions exist (Complex is actually a class): | 
					
						
							|  |  |  | # Complex([re [,im]) -> creates a complex number from a real and an imaginary part | 
					
						
							|  |  |  | # IsComplex(z) -> true iff z is a complex number (== has .re and .im attributes) | 
					
						
							|  |  |  | # ToComplex(z) -> a complex number equal to z; z itself if IsComplex(z) is true | 
					
						
							|  |  |  | #                 if z is a tuple(re, im) it will also be converted | 
					
						
							|  |  |  | # PolarToComplex([r [,phi [,fullcircle]]]) -> | 
					
						
							|  |  |  | #	the complex number z for which r == z.radius() and phi == z.angle(fullcircle) | 
					
						
							|  |  |  | #	(r and phi default to 0) | 
					
						
							| 
									
										
										
										
											1996-07-30 19:02:01 +00:00
										 |  |  | # exp(z) -> returns the complex exponential of z. Equivalent to pow(math.e,z). | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | # | 
					
						
							|  |  |  | # Complex numbers have the following methods: | 
					
						
							|  |  |  | # z.abs() -> absolute value of z | 
					
						
							|  |  |  | # z.radius() == z.abs() | 
					
						
							|  |  |  | # z.angle([fullcircle]) -> angle from positive X axis; fullcircle gives units | 
					
						
							|  |  |  | # z.phi([fullcircle]) == z.angle(fullcircle) | 
					
						
							|  |  |  | # | 
					
						
							|  |  |  | # These standard functions and unary operators accept complex arguments: | 
					
						
							|  |  |  | # abs(z) | 
					
						
							|  |  |  | # -z | 
					
						
							|  |  |  | # +z | 
					
						
							|  |  |  | # not z | 
					
						
							|  |  |  | # repr(z) == `z` | 
					
						
							|  |  |  | # str(z) | 
					
						
							|  |  |  | # hash(z) -> a combination of hash(z.re) and hash(z.im) such that if z.im is zero | 
					
						
							|  |  |  | #            the result equals hash(z.re) | 
					
						
							|  |  |  | # Note that hex(z) and oct(z) are not defined. | 
					
						
							|  |  |  | # | 
					
						
							|  |  |  | # These conversions accept complex arguments only if their imaginary part is zero: | 
					
						
							|  |  |  | # int(z) | 
					
						
							|  |  |  | # long(z) | 
					
						
							|  |  |  | # float(z) | 
					
						
							|  |  |  | # | 
					
						
							|  |  |  | # The following operators accept two complex numbers, or one complex number | 
					
						
							|  |  |  | # and one real number (int, long or float): | 
					
						
							|  |  |  | # z1 + z2 | 
					
						
							|  |  |  | # z1 - z2 | 
					
						
							|  |  |  | # z1 * z2 | 
					
						
							|  |  |  | # z1 / z2 | 
					
						
							|  |  |  | # pow(z1, z2) | 
					
						
							|  |  |  | # cmp(z1, z2) | 
					
						
							|  |  |  | # Note that z1 % z2 and divmod(z1, z2) are not defined, | 
					
						
							|  |  |  | # nor are shift and mask operations. | 
					
						
							|  |  |  | # | 
					
						
							|  |  |  | # The standard module math does not support complex numbers. | 
					
						
							|  |  |  | # (I suppose it would be easy to implement a cmath module.) | 
					
						
							|  |  |  | # | 
					
						
							|  |  |  | # Idea: | 
					
						
							|  |  |  | # add a class Polar(r, phi) and mixed-mode arithmetic which | 
					
						
							|  |  |  | # chooses the most appropriate type for the result: | 
					
						
							|  |  |  | # Complex for +,-,cmp | 
					
						
							|  |  |  | # Polar   for *,/,pow | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							|  |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | import types, math | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | twopi = math.pi*2.0 | 
					
						
							|  |  |  | halfpi = math.pi/2.0 | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | def IsComplex(obj): | 
					
						
							|  |  |  | 	return hasattr(obj, 're') and hasattr(obj, 'im') | 
					
						
							| 
									
										
										
										
											1993-12-17 14:23:52 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | def ToComplex(obj): | 
					
						
							|  |  |  | 	if IsComplex(obj): | 
					
						
							|  |  |  | 		return obj | 
					
						
							|  |  |  | 	elif type(obj) == types.TupleType: | 
					
						
							|  |  |  | 		return apply(Complex, obj) | 
					
						
							|  |  |  | 	else: | 
					
						
							|  |  |  | 		return Complex(obj) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | def PolarToComplex(r = 0, phi = 0, fullcircle = twopi): | 
					
						
							|  |  |  | 	phi = phi * (twopi / fullcircle) | 
					
						
							|  |  |  | 	return Complex(math.cos(phi)*r, math.sin(phi)*r) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | def Re(obj): | 
					
						
							|  |  |  | 	if IsComplex(obj): | 
					
						
							|  |  |  | 		return obj.re | 
					
						
							|  |  |  | 	else: | 
					
						
							|  |  |  | 		return obj | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | def Im(obj): | 
					
						
							|  |  |  | 	if IsComplex(obj): | 
					
						
							|  |  |  | 		return obj.im | 
					
						
							|  |  |  | 	else: | 
					
						
							|  |  |  | 		return obj | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | class Complex: | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __init__(self, re=0, im=0): | 
					
						
							|  |  |  | 		if IsComplex(re): | 
					
						
							|  |  |  | 			im = i + Complex(0, re.im) | 
					
						
							|  |  |  | 			re = re.re | 
					
						
							|  |  |  | 		if IsComplex(im): | 
					
						
							|  |  |  | 			re = re - im.im | 
					
						
							|  |  |  | 			im = im.re | 
					
						
							|  |  |  | 		self.__dict__['re'] = re | 
					
						
							|  |  |  | 		self.__dict__['im'] = im | 
					
						
							|  |  |  | 	 | 
					
						
							|  |  |  | 	def __setattr__(self, name, value): | 
					
						
							|  |  |  | 			raise TypeError, 'Complex numbers are immutable' | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __hash__(self): | 
					
						
							|  |  |  | 		if not self.im: return hash(self.re) | 
					
						
							|  |  |  | 		mod = sys.maxint + 1L | 
					
						
							|  |  |  | 		return int((hash(self.re) + 2L*hash(self.im) + mod) % (2L*mod) - mod) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							|  |  |  | 	def __repr__(self): | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 		if not self.im: | 
					
						
							|  |  |  | 			return 'Complex(%s)' % `self.re` | 
					
						
							|  |  |  | 		else: | 
					
						
							|  |  |  | 			return 'Complex(%s, %s)' % (`self.re`, `self.im`) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	def __str__(self): | 
					
						
							|  |  |  | 		if not self.im: | 
					
						
							|  |  |  | 			return `self.re` | 
					
						
							|  |  |  | 		else: | 
					
						
							|  |  |  | 			return 'Complex(%s, %s)' % (`self.re`, `self.im`) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __neg__(self): | 
					
						
							|  |  |  | 		return Complex(-self.re, -self.im) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __pos__(self): | 
					
						
							|  |  |  | 		return self | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __abs__(self): | 
					
						
							|  |  |  | 		# XXX could be done differently to avoid overflow! | 
					
						
							|  |  |  | 		return math.sqrt(self.re*self.re + self.im*self.im) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __int__(self): | 
					
						
							|  |  |  | 		if self.im: | 
					
						
							|  |  |  | 			raise ValueError, "can't convert Complex with nonzero im to int" | 
					
						
							|  |  |  | 		return int(self.re) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __long__(self): | 
					
						
							|  |  |  | 		if self.im: | 
					
						
							|  |  |  | 			raise ValueError, "can't convert Complex with nonzero im to long" | 
					
						
							|  |  |  | 		return long(self.re) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							|  |  |  | 	def __float__(self): | 
					
						
							|  |  |  | 		if self.im: | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 			raise ValueError, "can't convert Complex with nonzero im to float" | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 		return float(self.re) | 
					
						
							|  |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	def __cmp__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		return cmp((self.re, self.im), (other.re, other.im)) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	def __rcmp__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		return cmp(other, self) | 
					
						
							|  |  |  | 	 | 
					
						
							|  |  |  | 	def __nonzero__(self): | 
					
						
							|  |  |  | 		return not (self.re == self.im == 0) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	abs = radius = __abs__ | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	def angle(self, fullcircle = twopi): | 
					
						
							|  |  |  | 		return (fullcircle/twopi) * ((halfpi - math.atan2(self.re, self.im)) % twopi) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	phi = angle | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	def __add__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		return Complex(self.re + other.re, self.im + other.im) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	__radd__ = __add__ | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	def __sub__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		return Complex(self.re - other.re, self.im - other.im) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __rsub__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		return other - self | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __mul__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		return Complex(self.re*other.re - self.im*other.im, | 
					
						
							|  |  |  | 		               self.re*other.im + self.im*other.re) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	__rmul__ = __mul__ | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __div__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		d = float(other.re*other.re + other.im*other.im) | 
					
						
							|  |  |  | 		if not d: raise ZeroDivisionError, 'Complex division' | 
					
						
							|  |  |  | 		return Complex((self.re*other.re + self.im*other.im) / d, | 
					
						
							|  |  |  | 		               (self.im*other.re - self.re*other.im) / d) | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __rdiv__(self, other): | 
					
						
							|  |  |  | 		other = ToComplex(other) | 
					
						
							|  |  |  | 		return other / self | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | 	def __pow__(self, n, z=None): | 
					
						
							|  |  |  | 		if z is not None: | 
					
						
							|  |  |  | 			raise TypeError, 'Complex does not support ternary pow()' | 
					
						
							|  |  |  | 		if IsComplex(n): | 
					
						
							| 
									
										
										
										
											1996-07-30 19:02:01 +00:00
										 |  |  | 			if n.im:  | 
					
						
							|  |  |  | 			  if self.im: raise TypeError, 'Complex to the Complex power' | 
					
						
							|  |  |  | 			  else: return exp(math.log(self.re)*n) | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 			n = n.re | 
					
						
							|  |  |  | 		r = pow(self.abs(), n) | 
					
						
							|  |  |  | 		phi = n*self.angle() | 
					
						
							|  |  |  | 		return Complex(math.cos(phi)*r, math.sin(phi)*r) | 
					
						
							|  |  |  | 	 | 
					
						
							|  |  |  | 	def __rpow__(self, base): | 
					
						
							|  |  |  | 		base = ToComplex(base) | 
					
						
							|  |  |  | 		return pow(base, self) | 
					
						
							| 
									
										
										
										
											1996-07-30 19:02:01 +00:00
										 |  |  | 		 | 
					
						
							|  |  |  | def exp(z): | 
					
						
							|  |  |  | 	r = math.exp(z.re) | 
					
						
							|  |  |  | 	return Complex(math.cos(z.im)*r,math.sin(z.im)*r) | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 
 | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | def checkop(expr, a, b, value, fuzz = 1e-6): | 
					
						
							|  |  |  | 	import sys | 
					
						
							|  |  |  | 	print '       ', a, 'and', b, | 
					
						
							|  |  |  | 	try: | 
					
						
							|  |  |  | 		result = eval(expr) | 
					
						
							|  |  |  | 	except: | 
					
						
							|  |  |  | 		result = sys.exc_type | 
					
						
							|  |  |  | 	print '->', result | 
					
						
							|  |  |  | 	if (type(result) == type('') or type(value) == type('')): | 
					
						
							|  |  |  | 		ok = result == value | 
					
						
							|  |  |  | 	else: | 
					
						
							|  |  |  | 		ok = abs(result - value) <= fuzz | 
					
						
							|  |  |  | 	if not ok: | 
					
						
							|  |  |  | 		print '!!\t!!\t!! should be', value, 'diff', abs(result - value) | 
					
						
							| 
									
										
										
										
											1992-08-13 12:14:11 +00:00
										 |  |  | 
 | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | def test(): | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 	testsuite = { | 
					
						
							|  |  |  | 		'a+b': [ | 
					
						
							|  |  |  | 			(1, 10, 11), | 
					
						
							|  |  |  | 			(1, Complex(0,10), Complex(1,10)), | 
					
						
							|  |  |  | 			(Complex(0,10), 1, Complex(1,10)), | 
					
						
							|  |  |  | 			(Complex(0,10), Complex(1), Complex(1,10)), | 
					
						
							|  |  |  | 			(Complex(1), Complex(0,10), Complex(1,10)), | 
					
						
							|  |  |  | 		], | 
					
						
							|  |  |  | 		'a-b': [ | 
					
						
							|  |  |  | 			(1, 10, -9), | 
					
						
							|  |  |  | 			(1, Complex(0,10), Complex(1,-10)), | 
					
						
							|  |  |  | 			(Complex(0,10), 1, Complex(-1,10)), | 
					
						
							|  |  |  | 			(Complex(0,10), Complex(1), Complex(-1,10)), | 
					
						
							|  |  |  | 			(Complex(1), Complex(0,10), Complex(1,-10)), | 
					
						
							|  |  |  | 		], | 
					
						
							|  |  |  | 		'a*b': [ | 
					
						
							|  |  |  | 			(1, 10, 10), | 
					
						
							|  |  |  | 			(1, Complex(0,10), Complex(0, 10)), | 
					
						
							|  |  |  | 			(Complex(0,10), 1, Complex(0,10)), | 
					
						
							|  |  |  | 			(Complex(0,10), Complex(1), Complex(0,10)), | 
					
						
							|  |  |  | 			(Complex(1), Complex(0,10), Complex(0,10)), | 
					
						
							|  |  |  | 		], | 
					
						
							|  |  |  | 		'a/b': [ | 
					
						
							|  |  |  | 			(1., 10, 0.1), | 
					
						
							|  |  |  | 			(1, Complex(0,10), Complex(0, -0.1)), | 
					
						
							|  |  |  | 			(Complex(0, 10), 1, Complex(0, 10)), | 
					
						
							|  |  |  | 			(Complex(0, 10), Complex(1), Complex(0, 10)), | 
					
						
							|  |  |  | 			(Complex(1), Complex(0,10), Complex(0, -0.1)), | 
					
						
							|  |  |  | 		], | 
					
						
							|  |  |  | 		'pow(a,b)': [ | 
					
						
							|  |  |  | 			(1, 10, 1), | 
					
						
							| 
									
										
										
										
											1997-12-09 19:38:39 +00:00
										 |  |  | 			(1, Complex(0,10), 1), | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 			(Complex(0,10), 1, Complex(0,10)), | 
					
						
							|  |  |  | 			(Complex(0,10), Complex(1), Complex(0,10)), | 
					
						
							| 
									
										
										
										
											1997-12-09 19:38:39 +00:00
										 |  |  | 			(Complex(1), Complex(0,10), 1), | 
					
						
							| 
									
										
										
										
											1994-10-08 18:56:41 +00:00
										 |  |  | 			(2, Complex(4,0), 16), | 
					
						
							|  |  |  | 		], | 
					
						
							|  |  |  | 		'cmp(a,b)': [ | 
					
						
							|  |  |  | 			(1, 10, -1), | 
					
						
							|  |  |  | 			(1, Complex(0,10), 1), | 
					
						
							|  |  |  | 			(Complex(0,10), 1, -1), | 
					
						
							|  |  |  | 			(Complex(0,10), Complex(1), -1), | 
					
						
							|  |  |  | 			(Complex(1), Complex(0,10), 1), | 
					
						
							|  |  |  | 		], | 
					
						
							|  |  |  | 	} | 
					
						
							|  |  |  | 	exprs = testsuite.keys() | 
					
						
							|  |  |  | 	exprs.sort() | 
					
						
							|  |  |  | 	for expr in exprs: | 
					
						
							|  |  |  | 		print expr + ':' | 
					
						
							|  |  |  | 		t = (expr,) | 
					
						
							|  |  |  | 		for item in testsuite[expr]: | 
					
						
							|  |  |  | 			apply(checkop, t+item) | 
					
						
							|  |  |  | 	 | 
					
						
							|  |  |  | 
 | 
					
						
							|  |  |  | if __name__ == '__main__': | 
					
						
							|  |  |  | 	test() |