返回顶部

杜教筛学习笔记

$\quad $ 《一种可以在低于线性时间内计算出积性函数前缀和的筛法》

$\quad $ 其实说是筛法,倒不如说是技巧。先说一下大致过程吧。

$\quad $ 假设你要求 \(\sum _{i=1}^{n}f(i)\)\(f(x)\) 是一个积性函数,记 \(s(n)=\sum _{i=1}^{n}f(i)\)

$\quad $ 那么你要找出一个积性函数 \(g\) ,使得 \(g\)\(f\) 的狄利克雷卷积(记为 \(h\) )的前缀和方便计算。

主要过程如下:

\[\sum _{i=1}^{n}h(i)=\sum _{i=1}^{n}\sum _{d|i}g(d)f(\frac{i}{d}) \]

\[\sum _{i=1}^{n}h(i)=\sum _{d=1}^{n}g(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor}f(i) \]

\[\sum _{i=1}^{n}h(i)=\sum _{d=1}^{n}g(d)s(\lfloor\frac{n}{d}\rfloor) \]

然后提出含 \(s(n)\) 的项。

\[g(1)s(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)s(\lfloor\frac{n}{d}\rfloor) \]

然后递归求解+记忆化即可(吉大好像有个递推求解的方法,感兴趣的可以去看看),时间复杂度 \(O(n ^{\frac{2}{3}})\)

$\quad $ 写两个常用东西:

\(\varphi * I=id\)

\(I*\mu=\epsilon\)

$\quad $ 上面的乘号代表卷积运算,这两个东西其实就是我在莫反做题记开头给到的两个结论。这两个东西的应用似乎有人叫做欧拉反演和莫比乌斯反演?不是很清楚。

看看例题吧。

杜教筛模板

$\quad $ 即:

\[ans1=\sum _{i=1}^{n}\varphi(i) \]

\[ans2=\sum _{i=1}^{n}\mu(i) \]

$\quad $ 先去处理 \(ans1\) ,按照上面的方法,我们令 \(g\) 函数为 \(I\) ,则:

\[h(n)=\sum _{d|n}\varphi(d)=n \]

于是有:

\[s(n)=\sum _{i=1}^{n}i - \sum _{i=2}^{n}s(\lfloor\frac{n}{d}\rfloor)=\frac{n\times(n+1)}{2}-\sum _{d=2}^{n}s(\lfloor\frac{n}{d}\rfloor) \]

$\quad $ 然后处理 \(ans2\) ,我们仍令 \(g=I\) ,则:

\[h=\mu*I=\epsilon \]

$\quad $ 这个的前缀和显然就是 \(1\) 了。于是我们得到:

\[s(n)=1-\sum _{i=2}^{n}s(\lfloor\frac{n}{d}\rfloor) \]

$\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 $ 小清新题,求:

\[ans1=\sum _{i=1}^{n}\mu(i^2) \]

\[ans2=\sum _{i=1}^{n}\varphi(i^2) \]

$\quad $ \(ans1\) 一眼就是 \(1\) ,没什么好说的,我们来看 \(ans2\)

\[ans2=\sum _{i=1}^{n}\varphi(i^2)=\sum _{i=1}^{n}i\varphi(i) \]

$\quad $ 然后令 \(g=id\)
那么:

\[h(n)=\sum _{d|n}d\varphi(d)\frac{n}{d}=n^2 \]

$\quad $ 然后我们还需要知道:\(\sum _{i=1}^{n}i^2=\frac{n\times(n+1)\times(2n+1)}{6}\)

$\quad $ 然后就可以得出:

\[s(n)=\frac{n\times(n+1)\times(2n+1)}{6}-\sum _{i=2}^{n}is(\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 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\) ,那么:

\[h(n)=\sum _{d|n}d ^2 \varphi(d) (\frac{n}{d}) ^2 =n ^2 \sum _{d|n}\varphi(d)=n ^3 \]

$\quad $ 这里要用到一个等式:

\[\sum _{i=1}^{n}i ^3 =(\sum _{i=1}^{n}i) ^2 \]

证明

$\quad $ 然后我们把柿子摆出来:

\[s(n)=(\frac{n(n+1)}{2}) ^2 -\sum _{i=2}^{n}i ^2 s(\lfloor \frac{n}{d}\rfloor) \]

$\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]\) ,求:

\[\sum _{i=1}^{n}\sum _{j=1}^{m}\varphi(ij) \]

$\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;
}

\[To Be Continued…… \]

posted @ 2024-10-16 15:10  无敌の暗黑魔王  阅读(88)  评论(13编辑  收藏  举报