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