数字N的各位数字之和的某个幂等于N（运行速度太慢）

Question

我正在尝试编写一个Python脚本，查找所有整数(N)，使得数字N的各位数之和的某个次幂等于N。例如，N=81就符合条件，因为8 + 1 = 9，并且9的某个次幂(即2)等于81。

我选择的范围是任意的。我的脚本可以工作，但非常非常慢。理想情况下，我希望在大约6000毫秒内找到前30个这样的整数。

我的第一个解决方案:

def powerOfSum1():
    listOfN = []
    arange = [a for a in range(11, 1000000)] #range of potential Ns
    prange = [a for a in range(2, 6)] # a range for the powers to calculate
    for num in arange:
        sumOfDigits = sum(map(int, str(num)))
        powersOfSum = [sumOfDigits**p for p in prange]
        if num in powersOfSum:
            listOfN.append(num)
    return listOfN

在我的第二个解决方案中，我尝试存储每个sumOfDigits的所有幂次，但这并没有显著提高性能。

def powerOfSum2():
    listOfN = []
    powers= {}
    for num in range(11, 1000000):
        sumOfDigits = sum(map(int, str(num)))
        summ = str(sumOfDigits)
        if summ in powers:
            if num in powers[summ]:
                listOfN.append(num)
        else:
            powersOfSum = [sumOfDigits**p for p in range(2,6)]
            powers[summ] = powersOfSum
            if num in powers[summ]:
                listOfN.append(num)
    return listOfN

我还没有学习数据结构和算法，因此希望能够得到任何关于如何使这个脚本更加高效的指导。

- 5u2ie

---

请修正您代码的缩进。 - thefourtheye

arange = range(11, 1000000) 的作用类似于 arange = [a for a in range(11, 1000000)]。 - Jose Ricardo Bustos M.

@5u2ie 你使用的是哪个版本的Python？ - Anand S Kumar

选项2在我的电脑上需要2.1秒（使用Py3；Py2大约需要3秒）。这不是远远在您的6000毫秒限制范围内吗？ - paxdiablo

在计算整数序列时，应首先访问 Sloane's 并输入该序列。这是序列 [A023106 "a(n) is a power of the sum of its digits."] (https://oeis.org/A023106)。可以通过单击“列表”链接找到最多达 32 个这样的数字，最高可达 68719476736。通常还会提供算法（可能有效或无效），以及参考文献。 - ninjagecko

Answer 1

您的解决方案检查每个可能的整数，以确定它是否是一个解。只检查实际上是幂的整数是否为有效答案更有效，因为这样的整数较少。以下是一个可以实现此功能的代码示例。但是，要找到30可能需要比6秒更长的时间，因为这样的数字很快就会变得稀缺。

编辑 - 根据Padraic Cunningham和fjarri在评论中提出的更快的数字求和建议进行了更新，然后更新了一些其他调整，使其成为生成器，并使其适用于Python 3。

它仍然很快变慢，但该算法可能是可并行化的 - 也许您可以将数字求和放在单独的进程中。

编辑2 - 通过快速检查基数和结果值在模9意义下是否相等，节省了一些时间。可能还有其他数论技巧...

import heapq
import functools


def get_powers():
    heap = []
    push = functools.partial(heapq.heappush, heap)
    pop = functools.partial(heapq.heappop, heap)
    nextbase = 3
    nextbasesquared = nextbase ** 2
    push((2**2, 2, 2))
    while 1:
        value, base, power = pop()
        if base % 9 == value % 9:
            r = 0
            n = value
            while n:
                r, n = r + n % 10, n // 10
            if r == base:
                yield value, base, power
        power += 1
        value *= base
        push((value, base, power))
        if value > nextbasesquared:
            push((nextbasesquared, nextbase, 2))
            nextbase += 1
            nextbasesquared = nextbase ** 2


for i, info in enumerate(get_powers(), 1):
    print(i, info)

- Patrick Maupin

Answer 2

现在正是运用性能分析器查看代码耗时情况的绝佳时机。为了达到这一目的，我给你的代码添加了一个小cProfiler包装器：

#!/usr/bin/env python

import cProfile

import math


def powerOfSum1():
    listOfN = []
    arange = [a for a in range(11, 1000000)] #range of potential Ns
    prange = [a for a in range(2, 6)] # a range for the powers to calculate
    for num in arange:
        sumOfDigits = sum(map(int, str(num)))
        powersOfSum = [sumOfDigits**p for p in prange]
        if num in powersOfSum:
            listOfN.append(num)
    return listOfN


def main():
    cProfile.run('powerOfSum1()')

if __name__ == "__main__":
    main()

运行后，我得到了以下结果：

⌁ [alex:/tmp] 44s $ python powers.py
         1999993 function calls in 4.089 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.005    0.005    4.089    4.089 <string>:1(<module>)
        1    0.934    0.934    4.084    4.084 powers.py:7(powerOfSum1)
   999989    2.954    0.000    2.954    0.000 {map}
       10    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        2    0.017    0.009    0.017    0.009 {range}
   999989    0.178    0.000    0.178    0.000 {sum}

如果你看一下代码，似乎大部分时间都花在了map调用上，其中你将num转换为字符串，然后将每个数字转换为整数，最后进行求和。

这应该是程序的瓶颈所在。不仅你要频繁地执行这个操作，而且它本身就很慢：在这一行代码中，你进行了一个字符串解析操作，然后在字符串的每个字符上映射了一个整数转换函数，并将它们相加。

我敢打赌，如果你能在不先进行字符串转换的情况下计算出数字的总和，程序会快得多。

让我们试试吧。我还做了其他一些更改，例如删除开头多余的列表推导式。以下是我的修改结果：

#!/usr/bin/env python

#started at 47.56

import cProfile

import math

MAXNUM = 10000000

powersOf10 = [10 ** n for n in range(0, int(math.log10(MAXNUM)))]

def powerOfSum1():
    listOfN = []
    arange = range(11, MAXNUM) #range of potential Ns
    prange = range(2, 6) # a range for the powers to calculate
    for num in arange:
        sumOfDigits = 0
        for p in powersOf10:
            sumOfDigits += num / p % 10
        powersOfSum = []
        curr = sumOfDigits
        for p in prange:
            curr = curr * sumOfDigits
            if num < curr:
                break
            if num == curr:
                listOfN.append(num)
    return listOfN

def main():
    cProfile.run('powerOfSum1()')

if __name__ == "__main__":
    main()

是什么?它有什么作用?

⌁ [alex:/tmp] 3m42s $ python powers.py
         15 function calls in 0.959 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.006    0.006    0.959    0.959 <string>:1(<module>)
        1    0.936    0.936    0.953    0.953 powers.py:13(powerOfSum1)
       10    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        2    0.017    0.009    0.017    0.009 {range}

从4秒降到0.9秒？好多了。

如果你想真正看到效果，请将arange的上限额外增加一个零。在我的电脑上，原始代码需要约47秒，修改后的代码只需约10秒。

分析器是您的朋友，当您进行数十万次转换时，字符串转换并不是免费的 :)

- Alex

---

了解分析器的知识是至关重要的，而这是一个很好的教程。但另一方面，现在使用它可能还为时过早。我还没有运行过分析器，而在我的机器上的版本中，相同的前11个数字在不到1/100秒的时间内就显示出来了... - Patrick Maupin

在append()之后使用break可以避免多余的循环。 - AChampion

Answer 3

更新：我发现这是一个欧拉计划问题（#119），并且已经有基本相同的解决方案记录在这里：http://www.mathblog.dk/project-euler-119-sum-of-digits-raised-power/ 我不确定是否过于简化，但只需迭代一系列数字的幂似乎很快。您不能保证顺序，因此可以计算比所需更多的内容，然后进行排序并取前30个。我无法证明我已经得到了所有结果，但我已经测试了base高达500和exp高达50，并返回相同的前30个结果。应该注意OP仅测试指数高达5，这显着限制了结果数量。

def powerOfSum():
    listOfN = []
    for base in range(2, 100):
        num = base
        for _ in range(2, 10):
            num *= base
            if sum(map(int, str(num))) == base:
                listOfN.append(num)
    return sorted(listOfN)[:30]
powerOfSum()

输出

[81,
 512,
 2401,
 4913,
 5832,
 17576,
 19683,
 234256,
 390625,
 614656,
 1679616,
 17210368,
 34012224,
 52521875,
 60466176,
 205962976,
 612220032,
 8303765625,
 10460353203,
 24794911296,
 27512614111,
 52523350144,
 68719476736,
 271818611107,
 1174711139837,
 2207984167552,
 6722988818432,
 20047612231936,
 72301961339136,
 248155780267521]

在其中运行timeit（包括排序）后，我得到了以下结果：

100 loops, best of 3: 1.37 ms per loop

- AChampion

---

如果没有生成它们的顺序要求，或者事先不知道您是否拥有正确的常量进行迭代的要求，那么这种方法非常快速。 - Patrick Maupin

是的，如果您知道常量，速度很快，但即使有500、50，它只需要214毫秒就能找到197个结果。我真的很喜欢你的堆排序，用于按顺序幂，但它变慢得太快了（中间20多岁）。看起来很奇怪，与愚蠢的循环相比，它要慢得多。 - AChampion

@achampion，这是一个非常优雅的解决方案。然而，当我执行您的代码（即使使用常量500、50）时，最终排序列表中缺少第11和第12个元素（17210368和34012224），这让我感到疑惑。此外，我只得到了149个结果，有什么想法吗？ - 5u2ie

哎呀，我试图通过添加一个break语句来进行优化，但已经修复了...现在应该可以正常工作了。它们是同一底数的二次幂。 - AChampion

Answer 4

[编辑：由于我转录的算法存在错误，因此该方法相当慢（与其他方法一样快）。我将保留此处作为代码参考。然而，如果不诉诸于数论技巧，这似乎是最好的方法。]

在计算整数序列时，您应该首先访问 Sloane's 并输入该序列。这是序列 A023106 "a(n) 是其数字和的幂次."。点击“列表”链接可以找到前32个这样的数字，最高可达 68719476736。通常情况下，您还可以找到给出的算法（可能有效，也可能无效），以及参考文献。同时，还有前1137个这样的数字，直到[足以填满几段文字的某个大数字]。

关于高效算法：除非你有一种方法可以跳过一些数字范围而不查看它们，或者除非我们可以利用数字的某些数学属性，否则你只能使用O（N）算法。另一种方法是尝试计算所有幂次（这样可以跳过所有数字），然后测试每个幂次P=n^m，以查看“是否存在一个数字，其数字总和为P的某个幂次（或其他内容）”。

实际上，在上面的链接中已经为我们提供了此算法。上述链接中给出的算法是（在Mathematica中）：

fQ[n_] := Block[
  {b = Plus @@ IntegerDigits[n]}, 
  If[b > 1, IntegerQ[ Log[b, n]] ]
]; 
Take[
  Select[ 
    Union[ Flatten[ 
      Table[n^m, {n, 55}, {m, 9}]
    ]], fQ[ # ] &], 
  31
]
(* Robert G. Wilson v, Jan 28 2005 *)

松散的翻译如下：

def test(n):
     digitSum = sum of digits of n
     return n is a power of digitSum
candidates = set(n^m for n in range(55) for m in range(9))
matches = [c for c in candidates if test(c)]

完整的实现如下：

from math import *  # because math should never be in a module

def digitSum(n):
    return sum(int(x) for x in str(n))

def isPowerOf(a,b):
    # using log WILL FAIL due to floating-point errors
    # e.g. log_3{5832} = 3.0000..04
    if b<=1:
        return False
    # using https://dev59.com/AG855IYBdhLWcg3wSiQ7#4429063
    while a%b==0:
        a = a / b
    return a==1

def test(n):
    return isPowerOf(n, digitSum(n))

M = 723019613391360  # max number to check
candidates = set(n**m for n in xrange(int(sqrt(M)+1)) 
                       for m in xrange(int(log(M,max(n,2)))+1))
result = list(sorted([c for c in candidates if test(c)]))

输出：

>>> result
[2, 3, 4, 5, 6, 7, 8, 9, 81, 512, 2401, 4913, 5832, 17576, 19683, 234256, 390625, 614656, 1679616, 17210368, 34012224, 52521875, 60466176, 205962976, 612220032, 8303765625, 10460353203, 24794911296, 27512614111, 52523350144, 68719476736, 271818611107, 1174711139837, 2207984167552, 6722988818432, 20047612231936, 72301961339136, 248155780267521]

不幸的是，这需要相当长的时间。例如，在上面的例子中，我们必须检查53,863,062个候选项，这可能需要几分钟。

- ninjagecko

---

啊，是的，我正要发布一篇编辑文章，因为我同时发现了转录代码中的这个错误（修复方法是检查所有候选项，“candidates = set（...）”替换了不正确的“in range（55）... range（9）”）。将此答案保留在此处供其他人参考。 - ninjagecko