1. 程式人生 > >【BZOJ2510】弱題

【BZOJ2510】弱題

void 時間 1-n [] fixed 矩陣快速冪 1.0 for 要求

題目大意

\(M\)個球,一開始每個球均有一個初始標號,標號範圍為\(1-N\)且為整數,標號為i的球有\(a_i\)個,並保證\(\sum a_i=M\)
  每次操作等概率取出一個球(即取出每個球的概率均為\(1/M\)),若這個球標號為\(k(k < N)\),則將它重新標號為\(k + 1\);若這個球標號為\(N\),則將其重標號為\(1\)。(取出球後並不將其丟棄)
  現在你需要求出,經過K次這樣的操作後,每個標號的球的期望個數。

\(N ≤ 1000, M ≤ 100,000,000, K ≤ 2,147,483,647\)

題目分析

遞推方程很好想,但是需要優化。
由於\(K\)很大,考慮使用矩陣乘法優化\(dp\)

。但是\(n\)的範圍太大,不能直接化出矩陣來相乘。通過觀察矩陣,我們很容易發現矩陣的每一行之間是循環的,因此,我們可以只算一行即可,矩陣相乘的時間負責度瞬間降低為\(O(n^2)\),矩陣快速冪優化後總時間復雜為\(O(n^2log_2 n)\)

#include <bits/stdc++.h>
using namespace std;
int n,m,k;
namespace Task1{
    double p[2][1005];
    void solve(){
        for(int i=1;i<=n;i++)scanf("%lf",&p[0][i]);
        int cur=0,pre;
        for(int i=1;i<=k;i++){
            pre=cur;cur^=1;
            for(int j=1;j<=n;j++){
                if(j==1)p[cur][j]=p[pre][j]+p[pre][n]/m-p[pre][j]/m;
                else if(j==n)p[cur][j]=p[pre][j]+p[pre][j-1]/m-p[pre][j]/m;
                else p[cur][j]=p[pre][j]+p[pre][j-1]/m-p[pre][j]/m;
            }
        }
        for(int i=1;i<=n;i++)cout<<fixed<<setprecision(3)<<p[cur][i]<<"\n";
    }
}
namespace Task2{
    double tmp[1005][1005];
    struct node{
        double a[1005];
        node(){memset(a,0,sizeof(a));}
        double&operator[](int x){return a[x];}
        node operator*(node &b){
            node c;
            for(int i=1;i<=n;i++)
                for(int j=1;j<=n;j++){
                    int pos=i+j-1;
                    if(pos>n)pos-=n;
                    tmp[j][pos]=b[i];
                }
            for(int i=1;i<=n;i++)
                for(int j=1;j<=n;j++)
                    c[i]+=a[j]*tmp[j][i];
            return c;
        }
        node operator^(int cnt){
            node ret,mul=*this;
            bool flag=0;
            for(;cnt;cnt>>=1,mul=mul*mul)if(cnt&1){if(!flag)ret=mul,flag=1;else ret=ret*mul;}
            return ret;
        }
    }p;
    void solve(){
        for(int i=1;i<=n;i++)scanf("%lf",&p[i]);
        node mul;
        mul[1]=(m-1)/double(m);
        mul[2]=1.0/m;
        mul=mul^k;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++){
                int pos=i+j-1;
                if(pos>n)pos-=n;
                tmp[j][pos]=mul[i];
            }
        node ans;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                ans[i]+=p[j]*tmp[j][i];
        for(int i=1;i<=n;i++)printf("%.3lf\n",ans[i]);
    }
}
int main(){
    scanf("%d%d%d",&n,&m,&k);
    if(1ll*n*k<=5e7)Task1::solve();
    else Task2::solve();
}
/*
[(m-1)/m,1/m,0,0,0]
[0,(m-1)/m,1/m,0,0]
[0,0,(m-1)/m,1/m,0]
[0,0,0,(m-1)/m,1/m]
[1/m,0,0,0,(m-1)/m]

[a1,a2,a3] [a1,a2,a3] [a1*a1+a2*a3+a3*a2,a1*a2+a2*a1+a3*a3,a1*a3+a2*a2+a3*a1]
[a3,a1,a2]*[a3,a1,a2]=[a3*a1+a1*a3+a2*a2,......
[a2,a3,a1] [a2,a3,a1] [......

[a1,a2,a3] [b1,b2,b3] [a1*b1+a2*b3+a3*b2,a1*b2+a2*b1+a3*b3,a1*b3+a2*b2+a3*b1]
[a3,a1,a2]*[b3,b1,b2]=[a3*b1+a1*b3+a2*b2,......
[a2,a3,a1] [b2,b3,b1] [......

[a1,a2,a3]*[b1,b2,b3]=>[a1*b1+a2*b3+a3*b2,a1*b2+a2*b1+a3*b3,a1*b3+a2*b2+a3*b1]
*/

【BZOJ2510】弱題