Part0:符号约定
([p]):艾弗森记号.对于命题(p),当(p)成立时,([p])为(1),否则为(0).
(a|b):(b)被(a)整除.即(exist kin)使得(b=ka).
(aperp b):(a)与(b)互质,即(gcd(a,b)=1).特殊地,(1perp 1)也成立.
本文中的函数,未经特殊说明,均指数论函数,即定义在正整数集(N^+)上的一类函数.本文中(n,m,t,i,j,k)默认表示某个正整数,(p)默认为某个素数.
Part1:整除分块
我们的目标是求数论函数
的值.显然,我们可以在(O(n))时间内朴素地求出该函数的值.但是,画出(leftlfloorfrac ni ight floor)的图像可知,该函数有许多值是一样的,并且呈块状分布.并且更一般地,每一个块的最后一个数就是(leftlfloor n/(leftlfloor frac ni ight floor) ight floor)
于是,我们就有一种(O(sqrt n))的算法:
( ext{CALC_DIVIDE_SUM}(n):)
(resleftarrow 0)
(lleftarrow 1)
(rleftarrow 0)
(mathbf{while} lle n:)
(quad rleftarrow leftlfloor n/(leftlfloor frac nl
ight
floor)
ight
floor)
(quad resleftarrow res+(r-l+1)cdot leftlfloorfrac nl
ight
floorleftlfloorfrac ml
ight
floor)
(quad lleftarrow r+1)
(mathbf{return} res)
C++实现如下:
inline int divide_sum(int n)
{
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ans+=(r-l+1)*(n/l);
}
return ans;
}
当然,我们有时候是求函数的和,这时候就要将r-l+1替换为函数在这一段区间上的和.
如果我们要求
也可设计出以下(O(sqrt n+sqrt m))的算法:
( ext{CALC_DIVIDE_SUM2}(n,m):)
(resleftarrow 0)
(lleftarrow 1)
(rleftarrow 0)
(mathbf{while} lle min(n,m):)
(quad rleftarrow min(leftlfloor n/(leftlfloor frac nl
ight
floor)
ight
floor,leftlfloor m/(leftlfloor frac ml
ight
floor)
ight
floor))
(quad resleftarrow res+(r-l+1)cdot leftlfloorfrac nl
ight
floorleftlfloorfrac ml
ight
floor)
(quad lleftarrow r+1)
(mathbf{return} res)
C++实现如下:
inline int divide_sum(int n,int m)
{
for(int l=1,r;l<=std::min(n,m);l=r+1)
{
r=std::min(n/(n/l),m/(m/l));
ans+=(r-l+1)*(n/l)*(m/l);
}
return ans;
}
Part2:狄利克雷卷积
定义两个数论函数(f,g)的狄利克雷卷积(Dirichlet convolution)为:
狄氏卷积有以下性质:
(交换律)(1.f*g=g*f);
(结合律)(2.f*g*h=(f*g)*h=f*(g*h)),这是因为
(分配律)(3.(f+g)*h=f*h+g*h);
(数乘律)(4.(mf)*g=f*(mg)=m(f*g));
(单位元)(5.e*f=f),其中
(逆元)(6.)对于每个(f(1) e0)的函数(f),都存在一个函数(g)使得(f*g=e),记作(g=f^{-1}).
从上述结论可以看出,所有(f(1) e0)的函数对狄氏卷积构成一个群.
那么,如何求出一个函数的逆呢?
事实上,有
因为
Part3:积性函数
若一个数论函数(f)满足:当(nperp m)时,有(f(nm)=f(n)f(m)),则称之为积性函数(multiplicative function).
一些常见的积性函数有:(e(n)=[n=1]; ext{id}(n)=n; ext{id}^k(n)=n^k;1(n)= ext{id}^0(n)=1).事实上,这些函数是完全积性函数(completelty multiplicative function),即不论(n,m)是否互质,都有(f(nm)=f(n)f(m)).
另外,(n)的因数个数(sigma_0(n))和欧拉函数(varphi(n))((1,2,dots,n)中与(n)互质的数的个数)也是积性函数.证明如下:
引理 如果(nperp m),我们发现每个(nm)的约数(t)都可以分解成一个(n)的约数(gcd(n,t))和一个(m)的约数(gcd(m,t))的积,并且这种分解是一一对应的.
这是因为(nperp m,t|nmRightarrow t=gcd(t,nm)=gcd(t,n)gcd(t,m)),而且(a|n,b|mRightarrow ab|nm),并且,(a|n,b|nRightarrowgcd(ab,n)=a).
于是(sigma_0(nm)=sigma_0(n)sigma_0(m)).
同样的,(tperp nmLeftrightarrow tperp n,tperp mLeftrightarrow(tmod n)perp n,(tmod m)perp m),所以每个(1,2,dots,nm)之间的与(nm)互质的数(t)都可以对应到一个(1,2,dots,n)的与(n)互质的数(tmod n)和一个(1,2,dots,m)的与(m)互质的数(tmod m).
并且根据中国剩余定理,这种对应是一一对应的(即已知(aperp n,bperp m)后可以唯一确定一个(1,2,dots,nm)之间的(t)使得(tmod n=a, tmod m=b),且(tperp nm)).因此(varphi(nm)=varphi(n)varphi(m)).
下面我们来证明两个结论:
(1.)两个积性函数的狄氏卷积是积性函数.
由于上面的引理,又因为(nperp m,a|n,b|mRightarrow aperp b),于是若(nperp m),则
(2.)积性函数的逆是积性函数.
对(nm)的大小运用数学归纳法.设(f^{-1}=g),则
(1.)当(nm=1)时,(g(1)=1),结论成立;
(2.)假设(nm>1),当(n'm'<nm)时有(g(n'm')=g(n')g(m')),
Part4:积性函数的线性筛
根据算数基本定理,由于每个正整数(n)可以写成惟一的分解形式(n=prodlimits_{i=1}^{m}p_i^{alpha_i}),所以
所以,一个积性函数在素数幂处的取值决定了整个积性函数.并且,由于在线性筛素数时我们可以顺便求出每个数的最小质因数(p_1)以及最小质因数的次数(alpha_1)以及(n/p_1^{alpha_1}),就可以利用递推式(f_n=f(p_1^{alpha_1})f(n/p_1^{alpha_1}))直接计算(f)了.
另外,由于(sigma_0)和(varphi)在素数幂处的取值容易得到:(sigma_0(p^k)=k+1,varphi(p^k)=p^{k-1}(p-1)),所以
因(varphi*1)也是积性函数,又((varphi*1)(p^k)=p^k),故有( ext{id}=varphi*1).
Part5:莫比乌斯函数与莫比乌斯反演
我们定义(1)的逆是(mu),称为莫比乌斯函数(Mobius function).这样的话,若(g=f*1),则就有(f=g*mu),更一般地,
或等价地,
上述公式称为莫比乌斯反演公式(Mobius inversion formula).或者类似地,我们可以发现( ext{id}^k)的逆是(f(n)=mu(n)n^k),于是若
就有
比如,可得(varphi=mu* ext{id}),即
那么,如何求(mu)呢?
由于(1)是积性的,而(mu)是(1)的逆,所以(mu)也是积性的.容易计算出:
于是,
莫比乌斯函数有以下性质:
(1.)对于任意的正整数(n),有
(2.)对于任意的正整数(n),有
运用上面的性质(1),我们容易得出(mu)的线性筛代码:
( ext{LINEAR_SIEVE_MU}(n):)
(muleftarrow new array filled with 0)
(mathbf{for} ileftarrow 1 mathbf{to} n:)
(quad tarleftarrow [i=1]//tar ext{ is the sum of Mobius functions so far})
(quad deltaleftarrow tar-mu[i]//delta ext{ is the differences between the functions and the sum})
(quad mu[i]leftarrow delta)
(quad jleftarrow 2i)
(quad mathbf{while} jle n:)
(quadquad mu[j]leftarrow mu[j]+delta)
(quadquad jleftarrow j+i)
(mathbf{return} mu)
C++实现如下:
void sieve_mu(int n,int mu[]) //mu是储存结果的数组
{
for(int i=1;i<=n;++i)
{
int tar=(i==1?1:0);
int delta=tar-mu[i];
mu[i]=delta;
for(int j=2*i;j<=n;j+=i)
mu[j]+=delta;
}
}
我们也可以在筛素数的同时筛出(mu):
( ext{LINEAR_SIEVE_PRIME_AND_MU}(n):)
(mu,pleftarrow new array filled with 0)
(vis new array filled with mathbf{FALSE}//vis ext{ is used to label numbers whether has divisors or not})
(mu[1]leftarrow 1)
(mathbf{for} ileftarrow 2 mathbf{to} n:)
(quad mathbf{if} vis[i]=mathbf{FALSE}:)
(quadquad p.append(i))
(quadquad mu[i]leftarrow-1)
(quad mathbf{for} each prime p_0 mathbf{in} p:)
(quadquad vis[p_0cdot i]leftarrow mathbf{TRUE})
(quadquad mathbf{if} p|i:)
(quadquadquad mathbf{break})
(quadquad mathbf{else}:)
(quadquadquad mu[icdot p_0]leftarrow-mu[i])
(mathbf{return} mu,p)
C++实现如下:
void sieve_mu(int n,int mu[],int p[],int vis[],int& cnt)
{
memset(vis,0,sizeof(vis));
mu[1]=1;
cnt=0;
for(int i=2;i<=n;++i)
{
if(!vis[i])
{
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt&&p[j]*i<=n;++j)
{
vis[p[j]*i]=1;
if(i%p[j]==0)
break;
else
mu[i*p[j]]=-mu[i];
}
}
}
Part6:简单例题
LG P2257 YYY的GCD
对于(n,m),我们有(d|gcd(a,b)Leftrightarrow d|a,d|b).于是,对一个关于(gcd(a,b))的函数(f(gcd(a,b))),如果有(f(n)=sumlimits_{d|n}g(d))的话,就可以写成(sumlimits_{d|a,d|b}g(d))了.
考虑令(f(n)=[n是素数]).我们要找到一个(g)使得(f=1*g),那么(g=mu*f),
于是我们可以这样求出(g):
( ext{SIEVE_G}(n):)
(gleftarrow new array filled with 0)
(mu,pleftarrow ext{LINEAR_SIEVE_PRIME_AND_MU}(n))
(mathbf{for} each p_0 mathbf{in} p:)
(quad ileftarrow 1)
(quad mathbf{while} icdot p_0le n:)
(quadquad g[icdot p]leftarrow g[icdot p]+mu[i])
(mathbf{return} g)
复杂度为(O(sumlimits_{ple n}(n/p))=O(nloglog n)).
我们要求的是
先预处理(g)的前缀和,然后就可以用整除分块之后(O(sqrt n+sqrt m))回答每个询问.C++实现如下:
const int Maxn=1e7+5,MaxLogn=15;
int n,m,cnt,p[Maxn/MaxLogn],mu[Maxn],g[Maxn];
bool vis[Maxn];
void sieve()
{
mu[1]=vis[1]=1;
for(int i=2;i<Maxn;++i)
{
if(!vis[i])
{
p[++cnt]=i;
mu[i]=-1;
g[i]=1;
}
for(int j=1;j<=cnt&&i*p[j]<Maxn;++j)
{
int x=i*p[j];
vis[x]=1;
if(i%p[j]==0)
{
g[x]=mu[i];
mu[x]=0;
break;
}
else
{
g[x]=-g[i]+mu[i];
mu[x]=-mu[i];
}
}
g[i]+=g[i-1];
}
}
int main()
{
sieve();
int T;
scanf("%d",&T);
for(;T--;)
{
scanf("%d%d",&n,&m);
long long ans=0;
if(n>m)
swap(n,m);
for(int l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans+=1LL*(g[r]-g[l-1])*(n/l)*(m/l);
}
printf("%lld
",ans);
}
}