一下子,我不认为有一个巧妙的方法使这真正有效,但它很容易使它更快。如果将示例视为矩阵,则一次将它们相加一行。这需要对每个i求出i的所有除数并求它们的立方。总之,这需要与x**2成比例的操作数。在
您可以很容易地将其切割成与x成比例的操作,方法是用列求和矩阵。给定一个整数j,那么{}中有多少个整数可以被j整除?这很简单:在这个范围内有x//j的j的倍数,所以除数j贡献了j**3 * (x // j)的总和。在def better(x):
return sum(j**3 * (x // j) for j in range(1, x+1))
它运行得更快,但仍然需要与x成比例的时间。在
有一些低级的技巧可以反过来通过常数因子来加快速度,但是它们总体上仍然需要O(x)时间。例如,注意x // j == 1表示所有j,这样x // 2 < j <= x。因此,可以跳过sum中约一半的项,而将其替换为连续多维数据集总和的闭式表达式:
^{pr2}$
better2()大约是better()的两倍,但是要想比O(x)快,就需要更深入的洞察力。在
更快
闲暇时想想,我还是没有一个真正聪明的主意。但我给出的最后一个想法可以得出一个合乎逻辑的结论:不要把范围只有一个倍数的除数组合在一起,也不要把范围内只有两个倍数的除数,以及三个、四个,以及。。。这导致下面的better3(),它执行了许多与x的平方根大致成比例的操作:def better3(x):
result = 0
for i in range(1, x+1):
q1 = x // i
# value i has q1 multiples in range
result += i**3 * q1
# which values have i multiples?
q2 = x // (i+1) + 1
assert x // q1 == i == x // q2
if i < q2:
result += i * (sum3(q1) - sum3(q2 - 1))
if i+1 >= q2: # this becomes true when i reaches roughly sqrt(x)
break
return result
当然O(sqrt(x))比原来的O(x**2)有了很大的改进,但是对于非常大的参数来说,它仍然是不切实际的。例如,better3(10**6)似乎可以立即完成,但是better3(10**12)需要几秒钟时间,而{}则是休息时间;-)
注意:我使用的是python3。如果您使用的是python2,请使用xrange(),而不是{}。在
再来一个
better4()具有与O(sqrt(x))相同的O(sqrt(x))时间行为,但以不同的顺序进行求和,这样可以简化代码并减少对sum3()的调用。对于“大”参数,它比我的框上的better3()快大约50%。在def better4(x):
result = 0
for i in range(1, x+1):
d = x // i
if d >= i:
# d is the largest divisor that appears `i` times, and
# all divisors less than `d` also appear at least that
# often. Account for one occurence of each.
result += sum3(d)
else:
i -= 1
lastd = x // i
# We already accounted for i occurrences of all divisors
# < lastd, and all occurrences of divisors >= lastd.
# Account for the rest.
result += sum(j**3 * (x // j - i)
for j in range(1, lastd))
break
return result
通过在"A Successive Approximation Algorithm for Computing the Divisor Summatory Function"中扩展算法,可能会做得更好。这需要O(cube_root(x))时间来处理除数相加这个可能更简单的问题。但它涉及的范围更广,我对这个问题的关心程度不够,我自己也没有去追究;—)
微妙
数学中有一个很容易被忽略的微妙之处,所以我要把它拼出来,但只涉及better4()。在
在d = x // i之后,注释声称d是出现i次的最大除数。但这是真的吗?实际出现d的次数是x // d,这是我们计算的。我们怎么知道x // d实际上等于i?在
这就是if d >= i:保护该评论的目的。在^{之后我们知道x == d*i + r
对于某个满足r的整数0 <= r < i。这基本上就是楼层划分的意思。但是既然d >= i也是已知的(这就是if测试所保证的),那么{}也一定是这样。这就是我们知道x // d是{}。在
当d >= i为而不是时,这可能会发生故障,这就是为什么需要使用不同的方法。例如,如果x == 500和i == 51,d(x // i)是9,但不是9是最大的除数出现51次。实际上,9出现了500 // 9 == 55次。而对于正实数d == x/i
当且仅当i == x/d
对于楼层划分来说,这是不可能的。但是,如上所述,如果我们也知道{},第一个确实意味着第二个。在
只是f或者有趣
better5()重写better4(),获得大约10%的速度增益。真正的教学重点是说明预先计算所有的循环极限是很容易的。上面奇数代码结构的一个要点是,它神奇地为0输入返回0,而不需要对此进行测试。better5()放弃这一点:def isqrt(n):
"Return floor(sqrt(n)) for int n > 0."
g = 1 << ((n.bit_length() + 1) >> 1)
d = n // g
while d < g:
g = (d + g) >> 1
d = n // g
return g
def better5(x):
assert x > 0
u = isqrt(x)
v = x // u
return (sum(map(sum3, (x // d for d in range(1, u+1)))) +
sum(x // i * i**3 for i in range(1, v)) -
u * sum3(v-1))