Zolotarev's lemma

zolotarev.py

"""
Code discussed in the video:

https://youtu.be/0A3mAUGFUow

* Motivation:
    - Relation between Legendre symbol and quadratic equations on ZxZ
        + See https://en.wikipedia.org/wiki/Legendre_symbol#Table_of_values

* For any fixed number, validate
    - the equivalence of the Legendre symbol L for different definitions
        + legendre_sigma
        + legendre_gauss
        + See /wiki/Legendre_symbol
    - the quadratic reciprocity law
        + See https://en.wikipedia.org/wiki/Quadratic_reciprocity
    - the Zolotarev lemma
        + See https://en.wikipedia.org/wiki/Zolotarev%27s_lemma
        + And see https://en.wikipedia.org/wiki/Signature_of_a_permutation

* Compute the Z[i]-prime (in the first quadrant), corresponding to any N-prime.
"""

def run_motivation():
    p = 17
    a = 5

    for x in range(0, 300):
        for y in range(-500, 0):
            eq_sol = x**2 + p * y == a  # (x,y) \in Z^2
            if eq_sol:
                print(f"The grid point (x,y)=({x}, {y}) fulfills    x·x + {p}·y == {a}")


################################################################################


def is_uneven_prime(n):
    return n > 2 and all(n % factor > 0 for factor in range(2, n))


def is_prime(n):
    return n == 2 or is_uneven_prime(n)


def Zi_norm(z):
    return z.real ** 2 + z.imag ** 2


def Zi_factor(p: int):
    # Notes:
    # * The second factor of p is the conjugate of the return value of this function.
    # * And indeed all the primes Z[i], up to multiplication by a unit, derive from a N-prime p in this way.
    assert is_prime(p)
    if p % 4 == 3:
        return p

    if p == 2:
        return complex(1, 1)

    for a in range(p):  # Search for root in first quadrant (a > b > 0)
        assert a ** 2 < p  # We don't actually need to loop up to 'a'
        for b in range(1, a):
            z = complex(a, b)
            if Zi_norm(z) == p:
                return z

    assert False, "Sanity check failure"  # All primes in N must correspond to a Z[i]-prime with norm p


def Zi_prime_factoization(z):
    """
    TODO (implementation):
    Compute the prime factorization of ||z|| and use Zi_factor to find the other factors.
    """
    assert False, "Not impmlemented yet"


def naturals():
    n = 0
    while True:
        yield n
        n += 1


def primes():
    for p in filter(is_prime, naturals()):
        yield p


def smallest_prime_factor(n):
    for m in range(2, n + 1):
        if n % m == 0:
            return m


def prime_factors(n: int):
    assert n >= 0  # Allow for n=0 and n=1 here
    ps = []
    if n < 2:
        ps = [n]
    while n > 1:
        p = smallest_prime_factor(n)
        ps.append(p)
        n = int(n / p)
    return ps


def function_range_ordered(f, domain: list): # Ordered set
    return [f(x) for x in domain]


def function_range(f, domain: set):
    return function_range_ordered(f, list(domain))


def is_inversion(permutation, i, j):
    return i < j and permutation[i] > permutation[j]


def get_inversions(permutation):
    n = len(permutation)
    return [
        (i, j)
        for j in range(1, n)
        for i in range(1, n)
        if is_inversion(permutation, i, j)
    ]


def num_inversions(permutation):
    return len(get_inversions(permutation))


def permuation_signature(permutation):
    return (-1) ** num_inversions(permutation)


def is_permutation(ints: list):
    return sorted(ints) == list(range(len(ints)))


def exponent_gauss_tmp_aux(p):
    return int((p - 1) / 2)


class Zn:
    def __init__(self, n):
        self.size = n
        self.range = set(range(n))
        self.__standard_representative = lambda x: x % n

    def left_action(self, g, x):
        f = lambda x: g * x
        return self._eval(f, x)

    def legendre_sigma(self, x):
        assert is_prime(self.size)
        return 0 if self._equals_zero(x) else (1 if self._has_root(x) else -1)

    def legendre_gauss(self, x):
        assert is_prime(self.size)
        f = lambda y: y ** exponent_gauss_tmp_aux(self.size)
        r = self._eval(f, x)
        assert r in (0, 1, self.__standard_representative(-1)), f"{self.size}, {r}"
        return r if r in (0, 1) else -1

    def left_coset(self, g):
        f = lambda x: self.left_action(g, x)
        return function_range(f, self.range)

    def left_coset_ordered(self, g):
        if g == 0:
            return {0}
        f = lambda x: self.left_action(g, x)
        permutations = function_range_ordered(f, self.range)
        assert (g == 0 and permutations == {0}) or is_permutation(permutations)  # action by g trivializes or is an automorphism
        return permutations

    def _equals_zero(self, x):
        return self.__standard_representative(x) == 0

    def _eval(self, f, x):
        return self.__standard_representative(f(x))

    def _is_fun_zero(self, f, x):
        return self._equals_zero(self._eval(f, x))

    def _get_zeros(self, f):
        return [x for x in self.range if self._is_fun_zero(f, x)]

    def _has_root(self, x):
        return self.__get_square_roots(x) != []

    def __get_square_roots(self, x):
        f = lambda y: self.left_action(y, y) - x
        return self._get_zeros(f)


def assert_theorems(n):
    zn = Zn(n)

    for a in zn.range:
        if a == 0:
            continue

        legendre = zn.legendre_sigma(a)
        reciprocity_law_rhs = (-1) ** (exponent_gauss_tmp_aux(a) * exponent_gauss_tmp_aux(n)) * legendre
        assert not is_uneven_prime(a) or Zn(a).legendre_sigma(n) == reciprocity_law_rhs
        assert zn.legendre_gauss(a) == legendre
        assert permuation_signature(zn.left_coset_ordered(a)) == legendre
        # Note: Check what the theorems say to see where all relations are really guaranteed to work out

    for g in zn.range:
        if g < 2:
            print(f"(n={n}) action of g={g} ...")
            continue  # Uninteresting
        permuation = zn.left_coset_ordered(g)
        n_i = num_inversions(permuation)
        sign = permuation_signature(permuation)
        perm = " ".join(map(str, permuation))
        print(f"(n={n}) action of g={g}, sign={sign}, inversions={n_i} ({round(100 * n_i * (n-1)**-2, 2)}%), permutation = [{perm}]")
    print("")


if __name__=="__main__":
    for p in primes():
        z = Zi_factor(p)
        if z.imag == 0:
            print(f"Zi-prime {p}.\n")
        else:
            print(f"N-prime {p} = {int(z.real)}^2 + {int(z.imag)}^2, which has Z[i]-prime factor z = {int(z.real)} + {int(z.imag)} i.\n")

        assert_theorems(p)
        if p > 20:
            break