有这样一段两两求最大公约数的程序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) o (x+y,y)
\
(x,y) o (x,x+y)
]
不难发现这样的操作与取模操作的逆运算一一对应（一个是让(k=k+1)，另一个是换方向(k=1)）。
并且这样找出来的二元组不会重复。
考虑用二元组(x,y)表示一个数(ax+(a+b)y)。显然合法的(x,y)是所有(gcd(x,y)=1)。
并且每一个((x,y))只会作为最大值出现2次。
问题即求：
[egin{aligned}
sum _{x,y} [gcd(x,y)=1][ax+(a+b)yle m]
\
=sum _d mu(d) sum _{x,y} [ax+(a+b)y le lfloorfrac {m}{d}
floor]
\
end{aligned}
]
杜教筛 + 类欧几里得即可。
复杂度(O(Tsqrt mlog 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
",ans);
}
}