杜教筛学习笔记

本来是想和 Min_25 一起写的，但是太长了（为了避免莫反那篇一样的惨案）所以就分开了w.

0 - 前置

莫比乌斯反演

详见 莫比乌斯反演学习笔记 ，这里摘录一些式子。

Dirichlet卷积：

[(f*g)(x)=sum_{d|x}f(d)g(frac{x}{d}) ]

莫比乌斯反演：

[f=g*1=>g=f*mu\ f(x)=sum_{d|x}g(d)=>g(x)=sum_{d|x}f(d) imes mu(frac{x}{d}). ]

欧拉函数

性质：(sum_{d|n} varphi(d)=n)

表示成 Dirichlet卷积 的形式：(varphi*I=id) .

卷上 (mu) 得

[varphi*I=id\ =>varphi*I*mu=id*mu\ =>varphi * epsilon=id*mu ]

所以有

[varphi(n)=sum_{d|n}mu(d)cdotfrac{n}{d} ]

故

[frac{varphi(n)}{n}=sum_{d|n}frac{mu(d)}{d} ]

1 - 杜教筛

Method

用途：以低于线性的复杂度计算积性函数前缀和 ，即计算：(sum_{i=1}^n f(i)) .

推导：

构造两个积性函数 (h,g) ，使得 (h=f*g) .

记 (S(n)=sum_{i=1}^n f(i)) ，那么有

[egin{aligned} sum_{i=1}^nh(i) &=sum_{i=1}^nsum_{d|i}g(d)cdot f(frac{i}{d})\\ &=sum_{d=1}^ng(d)cdot sum_{i=1}^{lfloor n/d floor} f(i)\\ &=sum_{d=1}^ng(d)cdot S(lfloor frac{n}{d} floor )\\ &=g(1)cdot S(n)+sum_{d=2}^{n}g(d)cdot S(lfloorfrac{n}{d} floor)\\ g(1)S(n)&=sum_{i=1}^nh(i)-sum_{d=2}^ng(d)cdot S(lfloor frac{n}{d} floor ) end{aligned} ]

所以关键就在于找一个容易得到的 (h) .求得之后对后面整除分块就好了，复杂度 (mathcal{O}(n^{frac 23})) .

关于复杂度 （@George1123）

直接递归求是 (Theta(n^{frac{3}{4}})) ：

[egin{split} T(n)=&sqrt n+sumlimits_{i=2}^{sqrt n}left(T(i)cdot T(lfloorfrac{n}{i} floor) ight)\ ge&sqrt n+sumlimits_{i=2}^{sqrt n}left(sqrt icdot sqrt{lfloorfrac{n}{i} floor)} ight)\ ge&sqrt n+sumlimits_{i=2}^{sqrt n}2sqrt{sqrt n}\ =&n^{frac{3}{4}}\ end{split} ]

线性筛求出前 (n^{frac23}) 个，然后再杜教筛，是 (Theta(n^{frac{2}{3}})) .

[egin{split} T(n)=&m+sumlimits_{i=2}^{lfloorfrac nm floor}sqrt{lfloorfrac ni floor}\ =&m+frac{n}{sqrt m}\ ge&2n^{frac{2}{3}}(=: operatorname{while} m=n^{frac{2}{3}})\ end{split} ]

Example

莫比乌斯函数

求

[S(n)=sum_{i=1}^nmu(i). ]

由于 (mu*I=epsilon) ，所以 令 (I=>g,epsilon=>h) ,根据上面的式子

[g(1)S(n)=sum_{i=1}^nh(i)-sum_{d=2}^ng(d)cdot S(lfloor frac{n}{d} floor ) ]

有

[egin{aligned} 1 imes S(n)&=sum_{i=1}^nepsilon(i)-sum_{d=2}^nS(lfloor frac{n}{d} floor )\\ S(n) &= 1-sum_{d=2}^nS(lfloor frac nd floor) end{aligned} ]

欧拉函数

求

[S(n)=sum_{i=1}^nvarphi(i) ]

由于 (varphi *I=ID) ，所以

[egin{aligned} S(n) &=sum_{i=1}^n i-sum_{d=2}^n S(lfloor frac{n}{d} floor)\\ &=frac{n(n+1)}{2}-sum_{d=2}^nS(lfloor frac{n}{d} floor) end{aligned} ]

？函数

求

[S(n)=sum_{i=1}^nicdot varphi(i) ]

考虑 Dirichlet卷积的形式：

[sum_{d|n}(dcdot varphi(d))cdot g(frac{n}{d}) ]

令 (ID=>g) ，有

[sum_{d|n}(dcdot varphi(d))cdot frac{n}{d}=sum_{d|n}ncdot varphi(d)=nsum_{d|n}varphi(d)=n^2 ]

所以有

[S(n)=sum_{i=1}^ni^2-sum_{d=2}^ndcdot S(lfloor frac{n}{d} floor) ]

？？函数

求

[sum_{i=1}^nmu(i)i^2 ]

令 (g(i)=i^2) ，有

[g*f(n)=sum_{d|n}mu(d)d^2cdotBig(frac{n}{d}Big)^2=sum_{d|n}mu(d)=[n=1] ]

所以

[S(n)=sum_{i=1}^nh(i)-sum_{d=2}^nd^2cdot S(lfloor frac{n}{d} floor )=1-sum_{d=2^n}d^2cdot S(lfloor frac{n}{d} floor) ]

Code

大常数选手打的 模板

//Author: RingweEH
const int N=5e5+10;
int pri[N],tot=0,n;
ll mu[N],phi[N];
bool is[N];

void Sieve()
{
    mu[1]=phi[1]=1; is[1]=1;
    for ( int i=2; i<=N-10; i++ )
    {
        if ( !is[i] ) { pri[++tot]=i,mu[i]=-1,phi[i]=i-1; }
        for ( int j=1; j<=tot && (i*pri[j]<=(N-10)); j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j]==0 ) { mu[i*pri[j]]=0,phi[i*pri[j]]=phi[i]*pri[j]; break; }
            mu[i*pri[j]]=-mu[i]; phi[i*pri[j]]=phi[i]*phi[pri[j]];
        }
    }
    for ( int i=2; i<=(N-10); i++ )
        mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}

map<int,ll> mp_mu,mp_phi;

ll Sieve_Phi( int n )
{
    if ( n<=N-10 ) return phi[n];
    if ( mp_phi[n] ) return mp_phi[n];
    ll _res=0;
    for ( int l=2,r; l<=n; l=r+1 )
    {
        r=n/(n/l); _res+=(r-l+1)*Sieve_Phi(n/l);
        if ( r>=2147483647 ) break;
    }
    ll res=(ull)n*(n+1ll)/2ll-_res;
    mp_phi[n]=res; return res;
}

ll Sieve_Mu( int n )
{
    if ( n<=N-10 ) return mu[n];
    if ( mp_mu[n] ) return mp_mu[n];
    ll res=1;
    for ( int l=2,r; l<=n; l=r+1 )
    {
        r=n/(n/l); res-=(r-l+1)*Sieve_Mu(n/l);
        if ( r>=2147483647 ) break;
    }
    return mp_mu[n]=res;
}

int main()
{
    Sieve();
    int T=read();
    while ( T-- )
    {
        int n=read();
        printf("%lld %lld
",Sieve_Phi(n),Sieve_Mu(n) );
    }
    return 0;
}

下面是例题。

• 校内赛-合唱队形

例题：毒瘤的数学题

求

[left(sumlimits_{i=1}^nsumlimits_{j=1}^nijgcd(i,j) ight)mod p ]

(1le nle 10^{10},5 imes10^{8}le ple 1.1 imes 10^{9}) .

Solution

[egin{aligned} sum_{i=1}^nsum_{j=1}^nijgcd(i,j) &=sum_{d=1}^n d^3sum_{i=1}^{lfloor n/d floor}sum_{j=1}^{lfloor n/d floor} ij[gcd(i,j)=1]\\ &=sum_{d=1}^nd^3sum_{i=1}^{lfloor n/d floor}sum_{j=1}^{lfloor n/d floor}ijsum_{k|d}mu(k)\\ 令S(n)=frac{n*(n+1)}2,&=sum_{d=1}^nd^3sum_{k|d}mu(k)k^2S(lfloor frac{n}{dk} floor)^2\\ &=sum_{d=1}^nd^3sum_{k=1}^{lfloor n/d floor}mu(k)k^2S(lfloor frac{n}{dk} floor)^2 end{aligned} ]

用 (T=dk) 进行代换，得到：

[egin{aligned} sum_{d=1}^nd^3sum_{k=1}^{lfloor n/d floor}mu(k)k^2S(lfloor frac{n}{dk} floor)^2 &=sum_{d=1}^nd^3sum_{k=1}^{lfloor n/d floor}mu(k)k^2S(lfloor frac{n}{T} floor)^2\\ &=sum_{T=1}^nS(lfloor frac nT floor)^2sum_{d|T}d^3muBig(frac{T}{d}Big)Big(frac{T}{d}Big)^2\\ &=sum_{T=1}^nS(lfloor frac nT floor)^2T^2sum_{d|T}muBig(frac{T}{d}Big)cdot d\\ end{aligned} ]

根据莫反中的前置知识，有

[varphi(x)=(mu*ID)(x)=sumlimits_{d|x}dcdotmu(frac xd)\\ ]

所以：

[=sum_{T=1}^nS(lfloor frac{n}{T} floor )^2T^2varphi(T)\\ ]

然后就可以愉快地杜教筛了。把模板掏下来以便观察：

[g(1)S(n)=sum_{i=1}^nh(i)-sum_{d=2}^ng(d)cdot S(lfloor frac{n}{d} floor )（杜教筛）\\ (f*g)(x)=sum_{d|x}f(d)g(frac{x}{d})（Dirichlet卷积）\\ varphi*1=ID（和varphi 有关的前置式子） ]

令 (sum(n)=sum_{i=1}^ni^2varphi(i)) ，显然我们需要把 (sum(n)) 中这个 (i^2) 消掉。看着卷积式子，考虑令 (g(n)=n^2) .

于是有：

[(f*g)(n)=sum_{d|n}d^2varphi(d)Big(frac{n}{d}Big)^2=n^2sum_{d|n}varphi(d) ]

然后再把最后一个式子给摁上去：

[(f*g)(n)=n^3. ]

那么就有：

[sum(n)=sum_{i=1}^ni^3-sum_{d=2}^nd^2sum(lfloor frac{n}{d} floor)\\ ]

对于原式：

[sum_{T=1}^nS(lfloor frac{n}{T} floor )^2T^2varphi(T)\\ ]

就可以杜教筛解决了。

//Author: RingweEH
const int N=5e6+10;
int pri[N],tot=0,phi[N],sum[N],Mod;
ll n,inv6,inv2;
bool is[N];

int power( int a,int b )
{
	int res=1;
	for ( ; b; b>>=1,a=1ll*a*a%Mod )
		if ( b&1 ) res=1ll*res*a%Mod;
	return res;
}

void Sieve()
{
	phi[1]=1; is[1]=1;
	for ( int i=2; i<=(N-10); i++ )
	{
		if ( !is[i] ) pri[++tot]=i,phi[i]=i-1;
		for ( int j=1; j<=tot && (i*pri[j]<=(N-10)); j++ )
		{
			int pos=i*pri[j]; is[pos]=1;
			if ( i%pri[j]==0 ) { phi[pos]=1ll*phi[i]*pri[j]%Mod; break; }
			else phi[pos]=1ll*phi[i]*phi[pri[j]]%Mod;
		}
	}
	for ( int i=1; i<=(N-10); i++ )
		sum[i]=(1ll*i*i%Mod*phi[i]%Mod+sum[i-1])%Mod;
}

int pow1( ll x )
{
	x%=Mod; int res=x*(x+1)%Mod*inv2%Mod;
	return res;
}

int pow2( ll x )
{
	x%=Mod; int res=1ll*x*(x+1)%Mod*(2*x+1)%Mod*inv6%Mod;
	return res;
}

int pow3( ll x )
{
	int res=pow1(x); res=1ll*res*res%Mod;
	return res;
}

unordered_map<ll,int> mp;
int Sieve_Mu( ll x )
{
	if ( x<=(N-10) ) return sum[x];
	if ( mp[x] ) return mp[x];
	int res=pow3(x);
	for ( ll l=2,r; l<=x; l=r+1 )
	{
		r=x/(x/l); res=(res-1ll*(pow2(r)-pow2(l-1))*Sieve_Mu(x/l)%Mod)%Mod;
	}
	mp[x]=(res+Mod)%Mod; return mp[x];
}

void init()
{
	Sieve(); inv2=power( 2ll,Mod-2 ); inv6=power( 6ll,Mod-2 );
}

int main()
{
	Mod=read(),n=read(); init();

	ll ans=0;
	for ( ll l=1,r; l<=n; l=r+1 )
	{
		r=n/(n/l); ans=(ans+1ll*(Sieve_Mu(r)-Sieve_Mu(l-1))*pow3(n/l)%Mod)%Mod;
	}

	printf( "%lld
",(ans+Mod)%Mod );

	return 0;
}

例题：循环之美

求对于十进制下的 (1leq xleq n,1leq yleq m) ，有多少不相等的小数 (dfrac{x}{y}) 在 (k) 进制下是纯循环小数。

(1leq n,mleq 1e9,kleq 2000) .

Solution

以下过程默认 (gcd(x,y)=1) 。（这个很显然吧）

对于一个 (k) 进制下能成为纯循环小数的数，一定存在一个 (t) ，使得 (k^tcdot dfrac{x}{y}) 的小数部分等于 (dfrac xy) ，即

[k^tcdot dfrac{x}{y}-Biglfloor k^tcdot dfrac{x}{y}Big floor =dfrac xy =>k^tx-ycdot Biglfloor k^tcdot dfrac xy Big floor =x\\ ]

于是显然有

[k^txequiv x(mod y)=>k^tequiv 1(mod y)\\ k^tequiv 1^t(mod y),k>0,1>0=>kequiv 1(mod y)\\ ]

所以有：

[gcd(k,y)=1 ]

（这里一开始推到 (y|(k-1)) 去了……我有问题）

这样就可以写出式子了：

[sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1] ]

然后就照常推理：

[egin{aligned} sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1] &=sum_{j=1}^m[gcd(j,k)=1]sum_{i=1}^n[gcd(i,j)=1]\\ &=sum_{j=1}^m[gcd(j,k)=1]sum_{i=1}^nsum_{d|gcd(i,j)}mu(d)\\ &=sum_{d=1}^nmu(d)sum_{j=1}^{lfloor m/d floor}[gcd(jd,k)=1]sum_{i=1}^{lfloor n/d floor}\\ &=sum_{d=1}^nmu(d)Biglfloor frac{n}{d}Big floor sum_{j=1}^{lfloor m/d floor}[gcd(jd,k)=1]\\ end{aligned} ]

显然，这个式子应该尽量表示成两个不相关的式子乘积，注意到 (k) 是个常量，于是有

[sum_{d=1}^nmu(d)Biglfloor frac{n}{d}Big floor [gcd(d,k)=1]sum_{j=1}^{lfloor m/d floor}[gcd(j,k)=1] ]

令

[F(n)=sum_{d=1}^nmu(d)[gcd(d,k)=1],G(n)=sum_{j=1}^{n}[gcd(j,k)=1] ]

（本来这里的 (F) 是带了一个 $Biglfloor dfrac{n}{d}Big floor $ 的项……但是后来发现这是一个巨大的错误……总之很累赘，所以就没有带，反正最后是要整除分块的，这一项没什么用处）

然后就可以分开考虑了。

[egin{aligned} F(n)&=sum_{i=1}^nmu(i)[gcd(i,k)=1]\\ &=sum_{i=1}^nmu(i) sum_{d|gcd(i,k)}mu(d)\\ &=sum_{i=1}^nmu(i) sum_{d|k}mu(d)[d|i]\\ &=sum_{d|k}mu(d)sum_{i=1}^{lfloor n/d floor}mu(id)\\ end{aligned} ]

暴毙了，不会搞了 /kk

……哦，对哦，显然还可以再套一个 ([gcd(i,d)=1]) .不然这个 (mu(id)) 就没了。

这不是更复杂了吗 并不，因为 (mu) 可是积性函数呢 /xyx

[egin{aligned} F(n)&=sum_{d|k}mu(d)sum_{i=1}^{lfloor n/d floor}mu(id)\\ &=sum_{d|k}mu(d)sum_{i=1}^{lfloor n/d floor}mu(i)mu(d)[gcd(i,d)=1]\\ &=sum_{d|k}mu(d)^2sum_{i=1}^{lfloor n/d floor} mu(i)[gcd(i,d)=1]\\ end{aligned} ]

掉线了掉线了，又不会推了 /kel . 虽然这个 (gcd) 还可以拆但是 (d) 的范围也降不下来了啊。

正在重连……连接失败 . 呜呼！

。去看具体数学了。下午再来推这个。

我回来了！我发现上午人傻掉了。

让我们把最初的样子和最后的结果放在一起，带上一个 (k) 的参数：

[F(n,k)=sum_{d=1}^nmu(d)[gcd(d,k)=1]\\ F(n,k)=sum_{d|k}mu(d)^2sum_{i=1}^{lfloor n/d floor} mu(i)[gcd(i,d)=1]\\ ]

发现下面式子的后半部分就是 (F(lfloor frac{n}{d} floor,d)) ！也就是说可以直接递归计算，然后加个记忆化， (F(n)) 就做完了。

边界为 (F(n,1)) ，杜教筛就好了。筛它就完了!

什么？前面一半漏了？(kleq 2000) 啊，随便跑跑就行了。

下面来看 (G(n)) .

[G(n)=sum_{j=1}^{n}[gcd(j,k)=1]=Biglfloor frac{n}{k}Big floor varphi(k) ]

没了！这是好的！

那么答案就是：

[Ans=sum_{d=1}^nmu(d)[gcd(d,k)=1]G(lfloor frac md floor)cdot Biglfloor frac{n}{d}Big floor ]

整除分块就没了 awa.

注意分子和分母不一样，所以不能 swap(n,m) .

//Author: RingweEH
const int N=5e6+10,K=2010;
int pri[N],tot=0,mu[N],pre_mu[N];
int sum_G[N],n,m,k;
bool is[N];

int gcd( int a,int b )
{
	return (b==0) ? a : gcd(b,a%b);
}

void Sieve()
{
	is[1]=1; mu[1]=1;
	for ( int i=2; i<=(N-10); i++ )
	{
		if ( !is[i] ) pri[++tot]=i,mu[i]=-1;
		for ( int j=1; j<=tot && (i*pri[j]<=(N-10)); j++ )
		{
			int pos=i*pri[j]; is[pos]=1;
			if ( i%pri[j]==0 ) { mu[pos]=0; break; }
			mu[pos]=-mu[i];
		}
	}
	for ( int i=1; i<=(N-10); i++ )
		pre_mu[i]=pre_mu[i-1]+mu[i];
	for ( int i=1; i<=k; i++ )
		sum_G[i]=sum_G[i-1]+(gcd(i,k)==1);
}

unordered_map<int,int> sum_Mu,f[K];
int Sieve_Mu( int n )
{
	if ( n<=(N-10) ) return pre_mu[n];
	if ( sum_Mu[n] ) return sum_Mu[n];
	int res=1;
	for ( int l=2,r; l<=n; l=r+1 )
	{
		r=n/(n/l);
		res-=(r-l+1)*Sieve_Mu(n/l);
	}
	return sum_Mu[n]=res;
}

int get_G( int n )
{
	return sum_G[k]*(n/k)+sum_G[n%k];
}

int get_F( int n,int k )
{
	if ( !n ) return 0;
	if ( k==1 ) return f[k][n]=Sieve_Mu(n);
	if ( f[k][n] ) return f[k][n];
	int res=0;
	for ( int i=1; i*i<=k; i++ )
		if ( k%i==0 )
		{
			int x=k/i;
			if ( mu[i] ) res+=mu[i]*mu[i]*get_F(n/i,i);
			if ( mu[x] && (x*x!=k) ) res+=mu[x]*mu[x]*get_F(n/x,x);
		}
	return f[k][n]=res;
}

int main()
{
	n=read(); m=read(); k=read();

	Sieve(); int mn=min( n,m ); ll ans=0;
	for ( int l=1,r; l<=mn; l=r+1 )
	{
		r=min( n/(n/l),m/(m/l) );
		ans+=1ll*(n/l)*get_G(m/l)*( get_F(r,k)-get_F(l-1,k) );
	}
	printf( "%lld
",ans );

    return 0;
}