杜教筛学习笔记

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

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

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

那么你要找出一个积性函数 $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}})$ 。

写两个常用东西：

$\phi \ast I=id$

$I\ast \mu =ϵ$

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

看看例题吧。

杜教筛模板

即：

$$ans1=\sum _{i=1}^n\phi (i)$$

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

先去处理 $ans1$ ，按照上面的方法，我们令 $g$ 函数为 $I$ ,则：

$$h(n)=\sum _{d|n}\phi (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 )$$

然后处理 $ans2$ ，我们仍令 $g=I$ ，则：

$$h=\mu \ast I=ϵ$$

这个的前缀和显然就是 $1$ 了。于是我们得到：

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

然后分别计算 $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;
}

神犇和蒟蒻

小清新题，求：

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

$$ans2=\sum _{i=1}^n\phi (i^2)$$

$ans1$ 一眼就是 $1$ ，没什么好说的，我们来看 $ans2$ ：

$$ans2=\sum _{i=1}^n\phi (i^2)=\sum _{i=1}^{n}i\phi (i)$$

然后令 $g=id$
那么：

$$h(n)=\sum _{d|n}d\phi (d)\frac{n}{d}=n^2$$

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

然后就可以得出：

$$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{array}{rl}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\phi (T)(\frac{n\times (n+1)}{2}{)}^2\end{array}$$

我们吧前面那部分单独拿出来，记为 $f$ ，则 $f(x)=x^2\phi (x)$ ，我们可以令 $g=i{d}^2$ ，那么：

$$h(n)=\sum _{d|n}{d}^2\phi (d)(\frac{n}{d}{)}^2=n^2\sum _{d|n}\phi (d)=n^3$$

这里要用到一个等式：

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

证明

然后我们把柿子摆出来：

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

这里还要用到上一个题的方法求和🌚。

点击查看代码

#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的数论

有一个这道题的弱化版本，就挂那道题的链接吧

推导过程我在莫反做题记里面写了，这里就不赘述了，自己去看即可（就是太懒不想写柿子了）。

这两个不一样的地方就在于 $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\phi (ij)$$

掘势好题，这里也不挂链接了。

首先定义一个二元函数 $f(n,m)=\sum _{i=1}^m\phi (ni)$ ，然后答案就是 $\sum _{i=1}^{n}f(i,m)$

下面考虑如何转化 $\phi (ij)$ ，我们知道当两个数 $ab$ 互质时 $\phi (ab)=\phi (a)\phi (b)$ 、当 $a|b$ 时，$\phi (ab)=a\phi (b)$ 。

于是，我们可以给 $i$ 或者 $j$ 提出一个 $gcd(i,j)$ 这样就可以将其拆开了。

先尝试化简一下：

$$\begin{array}{rl}f(n,m)& =\sum _{i=1}^m\phi (ni)\\ & =\sum _{i=1}^m\phi (i)\phi (\frac{n}{gcd(n,i)})gcd(n,i)\\ & =\sum _{i=1}^m\phi (i)\phi (\frac{n}{gcd(n,i)})\sum _{k|gcd(n,i)}\phi (d)\\ & =\sum _{i=1}^m\phi (i)\sum _{k|gcd(n,i)}\phi (\frac{n}{d})\\ & =\sum _{d|n}\phi (\frac{n}{d})\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\phi (di)\\ & =\sum _{d|n}\phi (\frac{n}{d})s(d,\lfloor \frac{m}{d}\rfloor )\end{array}$$

但是呢，这种推导真的对吗？我们观察一下第三行到第四行的那步转化。我们利用了 $ab$ 互质时 $\phi (ab)=\phi (a)\phi (b)$ ，但是在这里，$\frac{n}{gcd(n,i)}$ 与枚举出的 $d$ 不一定互质。那有没有办法让它们互质呢？

我们可以定义一个东西 $t_n=\prod _{p|np\in prime}p$ ，并用它替换原柿子里的 $n$ ，这样就可以保证上面的转化正确了。我们重新推导：

$$\begin{array}{rl}f(n,m)& =\sum _{i=1}^m\phi (ni)\\ & =\frac{n}{t_n}\sum _{i=1}^m\phi (t_{n}i)\\ & =\frac{n}{t_n}\sum _{i=1}^m\phi (i)\phi (\frac{t_n}{gcd(t_n,i)})\sum _{d|gcd(t_n,i)}\phi (d)\\ & =\frac{n}{t_n}\sum _{i=1}^m\phi (i)\sum _{d|gcd(t_n,i)}\phi (\frac{n}{d})\\ & =\frac{n}{t_n}\sum _{d|n}\phi (\frac{n}{d})\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\phi (di)\\ & =\frac{n}{t_n}\sum _{d|n}\phi (\frac{n}{d})f(d,\lfloor \frac{m}{d}\rfloor )\end{array}$$

对于 $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{array}{rl}ans& =\sum _{a_1=L}^H\sum _{a_2=L}^H\cdots \sum _{a_n=L}^H[gc{d}_{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 }[gc{d}_{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|gc{d}_{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{array}$$

直接整除分块加杜教筛即可。

点击查看代码

#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;
}

$$ToBeContinued\dots \dots$$