题目

题目描述
OneInDark\sf OneInDarkOneInDark 讨厌数学，如同老鼠讨厌猫；并不是猫不好，只是老鼠不喜欢。

可悲的是，他不得不解决这样一个问题：求数列 {an}\{a_n\}{an} 的第 n0n_0n0 项。其中 ana_nan 的递推公式是
an=an−1+f(n)⋅a⌊n+bc⌋(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=0mvixif(x)=\sum_{i=0}^{m}v_ix^if(x)=∑i=0mvixi 是给定的多项式。初值 a0a_0a0 给定。保证 b+1<cb+1<cb+1<c 。

数据范围与提示
n⩽1018n\leqslant 10^{18}n⩽1018 且 m⩽20m\leqslant 20m⩽20，保证 a0,b,ca_0,b,ca0,b,c 是非负整数。

思路

不难发现
an=a0+∑i=1nf(i)a⌊i+bc⌋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+bc⌋w(i)=\lfloor{i+b\over c}\rfloorw(i)=⌊ci+b⌋ 。右边的求和怎么处理？最简单的想法是，枚举 w(i)w(i)w(i) 的取值，就像类欧几里得。记 r(i)=(i+1)c−b−1r(i)=(i{+}1)c{-}b{-}1r(i)=(i+1)c−b−1，即使得 w(x)=iw(x)=iw(x)=i 的最大 xxx 。记 g(i)g(i)g(i) 为 f(i)f(i)f(i) 的前缀和，那么有
∑i=1nf(i)aw(i)=∑i=1w(n)−1ai⋅{g[r(i)]−g[r(i−1)]}+a0⋅g[r(0)]+aw(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)] 还是关于 iii 的多项式；所以我们又只需要求出形如 ∑i=1naig(i)\sum_{i=1}^{n}a_ig(i)∑i=1naig(i) 的式子。这个还是用老方法，大力展开，记 h(i)h(i)h(i) 为 g(i)g(i)g(i) 的后缀和（到 nnn 为止）则
∑i=1naig(i)=a0h(1)+∑i=1naw(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) 次，看起来也没什么大不了。

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

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

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

∑i=1naw(i)g(i)=∑i=1ng(i)(a0+∑j=1w(i)aw(j)f(j))=a0∑i=1ng(i)+∑j=1w(n)aw(j)f(j)∑w(i)⩾jg(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) 前缀和，此地空余 a0h(n)+∑j=1w(n)aw(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(m2log⁡3n)\mathcal O(m^2\log^3 n)O(m2log3n) 。

而且加上 FFT\textit{FFT}FFT 就可以变成毒瘤 O[mlog⁡2n⋅log⁡2(mlog⁡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 onewhile(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 firstwhile(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) // bruteforceres.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 bPolynomial 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) // bruteforceres.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-valuePolynomial 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 negativeres = 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 definitionint res = 0; // with rangerep(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 factorialinv[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;
}

[ACNOI2022]高中数列题相关推荐

[硫化铂]高中数列题
高中数列题题解首先,我们的aaa的递推式显然是一个前缀和的形式,我们不妨将其全部展开,我们定义pi=i+bc,qi=ci−bp_i=\frac{i+b}{c},q_i=ci-bpi=ci+b, ...
高中数列题目的C++实现
/*高中裂项相消数列原题是: 求a(n)= 1/(1*2)+1/(2*3)...+1/[n*(n+1)] 数列前n项和为S(n),(0<n<2^31-1,且n为整数) * (1)求a(n) ...
【校招VIP】产品经理行测之数列题
考点介绍: 行测题是产品经理的必考题,对于数推题目,建议先通过作差.作和,将难度较低.相对容易的题目解决,再递推解决剩下的非特征数列.数列题有些题目本身并不重要,决定能否通过面试的是面试者的逻辑表达是 ...
html三角形坐标图怎么看,地理中三角图_高中地理题那种三角形的坐标图要怎么看啊_淘题吧...
Ⅰ 地理中三角形代表什么(地图中) 这个,要看那幅地图,一般来说,是指上峰.但如果你说的图上,是将那个作为一个标注,就要具体分析了 Ⅱ 地理的三角坐标图怎么看请赐教三角来坐标图是近几年地理源高考试题 ...
可近似看作直线的是_高中物理题根之一：《匀变速直线运动的规律》
"九层之台,起于垒土".高中物理开篇就是运动学.作为动力学的基础知识,在每年的高考中,或者单独命题,或者渗透在动力学问题中,都要对运动学的概念和规律进行考查. 理顺基础知识,溯清本 ...
PAT - 天梯赛 L3-013 非常弹的球（高中物理题）
L3-013. 非常弹的球时间限制 100 ms 内存限制 65536 kB 代码长度限制 8000 B 判题程序 Standard 作者俞勇(上海交通大学) 刚上高一的森森为了学好物理,买了一个 ...
Python解高考数学数列题
题目解读对于这种给定的数列是不确定的,但是任意一组满足要求的数列都可以得到同样结论的题我们直接代特值判断即可. 最容易想的就是令 a(n) = 1 数列生成我们可以使用Python编程语言来生成和 ...
一道发源于游戏中的数列题（求项数最大值）
之前上高中的时候打游戏时想到的一道题目扑克飞镖:每使用一次法力值消耗增加1点初始法力值为零费每次使用后回手吹气:获得五点法力值手牌中法力值最高的牌消耗归零假设魔术师现在手中仅有上述两张手牌 ...
青蛙跳台阶：我如何得知它是一道斐波那契数列题？——应用题破题“三板斧”
本文以C语言实现. 目录前言一.斐波那契阿数列基础知识二.引例:青蛙跳台阶三.破题分析:举例归纳 1. 三板斧的使用举例模拟(必要时画图) 找规律 2. 代码展示四.拓展用例:矩形覆盖问 ...

[ACNOI2022]高中数列题

题目

思路

代码

[ACNOI2022]高中数列题相关推荐

最新文章

热门文章