【GDOI2020模拟02.19】A （数位dp+范德蒙恒等式）

\(n<=12,d<=500,l,r<=d^{500}\)

题解：

解法1：

使r-=l，并把所有l相加。

这样只用考虑与r的关系。

设\(f[i][S][j=0-n]\)表示确定了高位的i位，与r的数位关系的状态是S，需要低位进j为上来。

转移可以枚举S，再枚举S'∈S，做n次前缀和确定这一位系数，再转移。

复杂度：\(O(d*3^n*n^2*d)\)

也可以用高维前缀和优化这个转移，设\(g[i][S][j][k]\)表示确定了S的前i位，现在是S，需要进j位，现在和是k，初值通过f得来。

复杂度：\(O(d*2^n*n^4*d)\)

解法2：

前面的都是考虑<=r，这也是数位dp的常规做法。

现在考虑容斥>l。

记\(s=\sum l\)，如果最后的和是\(v\)，那么系数就是挡板问题：\(C_{v-s-1}^{n-1}\)。

范德蒙恒等式：\(C_{n+m}^k=\sum_{i=0}^kC_{n}^i*C_{m}^{k-i}\)

当n、m、k是非负数的时候显然成立，其实n，m是负数时也成立，可以用生成函数解释。

那么把\(C_{v-s-1}^{n-1}\)拆开，\(=\sum_{i=0}^{n-1}C_{v}^i*C_{-s-1}^{n-1-i}\)

下标为负的组合数的求法：
\(C_{-n}^m=(-n)^{m↓}/m!\)

或者\((1+x)^{-n}[x^m]={1\over (1+x)^n}[x^m]\)

只需要求出\(C_v^i*[v d进制下每个数位都出现在T中]*[v>s]\)的和即可。

那么还是数位dp，设\(f[i][0/1][j]\)表示确定了高位的i位，与s的数位关系，\(C_{v}^i\)的和。

转移时相当于求\(C_{v+d^i*num}^i\)同样可以代入范德蒙恒等式拆开。

我把组合数展开成下降幂，然后维护0~n次和，是差不多的。

注意下前导0的处理，可以dp多开一维，这样太慢了，不如每次手工加回前导零的系数。

时间复杂度：\(O(2^n*d*n^2)\)

Code:

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;

const int mo = 1e9 + 7;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

const int N = 510;


int n, d, ky[N];

struct P {
	int a0, a[N];
	P() {
		a0 = 0; memset(a, 0, sizeof a);
	}
};

P operator + (P a, P b) {
	a.a0 = max(a.a0, b.a0);
	fo(i, 0, a.a0) {
		a.a[i] += b.a[i];
		if(a.a[i] >= d) a.a[i + 1] ++, a.a[i] -= d;
	}
	while(a.a[a.a0 + 1]) a.a0 ++;
	return a;
}

P operator - (P a, P b) {
	a.a0 = max(a.a0, b.a0);
	fo(i, 0, a.a0) {
		a.a[i] -= b.a[i];
		if(a.a[i] < 0) a.a[i] += d, a.a[i + 1] --;
	}
	while(a.a0 > 0 && !a.a[a.a0]) a.a0 --;
	return a;
}

namespace read {
char str[N * 3];
int b[N * 3];
void read(P &a) {
	memset(b, 0, sizeof b);
	int n;
	scanf("%s", str + 1);
	n = strlen(str + 1);
	fo(i, 1, n) b[i] = str[n - i + 1] - '0';
	a.a0 = -1;
	while(n > 1 || b[n] > 0) {
		int y = 0;
		fd(i, n, 1) {
			b[i] += y * 10;
			y = b[i] % d; b[i] /= d;
		}
		while(n > 0 && !b[n]) n --;
		a.a[++ a.a0] = y;
	}
}
}

void write(P a) {
	fo(i, 0, a.a0) pp("%d ", a.a[i]); hh;
}

P l[15], r[15], sum, val_1;

int mx;

ll c[15][15], fac[15], nf[15];

void build(int n) {
	fo(i, 0, n) {
		c[i][0] = 1;
		fo(j, 1, i) c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mo;
	}
	fac[0] = nf[0] = 1;
	fo(i, 1, n) fac[i] = fac[i - 1] * i % mo, nf[i] = ksm(fac[i], mo - 2);
}

struct Q {
	ll x, y;
};
ll ad[N]; Q ax[N][N][15];

ll f[2][2][13]; int o;

ll cnt[15], cnt2[15], cnt3[15];

ll p[15];

ll ans;

void dp(P sum, int xs) {
	memset(f, 0, sizeof f);
	f[o][1][0] = 1;
	int q0 = 1;
	fd(w, mx + 5, 0) {
		int num = sum.a[w];
		if(num) q0 = 0;
		o = !o; memset(f[o], 0, sizeof f[o]);
		ff(i, 0, n) ff(j, i, n) {
			f[o][0][j] += f[!o][0][i] * c[j][i] % mo * ax[w][d - 1][j - i].y % mo;
			f[o][0][j] += f[!o][1][i] * c[j][i] % mo * (ax[w][d - 1][j - i].y - ax[w][num][j - i].y) % mo;
			f[o][1][j] += f[!o][1][i] * c[j][i] % mo * ax[w][num][j - i].x % mo;
		}
		if(q0 && !ky[0]) f[o][1][0] ++;
		fo(i, 0, n) f[o][0][i] %= mo, f[o][1][i] %= mo;
	}
	ff(i, 0, n) cnt[i] = (f[o][0][i] + f[o][1][i]) % mo;
	memset(p, 0, sizeof p);
	p[0] = 1; cnt2[0] = cnt[0];
	ff(i, 1, n) {
		fd(j, i, 0) p[j] = (p[j] * -(i - 1) + (j ? p[j - 1] : 0)) % mo;
		cnt2[i] = 0; fo(j, 0, i) cnt2[i] = (cnt2[i] + cnt[j] * p[j]) % mo;
		cnt2[i] = cnt2[i] * nf[i] % mo;
	}
	ll s = 0;
	fo(i, 0, sum.a0) s = (s + ad[i] * sum.a[i]) % mo;
	cnt3[0] = 1;
	ff(i, 1, n) cnt3[i] = cnt3[i - 1] * (-s - (i - 1)) % mo;
	ff(i, 1, n) cnt3[i] = cnt3[i] * nf[i] % mo;
	fo(i, 0, n - 1) ans = (ans + cnt2[i] * cnt3[n - 1 - i] % mo * xs) % mo;
}

void dg(int x, P s, int xs) {
	if(x > n) {
		dp(s, xs);
		return;
	}
	dg(x + 1, s + l[x], xs);
	dg(x + 1, s + r[x], -xs);
}

int main() {
	freopen("A.in", "r", stdin);
	freopen("A.out", "w", stdout);
	scanf("%d %d", &n, &d);
	build(n);
	ff(i, 0, d) scanf("%d", &ky[i]);
	val_1.a[0] = 1;
	fo(i, 1, n) {
		read :: read(l[i]);
		read :: read(r[i]);
		mx = max(mx, r[i].a0);
		l[i] = l[i] - val_1;
	}
	ad[0] = 1; fo(i, 1, mx + 5) ad[i] = ad[i - 1] * d % mo;
	fo(i, 0, mx + 5) {
		fo(j, 0, d - 1) fo(k, 0, n) {
			ax[i][j][k].x = ksm(j * ad[i] % mo, k);
			if(!ky[j]) ax[i][j][k].x = 0;
			ax[i][j][k].y = ((j ? ax[i][j - 1][k].y : 0) + ax[i][j][k].x) % mo;
		}
	}
	dg(1, val_1, 1);
	ans = (ans % mo + mo) % mo;
	pp("%lld\n", ans);
}