[CF868G]El Toll Caves

Description

有n个洞穴，其中一个有宝藏。
你每天有k次机会去洞穴中找宝藏，如果你去到的洞穴中有宝藏则有1/2的概率找到。每次寻找的概率是独立计算的。
问找到宝藏的期望天数。
k<=n<=5e8

Solution

可以感受出来一个结论是你需要让每个点的访问次数尽量平均
考虑如下策略，将洞穴编号为0~n-1，第i天访问ik%n,ik+1%n,ik+2%n…ik+k-1%n这k个洞穴
设E[i]表示宝藏在i的期望天数
可以发现，对于i∈[0,n-k)，E[i+k]=E[i]+1，因为无论如何你访问i+k之前会访问i
对于i∈[n-k,n)，E[i+k-n]=E[i]/2+1，除去第一次访问i+k-n就找到了
这样我们就得到了一种O(n)做法，但不够，我们需要考虑优化
举个栗子：n=8，k=3
那么E[3]=E[0]+1,E[6]=E[3]+1=E[0]+2,E[4]=E[1]+1,E[7]=E[4]+1=E[1]+2,E[5]=E[2]+1
于是ans=∑E[0]=3E[0]+3E[1]+2E[2]+7
推广一下我们需要把(n,k)的问题变成(k,n%k)的问题
设$$ans=\sum _{i=0}^{k-1}S1(E[i])+\sum _{i=k}^{n-1}S2(E[i])$$
和转移，$E[i+k]=A(E[i]),E[i+k-n]=B(E[i])$
其中S1，S2，A，B是形如kx+b的一次函数
再约定$A^i=A(A^{i-1})$即为A的i次复合
那么我们可以得到
$$S1=S1+\sum _{i=0}^{\lfloor \frac{n}{k}\rfloor }S2(A^i)$$
$$S2=S1+\sum _{i=0}^{\lfloor \frac{n}{k}\rfloor -1}S2(A^i)$$
考虑如何维护A和B，注意到$$E[i+\lfloor \frac{n}{k}\rfloor \ast k-n]=B(A^{\lfloor \frac{n}{k}\rfloor -1}(E[i]))$$
即
$$E[i-n\%k]=B(A^{\lfloor \frac{n}{k}\rfloor -1}(E[i]))$$
两边求逆就是$$E[i]=A^{-(\lfloor \frac{n}{k}\rfloor -1)}(B^{-1}(E[i-n\%k]))$$
这个就是A’了
同理B’可以推出来是$$A^{-\lfloor \frac{n}{k}\rfloor }(B^{-1})$$
然后就做完了

Code

#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef long long ll;

const int Mo=1e9+7,inv2=5e8+4;

int ty,n,k,p;

int gcd(int x,int y) {return y?gcd(y,x%y):x;}

struct P{
	int k,b;
	P(int _k=0,int _b=0) {k=_k;b=_b;}
};

P operator * (P a,P b) {return P((ll)a.k*b.k%Mo,((ll)b.k*a.b%Mo+b.b)%Mo);}
P operator - (P a,P b) {return P((a.k-b.k+Mo)%Mo,(a.b-b.b+Mo)%Mo);}
P operator + (P a,P b) {return P((a.k+b.k)%Mo,(a.b+b.b)%Mo);}

int calc(P a,int x) {return ((ll)a.k*x%Mo+a.b)%Mo;}

int pwr(int x,int y) {
	int z=1;
	for(;y;y>>=1,x=(ll)x*x%Mo)
		if (y&1) z=(ll)z*x%Mo;
	return z;
}

P inv(P a) {
	int ia=pwr(a.k,Mo-2);
	return P(ia,Mo-(ll)a.b*ia%Mo);
}

P pw(P x,int y) {
	P z=P(1,0);
	for(;y;y>>=1,x=x*x)
		if (y&1) z=z*x;
	return z;
}

P ps(P x,int y) {
	if (!y) return P(0,0);
	if (x.k==1) return P(y,(ll)y*(y+1)%Mo*inv2%Mo*x.b%Mo);
	P a=pw(x,y+1)-x,z;
	z.k=(ll)a.k*pwr(x.k-1,Mo-2)%Mo;
	z.b=(ll)(z.k-y+Mo)%Mo*pwr(x.k-1,Mo-2)%Mo*x.b%Mo;
	return z;
}

int solve(int n,int k,P A,P B,P S1,P S2) {
	if (!k) return ((ll)S2.k*A.b%Mo*pwr(Mo+1-A.k,Mo-2)%Mo+S2.b)%Mo;
	P a,b,s1,s2;
	int k1=n%k,t=n/k-1;
	a=inv(B)*pw(inv(A),t);
	b=inv(B)*pw(inv(A),t+1);
	s2=S1+ps(A,t)*S2;
	s1=S1+ps(A,t+1)*S2;
	int sum=((ll)(n%k)*s1.b%Mo+(ll)(k-n%k)*s2.b%Mo)%Mo;
	s1.b=s2.b=0;
	return (sum+solve(k,n%k,a,b,s1,s2))%Mo;
}

int main() {
	for(scanf("%d",&ty);ty;ty--) {
		scanf("%d%d",&n,&k);//p=(Mo+1-p)%Mo;
		p=inv2;
		int g=gcd(n,k);n/=g;k/=g;
		printf("%d\n",(ll)solve(n,k,P(1,1),P(p,1),P(1,0),P(1,0))*pwr(n,Mo-2)%Mo);
	}
	return 0;
}