一直想写个好点的多项式模板,最近终于想清楚了怎么写,于是就写了一个!代码在本文最后。
众所周知,常见的多项式模板经常有这样两个问题:速度很慢、用起来很麻烦。
速度问题是老生常谈的。NTT我们可以找一个好的板子对着抄。形式幂级数运算我们可以采用 https://negiizhao.blog.uoj.ac/blog/4671 中的优化。在线卷积我们可以使用多叉降低复杂度,循环卷积和小范围暴力降低常数。那么本文考虑的主要问题是... 如何让你的模板封装得阳间一些?
我随机选取了若干份uoj上相关题目的AC代码:
一个字,惨。
我们先考虑形式幂级数和多项式的部分。这部分我的主要做法是,将每个多项式/形式幂级数记录精确度,即记录 $\bmod x^n$ 中的 $n$,在运算时保持该精确度不变。在牛顿迭代解 $G(f)=0$ 时,我们可以使用两个 lambda 传入 $G(f)$ 和 $G'(f)$,lambda 保持 $f$ 的精确度不变即可。我的模板中具体对应的是 shrink_len
属性,-1
表示为多项式,否则为精确次幂等于该值的集合幂级数。使用这个思路实现各运算即可。
以下各题均省略推导过程,代码均省略模板部分,模数均为 $998244353$。
loj150 挑战多项式
给定 $n$ 次多项式 $F(x)$,求 $G(x)$ 满足 $\displaystyle G(x) \equiv \left({\left({1+\ln\left({2+F(x)-F(0)-{\exp\left({\int_0^x\frac{1}{\sqrt{F(t)}}\textrm{d}t}\right)}}\right)}\right)^k}\right)^\prime \pmod {x^n}$。$n\le 10^5$
int n,m,f[123456];
int main()
{
scanf("%d%d",&n,&m);
default_shrink=++n;
for(int i=0;i<n;++i) scanf("%d",f+i);
poly p(vector<int>(f,f+n));
p=p-gexp(gint(ginv(gsqrt(p,mod_sqrt(f[0],MOD)))));
p=gcorner(gln(gcorner(p,{1})),{1});
gde(gexp(mi(m)*gln(p))).print_len(n-1,'\n');
}
Quite easily done! 跑了 sum 4269ms/max 228ms,在写作本文时排在ranklist第一页。
luogu6112 直接自然溢出啥事没有 加强版 (30%)
求出多项式 $F(x)$ 使得 $\left(x^{7}+x^{6}+x^{5}-2 x^{4}+x^{2}\right) F(x)^{2}+\left(x^{5}-x^{4}-2 x^{3}+2 x^{2}+x-1\right) F(x)+\left(x^{4}-2 x^{2}+1\right)\equiv 0 \pmod {x^n}$。$n\le 10^5$
int n;
int main()
{
cin>>n;
poly
p="x^7+x^6+x^5-2x^4+x^2"_p,
q="x^5-x^4-2x^3+2x^2+x-1"_p,
r="x^4-2x^2+1"_p,
a=gnewton(
[&](poly g) {
return (g*p+q)*g+r;
},
[&](poly g) {
return 2*g*p+q;
},
1,n+1);
cout<<a.get(n)<<"\n";
}
C++11就是好。
uoj50 链式法则
求出多项式 $F(x)$ 使得 $F'(x) \equiv \frac{1}{2}A(x)F^2(x) + 1 \pmod {x^n}$。$n\le 200000$
模板里咋能没有一阶微分方程呢。
int n; poly a(-1);
char str[288888];
int main()
{
scanf("%d%s",&n,str);
a.reserve(n+1);
for(int i=0;i<n;++i)
if(str[i]=='1') a[i]=rfac[i];
auto g=gnewton_d(
[&](poly x) {
x=mi((MOD+1)/2)*x*x*a;
++x[0]; return x;
},
[&](poly x) {
return x*a;
},0,n+1);
for(int i=1;i<=n;++i) {
mi ans=g.get(i)*fac[i];
printf("%d\n",(int)ans);
}
}
跑了 sum 2851ms/max 1464ms,在写作本文时(可以)排在ranklist第五页。
loj6684 有根无标号「奇树」计数
求出多项式 $F(x)$ 使得 $F(x) \equiv x\text{MSet}(\text{MSet}(F(x))-1) \pmod {x^n}$($\text{MSet}$ 的定义可参见 此处)。$n\le 100000$
这里我们进行牛顿迭代,注意到 $\text{MSet} (f(z))=\exp\left(\sum_{k\ge 1}f(z^k)/k\right)$,而倍增时 $\exp\left(\sum_{k\ge 2}f(z^k)/k\right)$ 有足够的精确度,因此我们可以将该部分视为一个常多项式,认为 $\text{MSet} (f(z))=\exp f(z)~\underbrace{\exp\left(\sum_{k\ge 2}f(z^k)/k\right)}_{\text{常多项式}}$。
由于模板中显然没有提供计算 $\sum_{k\ge 2}f(z^k)/k$ 的接口,我们需要手(kao)写(bei)。
poly MSet_(poly p,int st) {
p.fit_shrink();
assert(p.get(0)==0);
poly q(p.shrink_len);
q.fit_shrink();
for(int i=st;i<p.size();++i) {
mi iv=rfac[i]*fac[i-1];
for(int j=1;i*j<p.size();++j)
q[i*j]=q[i*j]+p[j]*iv;
}
return q;
}
int n;
int main()
{
scanf("%d",&n);
poly o=gnewton(
[&](poly u) {
return gshl(MSet(gcorner(MSet(u),{0})),1)-u;
},
[&](poly u) {
poly P=MSet_(u,2),
z=gcorner(gexp(P+u),{0}),
Q=MSet_(z,2);
//x*exp(Pexp(F)-1)Q-F=0
return gcorner(gshl(gexp(z+u+P+Q),1),{-1});
},
0,n+1);
for(int i=1;i<=n;++i)
printf("%d\n",(int)o.get(i));
}
在经过一些简单的常数优化(注意到 $G$ 和 $G'$ 都需要求出 $\text{MSet}$,可以进行记忆化)后,跑了 max 267ms,在写作本文时(可以)排在ranklist第一(忽略一位打表的同学)。
loj6538 烷基计数 加强版 加强版
求出多项式 $F(x)$ 使得 $F(z)\equiv z \frac{F(z)^{3}+3 F(z) F\left(z^{2}\right)+2 F\left(z^{3}\right)}{6}+1\pmod {x^n}$。$n\le 100000$
类似上题进行牛顿迭代。
int n;
int main() {
cin>>n;
cout<<gnewton(
[&](const poly& s) {
poly t=s;
t=inv(6)*(t*(t*t+3*gamp(t,2))+2*gamp(t,3));
t.insert(0); t=t+1-s; return t;
},
[&](poly t) {
t=inv(2)*(t*t+gamp(t,2));
t.insert(0); t=t-1; return t;
},
1,n+1).get(n)<<"\n";
}
跑了 sum 298ms/max 98ms,在写作本文时排在ranklist第一页。
此外,模板中还实现了快速线性递推、快速多点求值插值、幂和、简单函数复合逆等算法。
2021年了,牛顿法等朴素迭代方法已经不足以满足我们。接下来我们来考虑在线卷积算法。
常见(wo'jian'guo'de)的半在线卷积算法实现方法为,传入 $i$,更新到 $f[i]$ 位置时调用 $relax(i)$,但这种方法实现扩展性较差。
于是我们可以考虑把序列真的作为一个流来处理,主要定义如下:
struct ocpoly {
function<mi(int)> get_handle;
vector<int> vis; vector<mi> val;
mi get(int x) {
if((int)vis.size()>x&&vis[x]) return val[x];
if((int)vis.size()<=x)
vis.resize(x+1),val.resize(x+1);
val[x]=get_handle(x); vis[x]=1; return val[x];
}
poly await(int n) {
poly s(-1);
for(int i=0;i<n;++i)
s[i]=get(i);
return s;
}
void copy(ocpoly*b) {
get_handle=[=](int c) {
return b->get(c);
};
}
ocpoly() {}
ocpoly(ocpoly const&)=delete;
ocpoly& operator = (ocpoly const&)=delete;
}
即我们对外的接口为 get
函数,该函数实现了对lambda get_handle
返回值的记忆化(我们希望 $\text{get_handle}(i)$ 返回 $x^i$ 位置上的系数)。这种写法的好处在于扩展性较强,我们只需要改变lambda就可以派生出各种流。例如我们可以定义在一个流上进行积分:
struct ocint: public ocpoly {
ocint(ocpoly*p,mi c=0) {
get_handle=[=](int u) {
if(u==0) return c;
return p->get(u-1)*rfac[u]*fac[u-1];
};
}
}
我们可以定义更多的组合子,例如如果我们定义 $\text{ocde}$ 为流上的求导,$\text{ocquo}(A,B)$ 为 $A(x)B(x)^{-1}$ 形成的形式幂级数,我们可以组合成流上的 $\ln$:
struct ocln: public ocpoly {
ocln(ocpoly*a) {
copy(new ocint(new ocquo(new ocde(a),a)));
}
};
(在内存管理的角度上这种写法不是很优美,但是懒得用智能指针)
当然,实现最困难的部分在于全在线多项式卷积 $\text{ocmul}$。首先全在线转半在线只需要倍增最大长度,即我们对于当前考虑的位置 $i$,我们保证所有 $i=a+b$ 有 $\min(a,b)\le u$,即要求 $2(u+1)>i$,并在不满足要求时修改 $u$。这样我们就可以将在线卷积转换为 $A[0..u],B[u+1..]$ 和 $B[0..u],A[u+1..]$ 的半在线卷积。
为了实现非递归的 $4$-叉 半在线卷积,我们进行形如二项式分组的过程(也可以理解为记录递归栈),对于每一组我们记录 DFT 后的点值,每次维护分组并计算贡献。在小范围上我们进行循环展开的暴力。
这里DFT的总长是 $O(n\log n)$ 级别的,因此需要进行内存回收。这里的内存管理需要格外小心。
利用 ocmul
全在线地写完 $\exp$, 根号, 求逆, 快速幂... 大功告成!
举个例子,实现快速幂 $F^k$:注意到 $(F^k)'F=kF^kF'$,在线地计算 $\int kF^kF'F^{-1}$ 即可($F^k$ 即本身)。
如果你和我一样认为形如 copy(new ocint(new ocquo(new ocde(
的一大串不好看,我们可以使用一个包裹类包起来:
struct pipeline {
ocpoly*p;
pipeline() {p=new ocpoly();}
pipeline(ocpoly *q) {assert(q);p=q;}
pipeline(poly s) {p=new ocfixed(s);}
mi get(int n) {return p->get(n);}
poly await(int n) {return p->await(n);}
void set(pipeline q) {p->copy(q.p);}
};
pipeline operator + (pipeline a,pipeline b)
{return pipeline(new ocadd(a.p,b.p));}
这样就可以把这个包裹类 pipeline
当做对象愉快地进行操作了!我们甚至可以保留与多项式部分相同的接口。
loj150 挑战多项式
没错,我们可以全在线地实现本题!代码也长得差不多。
int n,m,f[123456];
int main()
{
scanf("%d%d",&n,&m); ++n;
for(int i=0;i<n;++i) scanf("%d",f+i);
pipeline p(vector<int>(f,f+n));
p=p-gexp(gint(ginv(gsqrt(p,mod_sqrt(f[0],MOD)))));
p=gcorner(gln(gcorner(p,{1})),{1});
gde(gexp(gln(p)*m)).await(n).print_len(n-1,'\n');
}
跑了 sum 8958ms/max 478ms,也就慢了一倍左右。
uoj50 链式法则
$F(x) \equiv \int\left(\frac{1}{2}A(x)F^2(x) + 1 \right)dx$,直接照抄式子!
int n; poly a(-1);
char str[288888];
int main()
{
scanf("%d%s",&n,str);
a.reserve(n+1);
for(int i=0;i<n;++i)
if(str[i]=='1') a[i]=rfac[i];
pipeline o;
o.set(gint(a*o*o/2+1_p)); //!
poly g=o.await(n+1);
for(int i=1;i<=n;++i) {
mi ans=g.get(i)*fac[i];
printf("%d\n",(int)ans);
}
}
主要就一行代码,sum 1042ms/max 432ms,(在写作本文时)喜提ranklist第三。
loj6684 有根无标号「奇树」计数
$F(x) \equiv x\text{MSet}(\text{MSet}(F(x))-1)$,直接照抄。
int n;
int main() {
scanf("%d",&n);
pipeline p;
p.set(gshl(MSet(MSet(p)-1_p),1));
poly o=p.await(n+1);
for(int i=1;i<=n;++i)
printf("%d\n",(int)o.get(i));
}
同样是一行代码,跑了 max 200ms,在写作本文时喜提ranklist第一(忽略一位打表的同学)。
loj6538 烷基计数 加强版 加强版
也是照抄一行式子。
int n;
int main() {
scanf("%d",&n);
pipeline T;
T.set(gcorner(gshl((gamp(T,3)*2+T*(gamp(T,2)*3+T*T))/6,1),{1}));
poly p=T.await(n+1);
printf("%d\n",p.get(n));
}
跑了 sum 582ms/max 167ms,在写作本文时(可以)排在ranklist第六页(虽然这次有点惨但好歹能过)。
我们也可以方便地实现更为复杂的全在线过程,例如 $F\equiv GH(1+x)+3x^2$,$G\equiv FH(1+x^2)+5x$,$H\equiv FG(1-x^3)+3x^4$。
int main() {
pipeline f,g,h;
f.set(gcorner(g*h*"1+x"_p+"3x^2"_p,{0}));
g.set(gcorner(f*h*"1+x^2"_p+"5x"_p,{0}));
h.set(gcorner(f*g*"1-x^3"_p+"3x^4"_p,{0}));
poly F=f.await(100000);
poly G=g.await(100000);
poly H=h.await(100000);
F.dump('\n');G.dump('\n');H.dump('\n');
}
可以看到全在线方法写起来比较优美,适用性广泛,效率也比较快,是我们熟悉的迭代方法的有力补充。
最后我们给出以上模板的实现,如果有更好的实现 / 更好的函数类型名 / 发现bug(这个比较重要)欢迎在评论区讨论。由于代码较长,建议自行截取需要的部分进行使用。
当然,发这份代码的本意并不是希望大家做题的时候上来就贴这份板子,主要还是希望给大家提供写板子的思路,以及普及一下在线卷积的作用?
代码更新:2021.08.07 (paste, gist) 修复了一处内存泄漏;修改了一些代码风格;修改了幂和的实现;换了个快点的ntt、细节上的常数优化(减少bitrev的调用、优化ocmul的小范围暴力)、加入了半在线卷积;加入了多项式欧几里得算法和有理分式重建