######################################################################################
# GAUSSINT.PY
# Basic Number Theory functions implemented in Python
# Note: Currently requires Python 2.x (uses +=, %= and other 2.x-isms)
# Author: Robert Campbell, <campbell@math.umbc.edu>
# Date: 28 February, 2005
# Version 1.0
# Requirements:
#	Requires at least Python 2.x (runs fine on Python 2.2)
# Bugs:
#	Division relies on floating point (complex) computations - fails for large values
#		(Need to reimplement, even at the cost of speed)
#   Need to clean up input/output routines, so a-bi does not print a+-bi and so
#       one can easily cut and paste to input previous output
######################################################################################
"""Gaussian Integer functions.
   Functions implemented are:
   		Arithmetic functions: +,*,/,%,**(exponentiation)
   		a.gcd(b) - Compute the greatest common divisor of a and b.
   		a.xgcd(b) - Extended gcd - return gcd(a,b) and x,y such that gcd(a,b)=xa+yb.
   		b.powmod(e,m) - Compute b**e mod m.
   A list of the functions implemented in GaussInt is printed by the command info().
   Usage: from gaussint import * """
   
   
Version = 'GAUSSINT.PY, version 0.1, 10 Jan, 2003, by Robert Campbell, campbell@math.umbc.edu'

import math  # Use floor
import types # Use IntType, LongType

class GaussInt:
	"Gaussian Integers"
	def __init__(self, realpart=0, imagpart=0):
		self.r = realpart
		self.i = imagpart
	def __str__(self):  # Overload string conversion used by print
		return "(" + str(self.r) + " + " + str(self.i) + "i)"
	def __repr__(self):  # Overload conversion used for output
		return "(" + str(self.r) + " + " + str(self.i) + "i)"
	def __eq__(self,other):  # Overload the "==" test operator
		return (self.r == other.r) and (self.i == other.i)
	def __ne__(self,other):  # Overload the "==" test operator
		return not (self == other)
	def norm(self):
		return self.r*self.r + self.i*self.i
	def add(self,summand):
		sum = GaussInt()
		sum.r = self.r + summand.r
		sum.i = self.i + summand.i
		return sum
	def __add__(self,summand):  # Overload the "+" operator
		if ((type(summand) == types.IntType) or (type(summand) == types.LongType)):
			summand = GaussInt(summand)  # Coerce if adding integer and GaussInt
		return self.add(summand)
	def __radd__(self,summand):  # Overload the "+" operator
		if ((type(summand) == types.IntType) or (type(summand) == types.LongType)):
			summand = GaussInt(summand)  # Coerce if adding integer and GaussInt
		return self.add(summand)
	def __iadd__(self,summand): # Overload the "+=" operator
		self = self + summand
		return self
	def __neg__(self):  # Overload the "-" unary operator 
		return GaussInt(-self.r,-self.i)
	def __sub__(self,summand):  # Overload the "-" binary operator 
		return self.__add__(-summand)
	def __isub__(self,summand): # Overload the "-=" operator
		self = self - summand
		return self
	def mult(self,multip):
		prod = GaussInt()
		prod.r = (self.r * multip.r) - (self.i * multip.i)
		prod.i = (self.i * multip.r) + (self.r * multip.i)
		return prod
	def __mul__(self,multip):  # Overload the "*" operator
		if ((type(multip) == types.IntType) or (type(multip) == types.LongType)):
			multip = GaussInt(multip)  # Coerce if multiplying integer and GaussInt
		return self.mult(multip)
	def __rmul__(self,multip):  # Overload the "*" operator
		if ((type(multip) == types.IntType) or (type(multip) == types.LongType)):
			multip = GaussInt(multip)  # Coerce if multiplying integer and GaussInt
		return self.mult(multip)
	def __imul__(self,multip): # Overload the "*=" operator
		self = self * multip
		return self
	def div(self,divisor):
		q = complex(self.r,self.i)/complex(divisor.r,divisor.i)
		return GaussInt(int(round(q.real)),int(round(q.imag)))
	def __div__(self,divisor):  # Overload the "/" operator 
		return self.div(divisor)
	def __idiv__(self,divisor): # Overload the "/=" operator
		self = self/divisor
		return self
	def mod(self,divisor):
		return self - divisor * (self/divisor)
	def __mod__(self,divisor):  # Overload the "%" operator 
		return self.mod(divisor)
	def __imod__(self,divisor): # Overload the "%=" operator
		self = self%divisor
		return self
	def divmod(self,divisor):
		q = self/divisor
		return q, self - divisor * q
	def xgcd(self,other):
		quot=GaussInt(); a1=GaussInt(1,0); b1=GaussInt(0,0); a2=GaussInt(0,0); 
		b2=GaussInt(1,0); a = self; b = other;
		if(b.norm() > a.norm()):
			a,b = b,a  # Swap a and b - need to start with a>b
			a1,b1,a2,b2 = a2,b2,a1,b1 # Swap (a1,b1) with (a2,b2)
		while (1):
			quot = a / b
			a %= b
			a1 -= quot*a2; b1 -= quot*b2
#			print a,b,a1,b1,a2,b2
			if (a == GaussInt(0,0)):
				return b, a2, b2
			quot = b / a
			b %= a
			a2 -= quot*a1; b2 -= quot*b1
#			print a,b,a1,b1,a2,b2
			if (b == GaussInt()):
				return a, a1, b1
	def gcd(self,other):
		return self.xgcd(other)[0]
	def powmod(self,exp,mod):
		accum = GaussInt(1,0); basepow2 = self; i=0
		while ((exp>>i)>0):
			if (((exp>>i) & 1) == 1):
				accum = (accum*basepow2) % mod
			basepow2 = (basepow2*basepow2) % mod
			i+=1
			print i,((exp>>(i-1)) & 1),basepow2,accum
		return accum
	def pow(self,exp):
		accum = GaussInt(1,0); basepow2 = self; i=0
		while ((exp>>i)>0):
			if (((exp>>i) & 1) == 1):
				accum = (accum*basepow2)
			basepow2 = (basepow2*basepow2)
			i+=1
			print i,((exp>>(i-1)) & 1),basepow2,accum
		return accum
	def __pow__(self,exponent):  # Overload the "**" operator 
		return self.pow(exponent)
#	def __coerce__(self,other):
#		return
GaussIntType = type(GaussInt(1,2))
	


# >>> from GaussInt import *
# >>> a = GaussInt(1,0)