一些引理
引理1
(forall a,b,c in Z,leftlfloor frac{a}{bc}
ight
floor=leftlfloor frac{leftlfloor frac{a}{b}
ight
floor}{c}
ight
floor) 。
证明:
(frac{a}{b}=leftlfloor frac{a}{b}
ight
floor+r) ((0 leq r < 1))
(leftlfloor frac{a}{bc}
ight
floor=leftlfloor frac{leftlfloor frac{a}{b}
ight
floor}{c} + frac{r}{c}
ight
floor=leftlfloor frac{leftlfloor frac{a}{b}
ight
floor}{c}
ight
floor) 。
引理2
(forall n in N_+,|{ leftlfloor frac{n}{d}
ight
floor | d in N_+ , d leq n }| leq leftlfloor 2sqrt{n}
ight
floor) 。
证明: (leftlfloor frac{n}{d}
ight
floor) 的取值最多有 (leftlfloor 2sqrt{n}
ight
floor) 种。
数论分块
基本形式:
(sum limits_{i=1}^n leftlfloor frac{n}{i}
ight
floor)
对于每一个 (i) ,可以发现: (lfloor frac{n}{i}
floor) 的每一种值呈块状分布。
所以,对于任意一个 (i) ((i leq n)) ,需要找到一个最大的 (j) ((i leq j leq n)) ,使得 (leftlfloor frac{n}{i}
ight
floor=leftlfloor frac{n}{j}
ight
floor) 。
此时,(j=leftlfloor frac{n}{leftlfloor frac{n}{i}
ight
floor}
ight
floor) 。
考虑证明这一结论。
首先,证明 (i leq j leq n) 。
显然, (j leq n) ,所以,只需证明 (i leq j) 。
(leftlfloor frac{n}{leftlfloor frac{n}{i}
ight
floor}
ight
floor geq leftlfloor frac{n}{frac{n}{i}}
ight
floor=leftlfloor i
ight
floor=i) ,
所以 (i leq j) 。
接下来,设 (k=leftlfloor frac{n}{i}
ight
floor) ,考虑证明 (j=leftlfloor frac{n}{k}
ight
floor) 为满足 (leftlfloor frac{n}{j}
ight
floor=k) 的最大值。
(leftlfloor frac{n}{j}
ight
floor =k Rightarrow k leq frac{n}{j} <k+1 Rightarrow frac{1}{k+1} < frac{j}{n} leq frac{1}{k} Rightarrow frac{n}{k+1} <j leq frac{n}{k}) ,
又因为 (j) 是整数,所以 (j_{max}=leftlfloor frac{n}{k}
ight
floor) 。
所以,每次以 ([i,j]) 为一块,直接算即可。
例题
(sum limits_{i=1}^n k mod i=sum limits_{i=1}^n k-i imes leftlfloor frac{k}{i}
ight
floor) ,
直接数论分块算即可。
ans=n*k;
while(l<=n)
{
l=r+1;
if(k/l) r=MIN(k/(k/l),n);
else r=n;
ans-=(k/l)*(l+r)*(r-l+1)/2;
}
二维数论分块
考虑一个二维的版本。
求 (sum limits_{i=1}^{min(n,m)} leftlfloor frac{n}{i}
ight
floor leftlfloor frac{m}{i}
ight
floor) 。
把代码中 r=n/(n/l)
替换为 r=min(n/(n/l),m/(m/l))
即可。
狄利克雷卷积
定义
定义两个数论函数 (f,g) 的狄利克雷卷积为:
((f*g)(n)=sum limits_{d|n}f(d)g(frac{n}{d})) 。
性质
1.交换律:(f*g=g*f) 。
2.结合律:((f*g)*h=f*(g*h)) 。
3.分配律:(f*(g+h)=f*g+f*h) 。
4.(f*varepsilon=f) , (varepsilon) 即狄利克雷卷积的单位元。(其中,(varepsilon (n)=[n=1]))
莫比乌斯函数
定义
设 (n=prod limits _{i=1}^k p_i^{c_i}) ((p_i) 是质数) ,
则 (mu(n)= egin{cases} 1&n=1\ (-1)^k&c_i=1\ 0& ext{otherwise} end{cases}) 。
性质
- (sum_{d|n} mu(d)=varepsilon(n))
- ([gcd(i,j)=1] Longleftrightarrow sum limits_{d|gcd(i,j)} mu(d))
莫比乌斯反演
式子
设 (f(n),g(n)) 为两个数论函数,则
(f(n)=sum_{d|n} g(d) Rightarrow g(n)=sum_{d|n} mu(d) f(frac{n}{d})) ,
(f(n)=sum_{n|d} g(d) Rightarrow g(n)=sum_{n|d} mu(frac{d}{n}) f(d)) 。
证明
以第一个式子为例。
(sum limits_{d|n} mu(d) f(frac{n}{d}) = sum limits_{d|n} mu(d) sum limits_{k|frac{n}{d}} g(k) = sum limits_{k|n} g(k) sum limits_{d|frac{n}{k}} mu(d) = sum limits_{k|n} g(k) [frac{n}{k}=1] = g(n)) 。
例题
1
求 (sum limits_{x=a}^b sum limits_{y=c}^d [gcd(x,y)]=k) 。
容斥后,问题转化为求 (sum limits_{x=1}^n sum limits_{y=1}^m [gcd(x,y)]=k) 。
同时除以 (k) ,得到 (sum limits_{x=1}^{leftlfloor frac{n}{k}
ight
floor} sum limits_{y=1}^{leftlfloor frac{m}{k}
ight
floor} [gcd(i,j)=1]) ,
根据前面莫比乌斯函数的性质,得到 (sum limits_{x=1}^{leftlfloor frac{n}{k}
ight
floor} sum limits_{y=1}^{leftlfloor frac{m}{k}
ight
floor} sum limits_{d|gcd(i,j)} mu(d)) ,
变换求和顺序, (sum limits_{d=1}^{min(n,m)} mu(d) sum limits_{x=1}^{leftlfloor frac{n}{k}
ight
floor} [d|x] sum limits_{y=1}^{leftlfloor frac{m}{k}
ight
floor} [d|y]) ,
(1) ~ (n) 中 (d) 的倍数有 (leftlfloor frac{n}{d}
ight
floor) 个,得到 (sum limits_{d=1}^{min(n,m)} mu(d) leftlfloor frac{n}{kd}
ight
floor leftlfloor frac{m}{kd}
ight
floor) 。
数论分块求解即可。
#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int N=5e4+10;
int cnt,p[N],mu[N],f[N];
int T,a,b,c,d,k;
bool flag[N];
int MIN(int x,int y)
{
return x<y?x:y;
}
void O()
{
mu[1]=1;
for(int i=2;i<N;i++)
{
if(!flag[i]) p[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*p[j]<N;j++)
{
flag[i*p[j]]=1;
if(i%p[j]==0)
{
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
}
for(int i=1;i<N;i++) f[i]=f[i-1]+mu[i];
}
int calc(int n,int m)
{
int l=0,r=0,ans=0;
while(l<=MIN(n,m))
{
l=r+1;
if(n/l==0||m/l==0) r=MIN(n,m);
else r=MIN(n/(n/l),m/(m/l));
ans+=(n/l)*(m/l)*(f[r]-f[l-1]);
}
return ans;
}
signed main()
{
O();
scanf("%lld",&T);
while(T--)
{
scanf("%lld%lld%lld%lld%lld",&a,&b,&c,&d,&k);
printf("%lld
",calc(b/k,d/k)-calc((a-1)/k,d/k)-calc(b/k,(c-1)/k)+calc((a-1)/k,(c-1)/k));
}
return 0;
}
2
求 (sum limits_{i=1}^n lcm(i,n)) 。
原式即 (sum limits_{i=1}^n frac{i imes n}{gcd(i,n)}) 。
将原式复制一遍,颠倒顺序,并且提出 (n) 一项,得 (frac{1}{2} imes (sum limits_{i=1}^{n-1} frac{i imes n}{gcd(i,n)}+sum limits_{i=n-1}^1 frac{i imes n}{gcd(i,n)})+n) ,
根据 (gcd(i,n)=gcd(n-i,n)) ,得 (frac{1}{2} imes (sum limits_{i=1}^{n-1} frac{i imes n}{gcd(i,n)}+sum limits_{i=n-1}^1 frac{i imes n}{gcd(n-i,n)})+n) ,
合并分母相同的项, (frac{1}{2} imes sum limits_{i=1}^{n-1} frac{n^2}{gcd(i,n)}+n) ,
把一个 (n) 合回来, (frac{1}{2} imes sum limits_{i=1}^{n} frac{n^2}{gcd(i,n)}+ frac{n}{2}) 。
考虑枚举 gcd ,要统计 (gcd(i,n)=d) 的个数。
当 (gcd(i,n)=d) 时, (gcd(frac{i}{d},frac{n}{d})=1) ,所以 (gcd(i,n)=d) 的个数为 (varphi(frac{n}{d})) 个。
得到 (frac{1}{2} imes sum limits_{d|n} frac{n^2*varphi(frac{n}{d})}{d}+frac{n}{2}) 。
从 (frac{n^2*varphi(frac{n}{d})}{d}) 中提个 (n) 出来,再设 (d'=frac{n}{d}) ,
得到 (frac{1}{2}n imes (sum limits_{d'|n} d' imes varphi(d')+1)) 。
发现 (sum limits_{d'|n} d' imes varphi(d')) 为积性函数,直接筛即可。
#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int N=1e6+10;
int T,n,f[N],p[N],cnt;
bool flag[N];
void O()
{
flag[0]=1,flag[1]=1,f[1]=1;
for(int i=2;i<N;i++)
{
if(!flag[i]) p[++cnt]=i,f[i]=i*(i-1)+1;
for(int j=1;j<=cnt&&i*p[j]<N;j++)
{
flag[i*p[j]]=1;
if(i%p[j]==0)
{
f[i*p[j]]=f[i]+(f[i]-f[i/p[j]])*p[j]*p[j];
break;
}
f[i*p[j]]=f[i]*f[p[j]];
}
}
}
signed main()
{
O();
scanf("%lld",&T);
while(T--)
{
scanf("%lld",&n);
printf("%lld
",n*(f[n]+1)/2);
}
return 0;
}
3
求 (sum limits_{i=1}^n sum limits_{j=1}^m lcm(i,j)) 。
显然,原式即 (sum limits_{i=1}^n sum limits_{j=1}^m frac{i cdot j}{gcd(i,j)}) 。
枚举 (gcd) ,两数除以 (gcd) 后得到的数互质,得 (sum limits_{k=1}^n k cdot sum limits_{i=1}^{leftlfloor frac{n}{k}
ight
floor} sum limits_{j=1}^{leftlfloor frac{m}{k}
ight
floor} [gcd(i,j)=1] i cdot j) 。
考虑把后面的式子拿出来计算,记 (sum(n,m)=sum limits_{i=1}^n sum limits_{j=1}^m [gcd(i,j)=1] i cdot j) 。
接下来化简这个式子。利用莫比乌斯函数的性质,同时枚举 (gcd) ,得到 (sum limits_{k=1}^n sum limits_{i=1}^n sum limits_{j=1}^m [k|i][k|j]mu(d) cdot i cdot j) 。
把后面的部分除以 (k) ,得到 (sum limits_{k=1}^n mu(k) cdot k^2 cdot sum limits_{i=1}^{leftlfloor frac{n}{k}
ight
floor} sum limits_{j=1}^{leftlfloor frac{m}{k}
ight
floor} i cdot j) 。
发现式子的前半段可以前缀和预处理,后半段可以直接算,合并的部分用数论分块即可。
把 (sum) 代回原式, (sum limits_{k=1}^n k cdot sum(leftlfloor frac{n}{k}
ight
floor,leftlfloor frac{m}{k}
ight
floor)) 。
发现这个式子依然可以用数论分块求解。
#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int N=1e7+10;
const int P=20101009;
int cnt,p[N],mu[N],f[N];
int n,m;
bool flag[N];
int MIN(int x,int y)
{
return x<y?x:y;
}
void O()
{
mu[1]=1;
for(int i=2;i<N;i++)
{
if(!flag[i]) p[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*p[j]<N;j++)
{
flag[i*p[j]]=1;
if(i%p[j]==0)
{
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
}
for(int i=1;i<N;i++) f[i]=(f[i-1]+(mu[i]+P)%P*i%P*i%P)%P;
}
int sum(int n,int m)
{
return (n*(n+1)/2%P)*(m*(m+1)/2%P)%P;
}
int calc(int n,int m)
{
int l=0,r=0,ans=0;
while(l<=MIN(n,m))
{
l=r+1;
if(n/l==0||m/l==0) r=MIN(n,m);
else r=MIN(n/(n/l),m/(m/l));
ans=(ans+(f[r]-f[l-1]+P)%P*sum(n/l,m/l)%P)%P;
}
return ans;
}
int CALC(int n,int m)
{
int l=0,r=0,ans=0;
while(l<=MIN(n,m))
{
l=r+1;
if(n/l==0||m/l==0) r=MIN(n,m);
else r=MIN(n/(n/l),m/(m/l));
ans=(ans+calc(n/l,m/l)*((l+r)*(r-l+1)/2%P)%P)%P;
}
return ans;
}
signed main()
{
O();
scanf("%lld%lld",&n,&m);
printf("%lld",CALC(n,m));
return 0;
}
4
求 (sum limits_{i=1}^n sum limits_{j=1}^m d(i cdot j)) ,其中 (d) 为约数个数函数。
首先, (d) 函数有一个性质, (d(i cdot j)=sum limits_{x|i} sum limits_{y|j} [gcd(x,y)=1]) 。
推一推,
(d(i cdot j)=sum limits_{x|i} sum limits_{y|j} [gcd(x,y)=1])
( =sum limits_{x|i} sum limits_{y|j} sum limits_{k|gcd(x,y)} mu(k))
( =sum limits_{k=1}^{min(i,j)} sum limits_{x|i} sum limits_{y|j} [k|gcd(x,y)] cdot mu(k))
( =sum limits_{k|i,k|j} mu(k) sum limits_{x|i} sum limits_{y|j} [k|gcd(x,y)])
( =sum limits_{k|i,k|j} mu(k) sum limits_{x|leftlfloor frac{i}{k}
ight
floor} sum limits_{y|leftlfloor frac{j}{k}
ight
floor} 1)
( =sum limits_{k|i,k|j} mu(k) d(leftlfloor frac{i}{k}
ight
floor) d(leftlfloor frac{j}{k}
ight
floor)) .
把这个式子代回原式,
( sum limits_{i=1}^n sum limits_{j=1}^m d(i cdot j))
(=sum limits_{i=1}^n sum limits_{j=1}^m sum limits_{k|i,k|j} mu(k) d(leftlfloor frac{i}{k}
ight
floor) d(leftlfloor frac{j}{k}
ight
floor))
(=sum limits_{k=1}^{min(n,m)} sum limits_{i=1}^n sum limits_{j=1}^m [k|i][k|j]mu(k) d(leftlfloor frac{i}{k}
ight
floor) d(leftlfloor frac{j}{k}
ight
floor))
(=sum limits_{k=1}^{min(n,m)} mu(k) sum limits_{i=1}^{leftlfloor frac{n}{k}
ight
floor} d(i) sum limits_{j=1}^{leftlfloor frac{m}{k}
ight
floor} d(j)) .
发现每一部分均可前缀和预处理,用数论分块合并答案即可。
#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int N=5e4+10;
int T,n,m,f[N],F[N];
int cnt,p[N],d[N],t[N],mu[N];
bool flag[N];
int MIN(int x,int y)
{
return x<y?x:y;
}
void O()
{
flag[0]=1,flag[1]=1;
mu[1]=1,d[1]=1;
for(int i=2;i<N;i++)
{
if(!flag[i]) p[++cnt]=i,mu[i]=-1,d[i]=2,t[i]=1;
for(int j=1;j<=cnt&&i*p[j]<N;j++)
{
flag[i*p[j]]=1;
if(i%p[j]==0)
{
mu[i*p[j]]=0;
d[i*p[j]]=d[i]/(t[i]+1)*(t[i]+2);
t[i*p[j]]=t[i]+1;
break;
}
mu[i*p[j]]=-mu[i],d[i*p[j]]=d[i]*2,t[i*p[j]]=1;
}
}
for(int i=1;i<N;i++) f[i]=f[i-1]+mu[i];
for(int i=1;i<N;i++) F[i]=F[i-1]+d[i];
}
int calc(int n,int m)
{
int l=0,r=0,ans=0;
while(l<=MIN(n,m))
{
l=r+1;
if(n/l==0||m/l==0) r=MIN(n,m);
else r=MIN(n/(n/l),m/(m/l));
ans+=(f[r]-f[l-1])*F[n/l]*F[m/l];
}
return ans;
}
signed main()
{
O();
scanf("%lld",&T);
while(T--)
{
scanf("%lld%lld",&n,&m);
printf("%lld
",calc(n,m));
}
return 0;
}