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