from IPython.display import Math, display

from collections import defaultdict
from itertools import count

"""
Explanation: It was proposed by Christian Goldbach that every odd composite number can be written as the sum of a prime and twice a square.
$$
9 = 7 + 2 \cdot 1^2 \
15 = 7 + 2 \cdot 2^2 \
21 = 3 + 2 \cdot 3^2 \
25 = 7 + 2 \cdot 3^2 \
27 = 19 + 2 \cdot 2^2 \
33 = 31 + 2 \cdot 1^2
$$
It turns out that the conjecture was false.
What is the smallest odd composite that cannot be written as the sum of a prime and twice a square?
Version 1: Brute force
Formally stated, Goldbach's Other Conjecture says that all odd composite numbers $n$ can be expressed as 
$$
n = 2k^2 + p
$$
for some integer $k$ and prime number $p$. Let 
$$
S_n = { n - 2k^2 : k = 1, 2, \dotsc, \lfloor \sqrt{\frac{n}{2}} \rfloor }
$$
If any element $n - 2k^2$ of $S_n$ is prime, then we call $k$ a witness to Goldbach's Other Conjecture for $n$. We are asked to find the smallest $n$ that has no witness, i.e. such that all elements of $S_n$ are composite.
Let $P_n$ be the set of all prime numbers stricly smaller than $n$, then our algorithm searches for the smallest $n$ such that $P_n \cap S_n = \emptyset$, providing a counterexample to Goldbach's Other Conjecture.
<!-- TEASER_END -->
End of explanation
"""


factors = defaultdict(list)
witness = {}
primes = set()
for n in count(2):
    if factors[n]:
        # n is composite
        for m in factors.pop(n):
            factors[n+m].append(m)
        if n % 2:
            # n is odd and composite
            for k in range(1, int((n/2)**.5)+1):
                p = n - 2*k**2
                if p in primes: 
                    # TODO: `in` is O(len(primes))
                    # could optimize by using set
                    witness[n] = k
                    break
            if not n in witness:
                break
    else:
        factors[n*n].append(n)
        primes.add(n)

n

"""
Explanation: We modify and augment a very space-efficient implementation (due to David Eppstein, UC Irvine) of the Sieve of Erastothenes to generate composite numbers and keep track of all prime numbers below it. In the outermost-loop below, primes will always be the list of all prime number strictly below $n$, i.e. it is equivalent to $P_n$. Note that if $n$ is odd and composite, then the largest prime below it is no greater than $n-2$, since $n-1$ is even. Since $\max S_n = n - 2$ searching through $P_n$ is sufficient (i.e. if $n-2$ is prime, then $n-2 \in P_n$.) The loop terminates when we encounter an $n$ with no witnesses.
End of explanation
"""


lines = [r'{0} &= {1} + 2 \cdot {2}^2 \\'.format(n, n - 2*witness[n]**2, witness[n]) for n in sorted(witness)]

Math(r"""
     \begin{{align}}
     {body}
     \end{{align}}
     """.format(body='\n'.join(lines[:100])))

"""
Explanation: The answer is 5777.
Note that not only is this implementation space-efficient, it is also very time efficient, since we incrementally build our list of primes, composites and witnesses incrementally from bottom-up. If we hadn't augmented the implementation of the prime sieve, we would have had to use the prime sieve to obtain all odd composites, and perform a primality test on all elements of $S_n$, which would (usually) require a prime factorization algorithm, which in turn (usually) requires a prime sieve, not to mention the fact that there would be many overlapping subproblems, i.e. many $n < m$ such that $S_n \cap S_m \ne \emptyset$, so we'd have to use dynamic programming, by, for example memoizing the primality testing function or some other optimization - just a whole mess of redundancies and inefficiencies that we happily avoided with this method :)
Let's list out the witnesses to the first 100 numbers.
End of explanation
"""


