一些闲话 & 免责声明:感觉好久没有更过学习笔记了,现在才来学这个不知道来不来得及呢。因为笔者是一个数学菜鸡,所以并不会很严谨(。
问题描述
LSB-first 算法
从线性递推式到分式第 nn 项
记 Q(x)=1−f1x−f2x2−⋯−fkxkQ(x)=1−f1x−f2x2−⋯−fkxk,F(x)F(x) 为递推数列的 OGF。我们需要求解 F(x) 的 n 次项系数。这里直接给出一个结论:
存在一个最高次项的次数小于 k 的多项式 P(x) 满足
F(x)=P(x)Q(x)
由于 P(x) 的度数 <k,而 F(x) 的前 k 项是已知的,所以可以仅用 F(x) 的前 k 项与 Q(x) 进行多项式乘法求出 P(x).
剩下的问题是,要怎么快速计算 P(x)Q(x) 的第 n 项。
Bostan-mori 算法
首先有
[xn]P(x)Q(x)=[xn]P(x)Q(−x)Q(x)Q(−x)
可以发现,Q(x)Q(−x) 是一个偶函数(考虑将 x,−x 分别代入),也就是说,它的奇次项均为零,于是可以构造多项式 v(x2)=Q(x)Q(−x),同时对分子奇偶分类可得
[xn]ceven(x2)+x⋅codd(x2)v(x2)={[xn/2]ceven(x)v(x),if n is even[xn/2]codd(x)v(x),if n is odd
对于后面的一步推导,考虑 v(x) 的逆也是偶函数即可。于是可以进行递推求解。复杂度 O(klogklogn).
代码实现
三个 998244353 相加又忘记开 long long 了,我自裁 🔫。
/* 先复习动人的墓志,再背诵浪漫的情诗。 */
# include <cstdio>
# include <cctype>
# define print(x,y) write(x), putchar(y)
template <class T>
inline T read(const T sample) {
T x=0; char s; bool f=0;
while(!isdigit(s=getchar())) f|=(s=='-');
for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
return f? -x: x;
}
template <class T>
inline void write(T x) {
static int writ[50], w_tp=0;
if(x<0) putchar('-'), x=-x;
do writ[++w_tp]=x-x/10*10, x/=10; while(x);
while(putchar(writ[w_tp--]^48), w_tp);
}
# include <cstring>
# include <iostream>
using namespace std;
const int MAXN = 170000;
const int mod = 998244353, G = 3;
void mul(int& x,int y) { x=1ll*x*y%mod; }
int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }
int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
int inv(int x,int y=mod-2,int r=1) {
for(; y; y>>=1, x=1ll*x*x%mod)
if(y&1) r=1ll*r*x%mod; return r;
} const int iG = inv(G);
int P[MAXN], rev[MAXN], g[MAXN];
int lim, bit, Inv, n, k, Q[MAXN];
void ntt(int* f,bool opt=0) {
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1) {
const int wn = inv(opt?iG:G,(mod-1)/(mid<<1));
for(int i=0; i<lim; i+=(mid<<1)) { int w=1;
for(int j=0; j<mid; ++j, w=1ll*w*wn%mod) {
const int tmp = 1ll*f[i|j|mid]*w%mod;
f[i|j|mid] = dec(f[i|j], tmp),
f[i|j] = inc(f[i|j], tmp);
}
}
} if(!opt) return;
for(int i=0;i<lim;++i) mul(f[i],Inv);
}
int main() {
n=read(9), k=read(9); Q[0]=lim=1;
for(int i=1;i<=k;++i) Q[i] = (0ll+(mod<<1)-read(9))%mod;
for(int i=0;i<k;++i) P[i] = (0ll+(mod<<1)+read(9))%mod;
while(lim<(k<<1|1)) lim<<=1, ++bit;
for(int i=0;i<lim;++i)
rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
ntt(P), ntt(Q); Inv = inv(lim);
for(int i=0;i<lim;++i) mul(P[i],Q[i]);
ntt(P,1), ntt(Q,1);
const int SIZE = sizeof(int)*(lim-k-1);
memset(P+k,0,sizeof(int)*(lim-k));
memset(Q+k+1,0,SIZE);
for(; n; n>>=1) {
for(int i=0;i<=k;++i)
g[i] = (i&1)?(mod-Q[i])%mod:Q[i];
memset(g+k+1,0,SIZE);
ntt(g), ntt(P), ntt(Q);
for(int i=0;i<lim;++i)
mul(P[i],g[i]), mul(Q[i],g[i]);
ntt(P,1), ntt(Q,1);
for(int i=0;i<=k;++i) Q[i]=Q[i<<1];
for(int i=0;i<=k;++i) P[i]=P[i<<1|(n&1)];
memset(P+k+1,0,SIZE), memset(Q+k+1,0,SIZE);
} print((mul(P[0],inv(Q[0])), P[0]), '\n');
return 0;
}
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek 解答了困扰我五年的技术问题
· 为什么说在企业级应用开发中,后端往往是效率杀手?
· 用 C# 插值字符串处理器写一个 sscanf
· Java 中堆内存和栈内存上的数据分布和特点
· 开发中对象命名的一点思考
· DeepSeek 解答了困扰我五年的技术问题。时代确实变了!
· PPT革命!DeepSeek+Kimi=N小时工作5分钟完成?
· What?废柴, 还在本地部署DeepSeek吗?Are you kidding?
· DeepSeek企业级部署实战指南:从服务器选型到Dify私有化落地
· 程序员转型AI:行业分析
2021-07-22 [杂题合集] One Last Blog