JZOJ 5616. 【NOI2018模拟3.31】沧海尘记

Description

Input

Output

Data Constraint

Solution

• 观察题目，就可以得出有变换：

  $$A_i=\sum_{j=1}^{n}{A_j*P_{i,j}}$$

• 那么用 $O(N^3)$ 的高斯消元就可以通过 52% 的数据。

• 要 AC 此题，就要用到 MCMC（Markov Chain-Monte Carlo，⻢尔科夫链蒙特卡洛）算法 。

• 天啊！什么高大上的算法？

• 随机得出⼀个概率向量 $A$，计算 $A*P^k$，输出即可。

• 其实就是随意使某个位置为 $1$（例如使 $A_1=1$），其它位置为 $0$ 。

• 之后不断执行 $O(N^2)$ 的变换，$A$ 始终会变成所求的行向量（别问我为什么）。

• 要使用这个MCMC算法呢，就要什么变换需没周期性、不可约，具体可参考zhiyong_will的博客。

• 具体实现的时候可以写⼀个 $while$ 循环，直到

  $$Max^{n}_{i=1} |Ai − Bi| < 10^{−15}$$

• ⼀般来说 $k_{max}=12$ ，则时间复杂度 $O(k*N^2)$ 。

Code

#include<cstdio>
#include<cmath>
#include<cctype>
using namespace std;
typedef double DB;
const int N=2505;
const DB eps=1e-15;
int n,a,b,c;
int q[N];
DB p[N][N],f[N],g[N];
int main()
{
    scanf("%d%d%d%d",&n,&a,&b,&c);
    for(int i=1;i<=n;i++)
    {
        int s=q[1]=(a*q[n]+b)%c;
        for(int j=2;j<=n;j++) s+=q[j]=(a*q[j-1]+b)%c;
        s+=n;
        for(int j=1;j<=n;j++) p[i][j]=((DB)q[j]+1)/s;
    }
    f[1]=1;
    while(true)
    {
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++) g[j]+=f[i]*p[i][j];
        bool pd=true;
        for(int i=1;i<=n;i++)
            if(fabs(f[i]-g[i])>eps)
            {
                pd=false;
                break;
            }
        if(pd) break;
        for(int i=1;i<=n;i++) f[i]=g[i],g[i]=0;
    }
    for(int i=1;i<=n;i++) printf("%.18lf ",g[i]);
    return 0;
}