Lucy_Hedgehog techniques

在project euler 的第\(10\)题的 \(forum\) 中 Lucy Hedgehog 提到的这种方法。

求 \(n\) 以内素数个数以及求 \(n\) 以内素数和的算法。

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

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

\(1.\) \(p\le1\) 。此时所有数都还没有被筛掉，所以 \(S(v,p)=\sum_{i=2}^{v}i=\frac{(2+v)(v-1)}{2}\)。

\(2.\) \(p\) 不是素数。因为筛法中 \(p\) 早已被别的数筛掉，所以在这步什么都不会做，所以此时 \(S(v,p)=S(v,p-1)\)。

\(3.\) \(p\) 是素数，但是 \(v。因为每个合数都一定有一个不超过其平方根的素因子，如果筛到 \(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\) 的定义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)\)。

这是求前 \(n\) 的素数和的方法。

至于求前 \(n\) 的素数个数的方法也差不多。

只需要把代码修改一下即可。

复杂度： \(O(n^{0.75})\)

C++代码：

#includeusing namespace std;
typedef long long ll;ll check(ll v, ll n, ll ndr, ll nv) {return v >= ndr ? (n / v - 1) : (nv - v);
}// ll S[10000000];
// ll V[10000000];
ll primenum(ll n) // O(n^(3/4))
{ll r = (ll)sqrt(n);ll ndr = n / r;assert(r*r <= n && (r+1)*(r+1) > n);ll nv = r + ndr - 1;std::vector S(nv+1);std::vector V(nv+1);for(ll i=0;i S[nv-p+1]) {ll sp = S[nv-p+1]; // sum of primes smaller than pll p2 = p*p;//    std::cout << "p=" << p << '\n'; // p is prime for(ll i=0;i= p2) {S[i] -= 1LL  * (S[check(V[i] / p, n, ndr, nv)] - sp);// //求素数个数}else break;}}}return S[0];
}
ll primesum(ll n) // O(n^(3/4))
{ll r = (ll)sqrt(n);ll ndr = n / r;assert(r*r <= n && (r+1)*(r+1) > n);ll nv = r + ndr - 1;std::vector S(nv+1);std::vector V(nv+1);for(ll i=0;i S[nv-p+1]) {ll sp = S[nv-p+1]; // sum of primes smaller than pll p2 = p*p; for(ll i=0;i= p2) {S[i] -= p* (S[check(V[i] / p, n, ndr, nv)] - sp); //求素数和}else break;}}}return S[0];
}
int main(int argc, char const *argv[]) {
//  std::cout << primesum(1e6) << '\n';std::cout << primenum(1e10) << '\n';std::cout << primesum(2e6) << '\n';cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.\n";return 0;
}

转载于:https://www.cnblogs.com/LzyRapx/p/8228009.html

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！