sbp · August 1, 2011 20:10
diff --git a/gaussian.py b/gaussian.py
 #!/usr/bin/env python

 import math

 def complex_eps(c): 
   return complex(round(c.real, 6), round(c.imag, 6))

 def floor(c): 
   real = math.floor(round(c.real, 6))
   imag = math.floor(round(c.imag, 6))
   return complex(real, imag)

 def complex_modulus(a, b): 
   return a - (b * floor(a / b))

 def gaussian_gcd(a, b): 
   while b != 0+0j: 
      a, b = b, complex_modulus(a, b)
   return a

 def modular_exponentiation(base, exponent, modulus):
   result = 1
   while exponent > 0:
      if (exponent & 1) == 1: 
         result = (result * base) % modulus
      exponent >>= 1
      base = (base ** 2) % modulus
   return result

 def root4(p): 
   # 4th root of 1 modulo p
   # Derived from http://bit.ly/p4ESdk by Robin Chapman
   base = 2
   pow = p/4
   while True:
      a = modular_exponentiation(base, pow, p)
      modulus = (a ** 2) % p
      if modulus == (p - 1): return a
      base += 1

 def sq2(p): 
   a = root4(p) # 4th root of 1 modulo p
   return gaussian_gcd(complex(p, 0), complex(a, 1))

 def first_quadrant_associate(c): 
   if c.real >= 0: 
      if c.imag >= 0: return c # 1st quadrant
      return c * 1j # 2nd quadrant
   elif c.imag < 0: return -c # 3rd quadrant
   return c * -1j # 4th quadrant

 def factor(n): 
   if n < 2: return []

   if not (n & 1): # faster than n % 2
      res = factor(n / 2)
      return [2] + res

   for i in xrange(3, int(n ** .5) + 1, 2): 
      if not (n % i): 
         res = factor(n / i)
         return [i] + res
   return [n]

 def norm(G): 
   return long(G.real ** 2) + long(G.imag ** 2)

 def factor_gaussian(G): 
  # Algorithm from http://bit.ly/pCn9Hm by Jim Lewis
   N = norm(G)
   primes = factor(N)
   if primes == [1]: 
      return [G], 1

   factors = []
   while primes: 
      p = primes[0]

      if p == 2: 
         u = 1 + 1j
         if not complex_modulus(G, u): q = u
         else: q = 1 - 1j
         primes.remove(p)
      elif (p % 4) == 3: 
         q = p
         primes.remove(p)
         primes.remove(p)
      else: 
         u = sq2(p)
         u = first_quadrant_associate(u) # u[0] + (u[1] * 1j))
         if not complex_modulus(G, u): q = u
         else: q = u.conjugate()
         primes.remove(p)

      factors.append(q)
      G = complex_eps(G/q)
   return factors, G

 def product(numbers): 
   result = 1
   for number in numbers: 
      result = result * number
   return result

 def divisor_cardinality(factors): 
   powers = []
   unique_factors = set(factors)
   for unique_factor in unique_factors: 
      power = factors.count(unique_factor)
      powers.append(power)
   return product(power + 1 for power in powers)

 def canvas(): 
   print '<!DOCTYPE html>'
   print '<title>Gaussian Plane</title>'
   print '<script>'
   print 'function draw() {'
   print '   var canvas = document.getElementById("plane")'
   print '   if (canvas.getContext) {'
   print '      var ctx = canvas.getContext("2d")'
   r = 256
   for y in xrange(0, r + 1): 
      for x in xrange(0, r + 1): 
         print '      // %s' % ([x, r - y])
         factors, G = factor_gaussian(complex(x, (r - y)))
         C = divisor_cardinality(factors)
         a = 255 - (C * 5)
         if C == 2: 
            b = 255 # 128 + 32
         else: b = a
         if a < 0: a = 0
         print '      ctx.fillStyle = "rgb(%s, %s, %s)"' % (b, b, a)
         print '      ctx.fillRect(%s, %s, 2, 2)' % (x * 2, y * 2)
   print '   }'
   print '}'
   print '</script>'
   print '<body onload="draw()">'
   print '<canvas id="plane" width="514" height="514">...</canvas>'

 if __name__ == '__main__': 
   canvas()
	#!/usr/bin/env python

	import math

	def complex_eps(c):
	return complex(round(c.real, 6), round(c.imag, 6))

	def floor(c):
	real = math.floor(round(c.real, 6))
	imag = math.floor(round(c.imag, 6))
	return complex(real, imag)

	def complex_modulus(a, b):
	return a - (b * floor(a / b))

	def gaussian_gcd(a, b):
	while b != 0+0j:
	a, b = b, complex_modulus(a, b)
	return a

	def modular_exponentiation(base, exponent, modulus):
	result = 1
	while exponent > 0:
	if (exponent & 1) == 1:
	result = (result * base) % modulus
	exponent >>= 1
	base = (base ** 2) % modulus
	return result

	def root4(p):
	# 4th root of 1 modulo p
	# Derived from http://bit.ly/p4ESdk by Robin Chapman
	base = 2
	pow = p/4
	while True:
	a = modular_exponentiation(base, pow, p)
	modulus = (a ** 2) % p
	if modulus == (p - 1): return a
	base += 1

	def sq2(p):
	a = root4(p) # 4th root of 1 modulo p
	return gaussian_gcd(complex(p, 0), complex(a, 1))

	def first_quadrant_associate(c):
	if c.real >= 0:
	if c.imag >= 0: return c # 1st quadrant
	return c * 1j # 2nd quadrant
	elif c.imag < 0: return -c # 3rd quadrant
	return c * -1j # 4th quadrant

	def factor(n):
	if n < 2: return []

	if not (n & 1): # faster than n % 2
	res = factor(n / 2)
	return [2] + res

	for i in xrange(3, int(n ** .5) + 1, 2):
	if not (n % i):
	res = factor(n / i)
	return [i] + res
	return [n]

	def norm(G):
	return long(G.real 2) + long(G.imag 2)

	def factor_gaussian(G):
	# Algorithm from http://bit.ly/pCn9Hm by Jim Lewis
	N = norm(G)
	primes = factor(N)
	if primes == [1]:
	return [G], 1

	factors = []
	while primes:
	p = primes[0]

	if p == 2:
	u = 1 + 1j
	if not complex_modulus(G, u): q = u
	else: q = 1 - 1j
	primes.remove(p)
	elif (p % 4) == 3:
	q = p
	primes.remove(p)
	primes.remove(p)
	else:
	u = sq2(p)
	u = first_quadrant_associate(u) # u[0] + (u[1] * 1j))
	if not complex_modulus(G, u): q = u
	else: q = u.conjugate()
	primes.remove(p)

	factors.append(q)
	G = complex_eps(G/q)
	return factors, G

	def product(numbers):
	result = 1
	for number in numbers:
	result = result * number
	return result

	def divisor_cardinality(factors):
	powers = []
	unique_factors = set(factors)
	for unique_factor in unique_factors:
	power = factors.count(unique_factor)
	powers.append(power)
	return product(power + 1 for power in powers)

	def canvas():
	print '<!DOCTYPE html>'
	print '<title>Gaussian Plane</title>'
	print '<script>'
	print 'function draw() {'
	print ' var canvas = document.getElementById("plane")'
	print ' if (canvas.getContext) {'
	print ' var ctx = canvas.getContext("2d")'
	r = 256
	for y in xrange(0, r + 1):
	for x in xrange(0, r + 1):
	print ' // %s' % ([x, r - y])
	factors, G = factor_gaussian(complex(x, (r - y)))
	C = divisor_cardinality(factors)
	a = 255 - (C * 5)
	if C == 2:
	b = 255 # 128 + 32
	else: b = a
	if a < 0: a = 0
	print ' ctx.fillStyle = "rgb(%s, %s, %s)"' % (b, b, a)
	print ' ctx.fillRect(%s, %s, 2, 2)' % (x * 2, y * 2)
	print ' }'
	print '}'
	print '</script>'
	print '<body onload="draw()">'
	print '<canvas id="plane" width="514" height="514">...</canvas>'

	if __name__ == '__main__':
	canvas()