Alatúlië,Heldo!

莫比乌斯反演学习笔记

Warning

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

0. 前置-整除分块

问题

给定一个 \(n\) ,求:

\[\sum_{i=1}^n\Big\lfloor \frac{n}{i}\Big\rfloor \]

思路&代码

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

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

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

\[n/l=k,n=r\times 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_1\sim a_n\) ,求最大的 \(d\) ,满足

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

展开取模,得到

\[\begin{aligned} \sum_{i=1}^nd-((a_i-1)-d\times \Big\lfloor\frac{a_i-1}{d}\Big\rfloor+1)\leq k &\iff \sum d\times\Big\lfloor\frac{a_i-1}{d}\Big\rfloor-(a_i-1)-1\leq k-n\times d\\\\ & \iff \sum d\times\Big\lfloor \frac{a_i-1}{d}\Big\rfloor \leq k-n\times d+\sum a_i\\\\ & \iff \sum \Big\lfloor \frac{a_i-1}{d}\Big\rfloor\leq \Big\lfloor\dfrac{k-n\times d+\sum a_i}{d}\Big\rfloor\\\\ & \iff \sum \Big\lfloor \frac{a_i-1}{d}\Big\rfloor\leq \Big\lfloor\dfrac{k+\sum a_i}{d}\Big\rfloor-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\n",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}) \]

常见的积性函数

  1. 单位函数。\(\varepsilon(x)=[x=1]\)
  2. 常数函数。\(1(x)=1\)
  3. 单位(?)函数。\(ID(x)=x\)
  4. 欧拉函数。\(\varphi(x)=\sum_{i=1}^x[\gcd(x,i)=1]\)
  5. 莫比乌斯函数

扩展:

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

Proof

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

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

\(m>1\) 时,假设只要 \(m_1\bot m_2\)\(m_1m_2<m\) ,这个性质都成立。

显然,当 \(m_1\bot 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_1\bot d_2 \]

根据归纳假设,只有在 \(d_1=m_1\)\(d_2=m_2\) 的情形下才可能不成立。那么把这一部分拆出来,得到

\[\begin{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)=\sum\limits_{d|x}d\\ \varepsilon(x)=(\mu*1)(x)=\sum\limits_{d|x}\mu(d)\\ \varphi(x)=(\mu*ID)(x)=\sum\limits_{d|x}d\cdot\mu(\frac xd) \]

1. 莫比乌斯函数

定义

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

另一种表述

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

  1. \(x=1\)\(\mu(x)=1\)
  2. \(\forall i,c_i=1\)\(\mu(x)=(-1)^k\)
  3. \(otherwise,\mu(x)=0\)

性质

  1. \(\mu*1=\varepsilon\) ,即 \(\sum_{d|x}\mu(d)=\varepsilon(x)\) ,也就是 \(\sum_{d|x}\mu(d)=1\) 当且仅当 \(x=1\) ,否则 \(\sum_{d|x} \mu(d)=0\) .

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

  3. \[[\gcd(i,j)=1]=\varepsilon(\gcd(i,j))=\sum_{d|\gcd(i,j)}\mu(d) \]

    反演结论。

  4. \[\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)\times \mu(\frac{x}{d}). \]

另一种表述方式(根据Dirichlet卷积):

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

Proof

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

3. 例题

注:以下默认 \(n\leq m\) .

Problem b

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

Solution

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

\[\begin{aligned} \sum_{x=a}^b\sum_{y=c}^d[\gcd(x,y)=k]&=\sum_{x=1}^b\sum_{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} \]

问题就转化为求

\[\begin{aligned} f(n,m)= \sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)=k] &=\sum_{x=1}^{\lfloor n/k\rfloor}\sum_{y=1}^{\lfloor m/k\rfloor}[\gcd(x,y)=1]\\ \because [\gcd(i,j)=1]=\sum_{d|\gcd(i,j)}\mu(d) \therefore & =\sum_{x=1}^{\lfloor n/k\rfloor}\sum_{y=1}^{\lfloor m/k\rfloor}\sum_{d|\gcd(x,y)}\mu(d)\\ &=\sum_{d=1}^{\lfloor n/k\rfloor} \mu(d)\sum_{x=1}^{\lfloor n/k\rfloor}\sum_{y=1}^{\lfloor m/k\rfloor}[d|\gcd(x,y)]\\ \because 1\sim n中d的倍数有 \lfloor n/d\rfloor 个\therefore &=\sum_{d=1}^{\lfloor n/k\rfloor} \mu(d)\Big\lfloor \frac{n}{kd}\Big\rfloor \Big\lfloor \frac{m}{kd}\Big\rfloor \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\n",ans );
	}
    return 0;
}

YY的GCD

给定 \(N,M(N,M\leq 1e7)\) ,求

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

多测,\(T\leq 1e4\) .

Solution

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

\(t=pd\) ,有

\[\begin{aligned} \sum_{p\in primes}\sum_{d=1}^{\lfloor N/p\rfloor}\mu(d)\Big\lfloor \frac{N}{pd}\Big\rfloor \Big\lfloor \frac{M}{pd}\Big\rfloor &=\sum_{t=1}^N\sum_{p\in primes,p|t}\mu(\frac{t}{p})\Big\lfloor \frac{N}{t}\Big\rfloor \Big\lfloor \frac{M}{t}\Big\rfloor\\ \end{aligned} \]

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

那么所求就是

\[\sum_{x=1}^nf(x)\Big\lfloor \frac{N}{x}\Big\rfloor \Big\lfloor \frac{M}{x}\Big\rfloor \]

\(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\n",ans );
    }
    return 0;
}

LCM Sum

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

Solution

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

\[\begin{aligned} \sum_{i=1}^n\text{lcm}(i,n) &=\sum_{i=1}^n \dfrac{i\times n}{\gcd(i,n)}\\ &=\sum_{d|n}\sum_{i=1}^n[\gcd(i,n)=d]\frac{i}{d}\times n\\ &=\sum_{d|n}\sum_{i=1}^{n/d}[\gcd(i,n/d)=1]\times i\times n\\ &=n\times\sum_{d|n}\sum_{i=1}^{n/d}[\gcd(i,n/d)=1]i\\ \end{aligned} \]

\(d\) 代换 \(n/d\)

\[\begin{aligned} & =n\times \sum_{d|n}\sum_{i=1}^d[\gcd(i,d)=1]\times i \end{aligned} \]

感觉推不下去了……吗?

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

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

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

特殊情况:当 \(d=1\) 时,和为 \(1\) . 因此,

\[\begin{aligned} f(n)=n\times \sum_{d|n}\sum_{i=1}^d[\gcd(i,d)=1]\times i &=n\times \sum_{d|n}\frac{\varphi(d)\times 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\n",ans );
	}
    return 0;
}

Crash的数字表格 / JZPTAB

\[\sum_{i=1}^n\sum_{j=1}^m\text{lcm}(i,j)\bmod 20101009 \]

Solution

\[\begin{split} \sum\limits_{i=1}^n\sum\limits_{j=1}^m\operatorname{lcm}(i,j)=&\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{ij}{\gcd(i,j)}\\ =&\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{ij}{d}[\gcd(i,j)=d]\\ =&\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor n/d\rfloor}\sum\limits_{j=1}^{\lfloor m/d\rfloor}ijd[\gcd(i,j)=1]\\ =&\sum\limits_{d=1}^n d\sum\limits_{i=1}^{\lfloor n/d\rfloor}i\sum\limits_{j=1}^{\lfloor m/d\rfloor}j\sum\limits_{k|\gcd(i,j)}\mu(k)\\ =&\sum\limits_{d=1}^n d\sum\limits_{k=1}^n\mu(k)\sum\limits_{i=1}^{\lfloor n/d\rfloor}i[k|i]\sum\limits_{j=1}^{\lfloor m/d\rfloor}j[k|j]\\ =&\sum\limits_{d=1}^n d\sum\limits_{k=1}^n\mu(k)\sum\limits_{i=1}^{\lfloor n/dk\rfloor}ik\sum\limits_{j=1}^{\lfloor m/dk\rfloor}jk\\ =&\sum\limits_{d=1}^n d\sum\limits_{k=1}^nk^2\mu(k)\frac{\lfloor n/dk\rfloor(\lfloor n/dk\rfloor+1)}{2}\cdot\frac{\lfloor m/dk\rfloor(\lfloor m/dk\rfloor+1)}{2}\\ \end{split} \]

\(x=dk\) 带入得

\[\sum\limits_{x=1}^nx\cdot\frac{\lfloor\frac{n}{x}\rfloor(\lfloor\frac{n}{x}\rfloor+1)}{2}\cdot\frac{\lfloor\frac{m}{x}\rfloor(\lfloor\frac{m}{x}\rfloor+1)}{2}\sum\limits_{k|x}k\mu(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\n",ans );
    return 0;
}

约数个数和

\(d(x)\)\(x\) 的约数个数,求

\[\sum_{i=1}^n\sum_{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\times j\) 就是:

\[i\times j=\prod p_k^{a_k+b_k} \]

考虑如何枚举 \(p_k\) 的所有幂次,设当前枚举的幂次为 \(c_k\) .

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

  • \(c_k\leq 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.

由引理,有:

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^md(ij) & =\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]\\ 交换求和顺序,&=\sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)=1]\sum_{i=1}^n\sum_{j=1}^m[x|i,y|j]\\ &=\sum_{x=1}^n\sum_{y=1}^m\sum_{d|\gcd(x,y)}\mu(d)\Big\lfloor \frac{n}{x}\Big\rfloor \Big\lfloor \frac{m}{y}\Big\rfloor\\ &=\sum_{d=1}^n\mu(d)\sum_{x=1}^n [d|x]\sum_{y=1}^m[d|y]\Big\lfloor \frac{n}{x}\Big\rfloor \Big\lfloor \frac{m}{y}\Big\rfloor\\ &=\sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor n/d\rfloor}\sum_{y=1}^{\lfloor m/d\rfloor}\Big\lfloor \frac{n}{xd}\Big\rfloor \Big\lfloor \frac{m}{yd}\Big\rfloor\\ &=\sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor n/d\rfloor}\Big\lfloor \frac{n}{xd}\Big\rfloor\sum_{y=1}^{\lfloor m/d\rfloor} \Big\lfloor \frac{m}{yd}\Big\rfloor\\ \end{aligned} \]

\(T_1=xd,T_2=yd\) ,有

\[\sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor n/d\rfloor}\Big\lfloor \frac{n}{xd}\Big\rfloor\sum_{y=1}^{\lfloor m/d\rfloor} \Big\lfloor \frac{m}{yd}\Big\rfloor=\sum_{d=1}^n\mu(d)\sum_{T_1=1}^n\Big\lfloor \frac{n}{T_1}\Big\rfloor\sum_{T_2=1}^m \Big\lfloor \frac{m}{T_2}\Big\rfloor \]

后面两项都是标准的整除分块式子,第一项就是 \(\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\n",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;
        }
    }
}

线性筛欧拉函数

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

性质:对于一个数 \(i\) ,如果 \(pri[j]|i\) ,那么有 \(phi[i\times pri[j]]=phi[i]\times pri[j]\) .

Proof

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

由此可得,\(x,x+i,\dots,x+i\times (pri[j]-1)\) 都与 \(i\) 互质,那么与 \(i\times pri[j]\) 互质的至少有 \(phi[i]\times 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\times 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\times pri[j])=f(i)\times 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]\times pri[j])=1\)

    \(f(i\times pri[j])=f(i/low[i])\times f(low[i]\times 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}^n\sum_{j=1}^m \gcd(i,j)^k\bmod 1e9+7 \]

思路:不妨设 \(n\leq m\) .

\[\sum_{i=1}^n\sum_{j=1}^m \gcd(i,j)^k\bmod 1e9+7 =\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^nd^k[\gcd(i,j)=1]\\ =\sum_{d=1}^nd^k\sum_{i=1}^{\lfloor n/d\rfloor}\sum_{j=1}^{\lfloor n/d\rfloor}\sum_{x|i,x|j}\mu(x)\\ =\sum_{d=1}^nd^k\sum_{x=1}^{\lfloor n/d\rfloor}\mu(x)\lfloor \frac{n}{dx}\rfloor\lfloor\frac{m}{dx}\rfloor \]

\(T=dx\) ,有:

\[\sum_{T=1}^n\lfloor n/T\rfloor \lfloor m/T\rfloor\sum_{d|T}d^k\mu(T/d) \]

\(g(T)=\sum_{d|T}d^k\mu(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\n",(ans+mod)%mod );
    }
}

学习材料

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

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

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

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

posted @ 2020-12-21 12:02  RingweEH  阅读(267)  评论(0编辑  收藏  举报