9

我正在尝试编写一个 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

我还没有研究过数据结构和算法，所以我很感激任何关于使这个脚本更高效的指针。

4

4 回答 4

4

您的解决方案会检查每个可能的整数以查看它是否可能是解决方案。仅检查实际上是幂的整数以查看它们是否为有效答案会更有效-因为它们的数量较少。这是可以做到这一点的东西。但是找到 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)
于 2015-09-14T01:47:30.247 回答
4

这是打破分析器并查看您的代码将所有时间花在哪里的好时机。为此，我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变成一个字符串，然后把每个数字变成一个 int，然后把它相加。

这将是程序的缓慢部分，这很有意义。你不仅做了很多，而且是一个缓慢的操作：在这一行，你做一个字符串解析操作，然后在字符串中的每个 char 上映射一个 int 转换函数，然后对它们求和。

我敢打赌，如果你可以在不先进行字符串转换的情况下计算数字的总和，它会快很多。

让我们试试看。我做了一些其他的改变，比如一开始就删除了多余的列表推导。这是我得到的：

#!/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()

有什么cProfile要说的？

⌁ [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 个。

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

于 2015-09-14T02:08:34.830 回答
3

更新：我发现这是一个 Project Euler 问题（#119），我发现基本上已经记录了相同的解决方案：http: //www.mathblog.dk/project-euler-119-sum-of-digits-升功率/

我不确定我是否过度简化，但只是迭代一系列数字的幂似乎很快。你不能保证订单，所以计算比你需要的更多，然后排序并获得前 30 个。我不能证明我已经得到了它们，但我已经测试base了 500 和exp50 并且它返回相同top 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
于 2015-09-14T04:17:19.267 回答
-1

[编辑：由于我正在转录的给定算法中的错误，这种方法相当慢（大约与其他方法一样快）。我将其保留在这里以供代码参考。然而，这似乎是在不诉诸数论技巧的情况下可以做到的最好的方法。]

在计算整数序列时，你应该先去 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 http://stackoverflow.com/a/4429063/711085
    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 名候选人，这可能需要几分钟。

于 2015-09-14T03:20:49.843 回答