可以发现两个矩阵可以凑成一个 $1\times1$ 的 $3$，因此后续只需在 $\bmod 3$ 意义下考虑。
另外大矩阵可以由小矩阵组成，因此只需考虑小矩阵即可。

\[\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\]

可以证明，充要条件是每行每列的和 $\bmod 3=0$。
必要性是注意到加上矩阵后，每行、列的和都不变。
充分性可以采用每次删掉一行来证明。

现在题意转化成了 $n\times m$ 的矩阵，其中 $a_{i,j}\in[0,\ell]$，求有多少个矩阵满足其每行每列的和 $\bmod 3=0$。

\[\sum_a \prod_{i=1}^n\left[\left(\sum_{j=1}^ma_{i,j}\right)\bmod 3=0\right] \prod_{j=1}^m\left[\left(\sum_{i=1}^na_{i,j}\right )\bmod 3=0\right]\]

我们想起来单位根反演。

\[[x\bmod p=0]=\frac1p\sum_{k=0}^{p-1}\omega_p^{kx}\] \[\begin{aligned} \text{answer}&=\sum_a \prod_{i=1}^n\left[\left(\sum_{j=1}^ma_{i,j}\right)\bmod 3=0\right] \prod_{j=1}^m\left[\left(\sum_{i=1}^na_{i,j}\right )\bmod 3=0\right]\\ &=\frac1{3^{n+m} }\sum_a \left(\prod_{i=1}^n\sum_{k=0}^2\omega_3^{k\sum\limits_{j=1}^ma_{i,j} }\right) \left(\prod_{j=1}^m\sum_{k=0}^2\omega_3^{k\sum\limits_{i=1}^na_{i,j} }\right)\\ &=\frac1{3^{n+m} }\sum_a \left(\prod_{i=1}^n\sum_{k=0}^2\prod_{j=1}^m\omega_3^{ka_{i,j} }\right) \left(\prod_{j=1}^m\sum_{k=0}^2\prod_{i=1}^n\omega_3^{ka_{i,j} }\right)\\ &=\frac1{3^{n+m} }\sum_a \left(\sum_{c_i\in[0,2]}\prod_{i=1}^n\prod_{j=1}^m\omega_3^{ka_{i,j} }\right) \left(\sum_{r_j\in[0,2]}\prod_{j=1}^m\prod_{i=1}^n\omega_3^{ka_{i,j} }\right)\\ &=\frac1{3^{n+m} }\sum_a\sum_{c_i,r_j\in[0,2]}\prod_{i=1}^n\prod_{j=1}^m\omega_3^{(c_i+r_j)a_{i,j} }\\ &=\frac1{3^{n+m} }\sum_{c_i,r_j\in[0,2]}\prod_{i=1}^n\prod_{j=1}^m\sum_{a_{i,j}=0}^\ell\omega_3^{(c_i+r_j)a_{i,j} }\\ \end{aligned}\]

这样，我们终于把 $a_{i,j}$ 拆出来了，只需令

\[f(t)=\sum_{i=0}^\ell w_3^{ti}\]

那么答案就是

\[\frac1{3^{n+m} }\sum_{c_i,r_j\in[0,2]}\prod_{i=1}^n\prod_{j=1}^mf(c_i+r_j)\\\]

注意到答案只和 $c_i,r_j$ 中 $0,1,2$ 的数量有关，暴力枚举，复杂度 $\mathcal O(n^4)$。

如果只枚举 $r$，令 $d_0,d_1,d_2$ 表示 $r_j$ 中 $0,1,2$ 的数量，那么

\[\begin{aligned} \text{answer}&=\frac1{3^{n+m} }\sum_{c_i,r_j\in[0,2]}\prod_{i=1}^n\prod_{j=1}^mf(c_i+r_j)\\ &=\frac1{3^{n+m} }\sum_{c_i\in[0,2]}\sum_{d_0+d_1+d_2=m}\binom m{d_0,d_1,d_2}\prod_{i=1}^n\prod_{j=0}^2f^{d_j}(c_i+j)\\ \end{aligned}\]

我们记

\[g(t,d_0,d_1,d_2)=\prod_{j=0}^2f^{d_j}(t+j)\] \[\begin{aligned} \text{answer}&=\frac1{3^{n+m} }\sum_{c_i\in[0,2]}\sum_{d_0+d_1+d_2=m}\binom m{d_0,d_1,d_2}\prod_{i=1}^n\prod_{j=0}^2f^{d_j}(c_i+j)\\ &=\frac1{3^{n+m} }\sum_{c_i\in[0,2]}\sum_{d_0+d_1+d_2=m}\binom m{d_0,d_1,d_2}\prod_{i=1}^ng(c_i,d_0,d_1,d_2)\\ &=\frac1{3^{n+m} }\sum_{d_0+d_1+d_2=m}\binom m{d_0,d_1,d_2}\prod_{i=1}^n\sum_{c_i\in[0,2]}g(c_i,d_0,d_1,d_2)\\ &=\frac1{3^{n+m} }\sum_{d_0+d_1+d_2=m}\binom m{d_0,d_1,d_2}\left(\sum_{c_i\in[0,2]}g(c_i,d_0,d_1,d_2)\right)^n\\ \end{aligned}\]

这样复杂度就降到了 $\mathcal O(n^2)$。

另外，本体取模会涉及到 $3$ 的单位根 $\omega_3$，模数 $P=10^9+9$，这要如何实现呢？

不妨先找个原根 $g$，有 $g^{P-1}=1$，那么，$\omega_3=g^{(P-1)/3}$，巧合的是，本体恰好满足 $3\mid P-1$。

#include<bits/stdc++.h>
#define siz(x) int((x).size())
#define all(x) std::begin(x),std::end(x)
#define fi first
#define se second
using namespace std;
using unt=unsigned;
using loli=long long;
using lolu=unsigned long long;
using pii=pair<int,int>;
mt19937_64 rng(random_device{}());
constexpr int N=1007,P=1e9+9;
struct mint{
	int d;
	mint()=default;
	mint(int x):d(x){}
	friend std::istream&operator>>(std::istream&x,mint&y){return x>>y.d;}
	friend std::ostream&operator<<(std::ostream&x,mint y){return x<<y.d;}
	friend mint operator+(mint x,mint y){return (x.d+=y.d)<P?x.d:x.d-P;}
	mint&operator+=(mint z){return (d+=z.d)<P?d:d-=P,*this;}
	friend mint operator-(mint x,mint y){return (x.d-=y.d)<0?x.d+P:x.d;}
	mint&operator-=(mint z){return (d-=z.d)<0?d+=P:d,*this;}
	friend mint operator*(mint x,mint y){return int(1ll*x.d*y.d%P);}
	mint&operator*=(mint z){return d=int(1ll*d*z.d%P),*this;}
	static mint qpow(int x,int y=P-2){int z=1;for(;y;y>>=1,x=int(1ll*x*x%P))if(y&1)z=int(1ll*x*z%P);return z;}
	friend mint operator/(mint x,mint y){return x*=qpow(y.d);}
	mint&operator/=(mint z){return (*this)*=qpow(z.d);}
	friend mint operator^(mint x,mint y){return qpow(x.d,y.d);}
	mint&operator^=(mint z){return *this=qpow(d,z.d);}
	mint operator()(mint z)const{return qpow(d,z.d);}
	mint&operator[](mint z){return *this=qpow(d,z.d);}
	mint inv()const{return qpow(d);}
	mint pow(mint z)const{return qpow(d,z.d);}
	int operator+()const{return d;}
	mint operator-()const{return d?P-d:0;}
	int operator~()const{return ~d;}
};
mint operator""_m(lolu x){return mint(int(x%P));}
int n,m,l;
mint fac[N],inv[N],w3=115381398;
signed main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
ios::sync_with_stdio(false);cin.tie(nullptr);
	cin>>n>>m>>l;
	fac[0]=inv[0]=1;
	for(int i=1;i<N;i++)fac[i]=fac[i-1]*i;
	inv[N-1]=fac[N-1].inv();
	for(int i=N-2;i;i--)inv[i]=inv[i+1]*(i+1);
	mint f0=l+1,f1=(1-w3(l+1))/(1-w3),f2=(1-(w3*w3)(l+1))/(1-w3*w3);
	mint ans=0;
	for(int i=0;i<=m;i++)for(int j=0;i+j<=m;j++){
		int k=m-i-j;
		mint g0=f0(i)*f1(j)*f2(k);
		mint g1=f1(i)*f2(j)*f0(k);
		mint g2=f2(i)*f0(j)*f1(k);
		ans+=(g0+g1+g2)(n)*fac[m]*inv[i]*inv[j]*inv[k];
	}
	cout<<ans/mint(3)(n+m);
	return 0;
}