先描述一下问题
已知m,n
求∑(i,j)=k1<=i<=m,1<=j<=n1 \sum_{(i,j)=k}^{1<=i<=m,1<=j<=n}1(i,j)=k∑1<=i<=m,1<=j<=n1
用文字描述即为已知整数n,m,
1<=i<=m,1<=j<=n,求有多少对(i,j)满足i,j的最大公约数为k
算法1
首先最简单的算法即为暴力枚举
枚举所有符合条件的数对,判断是否满足要求即可
时间复杂度o(mn)o(mn)o(mn)
算法1的优化
显然原式等价于∑(i,j)=11<=i<=m/k,1<=j<=n/k1 \sum_{(i,j)=1}^{1<=i<=m/k,1<=j<=n/k}1(i,j)=1∑1<=i<=m/k,1<=j<=n/k1
这样时间复杂度为o(mn/k2)o(mn/k^2)o(mn/k2)
算法2
显然这个时间复杂度是难以接受的
那么我们可以怎么做呢?
先考虑另一个问题,
有多少对(i,j)有公因子k
这个显然是[mk][nk][ {\frac m k} ][{\frac n k }][km][kn]
那么这跟正确答案相差多少?
显然这个数包含了最大公约数为2k,3k,4k…的情况,
那么我们减去这些数就好了
为保证无后效性,我们不得不采用倒推,即先推出n/k*k,再依次往下,直到枚举到k
由考虑最差情况,调和级数可知 时间复杂度为O(nlogn)
例题:
【NOI2010】能量采集
https://www.luogu.org/problemnew/show/P1447
在分析清楚思路之后
代码异常简洁
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
using namespace std;
#define ll long long
#define N 100007
ll n, m, ans, f[N];
int main()
{
cin >> n >> m;
if (n > m)
swap(m, n);
for (int i = n; i >= 1; i--)
{
f[i] = (m / i) * (n / i);
for (int j = i * 2; j <= n; j += i)
{
f[i] -= f[j];
}
ans += f[i] * i;
}
cout << 2 * ans - n * m;
return 0;
}
算法3
莫比乌斯反演
由算术基本定理知,任意正整数都可以拆分为多个质数的乘积
pip_ipi表示质数
则任意正整数n=p1a1p2a2...pqaqn = p_1^{a_1}p_2^{a_2}...p_q^{a_q}n=p1a1p2a2...pqaq
引入莫比乌斯函数μ(n)={1,n=1(−1)r,n=p1p2p3..pr,pi为互不相同的素数0,其他情况μ(n)= \begin {cases} 1,n = 1 \\(-1)^r, n=p_1p_2p_3..p_r,p_i为互不相同的素数\\0,其他情况\end{cases}μ(n)=⎩⎪⎨⎪⎧1,n=1(−1)r,n=p1p2p3..pr,pi为互不相同的素数0,其他情况
根据mobius反转定理
如果F(n)=∑d∣nf(n),则有f(n)=∑d∣nμ(d)F(nd)如果F(n)=\sum_{d|n}f(n),则有f(n)=\sum_{d|n}μ(d)F(\frac n d)如果F(n)=d∣n∑f(n),则有f(n)=d∣n∑μ(d)F(dn)
一种等价形式
如果F(n)=∑n∣df(d),则有f(n)=∑n∣dμ(d)F(dn)如果F(n)=\sum_{n|d}f(d),则有f(n)=\sum_{n|d}μ(d)F(\frac d n)如果F(n)=n∣d∑f(d),则有f(n)=n∣d∑μ(d)F(nd)
使用后一种形式
设F(x)为k|(i,j)的对数
那么F(d)=[nk][mk]=∑d∣if(i)F(d)=[\frac n k][\frac m k]=\sum_{d|i}f(i)F(d)=[kn][km]=∑d∣if(i)
反演后得f(d)=∑d∣iμ(i)F(id)f(d)=\sum_{d|i}μ(i)F(\frac i d)f(d)=∑d∣iμ(i)F(di)
于是时间复杂度就变成O(n)了
算法4
仔细观察发现(可以打表)
[nk][\frac n k][kn]只有n\sqrt nn种取值
再利用μ(i)的前缀和
可以优化时间复杂度至O(n)O(\sqrt n)O(n)
(预处理O(n))
例题
[HAOI2011]Problem b
https://www.luogu.org/problemnew/solution/P2522
题目描述
对于给出的n个询问,每次求有多少个数对(x,y),满足a≤x≤b,c≤y≤d,且gcd(x,y) = k,gcd(x,y)函数为x和y的最大公约数。
100%的数据满足:1≤n≤50000,1≤a≤b≤50000,1≤c≤d≤50000,1≤k≤50000
套用算法4,再利用容斥原理
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
using namespace std;
#define N 50000
#define ll long long
int t,k,a,b,c,d,mu[N+6],sumu[N+7],prime[N+6],tot;
bool notprime[N+6];
void init()
{
sumu[1]=mu[1]=1;
notprime[1]=1;
for(int i = 2; i <= N; i ++ )
{
if(!notprime[i])
{
prime[++tot]=i;
mu[i]=-1;
}
sumu[i]=sumu[i-1]+mu[i];
for(int j = 1;i*prime[j]<=N&&j <= tot; j ++)
{
notprime[i*prime[j]]=1;
if(i%prime[j]!=0)
{
mu[i*prime[j]]=-mu[i];
}
else break;
}
}
}
ll f(int m,int n)
{
ll ans = 0;
int last ;
if(m>n)swap(m,n);
for(int i = 1 ; i <= m; i =last +1)
{
last = min(n/(n/i),m/(m/i));
ans += ((ll)sumu[last]-sumu[i-1])*(n/i)*(m/i);
}
return ans;
}
int main()
{
init();
scanf("%d",&t);
while(t--)
{
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
a--;c--;
a/=k;b/=k;c/=k;d/=k;
printf("%lld\n",f(b,d)-f(a,d)-f(b,c)+f(a,c));
}
}