题目
题目背景
“人类一思考,上帝就发笑。” —— 犹太谚语
“ O n e I n D a r k \sf OneInDark OneInDark 研究数学, O U Y E \sf OUYE OUYE 就发笑。” —— 今日新闻
题目描述
记
f
1
(
n
)
f^1(n)
f1(n) 为将
n
n
n 用
m
m
m 的非负整数次幂的和来表示的方案数,如
m
=
2
m=2
m=2 时
n
=
4
n=4
n=4 的拆分数有
{
1
,
1
,
1
,
1
}
,
{
1
,
1
,
2
}
,
{
2
,
2
}
,
{
4
}
\{1,1,1,1\},\{1,1,2\},\{2,2\},\{4\}
{1,1,1,1},{1,1,2},{2,2},{4} 共
4
4
4 种。特别地,
f
1
(
0
)
=
1
f^1(0)=1
f1(0)=1 。
求 f k ( n ) f^k(n) fk(n) 。定义 f k ( n ) = ∑ i = 0 n f k − 1 ( i ) f 1 ( n − i ) ( k > 1 ) f^k(n)=\sum_{i=0}^{n}f^{k-1}(i)f^1(n{-}i)\;(k>1) fk(n)=∑i=0nfk−1(i)f1(n−i)(k>1) 。
数据范围与提示
n
,
m
⩽
1
0
18
n,m\leqslant 10^{18}
n,m⩽1018 但
1
⩽
k
⩽
20
1\leqslant k\leqslant 20
1⩽k⩽20 。显然
m
⩾
2
m\geqslant 2
m⩾2 。
思路
我连完全背包都没想到,但 O U Y E \sf OUYE OUYE 现场 A C \rm AC AC,令我自惭形秽啊!
神的做法
“数学的目标是不需要聪明的想法。” ——《具体数学》
我就是聪明的想法。—— O U Y E \sf OUYE OUYE
对于
f
1
(
n
)
f^1(n)
f1(n),问题就是完全背包。枚举
1
1
1 的使用量,不难发现剩余容量需为
m
m
m 的倍数,因此有
f
(
n
)
=
∑
i
=
0
⌊
n
m
⌋
f
(
i
)
f(n)=\sum_{i=0}^{\lfloor{n\over m}\rfloor}f(i)
f(n)=i=0∑⌊mn⌋f(i)
更进一步地,我们有
∑
i
=
0
n
g
(
n
−
i
)
f
(
i
)
=
∑
i
=
0
⌊
n
m
⌋
f
(
i
)
∑
j
=
i
m
n
g
(
n
−
j
)
\sum_{i=0}^{n}g(n{-}i)f(i)=\sum_{i=0}^{\lfloor{n\over m}\rfloor}f(i)\sum_{j=im}^{n}g(n{-}j)
i=0∑ng(n−i)f(i)=i=0∑⌊mn⌋f(i)j=im∑ng(n−j)
记为 g ( n − i ) g(n{-}i) g(n−i) 的含义就在于此:我们就可以对 g ( x ) g(x) g(x) 做前缀和了。
如何维护?可以利用点值。维护 g ( k ) ( 0 ⩽ k ⩽ deg g ) g(k)\;(0\leqslant k\leqslant\deg g) g(k)(0⩽k⩽degg) 的值,则前缀和是容易求解的。然后 h ( x ) h(x) h(x) 取了步长为 m m m 的若干位置,就是复合一次函数,拉格朗日插值可以得到新点值。
现在考虑求解
f
k
(
n
)
f^k(n)
fk(n) 。相当于每个数字做
k
k
k 次完全背包。因此
∑
i
=
0
n
g
(
n
−
i
)
f
k
(
i
)
=
∑
i
=
0
n
f
k
−
1
(
i
)
∑
j
=
i
n
g
(
n
−
j
)
\sum_{i=0}^{n}g(n{-}i)f^k(i)=\sum_{i=0}^{n}f^{k-1}(i)\sum_{j=i}^{n}g(n{-}j)
i=0∑ng(n−i)fk(i)=i=0∑nfk−1(i)j=i∑ng(n−j)
这里 f k − 1 ( n ) f^{k-1}(n) fk−1(n) 并不是题面中的含义,而是 1 1 1 还能用多少个,因此 f 0 ( n ) = [ m ∣ n ] f k ( n m ) f^0(n)=[m\mid n]f^k({n\over m}) f0(n)=[m∣n]fk(mn) 。
因此还是对 g g g 做前缀和,不过要做 k k k 次。然后复合一次函数 x m + ( n m o d m ) xm+(n\bmod m) xm+(nmodm) 并不复杂。
使用线性的拉格朗日插值,则每次复杂度 O [ ( k + deg g ) deg g ] \mathcal O[(k+\deg g)\deg g] O[(k+degg)degg],递归 O ( log n ) \mathcal O(\log n) O(logn) 层有 deg g ⩽ k log n \deg g\leqslant k\log n degg⩽klogn,因此复杂度 O ( k 2 log 3 n ) \mathcal O(k^2\log^3n) O(k2log3n) 。
由于 k ⩽ deg g k\leqslant\deg g k⩽degg,后略。若使用 O [ deg g log 2 deg g ] \mathcal O[\deg g\log^2\deg g] O[degglog2degg] 拉格朗日插值,则复杂度为 O [ k log 2 n log 2 ( k log n ) ] \mathcal O[k\log^2 n\log^2(k\log n)] O[klog2nlog2(klogn)] 。
事实上,在同等复杂度条件下,我们可以在 FFP \textrm{FFP} FFP 和点值和 OGF \textrm{OGF} OGF 间任意转换了,所以怎么搞都可以。
青蛙做法
这是 qwaszx \textsf{qwaszx} qwaszx 的做法。很不幸,他被 O U Y E \sf OUYE OUYE 吊打了。
答案即
[
x
n
]
1
(
1
−
x
)
k
+
1
(
1
−
x
m
)
k
(
1
−
x
m
2
)
k
⋯
[x^n]\frac{1}{(1-x)^{k+1}(1-x^m)^k(1-x^{m^2})^k\cdots}
[xn](1−x)k+1(1−xm)k(1−xm2)k⋯1
分子分母同时乘有限长多项式
(
1
−
x
m
)
k
+
1
(
1
−
x
)
k
+
1
\frac{(1-x^m)^{k+1}}{(1-x)^{k+1}}
(1−x)k+1(1−xm)k+1,则分母为关于
x
m
x^m
xm 的多项式,因此分子必须提取模
m
m
m 与
n
n
n 同余的项,往下递归就得到形式
[
x
n
]
P
(
x
)
(
1
−
x
)
c
(
1
−
x
m
)
k
(
1
−
x
m
2
)
k
⋯
[x^n]\frac{P(x)}{(1-x)^{c}(1-x^m)^k(1-x^{m^2})^k\cdots}
[xn](1−x)c(1−xm)k(1−xm2)k⋯P(x)
不断用这种方法做,有 deg P ( x ) = O ( c ) = O ( k log n ) \deg P(x)=\mathcal O(c)=\mathcal O(k\log n) degP(x)=O(c)=O(klogn) 。青蛙说提取 [ x i m + r ] P ( x ) ( 1 − x m ) c ( 1 − x ) c [x^{im+r}]P(x)\frac{(1-x^m)^c}{(1-x)^c} [xim+r]P(x)(1−x)c(1−xm)c 可以 O ( c 2 ) \mathcal O(c^2) O(c2) 。如果用 O G F \rm OGF OGF 形式维护 P ( x ) P(x) P(x),好像还挺困难,因为有 O ( c ) \mathcal O(c) O(c) 项要提取,枚举 P ( x ) P(x) P(x) 和 ( 1 − x m ) c (1-x^m)^c (1−xm)c 中的每一项算贡献,就 O ( c 3 ) \mathcal O(c^3) O(c3) 了……
从
OUYE
\textsf{OUYE}
OUYE 处获得一些启发:用点值维护
P
(
x
)
P(x)
P(x) 。那么
∑
i
x
i
[
z
i
m
+
r
]
P
(
z
)
(
1
−
z
m
)
c
(
1
−
z
)
c
=
[
z
r
]
[
P
(
z
)
(
1
−
z
m
)
c
(
1
−
z
)
c
m
o
d
(
z
m
−
x
)
]
\sum_{i} x^i[z^{im+r}]P(z)\frac{(1-z^m)^c}{(1-z)^c}=[z^r]\left[P(z)\frac{(1-z^m)^c}{(1-z)^c}\bmod(z^m-x)\right]
i∑xi[zim+r]P(z)(1−z)c(1−zm)c=[zr][P(z)(1−z)c(1−zm)cmod(zm−x)]
如果 m m m 次剩余数量够多且 m m m 次单位根存在,那还是可以计算的。但这要求太高了。
再从
O
U
Y
E
\sf OUYE
OUYE 处获得一些启发:维护形式
P
(
x
)
=
∑
j
=
0
c
−
1
f
j
(
1
−
x
)
j
P(x)=\sum_{j=0}^{c-1} f_j(1{-}x)^{j}
P(x)=j=0∑c−1fj(1−x)j
因为 P ( x ) ( 1 − x ) c P(x)\over(1-x)^c (1−x)cP(x) 就对应 “系数为上升幂多项式”。所以它显然是可维护的。说到底,这已经完全是在用数学给 O U Y E \sf OUYE OUYE 一个蹩脚的注解。毕竟,数学究其一生也就是想避免 O U Y E \sf OUYE OUYE 呀。
代码
我信仰什么,我便实现哪种方法。
#include <cstdio>
#include <algorithm> // Almighty XJX yyds!
#include <cstring> // crashed god of math.
#include <cctype> // rainybunny root of the evil
using llong = long long;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define rep0(i,a,b) for(int i=(a); i!=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
inline int readint(){
int a = 0, c = getchar(), f = 1;
for(; !isdigit(c); c=getchar()) if(c == '-') f = -f;
for(; isdigit(c); c=getchar()) a = a*10+(c^48);
return a*f;
}
const int MAXN = 1200, MOD = 1e9+7;
inline llong qkpow(llong b, int q){
llong a = 1;
for(; q; q>>=1,b=b*b%MOD) if(q&1) a = a*b%MOD;
return a;
}
int fac[MAXN], invfac[MAXN];
void prepare(const int n = MAXN-1){
rep(i,*fac=1,n) fac[i] = int(llong(i)*fac[i-1]%MOD);
invfac[n] = int(qkpow(fac[n],MOD-2));
drep(i,n,1) invfac[i-1] = int(llong(invfac[i])*i%MOD);
}
int y[MAXN];
int valueAt(const int &n, const int &x){
if(x <= n) return y[x]; // known value
static int tmp[MAXN]; *tmp = 1;
rep0(i,0,n) tmp[i+1] = int(llong(x-i)*tmp[i]%MOD);
int nowv = (n&1)?-1:1, res = 0;
drep(i,n,0){
res = int((res+llong((i&1)?-nowv:nowv)*tmp[i]
%MOD*invfac[i]%MOD*invfac[n-i]%MOD*y[i])%MOD);
nowv = int(llong(x-i)*nowv%MOD);
}
return res < 0 ? res+MOD : res;
}
int tmp[MAXN];
int main(){
prepare();
llong n, m; scanf("%lld %lld",&n,&m);
int k = readint(), len = 0;
y[0] = 1; // deg = 0
for(llong mm=m%MOD; n; n/=m){
rep(_,1,k){ // k times prefix-sum
++ len, y[len] = valueAt(len-1,len);
rep(i,1,len) y[i] = (y[i]+y[i-1])%MOD;
}
int r = int((n%m)%MOD);
rep(i,0,len) tmp[i] = valueAt(len,int((i*mm+r)%MOD));
memcpy(y,tmp,(len+1)<<2); // move
}
printf("%d\n",y[0]);
return 0;
}