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^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++代码:
#include<bits/stdc++.h>
using 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<ll> S(nv+1);
std::vector<ll> V(nv+1);
for(ll i=0;i<r;i++) {
V[i] = n / (i+1);
}
for(ll i=r;i<nv;i++) {
V[i] = V[i-1] - 1;
}
for(ll i = 0;i<nv;i++) {
S[i] = V[i] - 1; //求素数个数
}
for(ll p=2;p<=r;p++) {
if(S[nv-p] > S[nv-p+1]) {
ll sp = S[nv-p+1]; // sum of primes smaller than p
ll p2 = p*p;
// std::cout << "p=" << p << '\n'; // p is prime
for(ll i=0;i<nv;i++) {
if(V[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<ll> S(nv+1);
std::vector<ll> V(nv+1);
for(ll i=0;i<r;i++) {
V[i] = n / (i+1);
}
for(ll i=r;i<nv;i++) {
V[i] = V[i-1] - 1;
}
for(ll i = 0;i<nv;i++) {
S[i] = V[i] * ( V[i] + 1) / 2 - 1; //求素数和
}
for(ll p=2;p<=r;p++) { // p is prime
if(S[nv-p] > S[nv-p+1]) {
ll sp = S[nv-p+1]; // sum of primes smaller than p
ll p2 = p*p;
for(ll i=0;i<nv;i++) {
if(V[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;
}
最新文章
- ubuntu下安装lrzsz
- 葱类 Allium
- Iterator 迭代器(一)
- hdu 5441 travel 离线+带权并查集
- PHP学习笔记 - 进阶篇(4)
- (转)asp.net动态设置标题title 关键字keywords 描述descrtptions
- Oracle instr用法
- Setup Git Server in CentOS 6.3
- [国嵌笔记][024][ARM汇编编程概述]
- eclipse每次闪退后都提示查看\workspace\.metadata\.log
- 关闭sublime自动检测更新提示
- mac office2016
- React V16.x 生命周期调整
- Mac OS安装Windows各版本时注意事项(2014年后的Mac机相信会有这些问题)
- .net下WinDbg使用说明
- 九章面试题:Find first K frequency numbers 解题报告
- python dict与list
- (研)for循环的一个bug以及3个while循环的快排
- java并发回答
- Java50道经典习题-程序14 求日期