快速求素数和

记录一下第一回答

定义 $S(v,p)$ 为 $2$ 到 $v$ 所有整数中，在普通筛法中外层循环筛完 $p$ 时仍然幸存的数的和。因此这些数要不本身是素数，要不其最小的素因子也大于 $p$ 。因此我们需要求的是 $S(n,\lfloor\sqrt n\rfloor)$ ，其中 $n$ 是十亿。

为了计算 $S(v,p)$ ，先考虑几个特殊情况

$p\le1$ 。此时所有数都还没有被筛掉，所以 $S(v,p)=\sum_{i=2}^{v}i=\frac{(2+v)(v-1)}{2}$
$p$ 不是素数。因为筛法中 $p$ 早已被别的数筛掉，所以在这步什么都不会做，所以此时 $S(v,p)=S(v,p-1)$ ；
$p$ 是素数，但是 $v<p^2$ 。因为每个合数都一定有一个不超过其平方根的素因子，如果筛到 $p$ 时还没筛掉一个数，那么筛到 $p-1$ 时这个数也还在。所以此时也有 $S(v,p)=S(v,p-1)$ 。

现在考虑最后一种稍微麻烦些的情况： $p$ 是素数，且 $p^2\le v$ 。
此时，我们要用素数 $p$ 去筛掉剩下的那些数中 $p$ 的倍数。注意到现在还剩下的合数都没有小于 $p$ 的素因子。因此有：
$S(v,p)=S(v,p-1)-\sum_{\substack{2\le k \le v,\\ p\mbox{为}k\mbox{的最小素因子}}}k$

后面那项中提取公共因子 $p$ ，有：
$S(v,p)=S(v,p-1)-p\times\sum_{\substack{2\le k \le v,\\ p\mbox{为}k\mbox{的最小素因子}}}\frac{k}{p}$

因为 $p$ 整除 $k$ ，稍微变形一下，令 $t=\frac{k}{p}$ ，有：
$S(v,p)=S(v,p-1)-p\times\sum_{\substack{2\le t \le \lfloor\frac{v}{p}\rfloor,\\ t\mbox{的最小素因子}\ge p}}t$

注意到一开始提到的 $S$ 的定义（“这些数要不本身是素数，要不其最小的素因子也大于(注意!) $p$ ”），此时 $p$ 后面这项可以用 $S$ 来表达：
$S(v,p)=S(v,p-1)-p\times(S(\left\lfloor\frac{v}{p}\right\rfloor,p-1)-\{p-1\mbox{以内的所有素数和\}})$

再用 $S$ 替换素数和得到最终表达式：
$S(v,p)=S(v,p-1)-p\times(S(\left\lfloor\frac{v}{p}\right\rfloor,p-1)-S(p-1,p-1))$

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

 1 import time
 2 def P10(n):
 3     r = int(n ** 0.5)
 4     # assert r*r <= n and (r+1)**2 > n
 5     V = [n // i for i in range(1, r + 1)]
 6     #    print V
 7     V += list(range(V[-1] - 1, 0, -1))
 8     #print V
 9     S = {i: i * (i + 1) // 2 - 1 for i in V}
10     #print S
11     st = time.clock();
12     for p in range(2, r + 1):
13         if S[p] > S[p - 1]:  # p is prime
14             sp = S[p - 1]  # sum of primes smaller than p
15             p2 = p * p
16             for v in V:
17                 if v < p2: break
18                 S[v] -= p * (S[v // p] - sp)
19     end = time.clock();
20     print end - st
21     return S[n]
22 
23 while(True):
24     N = input()
25 
26     print P10(N)

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

————————————————————————更新——————————————————————————————

之后用C++实现了一遍，发现速度还没线性筛快，比较循环部分，发现python 的速度要比 C++快50+倍，如果去除 list 的影响，dict 与 map的效率可能相差两个数量级，去查找dict 与 map 的实现原理，原来， dict 是用 hash 复杂度O（1）而 map为了元素的有序性，采用红黑树复杂度为O（logn）。

 1 const int maxn = 1e6;
 2 map<ll, ll> mp;
 3 vector<ll> vr;
 4 ll arr[maxn];
 5 ll solve(ll n){
 6     ll r = sqrt(n);
 7     vr.clear();
 8     mp.clear();
 9     for (int i = 1; i <= r; ++i){
10         vr.pk(n / i);
11     }
12     for (int i = vr.back() - 1; i > 0; --i){
13         vr.pk(i);
14     }
15     int len = vr.size();
16     for (int i = 0; i < len; ++i){
17         arr[i] = vr[i];
18         mp[vr[i]] = vr[i] * (vr[i] + 1) / 2 - 1;
19     }
20     ll sp;
21     int p2;
22     double st = clock();
23     for (int i = 2; i <= r; ++i){
24         if (mp[i] > mp[i - 1]){
25             sp = mp[i - 1];
26             p2 = i * i;
27             for (int v = 0; v < len; ++v){
28                 if (arr[v] < p2) break;
29                 mp[arr[v]] -= i * (mp[arr[v] / i] - sp);
30             }
31         }
32     }
33     double ed = clock();
34     cout << ed - st << endl;
35     return mp[n];
36 }
37 
38 int main(){
39     //freopen("data.out", "w", stdout);
40     //freopen("data.in", "r", stdin);
41     //cin.sync_with_stdio(false);
42     int n;
43     while (cin >> n){
44         cout << solve(n) << endl;
45     }
46     return 0;
47 }

posted @ 2015-04-20 00:11 yeahpeng 阅读(1138) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

yeahpeng

快速求素数和

公告