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)


### Resources¶

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'}