MSE 3015816

mathse3015816.py

#!/usr/bin/python3

import math


def phi(n):
	# Naïve prime factorisation
	ps = set()
	m = n
	while m % 2 == 0:
		ps.add(2)
		m //= 2
	for p in range(3, int(math.sqrt(m)) + 1, 2):
		while m % p == 0:
			ps.add(p)
			m //= p
		if p > m:
			break
	if m > 1:
		ps.add(m)

	totient = n
	for p in ps:
		totient = totient * (p - 1) // p
	return totient


def ord(m, x):
	if m == 1:
		return 1

	xpow = 1
	for a in range(1, m + 1):
		xpow = (xpow * x) % m
		if xpow == 1:
			return a

	raise ValueError("ord: arguments %d,%d are not coprime" % (m, x))


def A(characteristic, root_set_size, degree):
	"""Compute a generating function for the building blocks of polynomials of interest and return it as a list.

	Arguments:
	characteristic: the characteristic of the field
	root_set_size: a function which takes integer n and gives the size of the sets of nth primitive roots of unity, c_K(n) * \rho_K(n)
	degree: the degree at which to truncate the generating function
	"""
	a = [0] * (degree + 1)
	# Exceptional factor: (x - 0)
	a[1] += 1
	# Primitive odd-power roots of unity
	# Since root_set_size(n) is a product of ord(n, 2), we want odd n for which ord(n, 2) <= degree.
	# That means that they're factors of (2^i - 1) for i <= degree.
	ns = set([1])
	for k in range(2, degree + 1):
		target = 2**k - 1
		for odd in range(3, target + 1, 2):
			if target % odd == 0:
				ns.add(odd)

	for n in ns:
		if characteristic > 0 and n % characteristic == 0:
			continue

		count = root_set_size(n)
		# Sanity check
		totient = phi(n)
		if totient % count != 0:
			raise ValueError("root_set_size(%d) is not a factor of the totient" % n)
		if count <= degree:
			a[count] += totient // count
	return a


def count_polynomials(p, root_set_size, degree):
	"""Compute a generating function for the polynomials of interest and return it as a list. See A."""
	b = [0] * (degree + 1)
	# General solution: Euler transform of building blocks
	a = A(p, root_set_size, degree)
	b[0] += 1
	for n in range(1, degree + 1):
		nbn = sum(b[n - k] * sum(d * a[d] for d in range(1, k + 1) if k % d == 0) for k in range(1, n + 1))
		if nbn % n != 0:
			raise Exception("Euler transform code is buggy")
		b[n] = nbn // n

	# Exceptional solution: P(x) = 0
	b[0] += 1
	return b


def root_set_size_GF(p, a):
	q = p**a
	def inner(n):
		k = ord(n, q)
		m = ord(n, 2)
		ii = set(q ** i % n for i in range(k))
		return k * min(j for j in range(1, m + 1) if 2 ** j % n in ii)
	return inner


if __name__ == "__main__":
	d = 12

	print("C")
	print(A(0, lambda n: ord(n, 2), d))
	print(count_polynomials(0, lambda n: ord(n, 2), d))
	print("")

	print("Q")
	print(A(0, phi, d))
	print(count_polynomials(0, phi, d))
	print("")

	print("R")
	def root_set_size_R(n):
		if n == 1:
			return 1

		pow2 = 1
		for a in range(1, n + 1):
			pow2 = (pow2 * 2) % n
			if pow2 == 1 or pow2 == n - 1:
				return 2 * a
		raise ValueError("Even n")
	print(A(0, root_set_size_R, d))
	print(count_polynomials(0, root_set_size_R, d))
	print("")

	for p in [3, 5, 7, 11, 13]:
		print("=========================================")
		print("Characteristic", p)
		for a in range(1, 100):
			print(a, count_polynomials(p, root_set_size_GF(p, a), d))