@51nod - 1196/1197/1198@ 字符串的数量

@description@
@solution@
@accepted code@
@details@

@description@

用 N 个不同的字符（编号1 - N），组成一个字符串，有如下要求：

(1) 对于编号为 i 的字符，如果2 * i > n，则该字符可以作为结尾字符。如果不作为结尾字符而是中间的字符，则该字符后面可以接任意字符。

(2) 对于编号为 i 的字符，如果2 * i <= n，则该字符不可以作为结尾字符。作为中间字符，那么后面接的字符编号一定要 >= 2 * i。

问有多少长度为M且符合条件的字符串，由于数据很大，只需要输出该数Mod 10^9 + 7的结果。

例如：N = 2，M = 3。则abb, bab, bbb是符合条件的字符串，剩下的均为不符合条件的字符串。

input

第1行：一个数T，表示后面用作输入测试的数的数量。（1 <= T <= 20)

第2 - T + 1行：每行2个数，N, M中间用空格分割，N为不同字符的数量，M为字符串的长度。(2 <= N, M <= 10^18)

output

共T行，对应符合条件的字符串的数量。由于数据很大，只需要输出该数Mod 10^9 + 7的结果。

sample input

1

6 3

sample output

73

@solution@

什么鬼题目。。。

@part - 1@

不难写出一个很 naive 的 dp：定义 dp[i][j] 长度为 i，末尾字符为 j 的方案数。

对 2*j > n 还是 <= n 进行分类讨论就可以得到转移式。

不难发现对于 2*j > n 的 dp 值是相同的，我们令 f[i] 表示这个值。

于是我们可以发现 dp[i][j] 总是可以表示为 f[i-1], f[i-2], ... 的线性组合。

考虑这个性质的实际意义：我们可以在一个字符后面添加若干的字符最终变为以 2*x > n 这一类字符结尾。

我们以 2*x > n 这一类字符作为隔断，将最终的字符串分成若干以这一类字符结尾的段。可以发现每一段的长度是 O(log) 的，因为后一个字符至少是前一个字符的两倍。

定义 g[i][x] 表示以字符 x 开头，产生长度 i 的段的方案数。于是根据定义就有 g[1][x] = 1，对于 2*x > n。

同时也有 f[i] = ∑g[j][x]*f[i-j]。因为段的长度不会超过 O(log)，所以可以预处理 g 然后矩乘。

@part - 2@

先来考虑矩乘这一部分，顺便把我们上述的内容整理一下。

定义 f[i] 表示问题所想我们求解的答案，再定义 g[i][x] 表示以字符 x 开头，最后一位且只有最后一位是 2*x > n 这类字符，长度为 i 的字符串方案数。

于是有转移 \(f[i] =\sum_{j, k}g[j][k]*f[i-j]\)。

注意到 g[i][x] 这种字符串，每多一位，末尾字符就会变为前一个的字符两倍。故它的长度是 O(log) 的。

所以我们在预处理好 g 之后，可以写一个 O(log^3 n*log m) 的矩乘求出 f[n]。

但是这样有些卡，我们可能需要进一步的优化。

使用一种名为“常系数齐次线性递推”的黑科技即可。（可以点击我的这篇博客查看详情）

不过这道题可以直接暴力多项式乘法与多项式取模，因为 log 也不大，说不定 fft 的常数反而会跑得慢些。

最后时间可以优化到 O(log^2 n*log m)

@part - 3@

考虑怎么预处理好 g。不难发现 g 有如下的转移式：

\[g[i][x]=\sum_{2*x\le j}g[i-1][j]
\]

然后开始玄幻起来了。我们使用万能的归纳法，尝试证明 g 始终满足以下形式：

\[g[i][x]=\begin{cases}
P1(x) & 1 \le x \le \lfloor{\frac{a}{2}}\rfloor \\
P2(x) & \lfloor{\frac{a}{2}}\rfloor + 1 \le x \le a
\end{cases}\]

即一个分段的多项式函数。

假设对于所有 <= i 的都成立，尝试证明对 i+1 成立。

如果 \(\lfloor{\frac{a}{2}}\rfloor + 1 \le 2*x \le a\)，则：

\[g[i+1][x] = \sum_{2*x\le j}^{j\le a}g[i][j] = \sum_{2*x\le j}^{j\le a}P2(j) = \sum_{j=0}^{a}P2(j) - \sum_{j=0}^{2*x-1}P2(j)
\]

前一项是一个常数，而后一项呢？我们不妨令 \(P2(x) = a_0+a_1x+...\)。

于是：

\[\sum_{j=0}^{2*x-1}P2(j)=a_0(\sum_{j=0}^{2*x-1}j^0)+a_1(\sum_{j=0}^{2*x-1}j^1)+...
\]

然后你发现那是个自然数幂和。然后根据你所了解到的知识，自然数幂和可以表示为多项式形式。然后得证。

如果 \(1 \le 2*x \le \lfloor{\frac{a}{2}}\rfloor\)，证明是类似的，不再赘述。

关于自然数幂和，这里有两篇（我认为可供参考的）博客：blog1, blog2。

找时间可以专门研究这玩意儿（咕咕咕）。

然后就没有然后了。预处理一下自然数幂和的系数之类的即可。

具体可以看代码。

@accepted code@

#include<cstdio>

typedef long long ll;

const int MOD = int(1E9) + 7;

const int MSIZE = 64;

int pow_mod(int b, int p) {

	int ret = 1;

	while( p ) {

		if( p & 1 ) ret = 1LL*ret*b%MOD;

		b = 1LL*b*b%MOD;

		p >>= 1;

	}

	return ret;

}

int len;

struct poly{

	int k[2*MSIZE + 5];

	poly() {for(int i=0;i<=2*MSIZE;i++) k[i] = 0;}

	int get_val(int x) {

		int ret = 0;

		for(int i=len;i>=0;i--)

			ret = (1LL*ret*x%MOD + k[i])%MOD;

		return ret;

	}

	friend poly operator +(poly A, poly B) {

		poly C;

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

			C.k[i] = (A.k[i] + B.k[i])%MOD;

		return C;

	}

	friend poly operator -(poly A, poly B) {

		poly C;

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

			C.k[i] = (A.k[i] + MOD - B.k[i])%MOD;

		return C;

	}

	friend poly operator *(int k, poly B) {

		poly C;

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

			C.k[i] = 1LL*k*B.k[i]%MOD;

		return C;

	}

	friend poly mul_mod(poly A, poly B, poly C) {

		poly D;

		for(int i=0;i<=2*len;i++)

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

				D.k[i] = (D.k[i] + 1LL*A.k[j]*B.k[i-j]%MOD)%MOD;

		for(int i=2*len;i>=len;i--) {

			for(int j=0;j<=len;j++)

				D.k[i-len+j] = (D.k[i-len+j] + MOD - 1LL*D.k[i]*C.k[j]%MOD)%MOD;

		}

		return D;

	}

	void debug() {

		for(int i=0;i<2*MSIZE;i++)

			printf("%d ", k[i]);

		puts("");

	}

}pw[MSIZE + 5], a[MSIZE + 5];

int comb[MSIZE + 5][MSIZE + 5];

void init() {

	for(int i=0;i<MSIZE;i++) {

		comb[i][0] = 1;

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

			comb[i][j] = (comb[i-1][j] + comb[i-1][j-1])%MOD;

	}

	for(int i=0;i<MSIZE;i++) {

		for(int j=0;j<=i+1;j++)

			pw[i].k[j] = comb[i+1][j];

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

			pw[i] = (pw[i] - comb[i+1][j]*pw[j]);

		pw[i] = 1LL*pow_mod(i+1, MOD-2)*pw[i];

	}

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

		for(int j=0;j<=i+1;j++)

			for(int k=0;k<=j;k++)

				a[i].k[k] = (a[i].k[k] + 1LL*pw[i].k[j]*(1LL*comb[j][k]*(1LL*pow_mod(2, k)*pow_mod(MOD-1, j-k)%MOD)%MOD)%MOD)%MOD;

}

int get_pw(ll x, int k) {return pw[k].get_val(x%MOD);}

int get_pw(ll l, ll r, int k) {return (get_pw(r, k) + MOD - get_pw(l - 1, k))%MOD;}

int get_val(poly A, ll l, ll r) {

	int ret = 0;

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

		ret = (ret + 1LL*get_pw(l, r, i)*A.k[i]%MOD)%MOD;

	return ret;

}

poly trans(poly A, ll r) {

	poly B; B.k[0] = get_val(A, 0, r);

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

		B = B - A.k[i]*a[i];

	return B;

}

poly get_M(ll n) {

	poly A, B, M; B.k[0] = M.k[len] = 1;

	for(int i=0;i<len;i++,n/=2){

		int sA = get_val(A, 1, n/2), sB = get_val(B, n/2 + 1, n);

		M.k[len-i-1] = (MOD - (sA + sB)%MOD)%MOD;

		A = trans(A, n/2), A.k[0] = (A.k[0] + sB)%MOD;

		B = trans(B, n);

	}

	return M;

}

poly p_pow(poly M, ll p) {

	poly ret, b; b.k[1] = 1, ret.k[0] = 1;

	while( p ) {

		if( p & 1 ) ret = mul_mod(ret, b, M);

		b = mul_mod(b, b, M);

		p >>= 1;

	}

	return ret;

}

int solve(ll n, ll m) {

	ll tmp = n; for(len = 0; tmp; len++, tmp /= 2);

	return p_pow(get_M(n), m+len-1).k[len-1];

}

int main() {

	init(); int T; scanf("%d", &T);

	for(int i=1;i<=T;i++) {

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

		printf("%d\n", solve(n, m));

	}

}

@details@

我可能永远也想不到把 g 表示成分段多项式函数的形式这一步。。。

写起来还是比较愉快的吧。。。毕竟没有毒瘤的多项式操作。。。所有多项式操作都可以直接暴力搞。。。

矩乘过不了是真的令人生气。。。

巴特西