卧薪尝胆，厚积薄发。
质数筛法
Date: Tue Jul 03 22:57:52 CST 2018
In Category:
线性筛
筛出一张质数表。
for(register int i = 2;i < MAXN;++i)isprime[i] = true;
for(register int i = 2;i < MAXN;++i)
{
if(isprime[i])prime[++tot] = i;
for(register int j = 1;j <= tot && i * prime[j] < MAXN;++j)
{
isprime[i * prime[j]] = false;
if(i % prime[j] == 0)break;
}
}
求得一些积性函数值。
for(register int i = 2;i < MAXN;++i)isprime[i] = true;
phi[1] = 1;/*欧拉函数*/mu[1] = 1;/*莫比乌斯函数*/cnt[1] = 0;/*最小质因子出现次数*/
mdiv[1] = 1;/*最小质因子*/mdivs[1] = 1;/*最小质因子出现次数次方*/
sigma[1] = 1;/*约数和 除数函数*/numdiv[1] = 1;/*约数个数*/
for(register int i = 2;i < MAXN;++i)
{
if(isprime[i])
{
prime[++tot] = i;
phi[i] = i - 1;mu[i] = -1;cnt[i] = 1;
mdiv[i] = i;mdivs[i] = i;
sigma[i] = i + 1;
numdiv[i] = 2;
}
for(register int j = 1;j <= tot && i * prime[j] < MAXN;++j)
{
int k = i * prime[j];
isprime[k] = false;
mdiv[k] = prime[j];
if(i % prime[j] == 0)
{
phi[k] = phi[i] * prime[j];mu[k] = 0;cnt[k] = cnt[i] + 1;
mdivs[k] = mdivs[i] * mdiv[i];
sigma[k] = sigma[i / mdivs[i]] * (sigma[i] / sigma[i / mdivs[i]] * mdiv[i] + 1);
numdiv[k] = (cnt[k] + 1) * numdiv[i] / (cnt[i] + 1);
break;
}
else
{
phi[k] = phi[i] * phi[prime[j]];mu[k] = mu[i] * mu[prime[j]];cnt[k] = 1;
mdivs[k] = prime[j];
sigma[k] = sigma[i] * sigma[prime[j]];
numdiv[k] = numdiv[i] * numdiv[prime[j]];
}
}
}
若
$f(p^k)$
可以在
$O(1)$
时间内算得，则
$f(x)$
可以在
$O(n)$
时间内筛得
$[1,n]$
中所有函数值。
杜教筛
在低于线性的时间复杂度内求得某个数论函数的前缀和。
原理
设要求
$\begin{align}F(n)=\sum_{i=1}^nf(i)\end{align}$
，
$f(x)$
是某个数论函数。
找到另外两个数论函数
$g$
和
$h$
使得
$f\times g=h$
，其中
$\times$
表示狄利克雷卷积，而
$g$
和
$h$
是两个方便求前缀和的函数。
$\begin{align}H(n)=\sum_{k=1}^n h(k)\end{align}$
$\begin{align}=\sum_{k=1}^n \sum_{d|k}g(d)f(\frac k d)\end{align}$
交换求和号，得：
$\begin{align}=\sum_{d=1}^n \sum_{d|k}g(d)f(\frac k d)\end{align}$
$\begin{align}=\sum_{d=1}^n \sum_{k=1}^{\lfloor\frac n d\rfloor}g(d)f(k)\end{align}$
$\begin{align}=\sum_{d=1}^ng(d)F(\lfloor\frac n d\rfloor)\end{align}$
于是
$\begin{align}g(1)F(n)=H(n)-\sum_{i=2}^ng(i)F(\lfloor\frac n i\rfloor)\end{align}$
复杂度证明
右边的
$F(\lfloor\frac n i\rfloor)$
可以用除法分块优化到
$O(\sqrt n)$
由于
$\lfloor\lfloor n/i\rfloor/j\rfloor=\lfloor n/(ij)\rfloor$
所以使用记忆化搜索递归实现：每一层使用的
$F(\lfloor\frac n i\rfloor)$
一定在第一层展开出现过，而第一层计算是
$O(\sqrt n)$
的。
由于每一次要递归计算的值都会在最外层循环被计算，所以只会展开一层每次展开复杂度
$\sqrt x$
，复杂度为：
$\begin{align}\sum_{i=1}^{\sqrt n}\sqrt{n / i}=O(n^{\frac 3 4})\end{align}$
优化
线性筛预处理出
$[1,n^{\frac 2 3}]$
范围内的
$F$
值，则复杂度降为
$O(n^{\frac 2 3})$
实现
筛
$\varphi$
和
$\mu$
的模板
筛
$\varphi$
时，
$g$
选
$1$
，
$h$
选
$Id$
。
筛
$\mu$
时，
$g$
选
$1$
，
$h$
选
$\varepsilon$
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<map>
#include<cstring>
using namespace std;
int testcases;
#define MAXN 5000010
typedef long long ll;
ll n;
map<int,ll> _phi,_mu;
ll phi[MAXN],mu[MAXN];
int prime[MAXN],tot = 0;
bool isprime[MAXN];
void init()
{
for(int i = 2;i < MAXN;++i)isprime[i] = true;
phi[1] = 1;
mu[1] = 1;
for(int i = 2;i < MAXN;++i)
{
if(isprime[i])
{
prime[++tot] = i;
phi[i] = i - 1;
mu[i] = -1;
}
for(int j = 1;j <= tot && i *prime[j] < MAXN;++j)
{
isprime[i * prime[j]] = false;
if(i % prime[j] == 0)
{
phi[i * prime[j]] = phi[i] * prime[j];
mu[i * prime[j]] = 0;
break;
}
else
{
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
mu[i * prime[j]] = -1 * mu[i];
}
}
}
for(int i = 2;i < MAXN;++i)
{
phi[i] += phi[i - 1];
mu[i] += mu[i - 1];
}
return;
}
ll calcphi(ll n)
{
map<int,ll>::iterator it;
if(n < MAXN)return phi[n];
if((it = _phi.find(n)) != _phi.end())return it -> second;
ll ans = (((ll)n * (n + 1)) >> 1);
for(ll l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ans -= (r - l + 1) * calcphi(n / l);
}
return _phi[n] = ans;
}
ll calcmu(ll n)
{
map<int,ll>::iterator it;
if(n < MAXN)return mu[n];
if((it = _mu.find(n)) != _mu.end())return it -> second;
ll ans = 1;
for(ll l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ans -= (r - l + 1) * calcmu(n / l);
}
return _mu[n] = ans;
}
int main()
{
init();
scanf("%d",&testcases);
for(int i = 1;i <= testcases;++i)
{
scanf("%lld",&n);
printf("%lld %lld\n",calcphi(n),calcmu(n));
}
return 0;
}
杜教筛的重点是找到一组
$g$
和
$h$
使得
$g\times f=h$
，一般不容易找到。
min
用于在低于线性的时间复杂度下求某个积性函数的前缀和。
具体来说时间复杂度是
$O(\frac{n^\frac 3 4}{\log n})$
。
要求
$f(p)$
是一个关于
$p$
的简单多项式，
$f(p^c)$
可以快速计算。
先不考虑求前缀和，考虑求一个函数：
$$
F(n)=\sum_{i=1}^n[i\ is\ prime]f(i)
$$
设：
$$
g(n,j)=\sum^n_{i=1}[i\ is\ prime\ or\ minp(i)>p_j]f(i)
$$
含义就是
$i$
是质数或者
$i$
的最小质因子大于
$p_j$
，把
$1\sim n$
内所有满足条件的
$f(i)$
加起来。
这个东西的实际意义就是埃氏筛法在筛掉前
$j$
个素数后没被筛掉的所有数的
$f(i)$
之和，包括所有的素数和还没有被筛掉的合数，因为这些合数
$minp(i)>p_j$
，所以
$\min\_25$
筛又被叫做埃筛的扩展。
考虑
$g(n,j)$
的转移：
若
$p_j^2>n$
，那显然这个质数不会筛掉任何其他数字，
$g(n,j)=g(n,j-1)$
。
若
$p_j^2\leqslant n$
，那么被筛掉的数一定是最小质因子
$=p_j$
的数，因为别的数要么还没被筛，要么已经筛掉了，被筛掉的数除掉
$p_j$
之后最小质因子一定大于等于
$p_j$
，因为原来最小质因子是
$p_j$
，而且又由于
$f$
是积性函数，因此我们可以提出来一个
$f(p_j)$
，但是这个时候我们又多减掉了原来的质数，因此再加回来，于是有：
$$
g(n,j)=g(n,j-1)-f(p_j)\bigl(g(\frac n{p_j},j-1)-\sum_{i=1}^{j-1}f(p_j)\bigr)
$$
但是我们的第一维还是
$n$
个，但我们只要求
$g(n,|p|)$
，而且又由于
$\lfloor\frac{\lfloor\frac n a\rfloor}{b}\rfloor=\lfloor\frac n{ab}\rfloor$
，因此实际有用的
$n$
由除法分块可知只有
$O(\sqrt n)$
个，而且当
$p_j^2>n$
的时候不转移，因此实际有用的质数也只有
$O(\sqrt{\frac n {\ln n}})$
个，因此减少了复杂度。
如果让
$f(i)=1$
，那么就有了一个
$\min\_25$
筛的经典应用：求区间质数个数。
代码：
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<cctype>
#include<cstring>
using namespace std;
typedef long long ll;
ll n;
#define MAXN 10000001
#define N 10000000
bool isprime[MAXN];
ll prime[MAXN],tot = 0;
void sieve()
{
for(ll i = 2;i <= N;++i)isprime[i] = true;
for(ll i = 2;i <= N;++i)
{
if(isprime[i])prime[++tot] = i;
for(ll j = 1;j <= tot && i * prime[j] <= N;++j)
{
isprime[i * prime[j]] = false;
if(i % prime[j] == 0)break;
}
}
return;
}
ll id1[MAXN],id2[MAXN];
ll id(ll k){return (k <= N ? id1[k] : id2[n / k]);}
ll w[MAXN],m = 0;
ll g[MAXN];
int main()
{
sieve();
scanf("%lld",&n);
for(ll l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
w[++m] = n / l;
if(w[m] <= N)id1[w[m]] = m;
else id2[n / w[m]] = m;
g[m] = w[m] - 1;
}
for(ll j = 1;j <= tot;++j)
{
for(ll i = 1;i <= m && prime[j] * prime[j] <= w[i];++i)
{
g[i] = g[i] - g[id(w[i] / prime[j])] + j - 1;
}
}
cout << g[1] << endl;
return 0;
}
然后我们继续考虑怎么用这个东西求积性函数前缀和，我们现在已经对所有
$x=\lfloor\frac n i\rfloor$
求出了：
$$
g(x,|p|)=\sum_{i=1}^x[i\ is\ prime]f(i)
$$
我们设：
$$
S(n,j)=\sum_{i=1}^n[minp(i)\geqslant p_j]f(i)
$$
那么我们要求的答案就是
$S(n,1)+f(1)$
。
现在考虑怎么求
$S(n,j)$
，首先是质数的部分：
$$
g(n,|p|)-\sum_{i=1}^{j-1}f(p_i)
$$
然后再考虑合数，枚举这个合数的最小质因子及其出现次数，然后由于积性把
$f$
提出来：
$$
S(n,j)=g(n,|p|)-\sum_{i=1}^{j-1}f(p_i)+\sum_{k=j}^{p_k^2\leqslant n}\sum_{e=1}^{p^{e+1}_k\leqslant n}\Bigl(S(\frac n{p^e_k},k+1)\times f(p^e_k)+f(p^{e+1}_k)\Bigr)
$$
这个的复杂度也是
$O(\frac{n^{\frac 3 4}}{\log n})$
。
局限性：在筛
$g$
的时候，要求
$f$
完全积性，满足这个的函数一般有
$id^k$
和
$1$
。
不算模板。
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<cctype>
#include<cstring>
using namespace std;
typedef long long ll;
ll n;
#define MAXN 1000010
int N;
#define MOD 1000000007
#define INV2 500000004
int prime[MAXN],tot = 0;
int ps[MAXN];
bool isprime[MAXN];
void sieve()
{
for(int i = 2;i <= N;++i)isprime[i] = true;
for(int i = 2;i <= N;++i)
{
if(isprime[i])
{
prime[++tot] = i;
ps[tot] = (ps[tot - 1] + i) % MOD;
}
for(int j = 1;j <= tot && i * prime[j] <= N;++j)
{
isprime[i * prime[j]] = false;
if(i % prime[j] == 0)break;
}
}
return;
}
ll w[MAXN],m = 0;
int id1[MAXN],id2[MAXN];
inline int id(ll k){return (k <= N ? id1[k] : id2[n / k]);}
int pc[MAXN];
int g[MAXN];
int S(ll n,int j)
{//cout << n << " " << j << endl;
if(n <= 1 || prime[j] > n)return 0;
int k = id(n);
int res = (((g[k] - pc[k] + MOD) % MOD - ps[j - 1] + MOD) % MOD + j - 1) % MOD;
if(j == 1)res = (res + 2) % MOD;
for(int i = j;i <= tot && 1ll * prime[i] * prime[i] <= n;++i)
{
ll p1 = prime[i],p2 = 1ll * prime[i] * prime[i];
for(int e = 1;p2 <= n;++e,p1 = p2,p2 *= prime[i])
{
res = ((res + 1ll * S(n / p1,i + 1) * (prime[i] ^ e) % MOD) % MOD + (prime[i] ^ (e + 1))) % MOD;
}
}
return res;
}
int main()
{
cin >> n;
N = sqrt(n);
sieve();
for(ll l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
w[++m] = n / l;
if(w[m] <= N)id1[w[m]] = m;
else id2[n / w[m]] = m;
pc[m] = (w[m] - 1) % MOD;
g[m] = ((w[m] + 2) % MOD) * ((w[m] - 1) % MOD) % MOD * INV2 % MOD;
}
for(int j = 1;j <= tot;++j)
{
for(int i = 1;i <= m && 1ll * prime[j] * prime[j] <= w[i];++i)
{
int k = id(w[i] / prime[j]);
pc[i] = ((pc[i] - pc[k] + MOD) % MOD + j - 1) % MOD;
g[i] = (g[i] - 1ll * prime[j] * ((g[k] - ps[j - 1] + MOD) % MOD) % MOD + MOD) % MOD;
}
}
cout << (S(n,1) + 1) % MOD << endl;
return 0;
}
In tag:
Timeline
About
Toolbox
Friends
Copyright © 2020
wjh15101051
ღゝ◡╹)ノ♡