影の魔女 (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); }