先描述一下問題
已知m,n
求 $\underset{(i,j)=}{\overset{}{\sum }}$

演算法1

首先最簡單的演算法即為暴力列舉
列舉所有符合條件的數對，判斷是否滿足要求即可
時間複雜度 $o(mn)$

演算法1的優化

顯然原式等價於 $\sum_{(i,j)=1}^{1<=i<=m/k,1<=j<=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 ++)
        {