快速求素数和-白红宇

快速求素数和

阅读量：4686 次

发布时间：2019-06-09

本文共 2350 字，大约阅读时间需要 7 分钟。

在知乎上看到一个问题：

记录一下第一回答

定义 $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
     
       mp; 3 vector
      
        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 }

转载于:https://www.cnblogs.com/yeahpeng/p/4440481.html

你可能感兴趣的文章