Skip to content

Instantly share code, notes, and snippets.

@pjt33

pjt33/mathse3015816.py

Last active Dec 5, 2018
Embed
What would you like to do?
MSE 3015816
#!/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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.