Relative primality

factor.py

from math import *

class Factor:

	def __init__(self, primes):
		self.primes = primes
	
	def factorize(self, N):
		factors = {}
		lim = int(floor(sqrt(N)))

		if N == 1:
			return {1:1}

		for f in self.primes:

			if f > lim:
				break
			
			while N % f == 0:
				if not f in factors:
					factors[f] = 0
					 
				factors[f] += 1
				N /= f
		
		if N != 1 and N in self.primes:
			factors[N] = 1
		
		return factors

image.py

class Image:

	def __init__(self, width, height, default):
		self.width = width
		self.height = height
		self.data = [int(default) for i in range(0, width * height * 3)]
	
	def save(self, filename):
		fp = open(filename, "w+")
		fp.write("P6\n" + str(self.width) + " " + str(self.height) + "\n255 ")
		fp.write("".join(map(chr, self.data)))
		fp.close()
	
	def set_rgb(self, x, y, r, g, b):
		r = max(0, min(255, int(r)))
		g = max(0, min(255, int(g)))
		b = max(0, min(255, int(b)))
		i = (x + y * self.width) * 3

		self.data[i] = r
		self.data[i + 1] = g
		self.data[i + 2] = b
	
	def set_gray(self, x, y, g):
		self.set_rgb(x, y, g, g, g)

main.py

import time

from image import *
from primes import *
from factor import *
from pgcd import *

N = 2
S = 2048

def relative_primality(A, B):
	fA = f.factorize(A)
	fB = f.factorize(B)
	pgcd = PGCD(fA, fB)
	fpgcd = f.factorize(pgcd)

	return 2.0 * sum(fpgcd.values()) / (sum(fA.values()) + sum(fB.values()))

img = Image(S, S, 255)

if N < 2:
	n = 2

p = Primes(N + S)

time_start = time.time()
p.compute()
time_primes = (time.time() - time_start) * 1000

print time_primes

f = Factor(p.get_primes())

time_start = time.time()
for y in range(N, N + S):
	for x in range(N, y):
		p = relative_primality(x, y)
		p = int(floor(255 * p))
		img.set_gray(x - N, y - N, p)
		img.set_gray(y - N, x - N, p)
time_factors = time.time() - time_start

print time_factors

img.save("relprim.ppm")

pgcd.py

def PGCD(a, b):
	p = 1
	for f in a:
		if f in b:
			p *= f ** min([a[f], b[f]])
	return p

primes.py

from math import *

class Primes:

	def __init__(self, n):
		self.n = n
		self.numbers = [-1 for i in range(0, n + 1)]
		self.numbers[1] = False
		self.numbers[2] = True
		self.numbers[3] = True
		self.numbers[4] = False
		self.computed = False
		self.primes = [2, 3]
	
	def compute(self):
		for i in range(5, self.n + 1):

			if i % 2 == 0:
				self.numbers[i] = False
				continue

			elif sum(map(int, str(i).split())) % 3 == 0:
				self.numbers[i] = False
				continue

			lim = int(floor(sqrt(i)))

			is_prime = True
			for d in self.primes:
				if d > lim:
					break

				if i % d == 0:
					is_prime = False
					break
			
			self.numbers[i] = is_prime

			if is_prime:
				self.primes.append(i)
		
		self.computed = True

	def is_prime(self, n):
		if not self.computed:
			self.compute()
		
		if n > self.n:
			return -1
		
		return self.numbers[n]
	
	def get_primes(self):
		return self.primes