CSU - 1917 There is no SSR （概率dp+矩阵优化）

In the ACMers of CSU, Zihao LI is a fan of yinyangshi. However, he has African blood in his veins. One day, he collected many charms to conjure Shikigami, which Chinese called shishen. Shishen has three levels:SSR, SR, R, and SSR is a rarity, R is ordinary, SR is between them. For most people, they can get one SSR, nineteen SR, and eighty R when he uses 100 charms. But Mr.Li is different. He can’t get SSR for his African blood...So, the score of luck for poor Mr.Li depend to the number of SR he gets in succession. Now, we know the probability of Mr.Li to conjure R, SR and the sum is equal to one. We want to know the probability of Mr.Li getting n SR in succession, which means he is lucky.

Input

There are many tests. Each case has two floats represent the probability of R, SR(0 < p1, p2 < 1)and two integer n (0 < n <= 100)as described above and m (0 < m <=1e7), the numbers of charms.

Output

For each test case, print a single float on its own line denoting the probability of Mr.Li getting n SR in succession. The answer should be around to the 1e-6.

Sample Input

0.0 1.0 1 2
1.0 0.0 1 2
0.5 0.5 1 2

Sample Output

1.000000
0.000000
0.750000

题意：阴阳师抽奖，得到R的概率为p，得到SR的概率为1-p，求抽m次奖连续获得n个SR的概率。

思路：我们可以反过来考虑，求出不会连续得到n个SR的概率，这个问题我们可以通过简单dp递推得到【dp[i]为抽i次奖不会连续得到n个SR的概率】：

但是由于m太大，因此我们没法线性递归，而由该递推式的特点，我们发现可以通过矩阵优化，并通过矩阵快速幂的方式在n^3logm的复杂度内解决。

#include<set>
#include<map>
#include<stack>
#include<queue>
#include<vector>
#include<string>
#include<math.h>
#include<stdio.h>
#include<string.h>
#include<iostream>
#include<algorithm>
#include<functional>
using namespace std;
#define ll long long
#define maxn 500050
struct mat{
	int n,m;
	double a[15][15];
	mat(int n,int m):n(n),m(m){}
	mat(){}
	void init(){
		for(int i=0;i<n;i++)
			for(int j=0;j<m;j++)
				a[i][j]=0;
	}
}a;
double po[15],p,q;
int n,m;
mat mul(mat a,mat b){
	mat c(a.n,b.m);
	c.init();
	for(int i=0;i<c.n;i++)
		for(int j=0;j<c.m;j++)
			for(int k=0;k<a.m;k++)
				c.a[i][j]+=a.a[i][k]*b.a[k][j];
	return c;
}
double matPow(int n){ 
	mat b(a.n,1); b.init();
	for(int i=0;i<a.n;i++)
		b.a[i][0]=1;
	while(n){
		if(n&1)
			b=mul(a,b);
		a=mul(a,a);
		n/=2;
	}
	return b.a[0][0];
}
int main(void){
	while(scanf("%lf%lf%d%d",&p,&q,&n,&m)!=EOF){
		if(q==0){
			printf("0.000000\n");
			continue;
		}
		po[0]=1;
		a=mat(n,n);
		a.init();
		for(int i=1;i<=n;i++)
			po[i]=po[i-1]*q;
		for(int i=n-1;i>=0;i--)
			a.a[0][i]=po[i]*p;
		for(int i=1;i<n;i++)
			a.a[i][i-1]=1.0;
		printf("%f\n",1.0-matPow(m-n+1));
	}
	return 0;
}