影の魔女 (AOJ2315)

nがk個・・・なんか正規分布に近づくとか・・・みたいにに騙されてはいけない(確信).
冷静に式に落とすと簡単.
でかいところではN-ボナッチみたいな遷移式(行列)になる(ので直近N個を縦ベクトルで見る)
小さいところははねかえるので線形代数です.

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;

typedef vector<double> vd;
typedef vector<vd> mat;
double eps=1e-9;
vd solve(mat a,vd b){
	int H=a.size(),W=a[0].size();
	bool used[101]={};
	rep(j,W){
		int i=0;
		while(i<H&&(used[i]||abs(a[i][j])<eps)) i++;
		if(i==H) continue;
		used[i]=1;
		double A=a[i][j];
		rep(l,W) a[i][l]/=A;
		b[i]/=A;
		rep(k,H){
			if(k==i) continue;
			double c=-a[k][j];
			rep(l,W) a[k][l]+=a[i][l]*c;
			b[k]+=b[i]*c;
		}
	}
	return b;
}
mat pro(mat a,mat b){
	int N=a.size();
	mat c(N,vd(N,0));
	rep(i,N) rep(j,N) rep(k,N) c[i][j]+=a[i][k]*b[k][j];
	return c;
}
mat ex(mat x,int p){
	int N=x.size();
	mat a(N,vd(N,0));
	rep(i,N) a[i][i]=1;
	while(p){
		if(p%2) a=pro(a,x);
		x=pro(x,x);
		p/=2;
	}
	return a;
}

double dp[11][101];
double p[101];
int main(){
	int S,n,k;
	cin>>S>>n>>k;
	if(S<0) S=-S;
	if(n==1&&S%k!=0){
		puts("-1");
		return 0;
	}
	dp[0][0]=1;
	rep(i,k) rep(j,n*k+1) rep1(l,n) if(dp[i][j]) dp[i+1][j+l]+=dp[i][j]/n;
	rep(i,n*k+1) p[i]=dp[k][i];
	int N=n*k;
	mat a(N,vd(N,0));
	vd b(N,0);
	rep1(i,N-1){
		rep1(j,N){
			a[i][abs(i-j)]+=p[j];
		}
		a[i][i]-=1;
		b[i]=-1;
	}
	a[0][0]=1,b[0]=0;
	b=solve(a,b);
	if(S<N){
		printf("%.12f\n",b[S]);
		return 0;
	}
	mat M(N+1,vd(N+1,0));
	rep(i,N) M[0][i]=p[i+1];
	M[0][N]=1;
	rep(i,N-1) M[i+1][i]=1;
	M[N][N]=1;
	reverse(all(b));
	b.pb(1);
	M=ex(M,S-N+1);
	double ans=0;
	rep(i,N+1) ans+=M[0][i]*b[i];
	printf("%.12f\n",ans);

}