【描述】
dolls终于找到了正确的路,来到了n个洞穴前,他得知其中一个有宝藏,dolls每天可以派k个机器人去一些洞穴找宝藏,如果机器人所在洞穴里有宝藏就有p的概率找到宝藏,求最优策略下找到宝藏的期望天数。
n
,
k
<
=
5
e
8
n,k<=5e8
n,k<=5e8
【思路】
期望神题。首先最优策略显然是优先探索探索次数少的洞穴。如果gcd(n,k)!=1,那么有gcd(n,k)块将永远被一起探索,所以我们可以先把n,k均除以gcd(n,k)并不影响答案。根据我们的探索顺序,我们把洞穴进行编号1-n。设
E
i
E_{i}
Ei表示宝藏在i号洞穴时找到的期望天数。显然:
a
n
s
=
∑
i
=
1
n
E
i
/
n
ans=\sum_{i=1}^{n} E_{i}/n
ans=∑i=1nEi/n
而显然E之间满足关系:
E
i
=
E
i
−
k
+
1
(
i
>
k
)
,
记
为
E
i
=
A
(
E
i
−
k
)
E_{i}=E_{i-k}+1(i>k),记为E_{i}=A(E_{i-k})
Ei=Ei−k+1(i>k),记为Ei=A(Ei−k)
E
i
=
(
1
−
p
)
(
E
i
−
k
+
n
+
1
)
+
p
∗
1
,
记
为
E
i
=
B
(
E
i
−
k
+
n
)
E_{i}=(1-p)(E_{i-k+n}+1)+p*1,记为E_{i}=B(E_{i-k+n})
Ei=(1−p)(Ei−k+n+1)+p∗1,记为Ei=B(Ei−k+n)
此时我们只需要把所有E用其中一个E表示出来即可实现O(n)求解。
考虑优化:
E的第一个公式很有趣,这意味着如果我们按k分为几块:
我们会把n分为这样的几块,若干部分大小为k,一部分大小为
n
m
o
d
k
n\mod k
nmodk。而我们只考虑最下面那一块,上面位置和它对应的E均可用最下面那一块的E表示出来,这样我们把对应位置求和可以把要求的东西转换为只与最下面那一块的E大小有关的一次函数S1(x),S2(x)。这样我们就只需要求出最下面的k个E了。这很有趣。但是这不够,我们可不可以考虑递归呢?注意到我们求S1,S2需要E的第一个公式,所以我们考虑维护新的子问题的E的第一个公式,记为
A
′
(
x
)
A'(x)
A′(x)。在新的子问题里面,我们定义A’表示
E
i
=
A
′
(
E
i
−
(
n
m
o
d
k
)
)
=
A
′
(
E
i
−
(
n
−
n
/
k
∗
k
)
)
=
A
−
n
/
k
+
1
(
B
−
1
(
E
i
−
(
n
m
o
d
k
)
)
)
E_{i}=A'(E_{i-(n\mod k)})=A'(E_{i-(n-n/k*k)})=A^{-n/k+1}(B^{-1}(E_{i-(n\mod k)}))
Ei=A′(Ei−(nmodk))=A′(Ei−(n−n/k∗k))=A−n/k+1(B−1(Ei−(nmodk)))。这样我们就找到了维护子问题公式A的办法。为了递归维护A’,我们需要维护B’,同理可得
E
i
=
B
′
(
E
i
−
(
n
m
o
d
k
)
+
n
)
=
A
′
(
E
i
−
(
n
−
n
/
k
∗
k
)
)
=
A
−
n
/
k
(
B
−
1
(
E
i
−
(
n
m
o
d
k
)
)
)
(
定
义
f
n
(
x
)
=
f
(
f
n
−
1
(
x
)
)
)
E_{i}=B'(E_{i-(n\mod k)}+n)=A'(E_{i-(n-n/k*k)})=A^{-n/k}(B^{-1}(E_{i-(n\mod k)}))(定义f^{n}(x)=f(f^{n-1}(x)))
Ei=B′(Ei−(nmodk)+n)=A′(Ei−(n−n/k∗k))=A−n/k(B−1(Ei−(nmodk)))(定义fn(x)=f(fn−1(x)))。这样我们就成功构造新的子问题了,可以递归求解了。递归边界时
E
=
A
(
E
)
E={A(E)}
E=A(E)即可解出E的答案。最初时为了和子问题统一形式,我们也定义S1(x)=x,S2(x)=x,A(x)=x+1,B(x)=(1-p)*x+1。
现在我们考虑维护这些一次函数:
1.求逆
在维护新的A’,B’时,我们需要对A,B求逆。这也很简单,即构造函数
A
−
1
A^{-1}
A−1满足
A
−
1
(
A
(
x
)
)
=
x
A^{-1}(A(x))=x
A−1(A(x))=x,代入A(x)即可求解
A
−
1
(
x
)
A^{-1}(x)
A−1(x)。
2.快速幂
这里我们需要快速对函数进行复合运算,函数的自我复合满足结合律,可以使用快速幂快速复合。
3.求和
在维护S1,S2时,我们需要求出
S
1
′
=
S
1
+
∑
i
=
1
n
/
k
−
1
S
2
(
A
i
(
x
)
)
,
S
2
′
=
S
1
+
∑
i
=
1
n
/
k
S
2
(
A
i
(
x
)
)
S1'=S1+\sum_{i=1}^{n/k-1} S2(A^{i}(x)),S2'=S1+\sum_{i=1}^{n/k} S2(A^{i}(x))
S1′=S1+∑i=1n/k−1S2(Ai(x)),S2′=S1+∑i=1n/kS2(Ai(x))。我们可以先求出
S
2
(
∑
i
=
1
m
A
i
(
x
)
)
S2(\sum_{i=1}^{m} A^{i}(x))
S2(∑i=1mAi(x)),再处理常数项的贡献。我们把里面相加,发现一次项系数是一个等比数列求和,常数项是一个差比数列求和,所以可以快速维护。
代码(有注释):
#include<bits/stdc++.h>
#define re register
using namespace std;
const int N=1e5+5;
int n,m,p;
inline int red()
{
int data=0;int w=1; char ch=0;
ch=getchar();
while(ch!='-' && (ch<'0' || ch>'9')) ch=getchar();
if(ch=='-') w=-1,ch=getchar();
while(ch>='0' && ch<='9') data=(data<<3)+(data<<1)+ch-'0',ch=getchar();
return data*w;
}const int mod=1e9+7;
inline int mul(const int&a,const int&b){return 1ll*a*b%mod;}
inline int add(const int&a,const int&b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(const int&a,const int&b){return a-b>=0?a-b:a-b+mod;}
struct node{//定义一次函数形式的变换的结构体
int a,b;
friend inline node operator*(const node&A,const node&B){//复合运算A(B(x))
return (node){mul(A.a,B.a),add(mul(A.a,B.b),A.b)};
}friend inline node operator+(const node&A,const node&B){return (node){add(A.a,B.a),add(A.b,B.b)};}//A(x)+B(x)
};
inline int ksm(int a,int b){
int ret=1;
while(b){
if(b&1)ret=mul(ret,a);
a=mul(a,a);b>>=1;
}return ret;
}
inline int inv(const int&a){return ksm(a,mod-2);}
inline node inv(const node&A){int c=inv(A.a);return (node){c,dec(0,mul(A.b,c))};}//一次函数变换求逆
inline node ksm(node a,int b){//一次函数变换的复合运算快速幂
node ret=(node){1,0};//一次函数变换的单位元
while(b){
if(b&1)ret=ret*a;
a=a*a;b>>=1;
}return ret;
}
inline node sum(const node&A,const int&n){//求和A^1+A^2+A^3......
if(!n)return (node){0,0};
if(A.a==1)return (node){n,mul(mul(mul(n+1,n),inv(2)),A.b)};//等比数列求和的特殊情况
return (node){mul(dec(ksm(A.a,n+1),A.a),inv(dec(A.a,1))),mul(A.b,mul(dec(add(n,ksm(A.a,n+1)),mul(n+1,A.a)),inv(ksm(dec(1,A.a),2))))};
}
inline int solve(int n,int k,node A,node B,node S1,node S2){
if(k==0)return mul(S2.a,mul(dec(0,A.b),inv(A.a-1)));//求解A(x)=x形式的方程
node a,b,s1,s2;int floor=n/k;//计算子问题的参数
a=ksm(inv(A),floor-1)*inv(B);
b=ksm(inv(A),floor)*inv(B);
s1=S1+(node){S2.a,0}*sum(A,floor)+(node){0,S2.b};
s2=S1+(node){S2.a,0}*sum(A,floor-1)+(node){0,S2.b};
int sum=add(mul(s1.b,n%k),mul(s2.b,k-n%k));//必须先处理常数项贡献,不能留到下一层,否则求和会出问题
return s1.b=0,s2.b=0,add(sum,solve(k,n%k,a,b,s1,s2));
}
inline int gcd(int a,int b){return b?gcd(b,a%b):a;}
int main(){
int T_T=red();
while(T_T--){
n=red();m=red();p=red();
int gd=gcd(n,m);n/=gd;m/=gd;
printf("%d\n",mul(solve(n,m,(node){1,1},(node){mod+1-p,1},(node){1,0},(node){1,0}),inv(n)));
}
}