Luogu4351 Frightful Formula

Arahc 11.0
  2022-03-09 08:00 2022-03-09 08:00 创建   2023-03-15 03:37:40 2023-03-15 03:37:40 更新      1.4k 字  8 分钟  

题面

简化题意

对于矩阵 ff，给定 n,A,B,C,f1∼n,1,f1,1∼nn,A,B,C,f_{1\sim n,1},f_{1,1\sim n}。有：

fi,j=Afi,j−1+Bfi−1,j+Cf_{i,j} = Af_{i,j-1}+Bf_{i-1,j} +C

求 fn,n mod (106+3)f_{n,n}\bmod (10^6+3)。

n≤2×105,A,B,C≤106n\leq 2\times 10^5,A,B,C\leq 10^6。

题解

把原贡献拆成两个部分：边界上的点（不考虑 CC 的贡献），还有非边界的点（考虑 CC 的贡献）。

对于前者，考虑 (1,x)(1,x) 的贡献是什么，显然要从 (2,x)(2,x) 开始转移，(2,x)→(n,n)(2,x)\to(n,n) 有 (2n−x−2n−2)\binom{2n-x-2}{n-2} 条路径。

这些路径横着走的长度和竖着走的长度分别一样。也就是算上了 n−1n-1 次 BB 和 n−xn-x 次 AA。于是 (1,x)(1,x) 的贡献为：

f1,x(2n−x−2n−2)An−xBn−1f_{1,x}\binom{2n-x-2}{n-2} A^{n-x}B^{n-1}

(x,1)(x,1) 贡献同理。

然后考虑中间的点即 CC 的贡献。考虑对于一个非边界的点 (x,y)→(n,n)(x,y)\to (n,n) 有 (2n−x−yn−x)\binom{2n-x-y}{n-x} 条路径。所有非边界的点的贡献就是：

C∑i=2n∑j=2n(2n−i−jn−i)An−jBn−iC\sum_{i=2}^n\sum_{j=2}^n \binom{2n-i-j}{n-i} A^{n-j}B^{n-i}

礼貌性拆一下组合数：

C∑i=2n∑j=2n(2n−i−j)!Bn−i(n−i)!An−j(n−j)!C\sum_{i=2}^n \sum_{j=2}^n (2n-i-j)!\frac{B^{n-i}}{(n-i)!}\frac{A^{n-j}}{(n-j)!}

这是一个卷积的形式，模数不好，要写 MTT。复杂度为 nlog⁡nn\log n。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#include<bits/stdc++.h>
#define int long long
#define double long double
using namespace std;
const int max_n=2000006,mod=1000003;
inline int read(){
    int x=0;bool w=0;char c=getchar();
    while(c<'0' || c>'9') w|=c=='-',c=getchar();
    while(c>='0' && c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return w?-x:x;
}
inline void write(int x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10^48);
}
inline int mi(int a,int p=mod-2){
    int res=1;
    while(p){
        if(p&1) res*=a,res%=mod;
        a*=a,a%=mod,p>>=1;
    }
    return res;
}
namespace DXS{
    struct com{
        double a,b; // a+bi
        com(double x=0,double y=0){a=x,b=y;}
        com operator + (const com b) const{
            com a=*this,res;
            res.a=a.a+b.a,res.b=a.b+b.b;
            return res;
        }
        com operator - (const com b) const{
            com a=*this,res;
            res.a=a.a-b.a,res.b=a.b-b.b;
            return res;
        }
        com operator * (const com b) const{
            com a=*this,res;
            res.a=a.a*b.a-a.b*b.b,
            res.b=a.a*b.b+a.b*b.a;
            return res;
        }
    };
    int rev[max_n*8];
    const double Pi=acos(-1.0);
    inline void FFT(com *g,int n,bool op){
        for(register int i=0;i<n;++i) if(i<rev[i])
            swap(g[i],g[rev[i]]);
        for(register int mid=1;mid<n;mid<<=1){
            com omega(cos(Pi/mid),sin(Pi/mid)*(op?1:-1));
            for(register int i=0;i<n;i+=(mid<<1)){
                com w(1,0);
                for(register int j=0;j<mid;++j,w=w*omega){
                    com x=g[i+j],y=w*g[i+j+mid];
                    g[i+j]=x+y,
                    g[i+j+mid]=x-y;
                }
            }
        }
    }
    const int K=32768,S=K*K;
    inline void mulMTT(int *f,int *g,int n,int m,int mod){
        static com P[max_n],Q[max_n],R[max_n];
        for(register int i=0;i<n;++i){
            P[i].a=f[i]/K,
            P[i].b=f[i]%K;
            R[i].a=f[i]/K,
            R[i].b=-f[i]%K;
        }
        for(register int i=0;i<m;++i){
            Q[i].a=g[i]/K,
            Q[i].b=g[i]%K;
        }
        int len=1;
        while(len<=n+m) len<<=1;
        for(register int i=0;i<len;++i)
            rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
        FFT(P,len,1),FFT(Q,len,1),FFT(R,len,1);
        for(register int i=0;i<len;++i){
            Q[i].a/=len,
            Q[i].b/=len;
        }
        for(register int i=0;i<len;++i){
            P[i]=P[i]*Q[i],
            R[i]=R[i]*Q[i];
        }
        FFT(P,len,0),FFT(R,len,0);
        for(register int i=0,a1b1=0,a1b2=0,a2b1=0,a2b2=0;i<len;++i){
            a1b1=(int)floor((P[i].a+R[i].a)/2+0.5)%mod,
            a1b2=(int)floor((P[i].b+R[i].b)/2+0.5)%mod,
            a2b1=((int)floor(P[i].b+0.5)-a1b2)%mod,
            a2b2=((int)floor(R[i].a+0.5)-a1b1)%mod;
            f[i]=(a1b1*S%mod+(a1b2+a2b1)%mod*K%mod+a2b2)%mod;
            f[i]=(f[i]%mod+mod)%mod;
        }
    }
}

int n,A,B,c,ans;
int fac[max_n*2],inv[max_n*2];
int f[max_n*8],g[max_n*8];

inline int C(int n,int m){
    if(n<0 || m<0 || n<m) return 0;
    return fac[n]*inv[m]%mod*inv[n-m]%mod;
}

signed main(){
    n=read(),A=read(),B=read(),c=read();
    fac[0]=inv[0]=1;
    for(register int i=1;i<=(n<<1);++i)
        fac[i]=fac[i-1]*i%mod;
    inv[n<<1]=mi(fac[n<<1]);
    for(register int i=(n<<1)-1;i;--i)
        inv[i]=inv[i+1]*(i+1)%mod;
    for(register int i=1;i<=n;++i){
        int x=read();
        if(i==1) continue;
        ans=(ans+C(2*n-i-2,n-2)*mi(A,n-1)%mod*mi(B,n-i)%mod*x%mod)%mod;
    }
    for(register int i=1;i<=n;++i){
        int x=read();
        if(i==1) continue;
        ans=(ans+C(2*n-i-2,n-2)*mi(A,n-i)%mod*mi(B,n-1)%mod*x%mod)%mod;
    }
    for(register int i=2;i<=n;++i){
        f[i]=mi(A,n-i)*inv[n-i]%mod,
        g[i]=mi(B,n-i)*inv[n-i]%mod;
    }
    DXS::mulMTT(f,g,n+1,n+1,mod);
    for(register int i=4;i<=(n<<1);++i)
        ans=(ans+c*fac[2*n-i]%mod*f[i]%mod)%mod;
    write(ans);
    return 0;
}
  • 标题: Luogu4351 Frightful Formula
  • 作者: Arahc
  • 创建于 : 2022-03-09 08:00:00
  • 更新于 : 2023-03-15 03:37:40
  • 链接: https://arahc.github.io/2022/03/09/【题解】Luogu4351-Frightful-Formula/
  • 版权声明: 本文章采用 CC BY-NC-SA 4.0 进行许可。
评论
目录
Luogu4351 Frightful Formula
  1. 题面
  2. 题解
  3. 代码