快速求素数和

标签：

在知乎上看到一个问题：求十亿内所有质数的和，怎么做最快?

记录一下第一回答

定义为到所有整数中，在普通筛法中外层循环筛完时仍然幸存的数的和。因此这些数要不本身是素数，要不其最小的素因子也大于。因此我们需要求的是，其中是十亿。

为了计算，先考虑几个特殊情况

1. 。此时所有数都还没有被筛掉，所以
2. 不是素数。因为筛法中早已被别的数筛掉，所以在这步什么都不会做，所以此时；
3. 是素数，但是。因为每个合数都一定有一个不超过其平方根的素因子，如果筛到时还没筛掉一个数，那么筛到时这个数也还在。所以此时也有。

现在考虑最后一种稍微麻烦些的情况：是素数，且。
此时，我们要用素数去筛掉剩下的那些数中的倍数。注意到现在还剩下的合数都没有小于的素因子。因此有：


后面那项中提取公共因子，有：


因为整除，稍微变形一下，令，有：


注意到一开始提到的的定义（“这些数要不本身是素数，要不其最小的素因子也大于(注意!)”），此时后面这项可以用来表达：


再用替换素数和得到最终表达式：



我们最终的结果是。计算时可以使用记忆化，也可以直接自底向上动态规划。
至于算法的复杂度就留作练习，是低于以上任何一种暴力方法的。

def P10(n):
    r = int(n ** 0.5)
    # assert r*r <= n and (r+1)**2 > n
    V = [n // i for i in range(1, r + 1)]
    #    print V
    V += list(range(V[-1] - 1, 0, -1))
    #    print V
    S = {i: i * (i + 1) // 2 - 1 for i in V}
    #    print S
    for p in range(2, r + 1):
        if S[p] > S[p - 1]:  # p is prime
            sp = S[p - 1]  # sum of primes smaller than p
            p2 = p * p
            for v in V:
                if v < p2: break
                S[v] -= p * (S[v // p] - sp)
    return S[n]


N = 1000000000
print P10(N)

View Code

1e9的数据能在1s能跑出来， 真是神一样的算法， 神一样的代码。

快速求素数和

标签：

原文地址：http://www.cnblogs.com/yeahpeng/p/4440481.html