import numpy as np from math import isqrt import time upto = 2**22 # 2**32 needs at least 32 Gb Ram t1 = time.time() A072103 = [] for m in range(2,isqrt(upto)+1): k = 2 while m**k < upto: A072103.append(m**k) k += 1 A001597 = sorted(set(A072103)) # eliminates multiples and sorting A362424 = np.zeros(upto+1, dtype="i4") A001597 = [1] + A001597 # we need a "1" in front of A001597 a = 0 while A001597[a] < upto // 2: b = a + 1 while b < len(A001597) and A001597[a] + A001597[b] < upto: A362424[A001597[a] + A001597[b]] += 1 b += 1 a += 1 A365295 = np.zeros(max(A362424)+1, dtype=object) for it in range (1, len(A362424)): if A365295[A362424[it]] == 0: A365295[A362424[it]] = it n = 0 while n < len(A365295) and A365295[n] != 0: n += 1 # trimming, because of possible zeros between A365295 = A365295[0:n] print(list(A365295)) t2 = time.time() print(t2-t1)