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

先介绍：把$m$拆成$n$个正整数之和，方案数为$C_{m-1}^{n-1}$。
证明见http://blog.csdn.net/xyz32768/article/details/78502117（[JSOI2011]分特产）中的相关部分。
首先可以发现，对于任意的$n1+1\leq i\leq n1+n2$，只需要把$m$减去$A[i]-1$，就可以不考虑$A[n1+1...n2]$的限制。
再考虑$A[1...n1]$的限制。看到$n1$只有 8 <script type="math/tex" id="MathJax-Element-20">8</script>，所以容易想到可以用容斥搞。
对于组合数，可以用扩展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;
}