python - N == N の桁数の合計の特定のべき乗 (実行が遅すぎる)

Question

N の数字の合計の特定の累乗が N に等しいすべての整数 (N) を見つける Python スクリプトを作成しようとしています。たとえば、8 + 1 = 9 であるため、N=81 が該当します。ある 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

2 番目の解決策では、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

私はまだデータ構造とアルゴリズムを勉強していないので、このスクリプトをより効率的にするためのヒントをいただければ幸いです。

Answer 1

これは、プロファイラーを分割して、コードがどこで時間を費やしているかを確認する絶好の機会です。そのために、コードの周りに小さな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.

これがプログラムの遅い部分であることは非常に理にかなっています。多くのことを行うだけでなく、操作が遅くなります。その 1 行で、文字列の解析操作を実行してから、文字列内の各文字に対して 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 かかります。

プロファイラーはあなたの友達です。文字列変換は、何十万回も行うと無料ではありません:)

Answer 2

あなたの解は、考えられるすべての整数を調べて、それが解になる可能性があることを確認します。実際に累乗である整数のみを調べて、それらが有効な答えであるかどうかを確認する方が効率的です。それらの数が少ないためです。これを行うものがあります。しかし、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)

Answer 3

更新: これが Project Euler の問題 (#119) であることを発見し、基本的に同じ解決策が既に文書化されていることを発見しました: http://www.mathblog.dk/project-euler-119-sum-of-digits-パワーアップ/

単純化しすぎているかどうかはわかりませんが、数値の範囲の累乗を反復するだけでかなり速いようです。順序を保証することはできないので、必要以上に計算してから、並べ替えて上位 30 を取得します。すべてを取得したことを証明することはできませんが、最大base500 とexp最大 50 をテストし、同じ結果を返します。 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

Answer 4

[編集: 私が転記していた特定のアルゴリズムのバグにより、この方法はかなり遅い (他の方法とほぼ同じ速さ)。コード参照用にこれをここに保持します。ただし、これは数論のトリックに頼らずに実行できる最善の方法のようです。]

整数シーケンスを計算するときは、最初に Sloane's に移動してシーケンスを入力する必要があります。これは、シーケンスA023106 「a(n) はその桁の合計の累乗です」です。. 「リスト」リンクをクリックすると、68719476736 までの最初の 32 個の番号を見つけることができます。多くの場合、アルゴリズム (効率的である場合とそうでない場合があります) と参考文献を見つけることができます。[いくつかの段落を埋めるのに十分な大きさの数] までの最初の 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 の候補をチェックする必要があり、数分かかる場合があります。