Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

Prime Decomposition of Sum of Three Cubes for 33

Project: Algorithms
Views: 260
Kernel: Python 3 (system-wide)
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))))
x=2371467137828918962011y=87784054428622391z=2751899171475457831\begin{align*} x &= 2^3\cdot7^1\cdot467^1\cdot378289^1\cdot896201^1 \\ y &= 8778405442862239^1 \\ z &= 2^7\cdot5^1\cdot89917^1\cdot47545783^1 \end{align*}
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'}