有这样一段两两求最大公约数的程序CoGcd，

int Gcd(int x, int y){
   if(y == 0)return x;
   return Gcd(y, x % y);
}

void CoGcd(int m){
   for(int i = 1; i <= m; i++)
       for(int j = 1; j <= m; j++)
           Gcd(i, j);
}

给出m的值，进行t次询问，每次询问包含一对xi，yi。针对每次询问，输出整个程序执行过程当中，Gcd(xi, yi)被执行了多少次。

例如：\(m = 20\)，

\(Gcd(8,5)\)会被执行4次，对应的x, y值是

(8,5) (5,8) (13,8) (8,13)，这4组x y，在调用Gcd时，会递归执行Gcd(8, 5)。

sol

考虑二元组\((a,b)\)，\(a>b\)，在一次操作之后是该二元组一定是\((ka+b,a)\)的形式。

考虑二元组\((a,a+b)\)，之后每次操作形如：
\[ （x,y) \to (x+y,y) \\ （x,y) \to (x,x+y) \]
不难发现这样的操作与取模操作的逆运算一一对应（一个是让\(k=k+1\)，另一个是换方向\(k=1\)）。

并且这样找出来的二元组不会重复。

考虑用二元组(x,y)表示一个数\(ax+(a+b)y\)。显然合法的\(x,y\)是所有\(\gcd(x,y)=1\)。

并且每一个\((x,y)\)只会作为最大值出现2次。

问题即求：
\[ \begin{aligned} \sum _{x,y} [\gcd(x,y)=1][ax+(a+b)y\le m] \\ =\sum _d \mu(d) \sum _{x,y} [ax+(a+b)y \le \lfloor\frac {m}{d}\rfloor] \\ \end{aligned} \]

杜教筛 + 类欧几里得即可。

复杂度\(O(T\sqrt m\log m + m^\frac 23)\)

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e6+6;
int mu[N], prime[N], pcnt, vis[N];
inline void sieve(int n = 5e6){
    mu[1] = 1;
    for(int i=2;i<=n;i++){
        if(!vis[i])prime[++pcnt]=i,mu[i]=-1;
        for(int j=1;j<=pcnt&&1ll*i*prime[j]<=n;j++){
            vis[i*prime[j]] = 1;
            if(i%prime[j])mu[i*prime[j]] = -mu[i];
            else {
                mu[i*prime[j]]=0;
                break;
            }
        }
    }
    for(int i=1;i<=n;i++)mu[i] += mu[i-1];
}
ll query(ll n, ll a, ll b, ll c) {
    ll x = a / c, y = b / c;
    if (!a) {
        return y * (n + 1);
    }else if (a >= c || b >= c) {
        return query(n, a % c, b % c, c) + n * (n + 1) / 2 * x + (n + 1) * y;
    } else {
        ll m = (1ll * a * n + b) / c;
        return n * m - query(m - 1, c, c - b - 1, a);
    }
}
map<int,int> MU;
inline int q1(int n){
    if(n<=5000000)return mu[n];
    if(MU.find(n)!=MU.end())return MU[n];
    int ans = 1;
    for(int u=2,v;u<=n;u=v+1){
        v=n/(n/u);
        ans -= (v-u+1)*q1(n/u);
    }
    return MU[n]=ans;
}
inline ll q(int a,int b,int n){
    int d = n/a;
    int t = n%a;
    return query(d-1,a,t,b);
}
ll ans = 0;
int t,m;
inline void find(int x,int y){
    if(x>m||y>m)return;
    ans+=2;
    find(x+y,y),find(x,y+x);
    
}

int main()
{
    cin >> t >> m;
    sieve();
    while(t--){
        ans = 0;
        int a,b;
        scanf("%d%d",&a,&b);
        if(a<=b){puts("1");continue;}
        b+=a;
        if(b>m){puts("2");continue;}
        ans = 4;
        for(int u=1,v;u<=m;u=v+1){
            v=m/(m/u);
            int f1 = q1(v) - q1(u-1);
            ans += 4ll*f1*q(a,b,m/u);
        }
        printf("%lld\n",ans);
    }
}