算法总结 给定范围内最大公约数为某一定值的数对个数的算法

先描述一下问题
已知m,n
求$\sum_{(i,j)=k}^{1<=i<=m,1<=j<=n}1$
用文字描述即为已知整数n,m，
1<=i<=m,1<=j<=n,求有多少对(i,j)满足i,j的最大公约数为k

算法1

首先最简单的算法即为暴力枚举
枚举所有符合条件的数对，判断是否满足要求即可
时间复杂度$o(mn)$

算法1的优化

显然原式等价于$\sum_{(i,j)=1}^{1&lt;=i&lt;=m/k,1&lt;=j&lt;=n/k}1$
这样时间复杂度为$o(mn/k^2)$

算法2

显然这个时间复杂度是难以接受的
那么我们可以怎么做呢？
先考虑另一个问题，
有多少对（i,j)有公因子k
这个显然是$[ {\frac m k} ][{\frac n k }]$
那么这跟正确答案相差多少？
显然这个数包含了最大公约数为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

莫比乌斯反演
由算术基本定理知，任意正整数都可以拆分为多个质数的乘积
$p_i$表示质数
则任意正整数$n = p_1^{a_1}p_2^{a_2}...p_q^{a_q}$
引入莫比乌斯函数$μ(n)= \begin {cases} 1,n = 1 \\(-1)^r, n=p_1p_2p_3..p_r,p_i为互不相同的素数\\0,其他情况\end{cases}$
根据mobius反转定理
$如果F(n)=\sum_{d|n}f(n)，则有f(n)=\sum_{d|n}μ(d)F(\frac n d)$
一种等价形式
$如果F(n)=\sum_{n|d}f(d)，则有f(n)=\sum_{n|d}μ(d)F(\frac d n)$
使用后一种形式
设F(x)为k|(i,j)的对数
那么$F(d)=[\frac n k][\frac m k]=\sum_{d|i}f(i)$
反演后得$f(d)=\sum_{d|i}μ(i)F(\frac i d)$
于是时间复杂度就变成O(n)了

算法4

仔细观察发现（可以打表）
$[\frac n k]$只有$\sqrt n$种取值
再利用μ(i)的前缀和
可以优化时间复杂度至$O(\sqrt 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));
    }
}