Open in CoCalc
In [ ]:
def gcd(m,n): k = m % n return n if k == 0 else gcd(n, k)
In [2]:
x = 8866128975287528 y = 8778405442862239 z = 2736111468807040 x**3 - y**3 - z**3
33
In [3]:
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.

In [4]:
Xdecomp = isPrime(x) Ydecomp = isPrime(y) Zdecomp = isPrime(z)
In [34]:
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))))
In [22]:
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

In [ ]:
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.

In [93]:
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
In [ ]:
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'])
In [ ]:
Xdecomp = E['x'] Ydecomp = E['y'] Zdecomp = E['z']

The following is just for a check

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