题面
对于一个整数序列 { a 1 , a 2 , . . . , a n } \{a_1,a_2,...,a_n\} {a1,a2,...,an},我们定义它的变换为 { b 1 , b 2 , . . . , b n } \{b_1,b_2,...,b_n\} {b1,b2,...,bn},其中 b i = a 1 ∣ a 2 ∣ . . . ∣ a i b_i=a_1|a_2|...|a_i bi=a1∣a2∣...∣ai,其中 ∣ | ∣ 表示二进制按位或运算。
现在求有多少个长为 n n n 的序列 a a a,满足 1 ≤ a i ≤ 2 k − 1 1\leq a_i\leq 2^k-1 1≤ai≤2k−1,使得它的变换 b b b 是严格单调递增的,对 1 0 9 + 7 10^9+7 109+7 取模。
1 ≤ n ≤ 1 0 18 1\leq n\leq10^{18} 1≤n≤1018, 1 ≤ k ≤ 3 × 1 0 4 1\leq k\leq3\times10^4 1≤k≤3×104。
题解
这题其实没有什么好说的,无非是几个步骤:推出式子 → 打出任意魔术NTT → 调试。我们一步一步来。
首先,这个变换其实是前缀按位或,前缀按位或如果要严格单调递增,那么每个数字都要至少存在一位 1 是前面的数字都没有的。所以 n > k n>k n>k 时一定无解,但这不重要。按照这个思路,容易发现,后一个数字的合法方案只取决于前面所有数字按位或后的二进制 1 个数。后一个数字的这些位——这些之前存在过 1 的位——都可以任意取 1 或 0,然后在这些位以外的地方再添加至少一个 1 位。这里我们先不考虑出现过 1 的位和没出现过 1 的位之间的相对位置,因为它可以最后考虑。
我们设
f
(
n
,
x
)
f(n,x)
f(n,x) 表示前
n
n
n 个数字,按位或起来存在
x
x
x 个 1 的方案数(忽视掉没出现过 1 的位)。那么考虑当前数字新添加的 1 的个数
x
−
i
x-i
x−i 、这
x
−
i
x-i
x−i 个 1 与之前的 1 的相对位置方案、以及之前的 1 位下可任意取 0 或 1 的方案,就可以得到下面的式子:
f
(
1
,
x
)
=
[
x
>
0
]
,
f
(
n
,
x
)
=
∑
i
=
1
x
−
1
f
(
n
−
1
,
i
)
⋅
(
x
i
)
⋅
2
i
f(1,x)=[x>0]~,\\ f(n,x)=\sum_{i=1}^{x-1} f(n-1,i) \cdot {x\choose i}\cdot2^{i}
f(1,x)=[x>0] ,f(n,x)=i=1∑x−1f(n−1,i)⋅(ix)⋅2i
我们把组合数拆开,化成卷积形式:
f
(
n
,
x
)
x
!
=
∑
i
=
1
x
−
1
f
(
n
−
1
,
i
)
i
!
⋅
2
i
⋅
1
(
x
−
i
)
!
\frac{f(n,x)}{x!}=\sum_{i=1}^{x-1}\frac{f(n-1,i)}{i!}\cdot 2^i\cdot\frac{1}{(x-i)!}
x!f(n,x)=i=1∑x−1i!f(n−1,i)⋅2i⋅(x−i)!1
我们发现,最终的多项式它完全没法化成
g
(
x
)
⋅
h
(
x
)
n
g(x)\cdot h(x)^n
g(x)⋅h(x)n 的形式,没法朴素地进行多项式幂。但是我们可以进行更灵活的倍增,把上式写成下面这个样子或许能启发思考:
f
(
n
,
x
)
x
!
=
∑
i
=
1
x
−
1
f
(
n
−
1
,
i
)
i
!
⋅
2
1
i
⋅
f
(
1
,
x
−
i
)
(
x
−
i
)
!
\frac{f(n,x)}{x!}=\sum_{i=1}^{x-1}\frac{f(n-1,i)}{i!}\cdot 2^{1i}\cdot\frac{f(1,x-i)}{(x-i)!}
x!f(n,x)=i=1∑x−1i!f(n−1,i)⋅21i⋅(x−i)!f(1,x−i)
我们可以试着不只转移一个数字,而是连续转移一段。我们插入一段数字,令这一段数字出现过 1 的所有位与前面的数字都不同,然后再在后面所有数字中,于前面出现 1 的位处任意取 0 和 1 ,可以最终得到一个无比类似的式子:
f
(
n
+
m
,
x
)
x
!
=
∑
i
=
1
x
−
1
f
(
n
,
i
)
i
!
⋅
2
i
m
⋅
f
(
m
,
x
−
i
)
(
x
−
i
)
!
\frac{f(n+m,x)}{x!}=\sum_{i=1}^{x-1}\frac{f(n,i)}{i!}\cdot 2^{im}\cdot\frac{f(m,x-i)}{(x-i)!}
x!f(n+m,x)=i=1∑x−1i!f(n,i)⋅2im⋅(x−i)!f(m,x−i)
这样我们就可以做类快速幂了。
任意魔术NTT呢,一般是取 3 个不同的可以做NTT的模数分别做NTT,然后用中国剩余定理做回去。适用于每一次单独的多项式乘法,乘了过后要取模题目要求的模数,再继续下一步乘。
在用中国剩余定理的时候,由于数字极大,所以我们得用一个方法,在过程中取模题目要求的模数,并保证最终结果正确。
假设
N
T
T
NTT
NTT 用的三个模数分别为
M
1
,
M
2
,
M
3
M_1,M_2,M_3
M1,M2,M3 ,对于一个数字
x
x
x 有三个式子已知:
{
x
≡
a
1
m
o
d
M
1
x
≡
a
2
m
o
d
M
2
x
≡
a
3
m
o
d
M
3
\begin{cases} x\equiv a_1\mod M_1\\ x\equiv a_2\mod M_2\\ x\equiv a_3\mod M_3 \end{cases}
⎩⎪⎨⎪⎧x≡a1modM1x≡a2modM2x≡a3modM3
那么结果的表达式为
a
1
⋅
M
2
M
3
⋅
i
n
v
1
(
M
1
M
2
)
+
a
2
⋅
M
1
M
3
⋅
i
n
v
2
(
M
1
M
3
)
+
a
3
⋅
M
1
M
2
⋅
i
n
v
3
(
M
1
M
2
)
m
o
d
M
1
M
2
M
3
a_1\cdot M_2M_3\cdot inv_1(M_1M_2)+a_2\cdot M_1M_3\cdot inv_2(M_1M_3) +a_3\cdot M_1M_2\cdot inv_3(M_1M_2)\mod M_1M_2M_3
a1⋅M2M3⋅inv1(M1M2)+a2⋅M1M3⋅inv2(M1M3)+a3⋅M1M2⋅inv3(M1M2)modM1M2M3
这三个逆元可以健康地求出来,然后……用 _ _ i n t 128 \rm\_\_int128 __int128 🙂🤮
开玩笑的,具体过程很繁杂,可以参见这位高职人员的详细解说:🔗 。
或者,你可以做真的猛士,去用 MTT 。
最后,调试的时候有必须注意的地方。如果你的NTT模数是 998244353 或者更小,那么对答案多项式进行一次正变换和逆变换后,与原来的多项式将不一定相同。原因就是,你用的是一个比 1 0 9 + 7 10^9+7 109+7 更小的模数。所以,每次进行 3 模NTT 时,要把转移多项式给存下来。
或者,你可以不用 469762049、998244353 和 1004535809 ,而是做个真的猛士,像我一样用上三个比 1 0 9 + 7 10^9+7 109+7 大的模数:
- 1004535809
- 1007681537(最多 2 的 20 次方,原根 3)
- 1012924417(最多 2 的 21 次方,原根 5)
时间复杂度 O ( k log 2 k ) O(k\log^2k) O(klog2k) ,常数可观 ↓。
CODE
我用的 _ _ i n t 128 \rm\_\_int128 __int128 ,危险行为请勿模仿。
//真的猛士,敢于将打到一半的代码提交
#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<stack>
#include<random>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 30005
#define LL long long
#define ULL unsigned long long
#define ENDL putchar('\n')
#define DB long double
#define lowbit(x) (-(x) & (x))
#define FI first
#define SE second
#define eps (1e-1)
#define BIG __int128
int xchar() {
static const int maxn = 1000000;
static char b[maxn];
static int pos = 0,len = 0;
if(pos == len) pos = 0,len = fread(b,1,maxn,stdin);
if(pos == len) return -1;
return b[pos ++];
}
//#define getchar() xchar()
LL read() {
LL f = 1,x = 0;int s = getchar();
while(s < '0' || s > '9') {if(s<0)return -1;if(s=='-')f=-f;s = getchar();}
while(s >= '0' && s <= '9') {x = (x<<1) + (x<<3) + (s^48);s = getchar();}
return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar((x%10)^48);}
void putnum(LL x) {
if(!x) {putchar('0');return ;}
if(x<0) putchar('-'),x = -x;
return putpos(x);
}
void AIput(LL x,int c) {putnum(x);putchar(c);}
const int MOD = 1000000007;
const int M1 = 1012924417,R1 = 5;
const int M2 = 1007681537,R2 = 3;
const int M3 = 1004535809,R3 = 3;
int n,m,s,o,k;
int fac[MAXN],inv[MAXN],invf[MAXN];
int p2[MAXN];
int qkpow(int a,int b) {
int res = 1;
while(b > 0) {
if(b & 1) res = res *1ll* a % MOD;
a = a *1ll* a % MOD; b >>= 1;
}return res;
}
int qkpow(int a,int b,int MOD) {
int res = 1;
while(b > 0) {
if(b & 1) res = res *1ll* a % MOD;
a = a *1ll* a % MOD; b >>= 1;
}return res;
}
int om,xm[MAXN<<2];
int rev[MAXN<<2];
void NTT(int *s,int n,int op,int MOD,int R) {
for(int i = 1;i < n;i ++) {
rev[i] = (rev[i>>1]>>1) | ((i&1) ? (n>>1):0);
if(rev[i] < i) swap(s[i],s[rev[i]]);
}
om = qkpow(R,(MOD-1)/n,MOD); xm[0] = 1;
if(op < 0) om = qkpow(om,MOD-2,MOD);
for(int i = 1;i < n;i ++) xm[i] = xm[i-1] *1ll* om % MOD;
for(int k = 2,t = n>>1;k <= n;k <<= 1,t >>= 1) {
for(int j = 0;j < n;j += k) {
for(int i = j,l = 0;i < j+(k>>1);i ++,l += t) {
int A = s[i],B = s[i+(k>>1)];
s[i] = (A + xm[l]*1ll*B) % MOD;
s[i+(k>>1)] = (A +MOD- xm[l]*1ll*B%MOD) % MOD;
}
}
}
if(op < 0) {
int invn = qkpow(n,MOD-2,MOD);
for(int i = 0;i < n;i ++) s[i] = s[i]*1ll*invn % MOD;
}return ;
}
int a[MAXN<<2],b[MAXN<<2],c[MAXN<<2];
void mult(int *A,int *B,int le) {
// printf("A:");for(int i = 0;i < le;i ++) printf("%lld ",A[i]*1ll*fac[i]%MOD); ENDL;
// printf("B:");for(int i = 0;i < le;i ++) printf("%lld ",B[i]*1ll*fac[i]%MOD); ENDL;
NTT(A,le,1,M1,R1); NTT(B,le,1,M1,R1);
for(int i = 0;i < le;i ++) a[i] = A[i]*1ll*B[i] % M1;
NTT(A,le,-1,M1,R1);NTT(B,le,-1,M1,R1);
NTT(a,le,-1,M1,R1);
NTT(A,le,1,M2,R2); NTT(B,le,1,M2,R2);
for(int i = 0;i < le;i ++) b[i] = A[i]*1ll*B[i] % M2;
NTT(A,le,-1,M2,R2);NTT(B,le,-1,M2,R2);
NTT(b,le,-1,M2,R2);
NTT(A,le,1,M3,R3); NTT(B,le,1,M3,R3);
for(int i = 0;i < le;i ++) c[i] = A[i]*1ll*B[i] % M3;
NTT(A,le,-1,M3,R3);NTT(B,le,-1,M3,R3);
NTT(c,le,-1,M3,R3);
BIG v1 = qkpow(M2*1ll*M3%M1,M1-2,M1);
BIG v2 = qkpow(M1*1ll*M3%M2,M2-2,M2);
BIG v3 = qkpow(M1*1ll*M2%M3,M3-2,M3);
BIG MD = M1*1ll*M2; MD *= M3;
for(int i = 0;i < le;i ++) {
BIG nm = (v1*(BIG)M2*M3%MD*a[i]%MD + v2*(BIG)M1*M3%MD*b[i]%MD + v3*(BIG)M1*M2%MD*c[i]%MD) % MD;
A[i] = nm % MOD;
}
// printf("A:");for(int i = 0;i < le;i ++) printf("%lld ",A[i]*1ll*fac[i]%MOD); ENDL;
// printf("B:");for(int i = 0;i < le;i ++) printf("%lld ",B[i]*1ll*fac[i]%MOD); ENDL;
return ;
}
int A[MAXN<<2],B[MAXN<<2],C[MAXN<<2];
int main() {
LL N = read();m = read();
if(N > m) {
printf("0\n");
return 0;
}
n = N;
A[0] = 1;
int le = 1; while(le <= m*2) le <<= 1;
fac[0]=fac[1]=inv[0]=inv[1]=invf[0]=invf[1]=1;
for(int i = 2;i <= m;i ++) {
fac[i] = fac[i-1] *1ll* i % MOD;
inv[i] = (MOD-inv[MOD%i]) *1ll* (MOD/i) % MOD;
invf[i] = invf[i-1] *1ll* inv[i] % MOD;
}
p2[0] = 1;
for(int i = 1;i <= m;i ++) {
B[i] = invf[i];
p2[i] = (p2[i-1]<<1) % MOD;
}
int ct = 0,pw = 1;
while(n > 0) {
if(n & 1) {
for(int i = 1;i <= m;i ++) A[i] = A[i] *1ll* p2[i] % MOD;
mult(A,B,le);
for(int i = m+1;i < le;i ++) A[i] = 0;
ct += pw;
// printf("%d:\n",ct);
// for(int i = 1;i <= m;i ++) {
// printf("%lld ",A[i]*1ll*fac[i]%MOD);
// }ENDL;
}
for(int i = 0;i < le;i ++) {
if(i>m) C[i] = 0;
else C[i] = B[i] *1ll* p2[i] % MOD;
}
mult(B,C,le);
for(int i = 0;i < le;i ++) {
if(i>m) B[i] = 0;
else p2[i] = p2[i]*1ll*p2[i]%MOD;
}pw <<= 1;
n >>= 1;
}
int ans = 0;
for(int i = N;i <= m;i ++) {
(ans += A[i]*1ll*fac[m]%MOD*invf[m-i]%MOD) %= MOD;
}
AIput(ans,'\n');
return 0;
}