快速求素数和

时间 2019-11-15

标签快速素数繁體版

原文原文链接

记录一下第一回答算法

定义为到全部整数中，在普通筛法中外层循环筛完时仍然幸存的数的和。所以这些数要不自己是素数，要不其最小的素因子也大于。所以咱们须要求的是 $S(n,\lfloor\sqrt n\rfloor)$ ，其中是十亿。

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

$p\le1$ 。此时全部数都尚未被筛掉，因此 $S(v,p)=\sum_{i=2}^{v}i=\frac{(2+v)(v-1)}{2}$
不是素数。由于筛法中早已被别的数筛掉，因此在这步什么都不会作，因此此时；
是素数，可是。由于每一个合数都必定有一个不超过其平方根的素因子，若是筛到时还没筛掉一个数，那么筛到时这个数也还在。因此此时也有。

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

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

由于整除，稍微变形一下，令 $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(v,p)=S(v,p-1)-p\times(S(\left\lfloor\frac{v}{p}\right\rfloor,p-1)-\{p-1\mbox{之内的全部素数和\}})$

再用替换素数和获得最终表达式：
$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)$ 。计算时可使用记忆化，也能够直接自底向上动态规划。
至于算法的复杂度就留做练习，是低于以上任何一种暴力方法的。code

 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能跑出来，真是神同样的算法，神同样的代码。blog

————————————————————————更新——————————————————————————————ci

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

 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 }