[ACNOI2022]高中数列题

题目

题目描述
$\mathsf{O}\mathsf{n}\mathsf{e}\mathsf{I}\mathsf{n}\mathsf{D}\mathsf{a}\mathsf{r}\mathsf{k}$ 讨厌数学，如同老鼠讨厌猫；并不是猫不好，只是老鼠不喜欢。

可悲的是，他不得不解决这样一个问题：求数列 { a n } \{a_n\} {an​} 的第 $n_0$ 项。其中 $a_n$ 的递推公式是
$$a_n=a_{n-1}+f(n)\cdot a_{\lfloor \frac{n+b}{c}\rfloor }(n⩾1)$$

其中 $f(x)={\sum }_{i=0}^m{v}_i{x}^i$ 是给定的多项式。初值 $a_0$ 给定。保证 $b+1<c$ 。

数据范围与提示
$n⩽10^{18}$ 且 $m⩽20$，保证 $a_0,b,c$ 是非负整数。

思路

不难发现
$$a_n=a_0+\sum _{i=1}^{n}f(i)a_{\lfloor \frac{i+b}{c}\rfloor }$$

为了方便，暂且记 $w(i)=\lfloor \frac{i+b}{c}\rfloor$ 。右边的求和怎么处理？最简单的想法是，枚举 $w(i)$ 的取值，就像类欧几里得。记 $r(i)=(i+1)c-b-1$，即使得 $w(x)=i$ 的最大 $x$ 。记 $g(i)$ 为 $f(i)$ 的前缀和，那么有
$$\begin{gathered} \sum _{i=1}^{n}f(i)a_{w(i)}=\sum _{i=1}^{w(n)-1}{a}_i\cdot \{g[r(i)]-g[r(i-1)]\} \\ +a_0\cdot g[r(0)]+a_{w(n)}\{g(n)-g[r(n-1)]\} \end{gathered}$$

$h(i)=g[r(i)]-g[r(i-1)]$ 还是关于 $i$ 的多项式；所以我们又只需要求出形如 ${\sum }_{i=1}^n{a}_{i}g(i)$ 的式子。这个还是用老方法，大力展开，记 $h(i)$ 为 $g(i)$ 的后缀和（到 $n$ 为止）则
$$\sum _{i=1}^n{a}_{i}g(i)=a_{0}h(1)+\sum _{i=1}^n{a}_{w(i)}h(i)f(i)$$

又和上面的问题一样了，枚举 $w(i)$，总是得到一个递归的问题。

每递归一次，系数 $g(i)$ 就变成 $h(i)f(i)$ 的区间和，次数会增加 $(m+1)$ 。递归 $\mathcal{O}(\log n)$ 次，看起来也没什么大不了。

但是我们还需要单点求出 $a_{w(n)}$ 的值！也需要递归！也就是说，我们要递归到两个子状态；复杂度直接完蛋！

但是，当 $c$ 较小时，$\frac{n}{c^i}$ 的重复概率较大；若将单点求值的结果记忆化，可以很大的提升运行效率。但是再怎么说也不算正解。

所以正解是什么呢？我觉得有点离谱，就是 再次展开；很难想象，出题人是怎么预见这一做法的可行性的？不是按照定义展开，而是按照推论展开。

$$\begin{array}{rl}\sum _{i=1}^n{a}_{w(i)}g(i)& =\sum _{i=1}^{n}g(i)(a_0+\sum _{j=1}^{w(i)}{a}_{w(j)}f(j))\\ & =a_0\sum _{i=1}^{n}g(i)+\sum _{j=1}^{w(n)}{a}_{w(j)}f(j)\sum _{w(i)⩾j}g(i)\end{array}$$

推到这里，大家应该已经一眼望穿了。我还是要吐槽一句：为什么，要试图 保留 $w(i)$ 型求和，而不是前缀求和？真是神了！

于是记 $h(i)$ 为 $g(i)$ 前缀和，此地空余 $a_{0}h(n)+{\sum }_{j=1}^{w(n)}{a}_{w(j)}f(j)\{h(n)-h[r(j-1)]\}$，只需要递归一次；每次复杂度是 $\mathcal{O}({\mathrm{deg}}^2)$，总复杂度 $\mathcal{O}(m^2{\log }^{3}n)$ 。

而且加上 $\text{FFT}$ 就可以变成毒瘤 $\mathcal{O}[m{\log }^{2}n\cdot {\log }^2(m\log n)]$ 了……

代码

由于本来实现的是糟糕的做法，这份代码的封装过度，比较冗长。

#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;
}