CoCalc -- Collaborative Calculation in the Cloud
SharedsumOfThreeCubes.ipynbOpen in CoCalc

Prime Decomposition of Sum of Three Cubes for 33

def gcd(m,n):
    k = m % n
    return n if k == 0 else gcd(n, k) 
x = 8866128975287528
y = 8778405442862239
z = 2736111468807040

x**3 - y**3 - z**3
33
def isPrime(n):
    i = 2
    primes = {}
    while i*i <= n:
        count = 0
        while n % i == 0:
            if i in primes:
                primes[i] = primes[i] + 1
            else:
                primes[i] = 1
            n = n//i
        i = i + 1
        
    if n != 1:
        primes[n] = 1
    return primes

This seems faster thant the following sieve

Just brute force output the list of primes that divide each number.

Xdecomp = isPrime(x)
Ydecomp = isPrime(y)
Zdecomp = isPrime(z)
from IPython.display import display, Math, Latex
Math(r"""\begin{align*}
     x &= %s \\
     y &= %s \\
     z &= %s
     \end{align*}"""%("\cdot".join(map(lambda n: "%d^%d"%(n,Xdecomp[n]), Xdecomp)), 
                "\cdot".join(map(lambda n: "%d^%d"%(n,Ydecomp[n]), Ydecomp)),
                "\cdot".join(map(lambda n: "%d^%d"%(n,Zdecomp[n]), Zdecomp))))
Xdecomp, Ydecomp, Zdecomp
({2: 3, 7: 1, 467: 1, 378289: 1, 896201: 1}, {8778405442862239: 1}, {2: 7, 5: 1, 89917: 1, 47545783: 1})

Prime sieve

This is the original code

def gen_primes():

    """ Generate an infinite sequence of prime numbers.
    """
    # Maps composites to primes witnessing their compositeness.
    # This is memory efficient, as the sieve is not "run forward"
    # indefinitely, but only as long as required by the current
    # number being tested.
    #
    D = {}
    
    # The running integer that's checked for primeness
    q = 2
    
    while True:
        if q not in D:
            # q is a new prime.
            # Yield it and mark its first multiple that isn't
            # already marked in previous iterations
            # 
            yield q,D
            
            D[q * q] = [q]
        else:
            # q is composite. D[q] is the list of primes that
            # divide it. Since we've reached q, we no longer
            # need it in the map, but we'll mark the next 
            # multiples of its witnesses to prepare for larger
            # numbers
            # 
            for p in D[q]:
                D.setdefault(p + q, []).append(p)
            del D[q]
        
        q += 1

Prime sieve

This is a variant of the above that only looks for prime factors of x, if x = 0 it runs as the above.

def gen_primes(x = 0):
    """ Generates list of prime divisors of x
    """
    # Maps composites to primes witnessing their compositeness.
    # This is memory efficient, as the sieve is not "run forward"
    # indefinitely, but only as long as required by the current
    # number being tested.
    #
    D = {}
    
    # The running integer that's checked for primeness
    q = 2
    
    while q <= x or x == 0: #x = 0 gives infinite list of primes 
        if q not in D:
            # q is a new prime.
            # Yield it and mark its first multiple that isn't
            # already marked in previous iterations
            # 
            if x == 0:
                yield q
            else:
                while x % q == 0:
                    x = x//q
                    yield q 
            
            D[q * q] = [q]
        else:
            # q is composite. D[q] is the list of primes that
            # divide it. Since we've reached q, we no longer
            # need it in the map, but we'll mark the next 
            # multiples of its witnesses to prepare for larger
            # numbers
            # 
            for p in D[q]:
                D.setdefault(p + q, []).append(p)
            del D[q]
        
        q += 1
g = gen_primes()

E = {'x':[], 'y':[], 'z':[]}
while q < max([x, y, z]):
    q = next(g)
    if x % q == 0:
         while x % q == 0:
            x = x//q
            E['x'].append(q)
            print('x', E['x'])
        
    if y % q == 0:
         while y % q == 0:
            x = y//q
            E['y'].append(q)
            print('y', E['y'])
    
    if z % q == 0:
         while z % q == 0:
            z = z//q
            E['z'].append(q)
            print('z', E['z'])
Xdecomp = E['x']
Ydecomp = E['y']
Zdecomp = E['z']

The following is just for a check

from functools import reduce

def prod(a):
    return reduce(lambda x, y: x * y, a, 1)
prod(Xdecomp)
g = gen_primes()
g = gen_primes()
h = (next(g) for i in range(1000))
m = map(isPrime, list(h))
a ={}
a['a']='a'
a.setdefault('a','b')
'a'
a
{'a': 'a'}