Description：

给定 $n,m,k$ ，求： $$ \sum_{i=0}^n\sum_{j=0}^m[2|i][2|j]\binom ij $$ $1\leqslant n\leqslant 10^9,1\leqslant m\leqslant 10^6$

## Solution：

假设 $2|n,2|m$ 。 $$ \begin{align} &\sum_{i=0}^n\sum_{j=0}^m[2|i][2|j]\binom ij\\ =&\sum_{i=0}^n\sum_{j=0}^m[2|i][2|j][x^j](x+1)^i\\ =&\sum_{j=0}^m[2|j][x^j](\sum_{i=0}^n[2|i](x+1)^i)\\ =&\sum_{i=0}^m[2|i][x^i]\frac{(x+1)^{n+2}-1}{(x+1)^2-1}\\ =&\sum_{i=0}^m[2|i][x^i]\frac{\begin{align}\sum_{j=1}^{n+2}\binom{n+2}{j}x^j\end{align}}{x(x+2)}\\ =&\sum_{i=0}^m[2|i][x^i]\frac{\begin{align}\sum_{j=0}^{n+1}\binom{n+2}{j+1}x^{j}\end{align}}{x+2} \end{align} $$ 这个时候我们发现需要做除法，但是这个除法很不好做，于是我们就用最简单的方法手动模拟除法，也就是假设： $$ \sum_{j=1}^{n+2}\binom{n+2}{j}x^{j-1}=A(x),x+2=B(x),C(x)=A(x)/B(x) $$ 那我们可以手工列一下除法的式子： $$ C(x)(x+2)=A(x)\\ [x^{i+1}]A(x)=[x^i]C(x)+2[x^{i+1}]C(x)\\ [x^0]A(x)=2[x^0]C(x) $$ 这个时候我们会发现可能会有除二操作，而题目又没有保证模数是质数，因此我们利用中国剩余定理，把模数分解成 $s\times 2^k$ 的形式，先考虑系数对 $s$ 取模： $$ [x^0]C(x)=\frac{[x^0]A(x)}2\\ [x^i]C(x)=\frac{[x^i]A(x)-[x^{i-1}]C(x)}2 $$ 于是顺着推下来就好了。

再考虑系数对 $2^k$ 取模，这个时候我们不能做除法： $$ \begin{align} [x^i]C(x)&=[x^{i+1}]A(x)-2[x^{i+1}]C(x)\\ &=[x^{i+1}]A(x)-2[x^{i+2}]A(x)+4[x^{i+2}]C(x)\\ &=\sum_{j=1}(-2)^{j-1}[x^{i+j}]A(x) \end{align} $$ 当 $j-1\geqslant k$ 的时候，再往后的就被模掉了，因此： $$ [x^i]C(x)=\sum_{j=1}^k(-2)^{j-1}[x^{i+j}]A(x) $$ 由于 $k$ 不会特别大因此暴力算就行了。

但是这样还是有一个问题，就是我们不能直接利用组合数 $\binom nm=\frac{n!}{m!(n-m)!}$ 来计算，因为要做除法但是可能没有逆元，因此用 $C(n,i)=\frac{C(n,i-1)\times(n-i+1)}{i}$ 来计算，并且就像扩展卢卡斯定理那样把模数中有的质因数拿出来记下他们的次数 $cnt[i]$ ，然后并不做除法，而是每次暴力计算一下每个质因数的幂次然后预处理他们的幂

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

typedef long long ll;

ll n,m,k;

#define MAXN 1000040

void exgcd(ll a,ll b,ll &x,ll &y)

{

	if(b == 0)

	{

		x = 1;y = 0;

		return;

	}

	exgcd(b,a % b,x,y);

	ll t = x;x = y;

	y = t - a / b * y;

	return; 

}

ll inv(ll a,ll mod)

{

	ll x,y;

	exgcd(a,mod,x,y);

	x = (x % mod + mod) % mod;

	return x;

}

ll cnt2,rem;

ll CRT(pair<ll,ll> a,pair<ll,ll> b)

{

	ll M = a.second * b.second;

	ll res1 = a.first * inv(b.second,a.second) % M * b.second % M;

	ll res2 = b.first * inv(a.second,b.second) % M * a.second % M;

	return (res1 + res2) % M;

}

ll A[MAXN];

ll pr[100],cnt[100];

ll power[100][35];

ll C[MAXN];

pair<ll,ll> solve1()// mod rem

{

	if(rem == 1)return make_pair(0,1);

	ll c = 1;

	for(int i = 1;i <= pr[0];++i)

	{

		power[i][0] = 1;

		for(int j = 1;j <= 34;++j)power[i][j] = power[i][j - 1] * pr[i] % rem;

	}

	A[0] = 1;

	for(int k = 1;k <= min(n + 2,m + 1);++k)

	{

		ll fz = n + 2 - k + 1,fm = k;

		for(int i = 1;i <= pr[0];++i)

		{

			while(fz % pr[i] == 0)++cnt[i],fz /= pr[i];

			while(fm % pr[i] == 0)--cnt[i],fm /= pr[i];

		}

		c = c * fz % rem * inv(fm,rem) % rem;

		ll res = c;

		for(int i = 1;i <= pr[0];++i)res = res * power[i][cnt[i]] % rem;

		A[k] = res;

	}

	ll res = 0;

	C[0] = A[1] * inv(2,rem) % rem;

	res = C[0];

	for(int i = 1;i <= m;++i)

	{

		C[i] = (A[i + 1] - C[i - 1] + rem) % rem * inv(2,rem) % rem;

		if(i % 2 == 0)res = (res + C[i]) % rem;

	}//cout << res << " " << rem << endl;

	return make_pair(res,rem);

}

pair<ll,ll> solve2()// mod 2^cnt2

{

	if(cnt2 == 0)return make_pair(0,1);

	ll c = 1;

	A[0] = 1;

	ll p = (1 << cnt2);

	int cnt = 0;

	for(int k = 1;k <= min(n + 2,m + cnt2 + 1);++k)

	{

		ll fz = n + 2 - k + 1,fm = k;

		while(fz % 2 == 0)++cnt,fz /= 2;

		while(fm % 2 == 0)--cnt,fm /= 2;

		c = c * fz % p * inv(fm,p) % p;

		ll res = c;

		res = res * (1 << cnt) % p;

		A[k] = res;

	}

	ll res = 0;

	for(int i = 0;i <= m;++i)

	{

		C[i] = 0;

		for(int j = 1,cur = 1;j <= cnt2;++j,cur = 1ll * cur * (p - 2) % p)

		{

			C[i] = (C[i] + 1ll * cur * A[i + j + 1] % p) % p;

		}

		if(i % 2 == 0)res = (res + C[i]) % p;

	}//cout << res << " " << p << endl;

	return make_pair(res,p);

}

int main()

{

	scanf("%lld%lld%lld",&n,&m,&k);

	if(k == 1){puts("0");return 0;}

	if(n % 2 == 1)--n;

	if(m % 2 == 1)--m;

	while(k % 2 == 0)++cnt2,k /= 2;

	rem = k;

	for(int i = 2;i * i <= rem;++i)

	{

		if(k % i == 0)

		{

			pr[++pr[0]] = i;

			while(k % i == 0)k /= i;

		}

	}

	if(k != 1)pr[++pr[0]] = k;

	cout << CRT(solve1(),solve2()) << endl;

	return 0;

}