bzoj 4451: [Cerc2015]Frightful Formula 数学+排列组合

题意

给你一个n*n矩阵的第一行和第一列，其余的数通过如下公式推出：
F[i,j]=a*f[i,j-1]+b*f[i-1,j]+c
求f[n][n]%(10^6+3)
2<=n<=200000
其余的数大于等于0小于等于10^6

分析

真是一道推式子的好题。。。

假设c=0，那么我们就考虑矩阵第一行和第一列对答案的贡献。
设某个点(i,j)走到(n,n)，要向上走x步，向右走y步，也就是x=n-i,y=n-j
那么显然其对答案的贡献就是

$$C_{x+y}^xb^xa^y*matrix[i][j]$$
那么我们就可以单独算出第一行和第一列的贡献。
接下来考虑常数项c对答案的贡献。
显然我们要求的是 $$c*\sum_{i=0}^{n-2}\sum_{j=0}^{n-2}C_{i+j}^ib^ia^j$$
将组合数拆开 $$c*\sum_{i=0}^{n-2}\sum_{j=0}^{n-2}\frac{(i+j)!}{i!j!}b^ia^j$$
$$=c*\sum_{i=0}^{n-2}\sum_{j=0}^{n-2}(i+j)!\frac{b^i}{i!}\frac{a^j}{j!}$$
推到这步就可以发现，这是两个指数生成函数的卷积的系数和，就可以用fft来做了，但是我们的梦想是O(n)!那么久接着推
我们考虑斜着算贡献，也就是先算第一斜行的总贡献，再转移到第二斜行的总贡献。
显然第一斜行的总贡献为1
显然在主对角线之前（包括主对角线）的任意斜行的总贡献为其上一斜行的总贡献*(a+b)
那么只要一次算出没斜行的总贡献即可。但考虑过了主对角线之后会有多余的贡献，那么就要减去这些多余的贡献。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;
const int N=400005;
const int MOD=1000003;

int n,a,b,c,jc[N],ny[N],ma[N],mb[N];

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

int ksm(int x,int y)
{
    int ans=1;
    while (y)
    {
        if (y&1) ans=(LL)ans*x%MOD;
        x=(LL)x*x%MOD;y>>=1;
    }
    return ans;
}

void prework()
{
    jc[0]=ny[0]=1;
    for (int i=1;i<=400000;i++)
    {
        jc[i]=(LL)jc[i-1]*i%MOD;
        ny[i]=ksm(jc[i],MOD-2);
    }
}

int main()
{
    prework();
    while (scanf("%d%d%d%d",&n,&a,&b,&c)!=EOF)
    {
        int ans=0;
        ma[0]=mb[0]=1;
        for (int i=1;i<=n;i++) ma[i]=(LL)ma[i-1]*a%MOD,mb[i]=(LL)mb[i-1]*b%MOD;
        for (int i=1;i<=n;i++)
        {
            int w=read();
            int x=n-i,y=n-1;
            if (i>1) ans=(ans+(LL)mb[x]*ma[y]%MOD*jc[x+y-1]%MOD*ny[x]%MOD*ny[y-1]%MOD*w%MOD)%MOD;
        }
        for (int i=1;i<=n;i++)
        {
            int w=read();
            int x=n-1,y=n-i;
            if (i>1) ans=(ans+(LL)mb[x]*ma[y]%MOD*jc[x+y-1]%MOD*ny[x-1]%MOD*ny[y]%MOD*w%MOD)%MOD;
        }
        n--;
        int now=1,w=(a+b)%MOD;
        ans=(ans+c)%MOD;
        for (int i=1;i<n;i++)
        {
            now=(LL)now*w%MOD;
            ans=(ans+(LL)now*c%MOD)%MOD;
        }
        int x1=0,y1=n-1,x2=n-1,y2=0;
        for (int i=1;i<n;i++)
        {
            int w1=(LL)jc[x1+y1]%MOD*ny[x1]%MOD*ny[y1]%MOD*mb[x1]%MOD*ma[y1]%MOD;
            int w2=(LL)jc[x2+y2]%MOD*ny[x2]%MOD*ny[y2]%MOD*mb[x2]%MOD*ma[y2]%MOD;
            now=(now-(LL)w1+MOD)%MOD;
            now=(now-(LL)w2+MOD)%MOD;
            now=(LL)now*w%MOD;
            now=(now+(LL)w1*b%MOD)%MOD;
            now=(now+(LL)w2*a%MOD)%MOD;
            ans=(ans+(LL)now*c%MOD)%MOD;
            x1++;y2++;
        }
        printf("%d\n",ans);
    }
    return 0;
}