There are seven prime numbers that meet the conditions. I need to find them all, but the program is just too slow, is there a better way to do it?
#edit Just found out that those numbers are called Narcissistic or Armstrong numbers
import math
def prime(n):
for i in range(2,int(math.sqrt(n))+1):
if n%i == 0: return False
return True
def sum_powers(n):
num_str = str(n)
sum = 0
for i in num_str:
sum += int(i)**len(num_str)
return sum
def findprimes(desired_amount):
amount = 0
number = 2
while amount < desired_amount:
if prime(number):
if number == sum_powers(number):
print(number)
amount += 1
number += 1
else:
number += 1
else:
number += 1
findprimes(7)
Is brute force possible?
First, it's clear that if you want to find all solutions (and be sure you have all possible solutions), generating primes will not work. Let's see why.
We want to find some upper bound for numbers we need to check. For a given number of digits n, the largest n-digit number is all 9's. This number also maximizes the sum over all digits d of d^n, with a sum of n*(9^n). So for the sum of digit powers to be as large as the original numbers, we at least need n*(9^n) >= 10^(n-1). This makes it clear that there's a finite bound, too, as the right side is growing exponentially faster than the left.
Using wolfram alpha, the largest zero of f(n) = n*(9^n) - 10^(n - 1) is at around n=60.8, so we can be sure we're done after checking all numbers with up to 60 digits. This is far, far larger than any consecutive list of primes ever computed (which appears to be all primes below around 10^20), so the naive approach isn't even brute-forceable by supercomputers.
Depth-first search over all digit multisets
However, using the strategy mentioned in my comment, there's a better way. We can just generate all of the ways to pick a multiset of any quantity of digits 0-9 up to 60: from combinatorics, there's (60 + 9 choose 9) total options, which is around 6 * 10^10. This is closer to being brute-forceable, especially with some clever filtering. To just find the seven answers, it turns out we only need to check up to 23 digits long.
Once you have a combination of m digits, you can add each digit to the m'th power and see whether you get the original digits back, and whether the sum-string has the same multiplicities of digits as your combination.
The full code is below. I've included my own Miller-Rabin primality test, but using a numerical library like gmpy2 for that will be faster. On my computer, it prints out the seven primes in around a minute. I compared this against @Alain T's solution, which is conceptually similar (but using a Miller-Rabin prime tester for both), and this runs about 3x faster despite not using their clever even/odd filtering.
import collections
from time import perf_counter
import math
from typing import List, Tuple
def main():
found_primes = 0
def solve(my_nums: List[Tuple[int, int]], actual_sum: int) -> None:
if not is_prime_miller(actual_sum):
return
sum_counts = collections.Counter(str(actual_sum))
for orig_dig, amount_to_use in my_nums:
if amount_to_use != sum_counts[str(orig_dig)]:
return
nonlocal found_primes
found_primes += 1
print(f"Found ans number {found_primes}: {actual_sum}")
def dfs(curr_digit: int, remaining: int, curr_sum: int, orig_num_digs: int, biggest: int,
used_dig_list: List[Tuple[int, int]]) -> None:
my_pow = POWS[curr_digit]
if curr_sum >= biggest or 10 * (remaining * my_pow + curr_sum) < biggest - 10:
return
if remaining == 0:
if len(str(curr_sum)) == orig_num_digs:
solve(used_dig_list, curr_sum)
return
assert remaining > 0
if curr_digit == 0:
# curr_sum += remaining*pow(9, orig_num_digs)
if len(str(curr_sum)) == orig_num_digs:
solve(used_dig_list + [(0, remaining)], curr_sum)
return
new_sum = curr_sum
for amount_to_use in range(remaining + 1):
dfs(curr_digit=curr_digit - 1, remaining=remaining - amount_to_use, curr_sum=new_sum,
orig_num_digs=orig_num_digs, biggest=biggest,
used_dig_list=used_dig_list + [(curr_digit, amount_to_use)])
new_sum += my_pow
if new_sum >= biggest:
return
for digits_to_use in range(1, 24):
print(f"Starting {digits_to_use}, {perf_counter() - start_time :.4f} sec. since start")
POWS = [pow(x, digits_to_use) for x in range(10)]
dfs(curr_digit=9, remaining=digits_to_use, curr_sum=0, orig_num_digs=digits_to_use,
biggest=pow(10, digits_to_use), used_dig_list=[])
def is_prime_miller(n: int) -> bool:
small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
if n < 2:
return False
for p in small_primes:
if n < p * p:
return True
if n % p == 0:
return False
r, s = 0, n - 1
while s % 2 == 0:
r += 1
s //= 2
for a in small_primes:
x = pow(a, s, n)
if x == 1 or x == n - 1:
continue
for _ in range(r - 1):
x = pow(x, 2, n)
if x == n - 1:
break
else:
return False
return True
if __name__ == '__main__':
start_time = perf_counter()
main()
which gives this output:
2
3
5
7
28116440335967
449177399146038697307
35452590104031691935943
Edit: I've corrected an earlier claim that 10^33 was a bound on these kinds of numbers, which the post mentions are called 'narcissistic numbers': the bound I derived should be 10^60, although there are indeed only 7 such primes, and the largest composite number of this kind has fewer than 40 digits.