[BZOJ4899]:记忆的轮廓(概率DP)
题目传送门
题目描述
通往贤者之塔的路上,有许多的危机。
我们可以把这个地形看做是一颗树,根节点编号为1,目标节点编号为n,其中1-n的简单路径上,编号依次递增,
在[1,n]中,一共有n个节点。我们把编号在[1,n]的叫做正确节点,[n+1,m]的叫做错误节点。一个叶子,如果是正
确节点则为正确叶子,否则称为错误叶子。莎缇拉要帮助昴到达贤者之塔,因此现在面临着存档位置设定的问题。
为了让昴成长为英雄,因此一共只有p次存档的机会,其中1和n必须存档。被莎缇拉设置为要存档的节点称为存档
位置。当然不能让昴陷入死循环,所以存档只能在正确节点上进行,而且同一个节点不能存多次档。因为通往贤者
之塔的路上有影响的瘴气,因此莎缇拉假设昴每次位于树上一个节点时,都会等概率选择一个儿子走下去。每当走
到一个错误叶子时,再走一步就会读档。具体的,每次昴到达一个新的存档位置,存档点便会更新为这个位置(假
如现在的存档点是i,现在走到了一个存档位置j>i,那么存档点便会更新为j)。读档的意思就是回到当前存档点
。初始昴位于1,当昴走到正确节点n时,便结束了路程。莎缇拉想知道,最优情况下,昴结束路程的期望步数是多
少?
输入格式
第一行一个正整数T表示数据组数。
接下来每组数据,首先读入三个正整数n,m,p。
接下来m-n行,描述树上所有的非正确边(正确边即连接两个正确节点的边)
用两个正整数j,k表示j与k之间有一条连边,j和k可以均为错误节点,也可以一个为正确节点另一个为错误节点。
输出格式
T行每行一个实数表示每组数据的答案。请保留四位小数。
样例
样例输入:
1
3 7 2
1 4
2 5
3 6
3 7
3 7 2
1 4
2 5
3 6
3 7
样例输出:
9.0000
数据范围与提示
数据保证j是k的父亲。
50≤p≤n≤700,m≤1500,T≤5。
数据保证每个正确节点均有至少2个儿子,至多3个儿子。
题解
大家可能是看了这两篇题解的其中一篇,或者都看了:
个人感觉这两篇题解都是不错的,看一篇就够了(毕竟第一篇是借鉴的第二篇)。
当然我也是看了第二篇题解才写出来的,还将那篇博客推荐给了我左边的大佬,但是他表示对DP式子的意义表示很疑惑,看不懂,想必很多人也有同样的疑惑,所以我就来解释一下,帮助大家的理解。
首先,我想对这道题的题意做一些解释:
1.对于“我们可以把这个地形看做是一颗树,根节点编号为1,目标节点编号为n,其中1-n的简单路径上,编号依次递增,在[1,n]中,一共有n个节点。我们把编号在[1,n]的叫做正确节点,[n+1,m]的叫做错误节点。”这句话的理解,很多同学表示疑惑,正确节点1-n一定是一条链,错误节点一定是这条链上的一些枝杈,为什么呢?
1-n递增,最短路径,有n个节点,如果一棵树的话,就不能同时满足这三条性质了。
2.对于“每当走到一个错误叶子时,再走一步就会读档。”这句话的理解,注意是“叶子”,意思就是说,你做到一个错误节点,会沿着这条枝杈“不撞南墙不回头”,然后“在走一步就会读档”你可以感性的理解为这个错误叶子向上一个读当点连了一条边。
来解释50%的算法:
n=p
注意这时候有一个性质,每一个正确节点都是存档点,那么走到一个错误叶子直接回到理他最近的那个正确节点,即为走出来的那个正确节点。
引用:首先设d[i]表示i的儿子数。设g[i]表示对于一个错误节点i,期望走多少步会读档。那么$ g[i]=1+ \frac{1}{d[i]× \sum g[j] }$其中j是i的儿子。
对于这个式子,你可以先将它化为$ g[i]= \frac{1}{d[i]× \sum g[j]+1 }$,g[j]表示走到这个点你期望走的步数,1表示你走到这个节点后,还需要花1个步数返回上一个存档点。
引用:对于每个正确节点i预处理s[i]表示i的错误儿子的g值和,那么$ s[i]= \sum g[j]$,j是i的错误儿子。
设f[i]表示正确节点i走到n的期望步数,显然f[n]=0,我们倒着递推。
$ f[i]=1+ \frac{1}{d[i]} ×f[i+1]+ \frac{1}{d[i]}× \sum g[j]+f[i] $[j是i的错误儿子]
设f[i]表示正确节点i走到n的期望步数,显然f[n]=0,我们倒着递推。
$ f[i]=1+ \frac{1}{d[i]} ×f[i+1]+ \frac{1}{d[i]}× \sum g[j]+f[i] $[j是i的错误儿子]
对于这个式子,你可以先将这个式子写成$ f[i]= \frac{1}{d[i]} + \frac{1}{d[i]} ×f[i+1]+ \frac{1}{d[i]}× \sum g[j]+1+f[i] $表示从i直接走向i+1所需要的一个步数,$\frac{1}{d[i]×f[i+1]}$表示你有$\frac{1}{d[i]}$的概率会走向i+1这个点,然后付出f[i+1]的步数走到n,$\frac{1}{d[i]× \sum g[j]+f[i]}$表示每一个儿子都有$\frac{1}{d[i]}$的概率会被走到,然后付出g[j]+1的步数(走到j需要一步),然后返回i这个点,再付出f[i]的步数。
引用:移项得f[i]=d[i]+f[i+1]+s[i]
引用:移项得f[i]=d[i]+f[i+1]+s[i]
至于式子的化简:
f[i]=1+1/d[i]*f[i+1]+1/d[i]*∑{g[j]+f[i]}
f[i]=1+1/d[i]*f[i+1]+1/d[i]*s[i]+1/d[i]*∑{f[i]}(s[i]=∑{g[j]})
f[i]=1+1/d[i]*f[i+1]+1/d[i]*s[i]+(d[i]-1)/d[i]*f[i](将f[i]提出,一共有d[i]个儿子,一个是正确儿子)
1/d[i]*f[i]=1+1/d[i]*f[i+1]+1/d[i]*s[i](移项)
f[i]=d[i]+f[i+1]+s[i](同乘d[i])
引用:复杂度线性。
引用:复杂度线性。
然后是70%算法:
引用:我们设dp,f[i,j]表示当前存档点为i,还剩j次存档机会。
首先我们需要预处理一个a[i,j],表示存档点为i,从i开始走到正确节点j的期望步数(中间不能存档)。
显然有边界条件a[i,i]=0。对于i<j,可以列出递推式:
a[i,j]=a[i,j-1]+1+1/d[j-1]*∑{g[k]+a[i,j]}[k是j-1的错误儿子]
首先我们需要预处理一个a[i,j],表示存档点为i,从i开始走到正确节点j的期望步数(中间不能存档)。
显然有边界条件a[i,i]=0。对于i<j,可以列出递推式:
a[i,j]=a[i,j-1]+1+1/d[j-1]*∑{g[k]+a[i,j]}[k是j-1的错误儿子]
对于这个式子,我还是先将它写成a[i,j]=a[i,j-1]+1/d[j-1]+1/d[j-1]*∑{g[k]+1+a[i,j]}a[i,j-1]表示它需要从j-1这个点走过来,所以肯定要加上它的步数;1/d[j-1]表示从j-1走到j有1/d[j-1]的概率直接走过来;1/d[j-1]表示从j-1走过来的时候会有概率走向错误节点;g[k]+1表示从错误节点k走到i要走的步数(走到这个节点需要1步);a[i,j]表示它会回到存档点i然后再把这部分走一遍;因为j-1有(d[j-1]-1)/d[j-1]个错误节点,然后我们把∑里的1提出来,再加上外面的1/d[j-1],就得了1。
引用:移项得a[i,j]=a[i,j-1]*d[j-1]+d[j-1]+s[j-1]
引用:移项得a[i,j]=a[i,j-1]*d[j-1]+d[j-1]+s[j-1]
至于式子的化简:
a[i,j]=a[i,j-1]+1+1/d[j-1]*∑{g[k]+a[i,j]}
a[i,j]=a[i,j-1]+1+1/d[j-1]*s[j-1]+1/d[i]*∑{a[i,j]}
a[i,j]=a[i,j-1]+1+1/d[j-1]*s[j-1]+(d[i]-1)/d[i]*a[i,j]
1/d[j-1]*a[i,j]=a[i,j-1]+1+1/d[j-1]*s[j-1]
a[i,j]=a[i,j-1]*d[j-1]+d[j-1]+s[j-1]
(对于式子化简的理解可以参考上面的50%算法,在此就不做过多的赘述)
引用:可以用n^2的时间预处理a,然后做dp就很好转移了。
枚举下一次的存档点k,那么f[i,j]可以由f[k,j-1]+a[i,k]转移而来。
引用:可以用n^2的时间预处理a,然后做dp就很好转移了。
枚举下一次的存档点k,那么f[i,j]可以由f[k,j-1]+a[i,k]转移而来。
相当与枚举存档点,然后转移。
引用:复杂度O(n2p)
引用:复杂度O(n2p)
最后是100%算法:
在此,我想将一下对于70%算法优化的算法。
引用:不过,有没有发现,70%和100%好像都没考虑一个问题。观察a数组,可以看到它是恐怖的增长的,我们最终答案会不会爆炸?
我们来估计答案的上界。考虑一种可行方案,每n/p个正确节点就设立一次存档位置,那么答案最大是多少呢?考虑最坏情况,观察a的转移,应该每变 换一次存档点,大约需要3(n/p)*s[i]+3(n/p-1)*s[i+1]+3(n/p-2)*s[i+2]+……
我们来估计答案的上界。考虑一种可行方案,每n/p个正确节点就设立一次存档位置,那么答案最大是多少呢?考虑最坏情况,观察a的转移,应该每变 换一次存档点,大约需要3(n/p)*s[i]+3(n/p-1)*s[i+1]+3(n/p-2)*s[i+2]+……
式子中3可以理解为最多有3个儿子(数据保证)
引用:因为最多m个节点,s的上限是1500(实际上也远远达不到),把所有s都视为这个上限,提取公因数,计算一下那个等比数列求和,由于p是有下界的, 因此n/p有上界14,发现最后也就是个12位数的样子,那么我们估计出答案最大也不会超过这个,可以放心做了。而至于a会爆炸的问题,double是可以 存很多位的,而且太大的a肯定不可能被用上。
那么其实,针对答案不会特别大,a的增长又很恐怖,我们还可以思考对70%的算法优化。那就是设定一个常数step,每次转移最多从距当前step步远的 位置转移过来。step取40多基本不会有问题了,因为a的下界已经是2^40了,而答案的上界远远没有达到,经过精确计算还可以再把step调小一点。
引用:因为最多m个节点,s的上限是1500(实际上也远远达不到),把所有s都视为这个上限,提取公因数,计算一下那个等比数列求和,由于p是有下界的, 因此n/p有上界14,发现最后也就是个12位数的样子,那么我们估计出答案最大也不会超过这个,可以放心做了。而至于a会爆炸的问题,double是可以 存很多位的,而且太大的a肯定不可能被用上。
那么其实,针对答案不会特别大,a的增长又很恐怖,我们还可以思考对70%的算法优化。那就是设定一个常数step,每次转移最多从距当前step步远的 位置转移过来。step取40多基本不会有问题了,因为a的下界已经是2^40了,而答案的上界远远没有达到,经过精确计算还可以再把step调小一点。
这部分我认为他说的还是挺清楚的,就不做过多的赘述。
引用:复杂度O(np log ans)
引用:复杂度O(np log ans)
注:以上引用全部来自第二篇博客,画横线部分为引用部分。
代码时刻
50%代码:
#include<bits/stdc++.h>
using namespace std;
struct rec
{
int nxt;
int to;
}e[5000];
int head[5000],cnt;
int n,m,p;
double g[5000],s[5000],dp[5000];
bool vis[5000];
int du[5000];
void pre_work()//多测不清空,爆零两行泪……
{
memset(head,0,sizeof(head));
memset(dp,0,sizeof(dp));
memset(g,0,sizeof(g));
memset(s,0,sizeof(s));
memset(du,0,sizeof(du));
cnt=0;
}
void add(int x,int y)
{
e[++cnt].nxt=head[x];
e[cnt].to=y;
head[x]=cnt;
}
void dfsgetG(int x)//计算g数组
{
vis[x]=1;
g[x]=1;
for(int i=head[x];i;i=e[i].nxt)
{
dfsgetG(e[i].to);
g[x]+=g[e[i].to]*1.0/(double)du[x];
}
}
int main()
{
int T;
scanf("%d",&T);
while(T--)
{
pre_work();
scanf("%d%d%d",&n,&m,&p);
for(int i=1;i<=n-1;i++)du[i]=1;//1-(n-1)中,每一个节点都要加上它到下一个正确节点的边
for(int i=n+1;i<=m;i++)
{
int a,b;
scanf("%d%d",&a,&b);
du[a]++;
add(a,b);
}
for(int i=n+1;i<=m;i++)
{
if(vis[i])continue;
dfsgetG(i);
}
for(int i=1;i<=n;i++)
for(int j=head[i];j;j=e[j].nxt)
s[i]+=g[e[j].to];//计算s数组
for(int i=n-1;i;i--)//倒推计算答案
dp[i]=dp[i+1]+du[i]+s[i];
cout<<fixed<<setprecision(4)<<dp[1]<<endl;//保留小数输出
}
return 0;
}
70%代码:
#include<bits/stdc++.h>
using namespace std;
struct rec
{
int nxt;
int to;
}e[2000];
int head[2000],cnt;
int n,m,p;
double g[2000],s[2000],dp[2000][2000],Map[2000][2000];
bool vis[2000];
int du[2000];
void pre_work()
{
for(int i=0;i<=1;i++)e[i].nxt=e[i].to=0;
memset(head,0,sizeof(head));
memset(dp,127,sizeof(dp));
memset(g,0,sizeof(g));
memset(s,0,sizeof(s));
memset(du,0,sizeof(du));
memset(Map,0,sizeof(Map));
memset(vis,0,sizeof(vis));
cnt=0;
}
void add(int x,int y)
{
e[++cnt].nxt=head[x];
e[cnt].to=y;
head[x]=cnt;
}
void dfsgetG(int x)
{
vis[x]=1;
g[x]=1;
for(int i=head[x];i;i=e[i].nxt)
{
dfsgetG(e[i].to);
g[x]+=g[e[i].to]*1.0/(double)du[x];
}
}
int main()
{
int T;
scanf("%d",&T);
while(T--)
{
pre_work();
scanf("%d%d%d",&n,&m,&p);
for(int i=1;i<=n-1;i++)du[i]=1;
for(int i=n+1;i<=m;i++)
{
int a,b;
scanf("%d%d",&a,&b);
du[a]++;
add(a,b);
}
for(int i=n+1;i<=m;i++)
{
if(vis[i])continue;
dfsgetG(i);
}
for(int i=1;i<=n;i++)
for(int j=head[i];j;j=e[j].nxt)
s[i]+=g[e[j].to];
for(int i=1;i<=n;i++)
for(int j=i+1;j<=n;j++)
Map[i][j]=(Map[i][j-1]+1)*(double)du[j-1]+s[j-1];//计算a数组
dp[n][1]=0;//dp初始
for(int j=2;j<=p;j++)
for(int i=1;i<=n;i++)
for(int k=i+1;k<=n;k++)
dp[i][j]=min(dp[i][j],dp[k][j-1]+Map[i][k]);
cout<<fixed<<setprecision(4)<<dp[1][p]<<endl;
}
return 0;
}
100%算法:
#include<bits/stdc++.h>
using namespace std;
struct rec
{
int nxt;
int to;
}e[2000];
int head[2000],cnt;
int n,m,p;
double g[2000],s[2000],dp[2000][2000],Map[2000][2000];
bool vis[2000];
int du[2000];
void pre_work()
{
for(int i=0;i<=1;i++)e[i].nxt=e[i].to=0;
memset(head,0,sizeof(head));
memset(dp,127,sizeof(dp));
memset(g,0,sizeof(g));
memset(s,0,sizeof(s));
memset(du,0,sizeof(du));
memset(Map,0,sizeof(Map));
memset(vis,0,sizeof(vis));
cnt=0;
}
void add(int x,int y)
{
e[++cnt].nxt=head[x];
e[cnt].to=y;
head[x]=cnt;
}
void dfsgetG(int x)
{
vis[x]=1;
g[x]=1;
for(int i=head[x];i;i=e[i].nxt)
{
dfsgetG(e[i].to);
g[x]+=g[e[i].to]*1.0/(double)du[x];
}
}
int main()
{
int T;
scanf("%d",&T);
while(T--)
{
pre_work();
scanf("%d%d%d",&n,&m,&p);
for(int i=1;i<=n-1;i++)du[i]=1;
for(int i=n+1;i<=m;i++)
{
int a,b;
scanf("%d%d",&a,&b);
du[a]++;
add(a,b);
}
for(int i=n+1;i<=m;i++)
{
if(vis[i])continue;
dfsgetG(i);
}
for(int i=1;i<=n;i++)
for(int j=head[i];j;j=e[j].nxt)
s[i]+=g[e[j].to];
for(int i=1;i<=n;i++)
for(int j=i+1;j<=n;j++)
Map[i][j]=(Map[i][j-1]+1)*(double)du[j-1]+s[j-1];
dp[n][1]=0;
for(int j=2;j<=p;j++)
for(int i=1;i<=n;i++)
for(int k=i+1;k<=min(n,i+12);k++)//优化上界
dp[i][j]=min(dp[i][j],dp[k][j-1]+Map[i][k]);
cout<<fixed<<setprecision(4)<<dp[1][p]<<endl;
}
return 0;
}
rp++