jzoj 5917.【NOIP2018模拟10.20】moon(矩阵乘法+概率dp)

Description

作为申国的学者，你需要严格遵守三大基本原则：
战争即和平
自由即奴役
无知即力量
你正在对一本书进行审核，其中片段写道：
“少焉，月出于东山之上，徘徊于斗牛之间。白露横江，水光接天。纵一苇之所如，凌万顷之茫然。浩浩乎如冯虚御风，而不知其所止；飘飘乎如遗世独立，羽化而登仙。”
这种行为明显不符合三大原则，比如“纵一苇之所如”中自由的意思已经在新话中杯删除了。
但是你在修改的同时，发现书中夹着一道问题：
酥室等人现在的位置是(x,y)，同时还有n个景点，坐标分别为(xi,yi)。
每次移动按照下面的顺序操作：
1、 选择一条直线，要求直线经过现在的位置和至少两个景点（如果现在在某个景点那里，也算一个）如果有多条直线满足要求，等概率选择一条。
2、 在选择的这条直线中，等概率选择一个直线覆盖了的景点移动过去，如果目前在景点上，也有可能停住不动。
酥室会进行若干次询问，第i次询问从一个你选的任意点出发（可以不是景点），然后连续移动mi步，最后到达ti的最大概率是多少。

---

Input

第一行一个整数n。
接下来n行每行两个整数，表示xi,yi。
接下来一行一个整数m，表示询问数。
接下来m行每行两个整数ti,mi，表示询问。

---

Output

对于每个询问输出一行，表示最大的概率，与标准答案误差在10^-6以内就算正确。

---

Sample Input5

0 0
1 3
2 2
3 1
4 4
10
1 1
2 1
3 1
4 1
5 1
3 2
3 3
3 4
3 5
3 6

---

Sample Output

0.50000000000000000000
0.50000000000000000000
0.33333333333333331483
0.50000000000000000000
0.50000000000000000000
0.18518518518518517491
0.15226337448559670862
0.14494741655235482414
0.14332164812274550414
0.14296036624949901017

---

Data Constraint

---

分析

这题面也没啥毛病啊
我们可以从终点开始走，这样就不需要去考虑起点。
设$f_{i,j}$表示第i步走到j的概率，转移为:
$$f_{i,k}+=\frac{f_{i-1,j}}{cnt\times num}$$
j,k表示从j走到k,cnt表示有cnt条直线经过k，num表示j,k所在这条直线的上的点数。
但这样时间是无法接受的，只能拿到七十分。
我们观察转移的式子，比较优美，可以用矩乘来优化，这样就可以拿到满分了。

---

code

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define db double
using namespace std;
const int N=210,M=1e4+10;
const double w=1e-8;
int n,m,bz[N][N],L[N*N][N],b[N],top;
struct matrix {
	double a[N][N];
	int n,m;
	matrix(int x=0,int y=0) {n=x,m=y,memset(a,0,sizeof(a));}
	matrix friend operator*(matrix x,matrix y) {
		matrix z(min(x.n,y.n),min(x.m,y.m));
		for (int i=1;i<=min(x.n,y.n);i++)
			for (int j=1;j<=min(x.m,y.m);j++)
				for (int k=1;k<=min(x.m,y.n);k++)
					z.a[i][j]+=x.a[i][k]*y.a[k][j];
		return z;
	}
}Ma[20];
struct P{
	int x,y;
}a[N];
bool pd(int i,int j,int k) {
	if (a[i].x==a[j].x) {
		if (a[i].x==a[k].x) return true;
	}
	else {
		if (a[i].x!=a[k].x && a[j].x!=a[k].x) {
			db k1=(db)(a[i].y-a[k].y)/(db)(a[i].x-a[k].x),k2=(db)(a[j].y-a[k].y)/(db)(a[j].x-a[k].x);
			if (fabs(k1-k2)<=w) return true;
		}
	}
	return false;
}
void work(int x) {
	L[top][++L[top][0]]=x;
	b[x]++;
}
int main() {
	freopen("moon.in","r",stdin);
	freopen("moon.out","w",stdout);
	int i,j,k;
	scanf("%d",&n);
	for (i=1;i<=n;i++) scanf("%d%d",&a[i].x,&a[i].y);
	for (i=1;i<=n;i++) {
		for (j=i+1;j<=n;j++) {
			if (bz[i][j]==0) {
				bz[j][i]=bz[i][j]=++top;	
				work(i);work(j);
				for (k=j+1;k<=n;k++)
					if (pd(i,j,k)) {
						for (int k1=1;k1<=L[top][0];k1++)
							bz[L[top][k1]][k]=bz[k][L[top][k1]]=top;
						work(k);
					}
			}
		}
	}
	for (i=1;i<=top;i++) {
		for (j=1;j<=L[i][0];j++) 
			for (k=1;k<=L[i][0];k++) {
				int x=L[i][j],y=L[i][k];
				Ma[0].a[x][y]+=1.0/(db)L[i][0]/(db)b[x];
			}
	}
	Ma[0].n=Ma[0].m=n;
	for (i=1;i<=log2(M);i++) Ma[i]=Ma[i-1]*Ma[i-1]; 
	scanf("%d",&m);
	for (int k1=1;k1<=m;k1++) {
		int z1,z2,x;scanf("%d%d",&z1,&z2);
		matrix t(n,1);
		t.a[z1][1]=1;
		for (z2--,x=floor(log2(z2));z2;z2-=(1<<x),x=floor(log2(z2))) t=Ma[x]*t;
		db ans=0.0,sum;
		for (i=1;i<=top;i++) {
			sum=0.0;
			for (j=1;j<=L[i][0];j++) sum+=t.a[L[i][j]][1];
			ans=max(ans,sum/(db)L[i][0]);
		}
		printf("%.10lf\n",ans);
	}
}