"""goldbach.py"""

import sys
from numpy import log2, ceil, maximum, array
from sympy import isprime, primepi
from matplotlib import pyplot as plt


def p1p2(n):
    """Returns minimal primes (p1, p2) summing to the even integer n >= 4"""
    p1 = 2
    while p1 < n:
        p2 = n - p1
        if isprime(p1) and isprime(p2):
            return (p1, p2)
        p1 += 1

    raise ValueError("Goldbach is false! ;-) Counterexample: %d" % n)


def makeplot(N):
    fourtoN = range(4, N + 1, 2)
    ps = [p1p2(n) for n in fourtoN]

    # First make the "encode everything" plot
    m1 = maximum.accumulate([ceil(log2(primepi(p1))) for p1, p2 in ps])
    m2 = maximum.accumulate([ceil(log2(primepi(p2))) for p1, p2 in ps])
    goldbach = m1 + m2
    conventional = [ceil(log2((n - 2) / 2)) for n in fourtoN]

    plt.plot(fourtoN, goldbach, label="Goldbach-Encoding")
    plt.plot(fourtoN, conventional, label="Conventional Encoding")
    plt.xlabel("Even integers (N)")
    plt.ylabel("Encoding length (bits)")
    plt.title("Bits to encode all even integers 4 <= n <= N")
    plt.legend(loc="lower right")
    plt.show()

    # Next make the "encode one integer at a time" plot
    log2p1 = array([log2(primepi(p1)) for p1, p2 in ps])
    log2p2 = array([log2(primepi(p2)) for p1, p2 in ps])
    conventional = [log2((n - 2) / 2) for n in fourtoN]

    plt.plot(fourtoN, log2p1 + log2p2, label="Goldbach")
    plt.plot(fourtoN, conventional, label="Ordinary")
    plt.plot(fourtoN, log2p1, label="log2(pi(p1))")
    plt.plot(fourtoN, log2p2, label="log2(pi(p2))")
    plt.xlabel("Even integers (n)")
    plt.ylabel("Encoding length (bits)")
    plt.legend(loc="upper left")
    plt.show()


if __name__ == "__main__":
    N = int(sys.argv[1]) if len(sys.argv) == 2 else 1000
    makeplot(N)