快速求素數和

時間 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 }

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。