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.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
9 Comments »

9 Responses to “Prime Power Triples”

  1. Zack said
    May 22, 2020 at 12:37 PM

    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

  2. Jan Van lent said
    May 22, 2020 at 3:28 PM 
    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)))
  3. Richard A. O'Keefe said
    May 26, 2020 at 10:45 AM

    Easy to solve by brute force, once you realise that there are not that many prime powers to look at.

  4. matthew said
    May 27, 2020 at 8:12 AM

    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))
  5. Alex B said
    May 28, 2020 at 1:36 PM

    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')
  6. r. clayton said
    June 7, 2020 at 2:26 AM

    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.

  7. Daniel said
    July 16, 2020 at 4:52 AM

    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))
  8. Daniel said
    July 16, 2020 at 5:01 AM

    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}))
  9. Daniel said
    July 16, 2020 at 5:06 AM

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

Leave a comment