Prime Power Triples
May 22, 2020
The easy brute force calculation uses three loops; the trick is that some numbers can be generated in more than one way, so they must be collected in a set instead of simply being counted:
(define (euler87 n)
(length (unique = (sort < (list-of s
(a in (primes (iroot 2 n)))
(b in (primes (iroot 3 n)))
(c in (primes (iroot 4 n)))
(s is (+ (* a a) (* b b b) (* c c c c)))
(< s n))))))
We collect solutions in a list, then sort, cast out duplicates, and count. Here’s an example:
> (euler87 50) 4
You can run the program at https://ideone.com/B8JkdH.
Nifty little drill. Here is my take on it using Julia 1.4.1: https://pastebin.com/DVjwfbAx
The original approach I took was way too sluggish, so I ended up using a slicker one. There is also a verbose option for viewing the specific combinations: main(n, true), where n is the maximum integer to explore (in this case 50 000 000). The default for this parameter is false.
Glad the blog is back! Cheers
from sympy import primerange def count(N): return len(set(i**2 + j**3 + k**4 for i in primerange(1, int(N**(1/2))+1) for j in primerange(1, int((N-i**2)**(1/3))+1) for k in primerange(1, int((N-i**2-j**3)**(1/4))+1)))Easy to solve by brute force, once you realise that there are not that many prime powers to look at.
Welcome back! Hope all is well with everybody.
Not sure I like those multiple calls to primerange etc. Here’s another Python solution that generates a single list of prime numbers. Takes a couple of seconds for the Euler problem (here we just count up to a million):
# Recursive sieve of Eratosthenes # Squares the size of the max prime each time def rsieve(primes): N = primes[-1]**2+1 sieve = [0]*N for p in primes: for q in range(p*p,N,p): sieve[q] = 1 return [i for i in range(2,N) if sieve[i] == 0] N = 1000000 primes = [2] while primes[-1]**2 < N: primes = rsieve(primes) s = set() for p in primes: ptotal = p**2 if ptotal > N: break for q in primes: qtotal = ptotal+q**3 if qtotal > N: break for r in primes: rtotal = qtotal+r**4 if rtotal > N: break s.add(rtotal) print(len(s))Very glad to hear you are well and that things are settling down somewhat.
Here’s my solution in Python, using a standard Sieve of Eratosthenes and precalculating the required prime powers. Runs in about half a second on my desktop.
#! /usr/bin/env python3 # Euler 87 # https://projecteuler.net/problem=87 # Also Programming Praxis # https://programmingpraxis.com/2020/05/22/prime-power-triples/ # Find count of numbers n < 50e6 where n = p^2 + q^3 + r^4; p, q, r all prime # Limits: # ceil (N ** 1/2) = 7072 # ceil (N ** 1/3) = 369 # ceil (N ** 1/4) = 85 import time import primes start = time.time() # Eratosthenes sieve primes up to largest limit ps = primes.sieve(7072) # Precalculate prime powers possible p2s = [p*p for p in ps] p3s = [p*p*p for p in ps if p <= 369] p4s = [p*p*p*p for p in ps if p <= 85] # add and collect unique answers sums = set() lim = 50_000_000 for p2 in p2s: for p3 in p3s: for p4 in p4s: tot = p2 + p3 + p4 if tot < lim: sums.add(tot) end = time.time() taken = end - start print(len(sums)) print(f'Calculated in {taken:.3} seconds')A solution in (not necessarily well-written) Icon. Also, note that prime-power triples aren’t unique; for example, 210 = 11^2 + 2^3 + 3^4 = 2^2 + 5^3 + 3^4. Code that doesn’t handle duplicates is producing incorrect counts.
Here’s a solution in Python.
import itertools def is_prime(n): x = 2 while x * x <= n: if n % x == 0: return False x += 1 + (x & 1) return n > 1 def prime_gen(): return (x for x in itertools.count(2) if is_prime(x)) pred = lambda x: x < 50_000_000 p2 = itertools.takewhile(pred, (x ** 2 for x in prime_gen())) p3 = itertools.takewhile(pred, (x ** 3 for x in prime_gen())) p4 = itertools.takewhile(pred, (x ** 4 for x in prime_gen())) pp_triples = set() for a, b, c in itertools.product(p2, p3, p4): sum = a + b + c if sum < 50_000_000: pp_triples.add(sum) print(len(pp_triples))Here’s a refactored version of my last solution.
import itertools def is_prime(n): x = 2 while x * x <= n: if n % x == 0: return False x += 1 + (x & 1) return n > 1 def prime_gen(pow, limit): return itertools.takewhile( lambda x: x < limit, (x ** pow for x in itertools.count(2) if is_prime(x))) N = 50_000_000 p2 = prime_gen(2, N) p3 = prime_gen(3, N) p4 = prime_gen(4, N) pp_triples = (a + b + c for a, b, c in itertools.product(p2, p3, p4)) print(len({x for x in pp_triples if x < N}))Here’s another refactor to make it shorter, at the cost of readability.
import itertools def is_prime(n): x = 2 while x * x <= n: if n % x == 0: return False x += 1 + (x & 1) return n > 1 def prime_gen(pow, limit): return itertools.takewhile( lambda x: x < limit, (x ** pow for x in itertools.count(2) if is_prime(x))) N = 50_000_000 pp_triples = (sum(x) for x in itertools.product(*(prime_gen(x, N) for x in (2, 3, 4)))) print(len({x for x in pp_triples if x < N}))