@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 的都成立，尝试证明对 i+1 成立。
如果 (lfloor{frac{a}{2}}
floor + 1 le 2*x le a)，则：
前一项是一个常数，而后一项呢？我们不妨令 (P2(x) = a_0+a_1x+...)。
于是：
然后你发现那是个自然数幂和。然后根据你所了解到的知识，自然数幂和可以表示为多项式形式。然后得证。
如果 (1 le 2*x le lfloor{frac{a}{2}}
floor)，证明是类似的，不再赘述。
关于自然数幂和，这里有两篇（我认为可供参考的）博客：blog1, blog2。
找时间可以专门研究这玩意儿（咕咕咕）。update in 2020/06/09：我没咕。真的。
然后就没有然后了。预处理一下自然数幂和的系数之类的即可。
具体可以看代码。
@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
", solve(n, m));
}
}
@details@
update in 2020/06/09：本题是可以用拉格朗日插值来做的。代码如下：
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MOD = int(1E9) + 7;
const int MSIZE = 64;
inline int add(int x, int y) {x += y; return x >= MOD ? x - MOD : x;}
inline int sub(int x, int y) {x -= y; return x < 0 ? x + MOD : x;}
inline int mul(int x, int y) {return (int)(1LL * x * y % MOD);}
int pow_mod(int b, int p) {
int ret = 1;
for(int i=p;i;i>>=1,b=mul(b,b))
if( i & 1 ) ret = mul(ret, b);
return ret;
}
int len;
struct poly{
int k[2*MSIZE + 5];
poly() {for(int i=0;i<=2*MSIZE;i++) k[i] = 0;}
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] = add(D.k[i], mul(A.k[j], B.k[i-j]));
for(int i=2*len;i>=len;i--)
for(int j=0;j<=len;j++)
D.k[i-len+j] = sub(D.k[i-len+j], mul(D.k[i], C.k[j]));
return D;
}
void debug() {
for(int i=0;i<2*MSIZE;i++)
printf("%d ", k[i]);
puts("");
}
}pw[MSIZE + 5], a[MSIZE + 5];
int ifct[MSIZE + 5];
void init() {
ifct[0] = 1; for(int i=1;i<=MSIZE;i++) ifct[i] = mul(ifct[i - 1], i);
for(int i=0;i<=MSIZE;i++) ifct[i] = pow_mod(ifct[i], MOD - 2);
}
int lf[MSIZE + 5], rf[MSIZE + 5];
int get(int n, int x, int *y) {
lf[0] = 1; for(int i=1;i<=n;i++) lf[i] = mul(lf[i - 1], sub(x, i));
rf[n + 1] = 1; for(int i=n;i>=1;i--) rf[i] = mul(rf[i + 1], sub(x, i));
int ans = 0;
for(int i=1;i<=n;i++) {
int del = mul(mul(lf[i - 1], rf[i + 1]), mul(ifct[i - 1], ifct[n - i]));
ans = ((n - i) & 1) ? sub(ans, mul(del, y[i])) : add(ans, mul(del, y[i]));
}
return ans;
}
int f[MSIZE + 5], df, g[MSIZE + 5], dg, t[MSIZE + 5], dt;
poly get_M(ll n) {
df = 0, g[dg = 1] = 1;
poly M; M.k[len] = 1;
for(int i=0;i<len;i++) {
ll l1 = 1, r1 = n / 2, l2 = n / 2 + 1, r2 = n;
if( dg + 1 <= r2 - l2 + 1 ) g[dg + 1] = get(dg, dg + 1, g), dg++;
for(int j=1;j<=dg;j++) g[j] = add(g[j], g[j - 1]);
if( df + 1 <= r1 - l1 + 1 ) f[df + 1] = get(df, df + 1, f), df++;
for(int j=1;j<=df;j++) f[j] = add(f[j], f[j - 1]);
int s = get(dg, (int)((r2 - l2 + 1) % MOD), g);
for(int j=1;j<=df;j++) f[j] = add(f[j], s);
if( i + 1 == len ) {
M.k[0] = sub(0, get(dg, 1, g));
break;
}
else M.k[len - i - 1] = sub(0, get(df, (int)(r1 % MOD), f));
n /= 2; ll l3 = 1, r3 = n / 2, l4 = n / 2 + 1, r4 = n;
dt = (int)min((ll)df, r3 - l3 + 1);
for(int j=1;j<=dt;j++) t[j] = get(df, (int)((r1 - 2*(r3 - j + 1) + 1) % MOD), f);
for(int j=1;j<=dt;j++) f[j] = t[j]; df = dt;
dt = (int)min((ll)dg, r4 - l4 + 1);
for(int j=1;j<=dt;j++) t[j] = get(dg, (int)((r2 - 2*(r4 - j + 1) + 1) % MOD), g);
for(int j=1;j<=dt;j++) g[j] = t[j]; dg = dt;
}
/*
for(int i=0;i<=len;i++)
printf("%d ", M.k[i]);
*/
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
", solve(n, m));
}
}
我可能永远也想不到把 g 表示成分段多项式函数的形式这一步。。。
写起来还是比较愉快的吧。。。毕竟没有毒瘤的多项式操作。。。所有多项式操作都可以直接暴力搞。。。
矩乘过不了是真的令人生气。。。