Sum of all number components powered to length of the number must equal the same number

Question

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)

Answer 1

You could probably skip a lot of numbers based on the fact that all primes, except 2 are odd numbers. This means that your candidate numbers must have an odd number of odd digits. Also, the order of digits in your candidate numbers will not matter for the resulting sum of powers so you should start with distinct combinations of digits and check if the result produces the same digits (in any permutations).

I believe that with this approach you should be able to use a prime checking algorithm on candidates that fit all of the above rather than start from a primes generator.

Here is my take on the approach:

• isPrime(N) is a simple prime check function
• oddDigits(N) is a generator that will produce distinct combinations of N digits that sum up to an odd number.
• match(digits) (also a generator), takes a combination of digits and determines if the sum of powers comes up to a number composed of the same digits (in any order). This is a generator because there could be more than one power to check when the combination of digits contain zeroes (that could be leading).
• The main loop goest through the digit combinations (there are only 5,007,002 for numbers up 20 digits) and prints out the matching numbers.

...

def isPrime(N):
    if N<3: return N==2
    return all(N%f for f in range(3,int(N**0.5)+1,2))

def oddDigits(size,minDigit=0,digitSum=0):
    for digit in range(minDigit,10):
        if size == 1:
            if (digitSum+digit)%2: yield [digit]
            continue
        for rest in oddDigits(size-1,digit,digitSum+digit):
            yield [digit]+rest
            
def match(digits):
    minLength = len(digits)-digits.count(0)
    for power in range(minLength,len(digits)+1):
        sp = sum(d**power for d in digits)
        if sorted(map(int,str(sp)))!=digits[-power:]:
            continue
        if isPrime(sp):
            yield sp
    
for digits in oddDigits(20):
    for N in match(digits):
        print(N)

The actual 7 primes must be quite large because this only finds 4 of them:

3
5
7
28116440335967

Hopefully, this can set you on the right track

Answer 2

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.