LOJ 2304 「NOI2017」泳池——思路+DP+常系数线性齐次递推

题目：https://loj.ac/problem/2304

看了各种题解……

\( dp[i][j] \) 表示有 i 列、第 j 行及以下默认合法，第 j+1 行至少有一个非法格子的概率，满足最大合法矩形面积 <= lm。其中第 j 行及以下的部分的贡献是 1 而不是 q 的几次方。

那么有 \( dp[i][j]=dp[i][j+1]*p^i + \sum\limits_{k=1}^{i}dp[k-1][j+1]*p^{k-1}*(1-p)*dp[i-k][j] \)

注意到当 i>k 的时候，最底下一行必然有至少一个位置是非法的。所以令 \(ans_i\) 表示 i 列的概率，有 \( ans_i = \sum\limits_{j=1}^{i}ans_{j-1}*(1-p)*dp[i-j][1]*p^{i-j} \)

\(ans_i\) 的初值就是 dp[i][0] 。注意 dp[0][*]=1 。然后可以用常系数线性齐次递推的知识优化。

注意清空数组。注意别把 n 的值真的改掉。

#include<cstdio>

#include<cstring>

#include<algorithm>

#define ll long long

using namespace std;

const int N=,M=N<<,mod=;

int upt(int x){while(x>=mod)x-=mod;while(x<)x+=mod;return x;}

int pw(int x,int k)

{int ret=;while(k){if(k&)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=;}return ret;}

int n,q,q2,f[N][N],bin[N],a[N],ans[M],b[M],c[M],lm;

void Mul(int *u,int *v)//(lm-1)'

{

  memset(c,,sizeof c);

  for(int i=;i<lm;i++)

    for(int j=;j<lm;j++)

      c[i+j]=(c[i+j]+(ll)u[i]*v[j])%mod;

  for(int i=*(lm-);i>=lm;i--)

    if(c[i])

      for(int j=;j<=lm;j++)

    c[i-j]=(c[i-j]+(ll)c[i]*a[j])%mod;

  memcpy(u,c,sizeof *lm);//0~lm-1

}

int solve(int tmp)

{

  lm=tmp; memset(f,,sizeof f);

  for(int j=;j<=lm+;j++)f[][j]=;//lm+1 not lm!!!

  for(int i=;i<=lm;i++)

    for(int j=lm/i;j>=;j--)

      {

    int tp=(ll)f[i][j+]*bin[i]%mod;

    for(int k=;k<=i;k++)

      {

        int ml=(ll)f[k-][j+]*f[i-k][j]%mod;

        ml=(ll)ml*q2%mod*bin[k-]%mod;

        tp=upt(tp+ml);

      }

    f[i][j]=tp;

      }

  if(n<=lm)return f[n][]; lm++;

  for(int i=;i<=lm;i++)

    {

      int tp=(ll)f[i-][]*bin[i-]%mod;

      a[i]=(ll)tp*q2%mod;//not lm-i

    }

  memset(ans,,sizeof ans);////

  memset(b,,sizeof b);////

  ans[]=b[]=; int tn=n;//////

  while(tn)

    {

      if(tn&)Mul(ans,b); Mul(b,b); tn>>=;

    }

  int ret=;

  for(int i=;i<lm;i++)

    ret=(ret+(ll)ans[i]*f[i][])%mod;

  return ret;

}

int main()

{

  int x,y,k;scanf("%d%d%d%d",&n,&k,&x,&y);

  q=(ll)x*pw(y,mod-)%mod; q2=upt(-q);

  bin[]=;

  for(int i=;i<=k;i++)bin[i]=(ll)bin[i-]*q%mod;

  printf("%d\n",upt(solve(k)-solve(k-)));

  return ;

}