[BZOJ3129][SDOI2013]方程（扩展Lucas定理+容斥）

先介绍：把m拆成n个正整数之和，方案数为Cn−1m−1。
证明见http://blog.csdn.net/xyz32768/article/details/78502117（[JSOI2011]分特产）中的相关部分。
首先可以发现，对于任意的n1+1≤i≤n1+n2，只需要把m减去A[i]−1，就可以不考虑A[n1+1...n2]的限制。
再考虑A[1...n1]的限制。看到n1只有8，所以容易想到可以用容斥搞。
对于组合数，可以用扩展Lucas定理求得。
扩展Lucas定理链接：http://blog.csdn.net/clove_unique/article/details/54571216
代码：

#include <cmath>#include <cstdio>#include <cstring>#include <iostream>#include <algorithm>using namespace std;inline int read() {    int res = 0; bool bo = 0; char c;    while (((c = getchar()) < '0' || c > '9') && c != '-');    if (c == '-') bo = 1; else res = c - 48;    while ((c = getchar()) >= '0' && c <= '9')        res = (res << 3) + (res << 1) + (c - 48);    return bo ? ~res + 1 : res;}const int N = 2e5 + 5, R = 4005, E = 14;typedef long long ll;int cnt, n1, n2, A1[E];ll bas[R], num[R], ex[R], PYZ, F[N];void exgcd(ll a, ll b, ll &x, ll &y) {    if (!b) return (void) (x = 1, y = 0);    exgcd(b, a % b, y, x); y -= a / b * x;}ll inv(ll a, ll LPF) {    ll x, y; exgcd(a, LPF, x, y);    x = (x % LPF + LPF) % LPF;    return x;}ll qpow(ll a, ll b, ll LPF) {    ll res = 1;    while (b) {        if (b & 1) (res *= a) %= LPF;        (a *= a) %= LPF;        b >>= 1;    }    return res;}ll crt() {    ll ans = 0; int i;    for (i = 1; i <= cnt; i++)        (ans += ex[i] * (PYZ / num[i]) * inv(PYZ / num[i], num[i])) %= PYZ;    return ans;}ll fac(ll n, ll x, ll LPF) {    if (n < x) return F[n]; ll ans = qpow(F[LPF - 1], n / LPF, LPF);    (ans *= F[n % LPF]) %= LPF; (ans *= fac(n / x, x, LPF)) %= LPF;    return ans;}ll comb(ll n, ll m, ll x, ll LPF) {    if (m > n) return 0; ll i; F[0] = 1; ll x1 = 0, x2 = 0, x3 = 0;    for (i = 1; i < LPF; i++) F[i] = i % x ? F[i - 1] * i % LPF : F[i - 1];    for (i = n; i;) x1 += (i /= x); for (i = m; i;) x2 += (i /= x);    for (i = n - m; i;) x3 += (i /= x); x1 -= x2 + x3;    ll tmp = qpow(x, x1, LPF);    x1 = fac(n, x, LPF); x2 = fac(m, x, LPF); x3 = fac(n - m, x, LPF);    return (tmp * x1 % LPF) * (inv(x2, LPF) * inv(x3, LPF) % LPF) % LPF;}void init() {    int i, x = PYZ, S = sqrt(PYZ); cnt = 0;    for (i = 2; i <= S; i++) if (x % i == 0) {        bas[++cnt] = i; num[cnt] = 1;        while (x % i == 0) num[cnt] *= i, x /= i;    }    if (x > 1) bas[++cnt] = num[cnt] = x;}ll C(ll n, ll m) {    int i; for (i = 1; i <= cnt; i++)        ex[i] = comb(n, m, bas[i], num[i]);    return crt();}void work() {    int i, j, n, m, x; n = read(); n1 = read(); n2 = read(); m = read();    for (i = 1; i <= n1; i++) A1[i] = read(); int Ans = 0;    for (i = 1; i <= n2; i++) if ((x = read()) >= 1) m -= x - 1;    if (n > m) return (void) puts("0");    for (i = 1; i <= n1; i++) if (A1[i] <= 0) return (void) puts("0");    for (i = 0; i < (1 << n1); i++) {        int cnt = m, tot = 0; for (j = 0; j < n1; j++)            if ((i >> j) & 1) cnt -= A1[j + 1], tot++;        int tmp = C(cnt - 1, n - 1);        if (tot & 1) Ans = (Ans - tmp + PYZ) % PYZ;        else (Ans += tmp) %= PYZ;    }    printf("%d\n", Ans);}int main() {    int T = read(); PYZ = read(); init();    while (T--) work();    return 0;}