我们发现，派机器人一定是尽量让每个洞穴被派的次数相同最优，那么，我们不妨以这样的策略去派机器人，设洞穴编号从 $0\sim n-1$，我们第 $i$ 天派机器人的洞穴为 $ik\bmod n\sim (i+1)k-1\bmod n$，这样就是最优策略，容易发现，每相邻 $\gcd(n,k)$ 个洞穴是完全等价的，所以我们可以令 $n,k$ 都除以 $\gcd(n,k)$。

设 $f_i$ 为宝藏在洞穴 $i$ 内，找到宝藏的期望天数，答案即为 $\frac{1}n\sum\limits^{n-1}_{i=0} f_i$。

对于 $k\le i< n$，我们有 $f_i=f_{i-k}+1$，因为在访问这些洞穴之前，我们一定会在前一天访问洞穴 $i-k$；对于 $0\le i<k$ ，我们有 $f_i=\frac 12+\frac{1}{2}(f_{i-k+n}+1)=\frac 12f_{i-k+n}+1$，即有 $\frac 12$ 的概率一次找到，有 $\frac 12$ 的概率第一次没找到，那么在再次访问洞穴 $i$ 的前一天，我们一定会访问 $i-k+n$。

由于 $\gcd(n,k)=1$，$ik$ 这个式子 $i$ 从 $0$ 到 $n-1$ 可以不重不漏的遍历 $\bmod n$ 的所有数，我们可以将 $f_{1\sim n-1}$ 表示为 $af_0+b$ 的形式，再根据 $f_{0}=\frac 12f_{n-k}+1$，就可以得到一个关于 $f_0$ 的一元一次方程，就可以解出 $f_0$，进而得到所有 $f_i$ 的值，时间复杂度 $O(n)$。

介绍一个东西：万能欧几里得算法，在平面直角坐标系内有一条不包含左端点的线段（$y=\frac{ax+b}{c},x\in(0,n],a,b\in\mathbb N,c,x\in\mathbb Z^+$），在其定义域内，从左到右，维护一个字符串，直线每碰到一条水平的整线就写下一个 U，每碰到一条竖直的整线就写下一个 R，碰到整点就先写 U 再写 R，给定一个字符串函数 $F$，求 $F(S)$。

要求 $F$ 满足，我们可以实现运算 $\cdot$：$F(S_1)\cdot F(S_2)=F(S_1+S_2)$，且 $\cdot$ 有结合律。假设 $\cdot$ 的时间复杂度为 $O(t)$，那么总时间复杂度为 $O(t\log \max(a,c))$，

万能欧几里得算法学习笔记

对于这题，考虑我们要维护的东西，一个是关于 $f_0$ 的一次函数，表示 $f_i$；一个是关于 $f_0$ 的一次函数，表示 $f_i$ 的前缀和。

那么我们有 $a=k,b=0,c=n$，我们令 $F$ 的取值为两个一次函数，表示为 $(A(x),B(x))$，那么有 $F(\text U)=(\frac 12x,0),F(\text R)=(x+1,x+1)$，$\cdot$ 定义如下 $(A(x),B(x))\cdot (C(x),D(x))=(C(A(x)),B(x)+D(A(x)))$。

这样我们的结果就是 $f_0$ 和 $\sum f_i$ 了。

我们得到了一个关于 $f_0$ 的一元一次方程，就可以解出 $f_0$，代入 $\sum f_i$ 中即可得到 $\sum f_i$ 的真实值，这题就做完了，时间复杂度 $O(T\log \max(n,k))$。

#include<iostream>
#include<numeric>
using std::cin;
using std::cout;
int const mod=1e9+7,inv2=(mod+1)/2;
int pow(int x,int y){
	int res=1;
	while(y){
		if(y&1) res=1ll*res*x%mod;
		x=1ll*x*x%mod,y>>=1;
	}
	return res;
}
struct node{
	int k,b,sumk,sumb;
	node():k(1),b(0),sumk(0),sumb(0){}
	node(int x,int y,int z,int w):k(x),b(y),sumk(z),sumb(w){}
	friend node operator *(node const &a,node const &b){
		return node(1ll*a.k*b.k%mod,(1ll*b.k*a.b+b.b)%mod,(a.sumk+1ll*a.k*b.sumk)%mod,(a.sumb+1ll*b.sumk*a.b+b.sumb)%mod);
	}
};
node pow(node x,int y){
	node res;
	while(y){
		if(y&1) res=res*x;
		x=x*x,y>>=1;
	}
	return res;
}
node solve(int a,int b,int c,int n,node u,node r){
	if(n==0) return node();
	if(b>=c) b%=c;
	if(a>=c) return solve(a%c,b,c,n,u,pow(u,a/c)*r);
	int m=(1ll*a*n+b)/c;
	if(m==0) return pow(r,n);
	return pow(r,(c-b-1)/a)*u*solve(c,c-b-1,a,m-1,r,u)*pow(r,n-(1ll*c*m-b-1)/a);
}
int main(){
	std::ios::sync_with_stdio(false),cin.tie(nullptr);
	int T;
	cin>>T;
	node u(inv2,0,0,0),r=node(1,1,1,1);
	while(T--){
		int n,k;
		cin>>n>>k;
		int g=std::gcd(n,k);
		n/=g,k/=g;
		node p=solve(k,0,n,n,u,r);
		int x=1ll*p.b*pow(1-p.k+mod,mod-2)%mod;
		cout<<(1ll*p.sumk*x+p.sumb)%mod*pow(n,mod-2)%mod<<'\n';
	}
	return 0;
}