UOJ #420. 【集训队作业2018】矩形
题目链接
题目大意
\[F(i,j)=
\begin{cases}
0 & i = 0 \\
g_i & i>0,j=0 \\
aF(i-1,j)+bF(i,j-1)+c & i>0,j>0
\end{cases}
\]
给定 \(n,m,p\)，\(a,b,c\)，\(g_1,g_2,..,g_n\)，求
\[\left(\sum_{i=1}^n\sum_{j=1}^m F(i,j)h^{(i-1)m+j-1}\right)\!\!\!\mod p
\]
\(n\leq 10^6,m\leq 10^9,p\in \text{prime}\)
思路
设 \(\displaystyle G(x)=\sum_{i=1}^n g_ix^i\)：
\[\begin{align}
f_{i,j} &= af_{i-1,j}+bf_{i,j-1}+c \\
\sum_{i\geq 1}\sum_{j\geq 1}f_{i,j}x^iy^j &= ax\sum_{i\geq 0}\sum_{j\geq 1}f_{i,j}x^iy^j + by\sum_{i\geq 1}\sum_{j\geq 0}f_{i,j}x^iy^j+c\sum_{i\geq 1}\sum_{j\geq 1}x^iy^j \\
F(x,y)-G(x) &= ax(F(x,y)-G(x))+byF(x,y)+c(\frac{1}{1-x}-1)(\frac{1}{1-y}-1) \\
F(x,y) &= (1-ax-by)^{-1}\left((1-ax)G(x)+\frac{cxy}{(1-x)(1-y)}\right) \\
&= \sum_{i\geq 0}\sum_{j \geq 0}\binom{-1}{i\;\;j}(-a)^i(-b)^jx^iy^j\left((1-ax)G(x)+\frac{cxy}{(1-x)(1-y)}\right) \\
&= \sum_{i\geq 0}\sum_{j \geq 0}\binom{i+j}{i}a^ib^jx^iy^j\left((1-ax)G(x)+\frac{cxy}{(1-x)(1-y)}\right) \\
[x^uy^v]F(x,y) &= \sum_{i=0}^u\sum_{j=0}^v \binom{i+j}{i}a^ib^j[x^{u-i}y^{v-j}]\left((1-ax)G(x)+\frac{cxy}{(1-x)(1-y)}\right) \\
&= b^v\sum_{i=0}^ua^i\left(\binom{i+v}{i}-\binom{i+v-1}{i-1}\right)g_{u-i}+c\sum_{i=0}^{u-1}\sum_{j=0}^{v-1}\binom{i+j}{i}a^ib^j \\
&= b^v\sum_{i=0}^ua^i\binom{i+v-1}{v-1}g_{u-i}+c\sum_{i=0}^{u-1}\sum_{j=0}^{v-1}\binom{i+j}{i}a^ib^j
\end{align}
\]
设 \(g_i\) 所贡献的系数为 \(p_i\)：
\[p_i=\sum_{x=i}^n\sum_{y=1}^ma^{x-i}\binom{x-i+y-1}{y-1}b^y\cdot h^{(x-1)m+y-1}
\]
另设 \(q(x)\) 来计算后面的和式
\[\begin{align}
q(x)=&\sum_{y=0}^{m-1}\binom{x+y}{y}(bh)^y\\
=& \sum_{y=0}^{m-1}\left(\binom{x+y-1}{y}+\binom{x+y-1}{y-1}\right)(bh)^y\\
=& q(x-1)+bh(q(x)-\binom{x+m-1}{m-1}(bh)^{m-1})
\end{align}
\]
\[\therefore q(x)=\begin{cases}
\;1 & x < 0 \\
\displaystyle\frac{1}{1-bh}\left(q(x-1)-\binom{x+m-1}{m-1}(bh)^m\right) & x\geq 0,bh\neq 1\\
\displaystyle\binom{x+m}{m-1} & x\geq 0, bh=1
\end{cases}
\]
则
\[\begin{align}
p_i&=\sum_{x=i}^na^{x-i}h^{(x-1)m-1}\sum_{y=1}^m\binom{x-i+y-1}{y-1}(bh)^y \\
&= bh\sum_{x=i}^n a^{x-i}h^{(x-1)m-1}q(x-i)\\
&= b\cdot h^{(i-1)m-1}\sum_{x=0}^{n-i}(ah^m)^xq(x)
\end{align}
\]
于是 \(O(n)\) 预处理前缀和后计算 \(\displaystyle \sum_{i=1}^n p_ig_i\) 即可。
考虑式子后半部分的贡献，即
\[\begin{align}
&\sum_{x=1}^n\sum_{y=1}^mh^{(x-1)m+y-1}\cdot c\sum_{i=0}^{x-1}\sum_{j=0}^{y-1}\binom{i+j}{i}a^ib^j \\
=& c\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}\binom{i+j}{i}a^ib^j\left(\sum_{x=i}^{n-1}h^{mx}\right)\left(\sum_{y=j}^{m-1}h^{y}\right) \\
=& c\sum_{i=0}^{n-1}a^i\sum_{j=0}^{m-1}\binom{i+j}{i}b^j\cdot\frac{h^{mn}-h^{mi}}{h^m-1}\cdot\frac{h^m-h^j}{h-1} \\
=& \frac{c}{(h-1)(h^m-1)}\sum_{i=0}^{n-1}a^i(h^{mn}-h^{mi})\sum_{j=0}^{m-1}\binom{i+j}{j}b^j(h^m-h^j)
\end{align}
\]
将第二个和式拆开：\(\displaystyle h^m\sum_{j=0}^{m-1}\binom{i+j}{j}b^j-\sum_{j=0}^{m-1}\binom{i+j}{j}(bh)^j\) ，两项都与 \(q(x)\) 的形式类似，可以使用相同方法计算。不过注意到推导只对 \(h\neq 1\) 适用，\(h=1\) 时是另外一个式子：
\[\begin{align}
& c\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}\binom{i+j}{i}a^ib^j\left(\sum_{x=i}^{n-1}h^{mx}\right)\left(\sum_{y=j}^{m-1}h^{y}\right) \\
= & c\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}\binom{i+j}{i}a^ib^j(n-i)(m-j)\\
=& c\sum_{i=0}^{n-1}a^i(n-i)\sum_{j=0}^{m-1}\binom{i+j}{j}b^j(m-j)
\end{align}
\]
拆开后我们得到一个 \(q(x)\) 形式的式子和一个 \(\displaystyle \sum_{j=0}^{m-1}\binom{i+j}{j}b^jj\) ，考虑类似的递推：
\[\begin{align}
s(x)&=\sum_{j=0}^{m-1}\binom{i+j}{j}b^jj \\
&= \sum_{i=0}^{m-1}\left(\binom{x+i-1}{i}+\binom{x+i-1}{i-1}\right)b^i\cdot i\\
&= s(x-1)+b\left(s(x)-\binom{x+m-1}{m-1}b^{m-1}m+q_b(x)\right) \\
&= \frac{1}{1-b}\left(s(x-1)+b\cdot q_b(x)-\binom{x+m-1}{m-1}b^mm\right)
\end{align}
\]
\(q_b(x)\) 表示将前面的 \((bh)^y\) 换成 \(b^y\) 对应的式子。而当 \(b=1\) 时
\[\begin{align}
s(x)&=\sum_{i=0}^{m-1}\binom{x+i}{i}i\\
&= \sum_{i=1}^{m-1}\sum_{j=i}^{m-1}\binom{x+j}{j}\\
&= \sum_{i=1}^{m-1}\left(\binom{x+m}{m-1}-\binom{x+i}{i-1}\right)\\
&= (m-1)\binom{x+m}{m-1}-\sum_{i=0}^{m-2}\binom{x+i+1}{i}\\
&= (m-1)\binom{x+m}{m-1}-\binom{x+m}{m-2}
\end{align}
\]
综上时间复杂度为 \(O(n+\log m)\) 。
Code
#include<iostream>
#define rep(i,a,b) for(int i = (a); i <= (b); i++)
#define per(i,b,a) for(int i = (b); i >= (a); i--)
#define ll long long
#define N 1001000
using namespace std;
int mod;
ll qpow(ll a, int b){
ll ret = 1;
for(; b; b >>= 1){ if(b&1) (ret *= a) %= mod; (a *= a) %= mod; }
return ret;
}
int T, n, m, h;
int a, b, c;
ll g[N];
ll fact[2][N], invf[2][N];
ll q0[N], q1[N];
ll pref[N], p[N], s[N];
ll ans;
void prework(){
fact[0][0] = 1;
rep(i,1,N-1) fact[0][i] = fact[0][i-1] * i % mod;
invf[0][N-1] = qpow(fact[0][N-1], mod-2);
per(i,N-2,0) invf[0][i] = invf[0][i+1] * (i+1) % mod;
fact[1][0] = max(m-2, 1);
rep(i,m-1,m+N-3) fact[1][i-m+2] = fact[1][i-m+1] * max(i, 1) % mod;
invf[1][N-1] = qpow(fact[1][N-1], mod-2);
per(i,m+N-4,m-2) invf[1][i-m+2] = invf[1][i-m+3] * (i+1) % mod;
}
ll C(int n, int m){
if(n < m || m < 0) return 0;
return fact[1][n-::m+2] * invf[1][m-::m+2] % mod * invf[0][n-m] % mod;
}
void cal(ll *q, ll k){
if(k == 0){ rep(x,0,n) q[x] = 1; return; }
if(k == 1){
rep(x,0,n) q[x] = C(x+m, m-1);
return;
}
ll div = qpow(1+mod-k, mod-2), pow = qpow(k, m);
rep(x,0,n) q[x] = div * ((x ? q[x-1] : 1) + mod - C(x+m-1, m-1) * pow % mod) % mod;
}
int main(){
ios::sync_with_stdio(false);
cin>>T>>n>>m>>h>>mod;
cin>>a>>b>>c;
rep(i,1,n) cin>>g[i];
if(h == 0){
cout<< (b*g[1] + c) % mod <<endl;
return 0;
}
prework();
cal(q0, b), cal(q1, (ll)b*h % mod);
ll coef = a * qpow(h, m) % mod, cur = 1;
rep(x,0,n){
pref[x] = ((x ? pref[x-1] : 0) + cur * q1[x] % mod * b % mod * h) % mod;
(cur *= coef) %= mod;
}
coef = qpow(h, m), cur = qpow(h, mod-2) % mod;
rep(i,1,n){
p[i] = cur * pref[n-i] % mod, (ans += p[i] * g[i]) %= mod;
(cur *= coef) %= mod;
}
if(h == 1){
if(b == 1) rep(x,0,n)
s[x] = (C(x+m, m-1) * (m-1) + mod - C(x+m, m-2)) % mod;
else{
ll div = qpow(1+mod-b, mod-2), pow = qpow(b, m);
rep(x,0,n) s[x] = div * ((x ? s[x-1] : 0) + b * q0[x] % mod + mod - C(x+m-1, m-1) * pow % mod * m % mod) % mod;
}
ll sum = 0;
for(ll i = 0, pow = 1; i < n; i++, (pow *= a) %= mod)
(sum += pow * (n-i) % mod * (m * q0[i] % mod + mod - s[i])) %= mod;
(ans += c * sum) %= mod;
} else{
ll pow = qpow(h, m), cur = 1, tot = qpow(h, (ll)n*m % (mod-1));
ll sum = 0;
for(ll i = 0, pa = 1; i < n; i++, (pa *= a) %= mod, (cur *= pow) %= mod)
(sum += pa * (tot + mod - cur) % mod * (pow * q0[i] % mod + mod - q1[i])) %= mod;
(ans += sum * c % mod * qpow(h-1, mod-2) % mod * qpow(pow-1, mod-2)) %= mod;
}
cout<< ans <<endl;
return 0;
}
