莫比乌斯反演学习笔记

Warning

由于式子较多，LaTeX 渲染可能会有点慢qwq

0. 前置-整除分块

问题

给定一个 (n) ，求：

[sum_{i=1}^nBiglfloor frac{n}{i}Big floor ]

思路&代码

这个式子显然可以 (mathcal{O}(n)) 计算，考虑优化。

对于小的数可能感觉不出来；但是对于较大的数，比如：(100/51,100/52,dots,100/100) ，你会发现它们全都是一样的！这说明其实刚才的过程中，有很大一部分的计算完全可以一次性解决！

现在来考虑如何计算左右端点。设当前左右端点为 (l,r) ，那么显然有：

[n/l=k,n=r imes k,=>r=n/(n/l). ]

于是就可以直接得到右端点了。容易得到这样的复杂度是 (mathcal{O}(sqrt n)) .

for(int l=1,r;l<=n;l=r+1)
{
	r=n/(n/l);
	ans+=(r-l+1)*(n/l);
}

习题

CF830-C Bamboo Partition

给定 (n) 个数 (a_1sim a_n) ，求最大的 (d) ，满足

[sum_{i=1}^n d-((a_i-1)mod d+1)leq k. ]

展开取模，得到

[egin{aligned} sum_{i=1}^nd-((a_i-1)-d imes Biglfloorfrac{a_i-1}{d}Big floor+1)leq k &iff sum d imesBiglfloorfrac{a_i-1}{d}Big floor-(a_i-1)-1leq k-n imes d\\ & iff sum d imesBiglfloor frac{a_i-1}{d}Big floor leq k-n imes d+sum a_i\\ & iff sum Biglfloor frac{a_i-1}{d}Big floorleq Biglfloordfrac{k-n imes d+sum a_i}{d}Big floor\\ & iff sum Biglfloor frac{a_i-1}{d}Big floorleq Biglfloordfrac{k+sum a_i}{d}Big floor-n\\ end{aligned} ]

于是发现了一个整除分块的基本式，可以直接求了。

//Author:RingweEH
const int N=110;
int n;
ll a[N],k,ans=0;
int main()
{
	n=read(); k=read(); ll mx=0;
	for ( int i=1; i<=n; i++ )
		a[i]=read()-1,k+=a[i]+1,mx=max( mx,a[i] );
	for ( ll l=1,r=1; r<=mx; l=r+1 )
	{
		ll sum=0; r=(1ll<<50);
		for ( int i=1; i<=n; i++ )
			if ( a[i]>=l ) r=min( r,a[i]/(a[i]/l) ),sum+=a[i]/l;
		ll tmp=k/(sum+n);
		if ( l<=tmp ) ans=max( ans,min( tmp,r ) );
	}
	if ( k/n>mx ) ans=max( ans,k/n );
	printf( "%lld
",ans );
	return 0;
}

0. 前置-积性函数

定义：如果 (gcd(x,y)=1) 且 (f(xy)=f(x)f(y)) ，那么 (f(n)) 为积性函数。

性质：如果 (f(x),g(x)) 均为积性函数，那么如下函数也是积性函数：

[h(x)=f(x^p)\ h(x)=f^p(x)\ h(x)=f(x)g(x)\ h(x)=sum_{d|x} f(d)g(frac{x}{d}) ]

常见的积性函数：

单位函数。(varepsilon(x)=[x=1])
常数函数。(1(x)=1)
单位（？）函数。(ID(x)=x)
欧拉函数。(varphi(x)=sum_{i=1}^x[gcd(x,i)=1])
莫比乌斯函数

扩展：

关于积性函数 (g(m)=sum_{d|m}f(d)) （如果满足 (g) 是积性那么 (f) 也是）性质的证明

Proof

采用数学归纳法证明这个性质。

首先，对于 (m=1) ，显然有 (g(1)=f(1)=1) .

当 (m>1) 时，假设只要 (m_1ot m_2) 且 (m_1m_2<m) ，这个性质都成立。

显然，当 (m_1ot m_2) 时，(m_1,m_2) 的因子也互质。

那么就有：

[g(m_1m_2)=sum_{d|m_1m_2}f(d)=sum_{d_1|m_1}sum_{d_2|m_2}f(d_1d_2),d_1ot d_2 ]
根据归纳假设，只有在 (d_1=m_1) 且 (d_2=m_2) 的情形下才可能不成立。那么把这一部分拆出来，得到

[egin{aligned} g(m_1m_2) &=sum_{d_1|m_1}f(d_1)sum_{d_2|m_2}f(d_2)-f(m_1)f(m_2)+f(m_1m_2)\\ &=g(m_1)g(m_2)-f(m_1)f(m_2)+f(m_1m_2) end{aligned} ]
所以有

[g(m_1m_2)=g(m_1)g(m_2) ]
即 (g) 是积性函数，从而 (f) 也是积性函数。

0. 前置-Dirichlet卷积

定义：两个数论函数 (f,g) 的Dirichlet卷积 ((f*g)) 为：

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

（会发现这个式子和上面积性函数的最后一个式子很像……）

性质：满足交换律和结合律，(varepsilon) 是该函数单位元。

举例：

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

1. 莫比乌斯函数

定义：

[mu(x)= egin{cases} 1 & x=1\ 0 & exists dinmathbb{Z}:d^2mid x\ (-1)^k & k 为 x 本质不同的的质因子个数\ end{cases} ]

另一种表述：

用算数基本定理把 (x) 写成质因数分解的形式得到 (x=prod_{i=1}^k p_i^{c_i}) ，那么有：

(x=1) ，(mu(x)=1)
(forall i,c_i=1) ，(mu(x)=(-1)^k)
(otherwise,mu(x)=0)

性质：

(mu*1=varepsilon) ，即 (sum_{d|x}mu(d)=varepsilon(x)) ，也就是 (sum_{d|x}mu(d)=1) 当且仅当 (x=1) ，否则 (sum_{d|x} mu(d)=0) .
[sum_{d|n}frac{mu(d)}{d}=frac{varphi(n)}{n} ]
[[gcd(i,j)=1]=varepsilon(gcd(i,j))=sum_{d|gcd(i,j)}mu(d) ]
反演结论。
[varphi*1=ID,sum_{d|nm}=sum_{x|n}sum_{y|m}[gcd(x,y)=1] ]
拓展。

int mu[N],pri[N],tot,is[N];
void sieve()
{
    is[1]=mu[1]=1;
    for ( int i=2; i<=n; i++ )
    {
        if ( !is[i] ) pri[++tot]=i,mu[i]=-1; //根据定义，质数的莫比乌斯函数显然是-1
        for ( int j=1; j<=tot && i*pri[j]<=n; j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j] ) mu[i*pri[j]]=-mu[i];  //多了一个质因子，取反
            else { mu[i*pri[j]]=0;break; }   //出现平方因子，置为0
        }
    }
}

2. 莫比乌斯反演

对于两个数论函数 (f,g) ，满足

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

那么有

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

另一种表述方式（根据Dirichlet卷积）：

[f=g*1=>g=f*mu ]

Proof

[egin{aligned} f*mu & =g*1*mu\ & =g*varepsilon（莫比乌斯函数性质1）\ & =g（单位函数） end{aligned} ]

3. 例题

注：以下默认 (nleq m) .

Problem b

求

[sum_{x=a}^bsum_{y=c}^d[gcd(x,y)=k] ]

Solution

原式有上下界限制，不好求，一个简单的想法是考虑容斥。

[egin{aligned} sum_{x=a}^bsum_{y=c}^d[gcd(x,y)=k]&=sum_{x=1}^bsum_{y=1}^d[gcd(x,y)=k]-sum_{x=1}^{a-1}sum_{y=1}^d[gcd(x,y)=k]-sum_{x=1}^{b}sum_{y=1}^{c-1}[gcd(x,y)=k]+sum_{x=1}^{a-1}sum_{y=1}^{c-1}[gcd(x,y)=k]\ end{aligned} ]

问题就转化为求

[egin{aligned} f(n,m)= sum_{x=1}^nsum_{y=1}^m[gcd(x,y)=k] &=sum_{x=1}^{lfloor n/k floor}sum_{y=1}^{lfloor m/k floor}[gcd(x,y)=1]\ ecause [gcd(i,j)=1]=sum_{d|gcd(i,j)}mu(d) herefore & =sum_{x=1}^{lfloor n/k floor}sum_{y=1}^{lfloor m/k floor}sum_{d|gcd(x,y)}mu(d)\ &=sum_{d=1}^{lfloor n/k floor} mu(d)sum_{x=1}^{lfloor n/k floor}sum_{y=1}^{lfloor m/k floor}[d|gcd(x,y)]\ ecause 1sim n中d的倍数有 lfloor n/d floor 个 herefore &=sum_{d=1}^{lfloor n/k floor} mu(d)Biglfloor frac{n}{kd}Big floor Biglfloor frac{m}{kd}Big floor end{aligned} ]

对 (mu(x)) 函数求前缀和，然后整除分块即可。

//Author: RingweEH
const int N=5e4+10;
int n,m,k,mu[N],sum[N],pri[N],tot=0;
bool is[N];
void init()
{
	is[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++ )
		{
			if ( i%pri[j] ) is[i*pri[j]]=1,mu[i*pri[j]]=-mu[i];
			else { mu[i*pri[j]]=0; is[i*pri[j]]=1; break; }
		}
	}
	sum[0]=0;
	for ( int i=1; i<=(N-10); i++ )
		sum[i]=sum[i-1]+mu[i];
}
int get_Fenk( int n,int m )
{
	int res=0;
	if ( n>m ) swap( n,m );
	for ( int l=1,r=0; l<=n; l=r+1 )
	{
		r=min( n/(n/l),m/(m/l) );
		res=res+(sum[r]-sum[l-1])*(n/l)*(m/l);
	}
	return res;
}
int main()
{
	int T=read(); init();
	while ( T-- )
	{
		int a=read(),b=read(),c=read(),d=read(),k=read();
		a--; c--; a/=k,b/=k,c/=k,d/=k;
		int ans=get_Fenk(b,d)-get_Fenk(a,d)-get_Fenk(b,c)+get_Fenk(a,c);
		printf( "%d
",ans );
	}
    return 0;
}

YY的GCD

给定 (N,M(N,Mleq 1e7)) ，求

[sum_{x=1}^Nsum_{y=1}^M[gcd(x,y)in primes] ]

多测，(Tleq 1e4) .

Solution

[egin{aligned} sum_{x=1}^Nsum_{y=1}^M[gcd(x,y)in primes] & =sum_{pin primes}sum_{x=1}^Nsum_{y=1}^M[gcd(x,y)=p]\ & =sum_{pin primes}sum_{x=1}^{lfloor N/p floor}sum_{y=1}^{lfloor M/p floor}[gcd(x,y)=1]\ & =sum_{pin primes}sum_{x=1}^{lfloor N/p floor}sum_{y=1}^{lfloor M/p floor}sum_{d|gcd(x,y)}mu(d)\ &=sum_{pin primes}sum_{d=1}^{lfloor N/p floor}mu(d)sum_{x=1}^{lfloor N/p floor}sum_{y=1}^{lfloor M/p floor}[d|gcd(x,y)]\ &=sum_{pin primes}sum_{d=1}^{lfloor N/p floor}mu(d)Biglfloor frac{N}{pd}Big floor Biglfloor frac{M}{pd}Big floor\ end{aligned} ]

令 (t=pd) ，有

[egin{aligned} sum_{pin primes}sum_{d=1}^{lfloor N/p floor}mu(d)Biglfloor frac{N}{pd}Big floor Biglfloor frac{M}{pd}Big floor &=sum_{t=1}^Nsum_{pin primes,p|t}mu(frac{t}{p})Biglfloor frac{N}{t}Big floor Biglfloor frac{M}{t}Big floor\ end{aligned} ]

令

[f(x)=sum_{pin primes,p|x}mu(frac{x}{p}) ]

那么所求就是

[sum_{x=1}^nf(x)Biglfloor frac{N}{x}Big floor Biglfloor frac{M}{x}Big floor ]

把 (f(x)) 及其前缀和预处理出来，然后整除分块即可。

~~怎么一模一样啊~~ 基本就是把上一题里面的 (k) 换成不确定的，多了一个求和而已。

//Author: RingweEH
const int N=1e7+10;
int n,m,mu[N],pri[N],tot=0,f[N];
ll sum[N];
bool is[N];
void init()
{
    is[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++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j] ) mu[i*pri[j]]=-mu[i];
            else { mu[i*pri[j]]=0; break; }
        }
    }
    for ( int i=1; i<=tot; i++ )
        for ( int j=pri[i]; j<=N-10; j+=pri[i] )
            f[j]+=mu[j/pri[i]];
    sum[0]=0;
    for ( int i=1; i<=N-10; i++ )
        sum[i]=sum[i-1]+f[i];
}
int main()
{
    int T=read(); init();
    while ( T-- )
    {
        n=read(); m=read();
        if ( n>m ) swap( n,m );
        ll ans=0;
        for ( int l=1,r=0; l<=n; l=r+1 )
        {
            r=min( n/(n/l),m/(m/l) );
            ll tmp=sum[r]-sum[l-1];
            ans+=tmp*(n/l)*(m/l);
        }
        printf( "%lld
",ans );
    }
    return 0;
}

LCM Sum

求

[sum_{i=1}^n ext{lcm}(i,n) ]

Solution

显然 ( ext{lcm}) 可以转化成 (gcd) .

[egin{aligned} sum_{i=1}^n ext{lcm}(i,n) &=sum_{i=1}^n dfrac{i imes n}{gcd(i,n)}\ &=sum_{d|n}sum_{i=1}^n[gcd(i,n)=d]frac{i}{d} imes n\ &=sum_{d|n}sum_{i=1}^{n/d}[gcd(i,n/d)=1] imes i imes n\ &=n imessum_{d|n}sum_{i=1}^{n/d}[gcd(i,n/d)=1]i\ end{aligned} ]

用 (d) 代换 (n/d) 得

[egin{aligned} & =n imes sum_{d|n}sum_{i=1}^d[gcd(i,d)=1] imes i end{aligned} ]

感觉推不下去了……吗？

会发现一个奇妙的东西：(sum_{i=1}^d[gcd(i,d)=1] imes i) 就是 ([1,d]) 中与 (d) 互质的数的和！~~但还是不会求~~

注意到如果 (gcd(i,d)=1) ，那么有 (gcd(d-i,d)=1) .也就是说，这样的数成对出现，且每一对的和为 (d) .所以有

[sum_{i=1}^d[gcd(i,d)=1] imes i=frac{varphi(d) imes d}{2} ]

特殊情况：当 (d=1) 时，和为 (1) . 因此，

[egin{aligned} f(n)=n imes sum_{d|n}sum_{i=1}^d[gcd(i,d)=1] imes i &=n imes sum_{d|n}frac{varphi(d) imes d}{2} end{aligned} ]

//Author: RingweEH
const int N=1e6+10;
int n,phi[N],pri[N],tot=0;
ll f[N];
bool is[N];
void init()
{
	is[1]=phi[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 && (pri[j]*i)<=(N-10); j++ )
		{
			if ( i%pri[j] ) is[i*pri[j]]=1,phi[i*pri[j]]=phi[i]*phi[pri[j]];
			else { is[i*pri[j]]=1; phi[i*pri[j]]=phi[i]*pri[j]; break; }
		}
	}
	f[1]=1;
	for ( int i=2; i<=N-10; i++ )
		f[i]=1ll*phi[i]*i/2;
	for ( int j=1; j<=tot; j++ )
		for ( int i=1; i*pri[j]<=N-10; i++ )
			f[i*pri[j]]+=f[i];
}
int main()
{
	int T=read(); init();
	while ( T-- )
	{
		int n=read();
		ll ans=f[n]*n;
		printf( "%lld
",ans );
	}
    return 0;
}

Crash的数字表格 / JZPTAB

求

[sum_{i=1}^nsum_{j=1}^m ext{lcm}(i,j)mod 20101009 ]

Solution

[egin{split} sumlimits_{i=1}^nsumlimits_{j=1}^moperatorname{lcm}(i,j)=&sumlimits_{i=1}^nsumlimits_{j=1}^mfrac{ij}{gcd(i,j)}\ =&sumlimits_{d=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^mfrac{ij}{d}[gcd(i,j)=d]\ =&sumlimits_{d=1}^nsumlimits_{i=1}^{lfloor n/d floor}sumlimits_{j=1}^{lfloor m/d floor}ijd[gcd(i,j)=1]\ =&sumlimits_{d=1}^n dsumlimits_{i=1}^{lfloor n/d floor}isumlimits_{j=1}^{lfloor m/d floor}jsumlimits_{k|gcd(i,j)}mu(k)\ =&sumlimits_{d=1}^n dsumlimits_{k=1}^nmu(k)sumlimits_{i=1}^{lfloor n/d floor}i[k|i]sumlimits_{j=1}^{lfloor m/d floor}j[k|j]\ =&sumlimits_{d=1}^n dsumlimits_{k=1}^nmu(k)sumlimits_{i=1}^{lfloor n/dk floor}iksumlimits_{j=1}^{lfloor m/dk floor}jk\ =&sumlimits_{d=1}^n dsumlimits_{k=1}^nk^2mu(k)frac{lfloor n/dk floor(lfloor n/dk floor+1)}{2}cdotfrac{lfloor m/dk floor(lfloor m/dk floor+1)}{2}\ end{split} ]

将 (x=dk) 带入得

[sumlimits_{x=1}^nxcdotfrac{lfloorfrac{n}{x} floor(lfloorfrac{n}{x} floor+1)}{2}cdotfrac{lfloorfrac{m}{x} floor(lfloorfrac{m}{x} floor+1)}{2}sumlimits_{k|x}kmu(k) ]

跟上题一样求即可。

//Author: RingweEH
const int N=1e7+10;
const ll Mod=20101009;
int n,m,mu[N],pri[N],tot=0;
ll f[N];
bool is[N];
void init()
{
    is[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++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j] ) mu[i*pri[j]]=-mu[i];
            else { mu[i*pri[j]]=0; break; }
        }
    }
    f[1]=1;
    for ( int i=2; i<=(N-10); i++ )
        f[i]=(1ll*mu[i]*i+Mod)%Mod;
    for ( int j=1; j<=tot; j++ )
        for ( int i=1; i*pri[j]<=(N-10); i++ )
            f[i*pri[j]]=(f[i*pri[j]]+f[i])%Mod;
}
ll func( int x )
{
    ll res=1ll*x*f[x]%Mod;
    (res*=1ll*(n/x+1)*(n/x)/2%Mod)%=Mod;
    (res*=1ll*(m/x+1)*(m/x)/2%Mod)%=Mod;
    return res;
}
int main()
{
    init();
    n=read(); m=read();
    if ( n>m ) swap( n,m );
    ll ans=0;
    for ( int i=1; i<=n; i++ )
        ans=(ans+func(i))%Mod;
    printf( "%lld
",ans );
    return 0;
}

约数个数和

设 (d(x)) 为 (x) 的约数个数，求

[sum_{i=1}^nsum_{j=1}^md(ij) ]

Solution

引理

[d(ij)=sum_{x|i}sum_{y|j}[gcd(x,y)=1] ]

Proof

一个很显然的想法是，枚举 (i,j) 中各自的因子，然后相乘得到 (ij) 的一个因子。但是这样会算重。

考虑如何计数。首先，用算数基本定理表示 (i,j) ：

[i=prod p_k^{a_k}\ j=prod p_k^{b_k}\ ]
那么 (i imes j) 就是：

[i imes j=prod p_k^{a_k+b_k} ]
考虑如何枚举 (p_k) 的所有幂次，设当前枚举的幂次为 (c_k) .

显然，我们可以 对无序的枚举钦定一个顺序 ，那么有如下情况：

(c_kleq a_k) 我们假定对于每个约数，都优先选择 (a) 中的部分。那么这里即是枚举 (a) 中所选的幂次。

(c_k>a_k) 这时候必须用到 (b) ，说明 (a) 中默认已经选满了，直接枚举的是 (b) 中的幂次。

显然，在这样的讨论下，无论何时 (a,b) 都不可能同时 被枚举 。

设当前 (i,j) 中枚举的因子分别为 (x|i,y|j) ，那么就相当于满足 (gcd(x,y)=1) .

因此，有等式：

[d(ij)=sum_{x|i}sum_{y|j}[gcd(x,y)=1] ]
Q.E.D.

由引理，有：

[egin{aligned} sum_{i=1}^nsum_{j=1}^md(ij) & =sum_{i=1}^nsum_{j=1}^msum_{x|i}sum_{y|j}[gcd(x,y)=1]\ 交换求和顺序，&=sum_{x=1}^nsum_{y=1}^m[gcd(x,y)=1]sum_{i=1}^nsum_{j=1}^m[x|i,y|j]\ &=sum_{x=1}^nsum_{y=1}^msum_{d|gcd(x,y)}mu(d)Biglfloor frac{n}{x}Big floor Biglfloor frac{m}{y}Big floor\ &=sum_{d=1}^nmu(d)sum_{x=1}^n [d|x]sum_{y=1}^m[d|y]Biglfloor frac{n}{x}Big floor Biglfloor frac{m}{y}Big floor\ &=sum_{d=1}^nmu(d)sum_{x=1}^{lfloor n/d floor}sum_{y=1}^{lfloor m/d floor}Biglfloor frac{n}{xd}Big floor Biglfloor frac{m}{yd}Big floor\ &=sum_{d=1}^nmu(d)sum_{x=1}^{lfloor n/d floor}Biglfloor frac{n}{xd}Big floorsum_{y=1}^{lfloor m/d floor} Biglfloor frac{m}{yd}Big floor\ end{aligned} ]

令 (T_1=xd,T_2=yd) ，有

[sum_{d=1}^nmu(d)sum_{x=1}^{lfloor n/d floor}Biglfloor frac{n}{xd}Big floorsum_{y=1}^{lfloor m/d floor} Biglfloor frac{m}{yd}Big floor=sum_{d=1}^nmu(d)sum_{T_1=1}^nBiglfloor frac{n}{T_1}Big floorsum_{T_2=1}^m Biglfloor frac{m}{T_2}Big floor ]

后面两项都是标准的整除分块式子，第一项就是 (mu) 函数的前缀和，直接求就好了。

//Author: RingweEH
const int N=5e4+10;
int n,m,mu[N],pri[N],tot=0,sum[N];
ll f[N];
bool is[N];

void init()
{
    is[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++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j] ) mu[i*pri[j]]=mu[i]*mu[pri[j]];
            else { mu[i*pri[j]]=0; break;}
        }
    }
    sum[0]=0;
    for ( int i=1; i<=(N-10); i++ )
        sum[i]=sum[i-1]+mu[i];
    memset( f,0,sizeof(f) );
    for ( int i=1; i<=(N-10); i++ )
    	for ( int l=1,r=0; l<=i; l=r+1 )
    	{
    		r=i/(i/l);
    		f[i]=f[i]+1ll*(r-l+1)*(i/l);
    	}
}

int main()
{
    int T=read(); init();
    while ( T-- )
    {
        n=read(); m=read();
        if ( n>m ) swap( n,m );
        ll ans=0;
        for ( int l=1,r=0; l<=n; l=r+1 )
        {
        	r=min( n/(n/l),m/(m/l) );
        	ans+=(sum[r]-sum[l-1])*f[n/l]*f[m/l];
        }
        printf( "%lld
",ans );
    }
}

4. 积性函数与线性筛

听说 所有积性函数都能线性筛 ，但是时间复杂度不能保证…… Link

线性筛质数

最基础的一种，保证每个质数只会被最小质因子筛掉。

int pri[N],tot,is[N];     //is[i]为1的表示不是质数
void sieve()
{
    is[1]=1;
    for ( int i=2; i<=n; i++ )
    {
        if ( !is[i] ) pri[++tot]=i;
        for ( int j=1; j<=tot && i*pri[j]<=n; j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j]==0 ) break;
        }
    }
}

线性筛欧拉函数

定义：(1sim n) 中与 (n) 互质的数的个数。

性质：对于一个数 (i) ，如果 (pri[j]|i) ，那么有 (phi[i imes pri[j]]=phi[i] imes pri[j]) .

Proof

如果一个数 (x) 与 (i) 互质，那么显然，(x+i) 依然与 (i) 互质。

由此可得，(x,x+i,dots,x+i imes (pri[j]-1)) 都与 (i) 互质，那么与 (i imes pri[j]) 互质的至少有 (phi[i] imes pre[j]) 个。而与 (i) 不互质的数 (y) ，加上 (i) 也不会与 (i) 互质，同理可得，不会出现新的与 (i) 互质的数，得证。

int phi[N],pri[N],tot,is[N];
void sieve()
{
    is[1]=phi[1]=1;
    for ( int i=2; i<=n; i++ )
    {
        if ( !is[i] ) pri[++tot]=i,phi[i]=i-1;  //质数除了本身都互质
        for ( int j=1; j<=tot && i*pri[j]<=n; j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j] ) phi[i*pri[j]]=phi[i]*phi[pri[j]];  //积性函数定义
            else { phi[i*pri[j]]=phi[i]*pri[j]; break; }  //上述性质
        }
    }
}

线性筛约数个数

令 (d[i]) 为 (i) 的约数个数。

将每个数 (x) 表示成算术基本定理的形式：(x=prod p_i^{a_i}) ，那么考虑每个质因子所取的幂次，有

[d[x]=prod_{i=1}^k(a_i+1) ]

维护每个数最小质因子的出现次数（(a_1) ）即可。

int d[N],a[N],pri[N],tot,is[N];
void sieve()
{
    is[1]=d[1]=1;
    for ( int i=2; i<=n; i++ )
    {
        if ( !is[i] ) pri[++tot]=i,d[i]=2,a[i]=1;   //质数的约数个数为2
        for ( int j=1; j<=tot && i*pri[j]<=n; j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j] ) d[i*pri[j]]=d[i]*d[pri[j]],a[i*pri[j]]=1;
            //新的质因数，i*pri[j]的最小质因数为pri[j]，幂次为1
            else { d[i*pri[j]]=d[i]/(a[i]+1)*(a[i]+2); a[i*pri[j]]=a[i]+1 ;break; }
            //最小质因数 pri[j]，也就是 a[i] 的幂次又多了1，除去之前乘的再乘上新的幂次
        }
    }
}

线性筛约数和

记 (sigma(i)) 表示 (i) 的约数和。考虑每个质因数取的幂次，同样有：（这个展开应该就能理解了吧）

[sigma(i)=prod_{i=1}^k(sum_{j=0}^{a_i}p_i^j) ]

令 (low(i)) 表示 (i) 的最小质因数的指数幂，即 (p_1^{a_1}) ；(sum(i)) 表示 (i) 的最小质因数对答案的贡献，(sum_{j=0}^{a_1}p_1^j) .

（小心爆 int /xyx）

int low[N],sum[N],sigma[N],pri[N],tot,is[N];
void sieve()
{
    is[1]=low[1]=sum[1]=sigma[1]=1;
    for ( int i=2; i<=n; i++ )
    {
        if ( !is[i] ) low[i]=pri[++tot]=i,sum[i]=sigma[i]=i+1;
        for ( int j=1; j<=tot && i*pri[j]<=n; j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j]==0 ) 
            {
                low[i*pri[j]]=low[i]*pri[j];
                sum[i*pri[j]]=sum[i]+low[i*pri[j]];
                sigma[i*pri[j]]=sigma[i]/sum[i]*sum[i*pri[j]];
                break;
            }
            //不能整除，也就是 pri[j] 是 i*pri[j] 的最小质因数
            low[i*pri[j]]=pri[j];
            sum[i*pri[j]]=pri[j]+1;
            sigma[i*pri[j]]=sigma[i]*sigma[pri[j]];
        }
    }
}

线性筛一般积性函数

前提：能快速计算出 (f(1),f(prime),f(prime^k)) .（也就是质因数个数小于等于 (1) 的所有数的函数值）

计算 :假设已经完成了上述值的计算。显然，含有至少两个质因数的数一定可以被分解成两个互质且不为 1 的数的乘积。根据积性函数的定义，可以用 (f(xy)=f(x)f(y)) 来得出相应的函数值。

现在来考虑具体怎么筛。

思考线性筛的之所以线性的原因，是因为 每个数都只被最小质因子筛了一次 ，也就是每个形如 (i imes pri[j]) 的数会在 (i) 的时候被 (i) 乘上 (pri[j]) 筛掉，且若将 (i) 唯一分解，有 (pri[j]leq p_1) .（显然，否则在 (p_1) 的时候就会 break 掉）

分类讨论：

(pri[j]<p_1)

由于 (gcd(pri[j],i)=1) ，直接计算即可 (f(i imes pri[j])=f(i) imes f(pri[j])) .
(pri[j]=p_1)

这时候需要记录 (i) 中最小质因子的幂次，也就是 (low[i]=p_1^{a_1}) .

显然，(i/low[i]) 中的最小质因数一定大于 (pri[j]) ，那么 (gcd(i/low[i],low[i] imes pri[j])=1) ，

有 (f(i imes pri[j])=f(i/low[i]) imes f(low[i] imes pri[j])) .

特别地，当 (low[i]=i) 时，(i) 只有一个质因数，要利用前面的前提特殊计算。

void sieve()
{
    is[1]=low[1]=1; f[1]=对1直接定义;
    for ( ll i=2; i<=n; i++ )
    {
        if ( !is[i] ) low[i]=pri[++tot]=i,f[i]=对质数直接定义;
        for ( ll j=1; j<=tot && i*pri[j]<=n; j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j]==0 )
            {
                low[i*pri[j]]=low[i]*pri[j];
                if ( low[i]==i ) f[i*pri[j]]=对质数的若干次幂进行定义(一般由f[i]递推);
                else f[i*pri[j]]=f[i/low[i]]*f[low[i]*pri[j]];
                break;
            }
            low[i*pri[j]]=pri[j]; f[i*pri[j]]=f[i]*f[pri[j]];
        }
    }
}

What's more

对于某些形如Dirichlet卷积形式，但 (f(x)) 或 (g(x)) 不是积性函数，数据范围较小可以暴力筛，枚举一个 (d) 计算对哪些 (x) 做贡献，较大的情况下可以考虑相关性质。

例题：P4449 于神之怒加强版

题意：求

[sum_{i=1}^nsum_{j=1}^m gcd(i,j)^kmod 1e9+7 ]

思路：不妨设 (nleq m) .

[sum_{i=1}^nsum_{j=1}^m gcd(i,j)^kmod 1e9+7 =sum_{i=1}^nsum_{j=1}^msum_{d=1}^nd^k[gcd(i,j)=1]\ =sum_{d=1}^nd^ksum_{i=1}^{lfloor n/d floor}sum_{j=1}^{lfloor n/d floor}sum_{x|i,x|j}mu(x)\ =sum_{d=1}^nd^ksum_{x=1}^{lfloor n/d floor}mu(x)lfloor frac{n}{dx} floorlfloorfrac{m}{dx} floor ]

令 (T=dx) ，有：

[sum_{T=1}^nlfloor n/T floor lfloor m/T floorsum_{d|T}d^kmu(T/d) ]

令 (g(T)=sum_{d|T}d^kmu(T/d)) ,预处理 (g(T)) 和前缀和（线筛），分块即可。

const int inf=0x7f7f7f7f,N=5e6+10,mod=1e9+7,M=5e6;
int pri[N],u[N],g[N],f[N],T,k,tot,n,m;
bool vis[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 init()
{
    f[1]=1;
    for ( int i=2; i<=M; i++ )
    {
        if ( !vis[i] )
        {
            pri[++tot]=i; g[tot]=power( i,k ); f[i]=(g[tot]-1+mod)%mod;
        }
        for ( int j=1; j<=tot && i*pri[j]<=M; j++ )
        {
            int tmp=i*pri[j];
            vis[tmp]=1;
            if ( i%pri[j]==0 )
            {
                f[tmp]=1ll*f[i]*g[j]%mod; break;
            }
            f[tmp]=1ll*f[i]*f[pri[j]]%mod;
        }
    }
    for ( int i=1; i<=M; i++ )
        f[i]=(f[i]+f[i-1])%mod;
}
int main()
{
    scanf( "%d%d",&T,&k );
    init();
    while ( T-- )
    {
        scanf( "%d%d",&n,&m ); 
        int maxx=min( n,m ),pos,ans=0;
        for ( int i=1; i<=maxx; i=pos+1 )
        {
            pos=min( n/(n/i),m/(m/i) );
            ans=(ans+1ll*(f[pos]-f[i-1]+mod)%mod*(n/i)%mod*(m/i)%mod)%mod;
        }
        printf( "%d
",(ans+mod)%mod );
    }
}

学习材料

窝tcl，年底才学 神George1123 年初就会的东西qaq

顺手丢一个 Baidu文库的莫反PPT。

致谢以下所有良心博主。感觉是全网的精华，而且几乎没什么好补充的了。

租酥雨-积性函数与线性筛 George1123-莫比乌斯反演学习笔记