杜教筛学习笔记
$\quad $ 《一种可以在低于线性时间内计算出积性函数前缀和的筛法》
$\quad $ 其实说是筛法，倒不如说是技巧。先说一下大致过程吧。
$\quad $ 假设你要求 \(\sum _{i=1}^{n}f(i)\) ，\(f(x)\) 是一个积性函数，记 \(s(n)=\sum _{i=1}^{n}f(i)\) 。
$\quad $ 那么你要找出一个积性函数 \(g\) ，使得 \(g\) 与 \(f\) 的狄利克雷卷积（记为 \(h\) ）的前缀和方便计算。
主要过程如下：
然后提出含 \(s(n)\) 的项。
然后递归求解+记忆化即可（吉大好像有个递推求解的方法，感兴趣的可以去看看），时间复杂度 \(O(n ^{\frac{2}{3}})\) 。
$\quad $ 写两个常用东西：
\(\varphi * I=id\)
\(I*\mu=\epsilon\)
$\quad $ 上面的乘号代表卷积运算，这两个东西其实就是我在莫反做题记开头给到的两个结论。这两个东西的应用似乎有人叫做欧拉反演和莫比乌斯反演？不是很清楚。
看看例题吧。
杜教筛模板
$\quad $ 即：
$\quad $ 先去处理 \(ans1\) ，按照上面的方法，我们令 \(g\) 函数为 \(I\) ,则：
于是有：
$\quad $ 然后处理 \(ans2\) ，我们仍令 \(g=I\) ，则：
$\quad $ 这个的前缀和显然就是 \(1\) 了。于是我们得到：
$\quad $ 然后分别计算 \(ans1\) 和 \(ans2\) 即可。
点击查看代码
#define yhl 0
#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/hash_policy.hpp>
using namespace __gnu_pbds;
#define ll long long
const int N=6e6+10,M=6e6;
int t,prime[N],mu[N],n,tot;
ll phi[N];
gp_hash_table<ll,ll>ph,m;
bool vis[N];
inline void init(){
mu[1]=phi[1]=1;vis[1]=1;
for(int i=1;i<=M;i++){
if(!vis[i]){
prime[++tot]=i;
mu[i]=-1;
phi[i]=i-1;
}
for(int j=1;j<=tot&&i*prime[j]<=M;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j])){
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
mu[i*prime[j]]=-mu[i];
phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}
for(int i=1;i<=M;i++){
mu[i]+=mu[i-1];
phi[i]+=phi[i-1];
}
}
inline ll s(ll x){
if(x<=M)return mu[x];
if(m[x])return m[x];
ll ans=1;
ll r=yhl;
for(ll l=2;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*s(x/l);
}
return m[x]=ans;
}
inline ll sl(ll x){
if(x<=M)return phi[x];
if(ph[x])return ph[x];
ll ans=x*(x+1)>>1;
ll r=yhl;
for(ll l=2;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*sl(x/l);
}
return ph[x]=ans;
}
signed main(){
init();
scanf("%d",&t);
while(t--){
scanf("%d",&n);
printf("%lld %lld\n",sl(n),s(n));
}
return yhl;
}
神犇和蒟蒻
$\quad $ 小清新题，求：
$\quad $ \(ans1\) 一眼就是 \(1\) ，没什么好说的，我们来看 \(ans2\) ：
$\quad $ 然后令 \(g=id\)
那么：
$\quad $ 然后我们还需要知道：\(\sum _{i=1}^{n}i^2=\frac{n\times(n+1)\times(2n+1)}{6}\)
$\quad $ 然后就可以得出：
点击查看代码
#define yhl 0
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e6,p=1e9+7;
int prime[N],tot,ph[N+10],n,inv;
bool vis[N+10];
unordered_map<int,int>f;
int qum(int a,int b){
int ans=1;
while(b){
(b&1)&&(ans=ans*a%p);
a=a*a%p;
b>>=1;
}
return ans;
}
void init(){
ph[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++tot]=i;
ph[i]=i-1;
}
for(int j=1;j<=tot&&i*prime[j]<=N;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j])){
ph[i*prime[j]]=ph[i]*prime[j]%p;
break;
}
ph[i*prime[j]]=ph[i]*ph[prime[j]]%p;
}
}
for(int i=1;i<=N;i++)ph[i]=ph[i]*i%p;
for(int i=1;i<=N;i++)ph[i]=(ph[i]+ph[i-1])%p;
}
inline int op(int x){
return x*(x+1)%p*(2*x%p+1)%p*inv%p;
}
inline int up(int x){
return x*(x+1)/2%p;
}
int get(int n){
if(n<=N)return ph[n];
if(f[n])return f[n];
int ans=op(n);
for(int l=2,r=0;l<=n;l=r+1){
r=n/(n/l);
ans=(ans-(up(r)-up(l-1)+p)*get(n/l)%p+p)%p;
}
return f[n]=ans;
}
signed main(){
init();inv=qum(6,p-2);
scanf("%lld",&n);
printf("1\n%lld",get(n));
return yhl;
}
简单的数学题
\begin{aligned}
ans&=\sum _{i=1}^{n}\sum _{j=1}^{n}ijgcd(i,j)\\
&=\sum _{d=1}^{n}d\sum _{i=1}^{n}\sum _{j=1}^{n}ij[gcd(i,j)=d]\\
&=\sum _{d=1}^{n} d ^{3} \sum _{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum _{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum _{k|gcd(i,j)}\mu(k)\\
&=\sum _{d=1}^{n}d ^3 \sum _{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)k ^2 (\sum _{i=1}^{\lfloor\frac{n}{kd}\rfloor}i) ^2\\
&=\sum _{T=1}^{n}T ^2 \varphi(T) (\frac{n\times(n+1)}{2}) ^2\\
\end{aligned}
$\quad $ 我们吧前面那部分单独拿出来，记为 \(f\) ，则 \(f(x)=x ^2 \varphi(x)\) ，我们可以令 \(g=id ^2\) ，那么：
$\quad $ 这里要用到一个等式：
$\quad $ 然后我们把柿子摆出来：
$\quad $ 这里还要用到上一个题的方法求和🌚。
点击查看代码
#define yhl 0
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=3e6;
int phi[N+10],prime[N+10],tot,n,p,inv,iv;
bool vis[N+10];
unordered_map<int,int>f;
void init(){
phi[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++tot]=i;
phi[i]=i-1;
}
for(int j=1;j<=tot&&i*prime[j]<=N;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j])){phi[i*prime[j]]=phi[i]*prime[j]%p;break;}
phi[i*prime[j]]=phi[i]*phi[prime[j]]%p;
}
}
for(int i=1;i<=N;i++)phi[i]=i*i%p*phi[i]%p;
for(int i=1;i<=N;i++)phi[i]=(phi[i]+phi[i-1])%p;
}
inline int qum(int a,int b){
int ans=1;
while(b){
(b&1)&&(ans=ans*a%p);
a=a*a%p;
b>>=1;
}
return ans;
}
inline int op(int x){
x%=p;
return x*(x+1)%p*iv%p;
}
inline int func(int x){
x%=p;
return x*(x+1)%p*(2*x%p+1)%p*inv%p;
}
int get(int x){
if(x<=N)return phi[x];
if(f[x])return f[x];
int r=op(x);
int ans=(r*r)%p;
for(int i=2;i<=x;i=r+1){
r=x/(x/i);
ans=(ans-(func(r)-func(i-1)+p)%p*get(x/i)%p+p)%p;
}
return f[x]=ans;
}
signed main(){
scanf("%lld%lld",&p,&n);
init();inv=qum(6,p-2);iv=qum(2,p-2);
int ans=yhl;
for(int i=1,r=yhl;i<=n;i=r+1){
r=n/(n/i);
ans=(ans+(get(r)-get(i-1)+p)%p*op(n/i)%p*op(n/i)%p+p)%p;
}
printf("%lld",ans);
return yhl;
}
[bzoj4176] Lucas的数论
有一个这道题的弱化版本，就挂那道题的链接吧
推导过程我在莫反做题记里面写了，这里就不赘述了，自己去看即可（就是太懒不想写柿子了）。
$\quad $ 这两个不一样的地方就在于 \(n\) 的取值极大，我们不能全部预处理出来，处理出 \(\mu\) 的前 \(n ^{\frac{2}{3}}\) 项即可，对于 \(\sum _{i=1}^{n}\lfloor\frac{n}{i}\rfloor\) 我们只能采用根号分治了。
点击查看代码
#define yhl 0
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e6,p=1e9+7;
int mu[N+10],ans,n,prime[N],tot;
bool vis[N+10];
unordered_map<int,int>f;
void init(){
mu[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1;j<=tot&&i*prime[j]<=N;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j]))break;
mu[i*prime[j]]=mu[i]*mu[prime[j]];
}
}
for(int i=1;i<=N;i++)mu[i]=(mu[i]+mu[i-1]+p)%p;
}
int get_f(int n){
int ans=yhl;
for(int i=1,r;i<=n;i=r+1){
r=n/(n/i);
ans=(ans+(r-i+1)*(n/i)%p)%p;
}
return ans*ans%p;
}
int get_mu(int x){
if(x<=N)return mu[x];
if(f[x])return f[x];
int ans=1;
for(int i=2,r;i<=x;i=r+1){
r=x/(x/i);
ans=(ans-(r-i+1)*get_mu(x/i)%p+p)%p;
}
return f[x]=ans;
}
signed main(){
init();
scanf("%lld",&n);
for(int i=1,r;i<=n;i=r+1){
r=n/(n/i);
ans=(ans+(get_mu(r)-get_mu(i-1)+p)*get_f(n/i)%p)%p;
}
printf("%lld",ans);
return yhl;
}
DZY Loves Math IV
题意
对于给定的 \(n\in [1,1e5]\)、 \(m\in [1,1e9]\) ，求：
$\quad $ 掘势好题，这里也不挂链接了。
$\quad $ 首先定义一个二元函数 \(f(n,m)=\sum _{i=1}^{m}\varphi(ni)\) ，然后答案就是 \(\sum _{i=1}^{n}f(i,m)\)
$\quad $ 下面考虑如何转化 \(\varphi(ij)\) ，我们知道当两个数 \(a b\) 互质时 \(\varphi(ab)=\varphi(a)\varphi(b)\) 、当 \(a|b\) 时，\(\varphi(ab)=a\varphi(b)\) 。
$\quad $ 于是，我们可以给 \(i\) 或者 \(j\) 提出一个 \(gcd(i,j)\) 这样就可以将其拆开了。
先尝试化简一下：
\begin{aligned}
f(n,m)&=\sum _{i=1}^{m}\varphi(ni)\\
&=\sum _{i=1}^{m}\varphi(i)\varphi(\frac{n}{gcd(n,i)})gcd(n,i)\\
&=\sum _{i=1}^{m}\varphi(i)\varphi(\frac{n}{gcd(n,i)})\sum _{k|gcd(n,i)}\varphi(d)\\
&=\sum _{i=1}^{m}\varphi(i)\sum _{k|gcd(n,i)}\varphi(\frac{n}{d})\\
&=\sum _{d|n}\varphi(\frac{n}{d})\sum _{i=1}^{\lfloor\frac{m}{d}\rfloor}\varphi(di)\\
&=\sum _{d|n}\varphi(\frac{n}{d})s(d,\lfloor\frac{m}{d}\rfloor)
\end{aligned}
$\quad $ 但是呢，这种推导真的对吗？我们观察一下第三行到第四行的那步转化。我们利用了 \(a b\) 互质时 \(\varphi(ab)=\varphi(a)\varphi(b)\) ，但是在这里，\(\frac{n}{gcd(n,i)}\) 与枚举出的 \(d\) 不一定互质。那有没有办法让它们互质呢？
$\quad $ 我们可以定义一个东西 \(t _n = \prod _{p|n p\in prime}p\) ，并用它替换原柿子里的 \(n\) ，这样就可以保证上面的转化正确了。我们重新推导：
\begin{aligned}
f(n,m)&=\sum _{i=1}^{m}\varphi(ni)\\
&=\frac{n}{t_n}\sum _{i=1}^{m}\varphi(t _{n}i)\\
&=\frac{n}{t_n}\sum _{i=1}^{m}\varphi(i)\varphi(\frac{t_n}{gcd(t _{n},i)})\sum _{d|gcd(t _{n},i)}\varphi(d)\\
&=\frac{n}{t_n}\sum _{i=1}^{m}\varphi(i)\sum _{d|gcd(t _{n},i)}\varphi(\frac{n}{d})\\
&=\frac{n}{t_n}\sum _{d|n}\varphi(\frac{n}{d})\sum _{i=1}^{\lfloor\frac{m}{d}\rfloor}\varphi(di)\\
&=\frac{n}{t_n}\sum _{d|n}\varphi(\frac{n}{d})f(d,\lfloor\frac{m}{d}\rfloor)
\end{aligned}
$\quad $ 对于 \(t_n\) ，这东西是个积性函数，直接线性筛即可。然后就没什么了。
点击查看代码
#define yhl 0
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e6,p=1e9+7;
int phi[N+1],prime[N+1],tot,t[N];
bool vis[N+1];
unordered_map<int,int>ph;
unordered_map<int,int>f[N];
void init(){
phi[1]=t[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++tot]=i;
phi[i]=i-1;
t[i]=i;
}
for(int j=1;j<=tot&&i*prime[j]<=N;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j])){
phi[i*prime[j]]=phi[i]*prime[j]%p;
t[i*prime[j]]=t[i];
break;
}
phi[i*prime[j]]=phi[i]*phi[prime[j]]%p;
t[i*prime[j]]=t[i]*prime[j];
}
}
for(int i=1;i<=N;i++){
phi[i]=(phi[i]+phi[i-1])%p;
}
}
int get_phi(int x){
if(x<=N)return phi[x];
if(ph[x])return ph[x];
int ans=x*(x+1)/2%p;;
for(int i=2,r;i<=x;i=r+1){
r=x/(x/i);
ans=(ans-(r-i+1)*get_phi(x/i)%p+p)%p;
}
return ph[x]=ans;
}
int get_s(int n,int m){
if(n==1)return f[n][m]=get_phi(m);
if(m==1)return (phi[n]-phi[n-1]+p)%p;
if(!m)return yhl;
if(f[n][m])return f[n][m];
int ans=yhl,op=n/t[n];
for(int i=1;i*i<=t[n];i++){
if(!(t[n]%i)){
ans=(ans+(phi[t[n]/i]-phi[t[n]/i-1]+p)*get_s(i,m/i)%p)%p;
if(i*i!=t[n]){
ans=(ans+(phi[i]-phi[i-1]+p)*get_s(t[n]/i,m/(t[n]/i))%p)%p;
}
}
}
return f[n][m]=ans*op%p;
}
signed main(){
init();int n,m;
scanf("%lld%lld",&n,&m);
int ans=yhl;
for(int i=1;i<=n;i++){
ans=(ans+get_s(i,m))%p;
}
printf("%lld",ans);
return yhl;
}
[CQOI2015]选数
题面
即：
\begin{aligned}
ans&=\sum _{a _1 =L}^{H}\sum _{a _2 =L}^{H}\cdots \sum _{a _n =L}^{H}[gcd _{i=1}^{n}a _{i} =K]\\
&=\sum _{a _1 =\lceil\frac{L}{K}\rceil}^{\lfloor\frac{H}{K}\rfloor}\sum _{a _2 =\lceil\frac{L}{K}\rceil}^{\lfloor\frac{H}{K}\rfloor}\cdots\sum _{a _n =\lceil\frac{L}{K}\rceil}^{\lfloor\frac{H}{K}\rfloor}[gcd _{i=1}^{n} =1]\\
&=\sum _{a _1 =\lceil\frac{L}{K}\rceil}^{\lfloor\frac{H}{K}\rfloor}\sum _{a _2 =\lceil\frac{L}{K}\rceil}^{\lfloor\frac{H}{K}\rfloor}\cdots\sum _{a _n =\lceil\frac{L}{K}\rceil}^{\lfloor\frac{H}{K}\rfloor}\sum _{d|gcd _{i=1}^{n}}\mu(d)\\
&=\sum _{d=1}^{\lfloor\frac{H}{K}\rfloor}\mu(d)(\lfloor\frac{H}{K}\rfloor-\lceil\frac{L}{K}\rceil+1) ^{n}\\
&=\sum _{d=1}^{\lfloor\frac{H}{K}\rfloor}\mu(d)(\lfloor\frac{H}{K}\rfloor-\lfloor\frac{L-1}{K}\rfloor) ^{n}
\end{aligned}
$\quad $ 直接整除分块加杜教筛即可。
点击查看代码
#define yhl 0
#include<bits/stdc++.h?
using namespace std;
#define int long long
const int N=1e6,p=1e9+7;
int prime[N],mu[N+10],tot,l,h,k,ans,n;
unordered_map<int,int> f;
bool vis[N+10];
inline int qum(int a,int b){
int ans=1;
while(b){
(b&1)&&(ans=ans*a%p);
a=a*a%p;
b>>=1;
}
return ans;
}
inline void init(){
vis[1]=mu[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1;j<=tot&&i*prime[j]<=N;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j]))break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=N;i++)mu[i]=(mu[i-1]+mu[i]+p)%p;
}
inline int func(int x){
if(x<=N)return mu[x];
if(f[x])return f[x];
int ans=1;
for(int i=2,r;i<=x;i=r+1){
r=x/(x/i);
ans=(ans-(r-i+1)*func(x/i)%p+p)%p;
}
return f[x]=ans;
}
signed main(){
scanf("%lld%lld%lld%lld",&n,&k,&l,&h);
init();
l=(l-1)/k,h=h/k;
for(int i=1,r;i<=h;i=r+1){
if(l/i)r=min(h/(h/i),l/(l/i));
else r=h/(h/i);
ans=(ans+(func(r)-func(i-1)+p)*qum(h/i-l/i,n)%p)%p;
}
printf("%lld",ans);
return yhl;
}
