如果你将其实现为一个递归生成器，你将不需要存储大量的值（甚至不需要存储结果）：

def powersum(n,k,b=1):
    bk = b**k
    while bk < n:
        for bases in powersum(n-bk,k,b+1):
            yield [b]+bases
        b += 1
        bk = b**k
    if bk == n :
        yield [b]
        
print(*powersum(100,2)) # [1, 3, 4, 5, 7] [6, 8] [10]
print(*powersum(100,3)) # [1, 2, 3, 4]
print(sum(1 for _ in powersum(1000,2))) # 1269 solutions
print(sum(1 for _ in powersum(2000,2))) # 27526 solutions (6 seconds)

• 注意，这仍然具有指数时间复杂度，因此对于n的稍大值，它将慢得多 *

print(sum(1 for _ in powersum(2200,2))) # 44930 solutions (12 seconds)
print(sum(1 for _ in powersum(2500,2))) # 91021 solutions (25 seconds)
print(sum(1 for _ in powersum(2800,2))) # 175625 solutions (55 seconds)
print(sum(1 for _ in powersum(3100,2))) # 325067 solutions (110 seconds)

你的问题是你在subsets中保存了太多不相关的组合，我在这里写了一些代码，只保存相关的列表，所以应该可以做到这一点，（我在这里写了多行，这样更容易理解）

def my_solutions(n,k):
    subsets = []
    kth_root = int(n**(1/k))# kth_root of n is the upper bound
    relevant_list = range(1,kth_root+1)
    for comb_size in range(1,kth_root+1):
        for comb in combinations(relevant_list, comb_size):
            sum_comb = sum([x**k for x in comb])
            if (sum_comb == n):
                subsets.append(list(comb))
    return subsets

下面是相同的代码，但用一行程序编写，以最大限度地提高效率

def my_solutions_one_linner(n,k):
    subsets = []
    kth_root = int(n**(1/k))# kth_root of n is the upprt bound
    relevant_list = range(1,kth_root+1)
    ans =  [[list(comb) for comb in combinations(relevant_list, comb_size) if sum([x**k for x in comb]) == n] for comb_size in range(1,kth_root+1)]
    return [lst for lst in ans if lst!=[]]

@Reggie-florade和@Ohad-sharet，可以让代码运行得更快。请看下面的代码。我的解决方案与@Ohad-sharet非常相似，但运行速度要快得多。我的函数名是“kth_power_and_n”。

from itertools import combinations
import time
import math
import pandas as pd

def my_solutions_one_linner(n,k):
    subsets = []
    kth_root = int(n**(1/k))# kth_root of n is the upprt bound
    relevant_list = range(1,kth_root+1)
    ans =  [[list(comb) for comb in combinations(relevant_list, comb_size) if sum([x**k for x in comb]) == n] for comb_size in range(1,kth_root+1)]
    return [lst for lst in ans if lst!=[]]  

def kth_power_and_n(n, k):
    UpperLimit = int(math.floor(n**(1/k)))
    Matrix = list(map(lambda x: x**k,range(1,UpperLimit+1)))
    Final_Solution_List = []
    for i in range(1,len(Matrix)+1):
        Generated_Combinations = combinations(Matrix,i)
        for eachCombination in Generated_Combinations:
            if sum(eachCombination) == n:
                Final_Solution_List.append(sorted(list(map(lambda x: int(round(x**(1/k))), eachCombination))))
    return Final_Solution_List

Solution_Dict = {"n":[],"k":[],"time taken 1":[],"time taken 2":[],"Solution Found 1":[],"Solution Found 2":[]}

                
for n in [50,100,250,500]:
    for k in [2,3,4]:
        ## Previous Solution time taken
        Start = time.time()
        Solution = my_solutions_one_linner(n,k)
        Time_Taken = time.time()-Start
        Solution_Dict["n"].append(n)
        Solution_Dict["k"].append(k)
        Solution_Dict["time taken 1"].append(Time_Taken)        
        if Solution:
            Solution_Dict["Solution Found 1"].append(True)
        else:
            Solution_Dict["Solution Found 1"].append(False)

        ## Significantly faster code
        Start = time.time()
        Solution = kth_power_and_n(n,k)
        Time_Taken = time.time()-Start
        Solution_Dict["time taken 2"].append(Time_Taken)        
        if Solution:
            Solution_Dict["Solution Found 2"].append(True)
        else:
            Solution_Dict["Solution Found 2"].append(False)

DF_One_Liner = pd.DataFrame(data=Solution_Dict)
DF_One_Liner["Time Gain Ratio"] = DF_One_Liner["time taken 1"]/DF_One_Liner["time taken 2"]
print(DF_One_Liner)

下面是两种方法的时间比较。

n  k  time taken 1  time taken 2  Solution Found 1  Solution Found 2  Time Gain Ratio
0    50  2      0.000170      0.000047              True              True         3.627551
1    50  3      0.000012      0.000007             False             False         1.758621
2    50  4      0.000007      0.000005             False             False         1.333333
3   100  2      0.001503      0.000185              True              True         8.142119
4   100  3      0.000019      0.000011              True              True         1.723404
5   100  4      0.000010      0.000006             False             False         1.833333
6   250  2      0.067290      0.006028              True              True        11.163476
7   250  3      0.000073      0.000018             False             False         4.135135
8   250  4      0.000010      0.000006             False             False         1.720000
9   500  2     11.933419      0.870651              True              True        13.706314
10  500  3      0.000185      0.000036             False             False         5.132450
11  500  4      0.000019      0.000009             False             False         2.135135

我的解决方案是基于这个问题Finding all possible combinations of numbers to reach a given sum的答案。链接问题中的目标是通过添加列表numbers中的元素来达到值target。在我们的例子中，numbers将是从1到target**(1/k)的整数提升到k。
这是另一个答案的函数。唯一的修改是我们知道numbers是按升序排列的。这意味着我们可以更早地停止迭代。如果numbers[0]太大，则以下幂的其余部分也会太大。

def subset_sum(numbers: list[int], target: int, partial: list[int]=[], partial_sum: int=0):
    if partial_sum == target:
        yield partial
        return
    if numbers and (partial_sum + numbers[0]) > target:
        return
    for i, n in enumerate(numbers):
        remaining = numbers[i + 1:]
        yield from subset_sum(remaining, target, partial + [n], partial_sum + n)

现在让我们正确地调用这个函数。这里重要的一点是subset_sum需要幂，但我们的解决方案将使用这些数字而不增加它们。要做到这一点，让我们定义一个字典，其中键是幂，值是答案中的数字。在此之后，我们只需要给予subset_sum一个包含字典键的列表。

def solutions(n, k):
    kth_root = math.ceil(n**(1/k))
    powers_dict = {i**k: i for i in range(1, kth_root+1)}
    powers = list(powers_dict.keys())
    return [[powers_dict[i] for i in answer] for answer in subset_sum(powers, n)]

在我的机器上，这段代码花了9秒才完成solutions(2000, 2)。

好吧，只是为了好玩，我试着写一个递归函数来做这个;下面是结果（power_sum（2000，2）在我的机器上大约10秒内就能得到答案，所以尽管它可能不是最优的，但它仍然不算太坏）。

from math import inf
def power_sums(total, p, max_allowed=inf):
    nmax = min(int(total**(1/p)), max_allowed)
    for i in range(nmax, 0, -1):
        if (ip := i**p) == total:
            yield [i]
        else:
            for x in power_sums(total - ip, p, i-1):
                yield [i] + x