N == N 的数字和的某些幂(运行 太慢了)

Certain Power of Sum of Digits of N == N (running too slowly)

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

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

您的解决方案会检查每个可能的整数,看看它是否是一个解决方案。只检查实际上是幂的整数以查看它们是否是有效答案会更有效——因为这样的整数更少。这是可以做到这一点的东西。但是找到 30 个可能需要 6 秒以上的时间——它们很快就会变得稀缺。

EDIT -- 更新以执行 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)

现在是打开探查器并查看您的代码将所有时间花在哪里的好时机。为此,我对您的代码进行了一些 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,即然后求和。

这将是程序中较慢的部分,这很有意义。你不仅做了很多,而且这是一个缓慢的操作:在那一行上,你做了一个字符串解析操作,然后在字符串中的每个字符上映射一个 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。

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

[编辑:由于我正在转录的给定算法中存在错误,此方法相当慢(与其他方法一样快)。我将其保留在这里以供代码参考。然而,这似乎是不求助于数论技巧的最佳方法。]

在计算整数序列时,您应该先去Sloane's 输入序列。这是序列 A023106 "a(n) is a power of the sum of its digits.". The first 32 such numbers up to 68719476736 can be found by clicking on the 'list' link. Often one can find algorithms (which may or may not be efficient) also given, as well as references. There is also linked the first such 1137 numbers up to [some number large enough to fill up a few paragraphs]

现在关于一个有效的算法:除非你有某种方法跳过 范围 的数字而不看它们,或者除非我们可以利用一些数学 属性数字,你被 O(N) 算法困住了。您可以采用的另一种方法是尝试计算所有幂(这可以让您跳过所有数字),然后测试每个幂 P=n^m 以查看 "is there a number whose digits sum to some power of P (or whatever)".

其实这个算法在上面link中已经提供给我们了。上面link给出的算法是(在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 
    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 个候选人,这可能需要几分钟。

更新:我发现这是一个欧拉项目问题(#119),我发现基本上相同的解决方案已经记录在案:http://www.mathblog.dk/project-euler-119-sum-of-digits-raised-power/

我不确定我是否过度简化了,但只是迭代一系列数字的幂似乎很快。你不能保证顺序,所以计算的比你需要的多,然后排序并取前 30 名。我无法证明我已经得到了所有但我已经测试了 base 高达 500 和 exp 最多 50 并且它 returns 相同的前 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