JZOJ5616. 【NOI2018模拟3.31】沧海尘记

问题描述

南阳渡旁有一个饿得发昏的道士，详聊过后，才知道他想要兜售身上的道书以求温 饱。但是他所开出的价码，实在让人不敢领教……
买下道书后，第⼀⻚上的内容就令你⼤吃⼀惊—— 如果把⼈按照拥有财富的多少分为 n 等，那么对于第 i 等的⼈，他的下⼀代就会有 Pi,j 的 概率成为第 j 等的⼈。对于⼀个居⺠点，如果某⼀代的贫富分布和下⼀代的贫富分布相同，那 么就称之为概率 P 下平衡的分布。 给定概率 P，你需要给出某⼀代的贫富分布 A（其中 Ai 表示：令这⼀代的总⼈数为 1，则 有 Ai ⼈是第 i 个财富等级），使得 A 是概率 P 下平衡的分布。 具体地，给定⼀个 n 阶转移概率矩阵 P，你需要求出⼀个 1×n 的⾏向量 A，使得 AP = A。

题解

如果知道：Markov Chain-Monte Carlo，⻢尔科夫链蒙特卡洛
那么答案很自然就会想到。
先随机一个A，最后就输出$AP^k$
时间复杂度是$O(kn^2)$

至于k到底是多少，
就直接循环，直到$AP^k$合法为止。
要注意精度误差。

code

#include <queue>
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <string.h>
#include <cmath>
#include <math.h>
#include <time.h>
#define ll long long
#define N 2503
#define M 103
#define db double
#define P putchar
#define G getchar
#define inf 998244353
#define pi 3.1415926535897932384626433832795
#define jd 0.000000000000001
using namespace std;
char ch;
void read(int &n)
{
    n=0;
    ch=G();
    while((ch<'0' || ch>'9') && ch!='-')ch=G();
    ll w=1;
    if(ch=='-')w=-1,ch=G();
    while('0'<=ch && ch<='9')n=(n<<3)+(n<<1)+ch-'0',ch=G();
    n*=w;
}

int max(int a,int b){return a>b?a:b;}
int min(int a,int b){return a<b?a:b;}
ll abs(ll x){return x<0?-x:x;}
ll sqr(ll x){return x*x;}
void write(ll x){if(x>9) write(x/10);P(x%10+'0');}

long db a[N],p[N][N],b[N];
int n,q[N][N],s[N],s1,s2,s3;

bool pd()
{
    long db s=0;
    for(int i=1;i<=n;i++)
    {
        if(fabs(a[i]-b[i])>jd || a[i]<0 || a[i]>1)return 0;
        s+=a[i];
    }
    if(fabs(s-1)>jd)return 0;
        else return 1;
}

void work()
{
    memset(b,0,sizeof(b));
    for(int i=1;i<=n;i++) 
        for(int j=1;j<=n;j++)
            b[j]+=a[i]*p[i][j];
}

int main()
{
    freopen("dust.in","r",stdin);
    freopen("dust.out","w",stdout);

    read(n);read(s1);read(s2);read(s3);
    for(int i=1;i<=n;i++)
    {
        q[i][1]=(s1*q[i-1][n]+s2)%s3;
        s[i]=q[i][1]+1;
        for(int j=2;j<=n;j++)
            q[i][j]=(q[i][j-1]*s1+s2)%s3,s[i]+=q[i][j]+1;
        for(int j=1;j<=n;j++)
            p[i][j]=(db)(q[i][j]+1)/s[i];
    }

    for(b[1]=1;!pd();work())
        swap(a,b);

    for(int i=1;i<=n;i++)
        printf("%.20lf ",(db)a[i]);

    return 0;
}