[51nod2583]數論只會Gcd

假設詢問的是\((x,y),x>y\)。

考慮建一棵樹，每個點代表一個二元組\((a,b)(a>b,a \le m)\)。如果\((a,b)\)進行一次操作后會變成\((a',b')\)，那么\((a',b')\)是\((a,b)\)的父節點。要求的就是\((x,y)\)的子樹大小\(\times 2\)。

容易發現\((a,b)\)的直接兒子是\((a+b,a),(2a+b,a),(3a+b,a)\dots\)。如下圖：

考慮改變連邊方式，變成這樣：

再將奇數層（假設\((a,b)\)在第0層）的點編號兩維取反，變成這樣：

問題變成了，一開始有一個二元組\((x,x+y)\)，每次可以進行以下兩種操作之一：

\((a,b)->(a+b,b)\)

\((a,b)->(a,a+b)\)

問能到達多少節點。

用二元組\((p,q)\)表示一個數\(px+q(x+y)\)，那么就是一開始有一個四元組\(((1,0),(0,1))\)，每次可以進行以下兩種操作之一：

\(((a,b),(c,d))->((a+c,b+d),(b,d))\)

\(((a,b),(c,d))->((a,b),(a+c,b+d))\)

可以發現這類似於Stern-Brocot Tree的構造過程，在不考慮\(\le m\)限制的前提下，SBTree上任意兩個相鄰節點代表一個合法的四元組。如下圖：

現在有了\(\le m\)的限制，一對合法四元組是由兩個二元組組成的，考慮在較大的那個二元組處統計它。（較大的二元組指兩維都較大的二元組，這里因為題目性質，兩維都較大的說法是不存在歧義的。）根據SBTree的性質，一對二元組\((i,j)\)出現過當且僅當\(\gcd(i,j)=1\)，且每個二元組會作為較大值貢獻答案\(2\)次（如上圖，在且僅在該二元組第一次出現的那一層，它比相鄰兩個二元組大）。\(\le m\)的限制就是額外要求\(xi+(x+y)j \le m\)。

那么要求的子樹大小就是\(2\times \sum_{i,j}[\gcd(i,j)=1][xi+(x+y)j \le m]\)。隨便莫比烏斯反演以后杜教篩+類歐搞搞就完事了！！1

#include<bits/stdc++.h>
using namespace std;
const int N=2e6+10;
typedef long long ll;

int gi() {
	int x=0,o=1;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
	if(ch=='-') o=-1,ch=getchar();
	while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
	return x*o;
}

int mu[N],smu[N],m;
map<int,int> M;
ll ans;

int getsmu(int n) {
	if(n<N) return smu[n];
	if(M.count(n)) return M[n];
	int ret=1;
	for(int i=2,j;i<=n;i=j+1) {
		j=n/(n/i);ret-=(j-i+1)*getsmu(n/i);
	}
	return M[n]=ret;
}

ll cal(ll n,ll a,ll b,ll c) {
	if(n<0) return 0;
	if(!a) return (b/c)*(n+1);
	if(a>=c||b>=c) return cal(n,a%c,b%c,c)+n*(n+1)/2*(a/c)+(n+1)*(b/c);
	ll m=(a*n+b)/c;
	return n*m-cal(m-1,c,c-b-1,a);
}

ll cal(int x,int y,int n) {
	return cal(n/x-1,x,n%x,y);
}

int main() {
	mu[1]=1;
	for(int i=1;i<N;i++)
		for(int j=i+i;j<N;j+=i) mu[j]-=mu[i];
	for(int i=1;i<N;i++) smu[i]=smu[i-1]+mu[i];
	int T=gi();m=gi();
	while(T--) {
		int x=gi(),y=gi();
		if(x>m||y>m) cout<<"0\n";
		else if(x<=y) cout<<"1\n";
		else {
			y+=x;
			if(y>m) cout<<"2\n";
			else {
				ans=4;
				for(int i=1,j;i<=m;i=j+1) {
					j=m/(m/i);ans+=4ll*(getsmu(j)-getsmu(i-1))*cal(x,y,m/i);
				}
				cout<<ans<<'\n';
			}
		}
	}
	return 0;
}