乘以数字长度的所有数字分量的总和必须等于相同的数字

Question

符合条件的素数有七个。我需要全部找到，但是程序太慢了，有没有更好的方法？

#edit 刚刚发现那些数字被称为自恋或阿姆斯特朗数字

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

你可能会跳过很多数字，因为除 2 之外的所有素数都是奇数。这意味着您的候选号码必须有奇数个奇数位。此外，候选数字中的数字顺序对得出的幂和并不重要，因此您应该从不同的数字组合开始，并检查结果是否产生相同的数字（在任何排列中）。

我相信通过这种方法，您应该能够对符合上述所有条件的候选人使用素数检查算法，而不是从素数生成器开始。

这是我对这种方法的看法：

isPrime(N) 是一个简单的质数检查函数
oddDigits(N) 是一个生成器，它将生成 N 个数字的不同组合，这些数字的总和为奇数。
match(digits)（也是一个生成器），采用数字组合并确定幂的总和是否达到由相同数字（以任何顺序）组成的数字。这是一个生成器，因为可能有不止一个幂来检查数字组合何时包含零（可能是前导）。
主循环遍历数字组合（最多 20 位的数字只有 5,007,002）并打印出匹配的数字。

...

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)

实际的 7 个素数一定很大，因为这只能找到其中的 4 个：

3
5
7
28116440335967

希望这能让你走上正轨

Answer 2

暴力破解可行吗？

首先，很明显，如果你想找到所有的解（并确保你有所有可能的解），生成素数是行不通的。让我们看看为什么。

我们想为需要检查的数字找到一些上限。对于给定的位数 n，最大的 n 位数是全 9。该数字还最大化了 d^n 的所有数字 d 的总和，总和为 n*(9^n)。所以要使数次幂之和和原数一样大，我们至少需要n*(9^n) >= 10^(n-1)。这清楚地表明也存在有限界限，因为右侧的增长速度比左侧快得多。

使用 wolfram alpha，f(n) = n*(9^n) - 10^(n - 1) 的最大零在 n=60.8 左右，因此我们可以确定在检查了最多 60 位的所有数字后我们就完成了。这 远远大于 比任何计算过的连续素数列表（似乎都是低于 10^20 左右的素数），所以天真的方法甚至不能被暴力破解超级计算机。

对所有数字多重集进行深度优先搜索

但是，使用我评论中提到的策略，还有更好的方法。我们可以生成所有方法来选择任意数量的数字 0-9 到 60 的多重集：从组合学来看，总共有（60 + 9 选择 9）个选项，大约 6 * 10^10。这更接近于暴力破解，尤其是通过一些巧妙的过滤。仅仅为了找到这七个答案，原来我们只需要检查最多 23 位数字。

一旦你有了 m 个数字的组合，你可以将每个数字的 m 次方相加，看看你是否能得到原来的数字，以及 sum-string 是否具有相同的重数数字作为您的组合。

完整代码如下。我已经包含了我自己的 Miller-Rabin 素数测试，但是使用像 gmpy2 这样的数字库会更快。在我的电脑上，它在大约一分钟内打印出七个素数。我将其与 @Alain T 的解决方案进行了比较，后者在概念上是相似的（但两者都使用 Miller-Rabin 素数测试仪），尽管没有使用他们聪明的 even/odd 过滤，但它的运行速度快了大约 3 倍。

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

给出了这个输出：

2
3
5
7
28116440335967
449177399146038697307
35452590104031691935943

编辑：我已经更正了之前的说法，即 10^33 是这些数字的界限，post 提到的称为 'narcissistic numbers'：我得出的界限应该是10^60，虽然这样的质数确实只有7个，而且这种最大的合数也不足40位。

乘以数字长度的所有数字分量的总和必须等于相同的数字

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

python

optimization

primes