题目
题目描述
O
n
e
I
n
D
a
r
k
\sf OneInDark
OneInDark 讨厌数学,如同老鼠讨厌猫;并不是猫不好,只是老鼠不喜欢。
可悲的是,他不得不解决这样一个问题:求数列
{
a
n
}
\{a_n\}
{an} 的第
n
0
n_0
n0 项。其中
a
n
a_n
an 的递推公式是
a
n
=
a
n
−
1
+
f
(
n
)
⋅
a
⌊
n
+
b
c
⌋
(
n
⩾
1
)
a_n=a_{n-1}+f(n)\cdot a_{\lfloor{n+b\over c}\rfloor}\quad(n\geqslant 1)
an=an−1+f(n)⋅a⌊cn+b⌋(n⩾1)
其中 f ( x ) = ∑ i = 0 m v i x i f(x)=\sum_{i=0}^{m}v_ix^i f(x)=∑i=0mvixi 是给定的多项式。初值 a 0 a_0 a0 给定。保证 b + 1 < c b+1<c b+1<c 。
数据范围与提示
n
⩽
1
0
18
n\leqslant 10^{18}
n⩽1018 且
m
⩽
20
m\leqslant 20
m⩽20,保证
a
0
,
b
,
c
a_0,b,c
a0,b,c 是非负整数。
思路
不难发现
a
n
=
a
0
+
∑
i
=
1
n
f
(
i
)
a
⌊
i
+
b
c
⌋
a_n=a_0+\sum_{i=1}^{n}f(i)a_{\lfloor{i+b\over c}\rfloor}
an=a0+i=1∑nf(i)a⌊ci+b⌋
为了方便,暂且记
w
(
i
)
=
⌊
i
+
b
c
⌋
w(i)=\lfloor{i+b\over c}\rfloor
w(i)=⌊ci+b⌋ 。右边的求和怎么处理?最简单的想法是,枚举
w
(
i
)
w(i)
w(i) 的取值,就像类欧几里得。记
r
(
i
)
=
(
i
+
1
)
c
−
b
−
1
r(i)=(i{+}1)c{-}b{-}1
r(i)=(i+1)c−b−1,即使得
w
(
x
)
=
i
w(x)=i
w(x)=i 的最大
x
x
x 。记
g
(
i
)
g(i)
g(i) 为
f
(
i
)
f(i)
f(i) 的前缀和,那么有
∑
i
=
1
n
f
(
i
)
a
w
(
i
)
=
∑
i
=
1
w
(
n
)
−
1
a
i
⋅
{
g
[
r
(
i
)
]
−
g
[
r
(
i
−
1
)
]
}
+
a
0
⋅
g
[
r
(
0
)
]
+
a
w
(
n
)
{
g
(
n
)
−
g
[
r
(
n
−
1
)
]
}
\sum_{i=1}^{n}f(i)a_{w(i)}= \sum_{i=1}^{w(n)-1}a_i\cdot\big\{g[r(i)]-g[r(i{-}1)]\big\}\\ +a_0\cdot g[r(0)]+a_{w(n)} \big\{g(n)-g[r(n{-}1)]\big\}
i=1∑nf(i)aw(i)=i=1∑w(n)−1ai⋅{g[r(i)]−g[r(i−1)]}+a0⋅g[r(0)]+aw(n){g(n)−g[r(n−1)]}
h
(
i
)
=
g
[
r
(
i
)
]
−
g
[
r
(
i
−
1
)
]
h(i)=g[r(i)]-g[r(i{-}1)]
h(i)=g[r(i)]−g[r(i−1)] 还是关于
i
i
i 的多项式;所以我们又只需要求出形如
∑
i
=
1
n
a
i
g
(
i
)
\sum_{i=1}^{n}a_ig(i)
∑i=1naig(i) 的式子。这个还是用老方法,大力展开,记
h
(
i
)
h(i)
h(i) 为
g
(
i
)
g(i)
g(i) 的后缀和(到
n
n
n 为止)则
∑
i
=
1
n
a
i
g
(
i
)
=
a
0
h
(
1
)
+
∑
i
=
1
n
a
w
(
i
)
h
(
i
)
f
(
i
)
\sum_{i=1}^{n} a_ig(i) = a_0 h(1)+ \sum_{i=1}^{n} a_{w(i)} h(i)f(i)
i=1∑naig(i)=a0h(1)+i=1∑naw(i)h(i)f(i)
又和上面的问题一样了,枚举 w ( i ) w(i) w(i),总是得到一个递归的问题。
每递归一次,系数 g ( i ) g(i) g(i) 就变成 h ( i ) f ( i ) h(i)f(i) h(i)f(i) 的区间和,次数会增加 ( m + 1 ) (m{+}1) (m+1) 。递归 O ( log n ) \mathcal O(\log n) O(logn) 次,看起来也没什么大不了。
但是我们还需要单点求出 a w ( n ) a_{w(n)} aw(n) 的值!也需要递归!也就是说,我们要递归到两个子状态;复杂度直接完蛋!
但是,当
c
c
c 较小时,
n
c
i
n\over c^i
cin 的重复概率较大;若将单点求值的结果记忆化,可以很大的提升运行效率。但是再怎么说也不算正解。
所以正解是什么呢?我觉得有点离谱,就是 再次展开;很难想象,出题人是怎么预见这一做法的可行性的?不是按照定义展开,而是按照推论展开。
∑ i = 1 n a w ( i ) g ( i ) = ∑ i = 1 n g ( i ) ( a 0 + ∑ j = 1 w ( i ) a w ( j ) f ( j ) ) = a 0 ∑ i = 1 n g ( i ) + ∑ j = 1 w ( n ) a w ( j ) f ( j ) ∑ w ( i ) ⩾ j g ( i ) \begin{aligned} \sum_{i=1}^{n}a_{w(i)} g(i) &=\sum_{i=1}^{n}g(i)\big(a_0+\sum_{j=1}^{w(i)}a_{w(j)}f(j)\big)\\ &=a_0\sum_{i=1}^{n}g(i)+\sum_{j=1}^{w(n)}a_{w(j)}f(j)\sum_{w(i)\geqslant j}g(i) \end{aligned} i=1∑naw(i)g(i)=i=1∑ng(i)(a0+j=1∑w(i)aw(j)f(j))=a0i=1∑ng(i)+j=1∑w(n)aw(j)f(j)w(i)⩾j∑g(i)
推到这里,大家应该已经一眼望穿了。我还是要吐槽一句:为什么,要试图 保留 w ( i ) w(i) w(i) 型求和,而不是前缀求和?真是神了!
于是记 h ( i ) h(i) h(i) 为 g ( i ) g(i) g(i) 前缀和,此地空余 a 0 h ( n ) + ∑ j = 1 w ( n ) a w ( j ) f ( j ) { h ( n ) − h [ r ( j − 1 ) ] } a_0h(n)+\sum_{j=1}^{w(n)}a_{w(j)}f(j)\big\{h(n)-h[r(j{-}1)]\big\} a0h(n)+∑j=1w(n)aw(j)f(j){h(n)−h[r(j−1)]},只需要递归一次;每次复杂度是 O ( deg 2 ) \mathcal O(\deg^2) O(deg2),总复杂度 O ( m 2 log 3 n ) \mathcal O(m^2\log^3 n) O(m2log3n) 。
而且加上 FFT \textit{FFT} FFT 就可以变成毒瘤 O [ m log 2 n ⋅ log 2 ( m log n ) ] \mathcal O[m\log^2 n\cdot \log^2(m\log n)] O[mlog2n⋅log2(mlogn)] 了……
代码
由于本来实现的是糟糕的做法,这份代码的封装过度,比较冗长。
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cctype>
using namespace std;
typedef long long llong;
# define rep(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<<3)+(a<<1)+(c^48);
return a*f;
}
const int MOD = 1004535809;
inline int modMinus(int x,const int y){
return ((x -= y) < 0) ? (x+MOD) : x;
}
inline int modAdd(int x,const int y){
return ((x += y) >= MOD) ? (x-MOD) : x;
}
const int MAXN = 1265;
int c[MAXN][MAXN], inv[MAXN];
struct Polynomial{
int x[MAXN], len;
void negate(){ rep(i,0,len) x[i] = MOD-x[i]; }
Polynomial operator + (const Polynomial &b) const {
Polynomial v;
if(len > b.len){
rep(i,0,b.len) v.x[i] = modAdd(x[i],b.x[i]);
memcpy(v.x+b.len+1,x+b.len+1,(len-b.len)<<2);
v.len = len; // no exception
}
else{
rep(i,0,len) v.x[i] = modAdd(x[i],b.x[i]);
memcpy(v.x+len+1,b.x+len+1,(b.len-len)<<2);
v.len = b.len; // become longer one
while(v.len && !v.x[v.len]) -- v.len;
}
return v;
}
Polynomial operator - (const Polynomial &b) const {
Polynomial v;
if(len > b.len){
rep(i,0,b.len) v.x[i] = modMinus(x[i],b.x[i]);
memcpy(v.x+b.len+1,x+b.len+1,(len-b.len)<<2);
v.len = len; // no exception
}
else{
rep(i,0,len) v.x[i] = modMinus(x[i],b.x[i]);
rep(i,len+1,b.len) v.x[i] = MOD-b.x[i];
v.len = b.len; // bigger first
while(v.len && !v.x[v.len]) -- v.len;
}
return v;
}
Polynomial operator * (const Polynomial &b) const {
Polynomial res; res.len = len+b.len;
memset(res.x,0,(res.len+1)<<2);
rep(i,0,len) rep(j,0,b.len) // bruteforce
res.x[i+j] = int((res.x[i+j]+llong(x[i])*b.x[j])%MOD);
return res; // no exception
}
Polynomial operator * (const llong &v) const {
Polynomial res; res.len = len;
rep(i,0,len) res.x[i] = int(v*x[i]%MOD);
return res; // like a transform
}
/// the @c Polynomial when variable is @p k * x + @p b
Polynomial variable(int k,int b) const {
static int powk[MAXN], powb[MAXN];
Polynomial res; res.len = len;
memset(res.x,0,(len+1)<<2);
rep(i,powk[0]=powb[0]=1,len){
powk[i] = int(llong(powk[i-1])*k%MOD);
powb[i] = int(llong(powb[i-1])*b%MOD);
}
rep(i,0,len) rep(j,0,i) // bruteforce
res.x[j] = int((res.x[j]+llong(x[i])*c[i][j]
%MOD*powk[j]%MOD*powb[i-j])%MOD);
return res;
}
int valueAt(int v) const {
int res = 0, t = 1;
rep(i,0,len){
res = int((res+llong(t)*x[i])%MOD);
t = int(llong(t)*v%MOD); // 0^0 = 1
}
return res;
}
};
Polynomial divide_monomial(Polynomial &a,int b){
Polynomial v; v.len = a.len-1;
for(int i=a.len-1,r=a.x[a.len]; ~i; --i){
v.x[i] = r, r = (a.x[i]-llong(b)*v.x[i]%MOD);
if(r < 0) r += MOD; // prefer positive
}
return v;
}
Polynomial prefix_sum(const Polynomial &a){
static int y[MAXN]; // point-value
Polynomial res, all; res.len = -1;
all.len = 0, all.x[0] = 1;
rep(i,0,a.len+1){
static Polynomial jb; jb.len = 1;
jb.x[0] = MOD-i, jb.x[1] = 1; // (x-i)
all = all*jb; y[i] = a.valueAt(i);
}
rep(i,1,a.len+1) y[i] = modAdd(y[i],y[i-1]);
rep(i,0,a.len+1){
if((a.len^i)&1) // minus bigger get negative
res = res+divide_monomial(all,MOD-i)
*inv[i]*inv[a.len+1-i]*y[i];
else res = res-divide_monomial(all,MOD-i)
*inv[i]*inv[a.len+1-i]*y[i];
}
return res;
}
const int MAXM = 16;
int a0, b0, c0, known[MAXM];
Polynomial trans;
llong solve(llong n,const Polynomial &g){
if(n < MAXM){ // calculate by definition
int res = 0; // with range
rep(i,1,n) res = int((res+g.valueAt(i)
*llong(known[(i+b0)/c0]))%MOD);
return res; // bruteforce
}
Polynomial h = prefix_sum(g); h.x[0] = 0;
const int hn = h.valueAt(int(n%MOD));
h = h.variable(c0,MOD-b0-1);
h.negate(); h.x[0] = modAdd(hn,h.x[0]);
return (solve((n+b0)/c0,h*trans)+llong(hn)*a0)%MOD;
}
int logint(llong n,int v){
int res = 0;
while(n != 1) n = (n+v-1)/v, ++ res;
return res;
}
int main(){
freopen("seq.in","r",stdin);
freopen("seq.out","w",stdout);
llong n; scanf("%lld",&n);
a0 = readint(), b0 = readint(), c0 = readint();
int m = readint(); trans.len = m;
rep(i,0,m) trans.x[i] = readint();
const int LEN = logint(n,c0)*(m+1);
rep(i,0,LEN+1) rep(j,c[i][0]=1,i)
c[i][j] = modAdd(c[i-1][j-1],c[i-1][j]);
rep(i,(inv[1]=1)<<1,LEN+1)
inv[i] = int(llong(MOD-MOD/i)*inv[MOD%i]%MOD);
rep(i,inv[0]=1,LEN+1) // inverse of factorial
inv[i] = int(llong(inv[i-1])*inv[i]%MOD);
known[0] = a0;
for(int i=1; i!=MAXM; ++i)
known[i] = int((known[i-1]+known[(i+b0)/c0]
*llong(trans.valueAt(i)))%MOD);
printf("%lld\n",(solve(n,trans)+a0)%MOD);
return 0;
}