写这玩意儿的感想就是:一开始常数比你小了一点的,到后来随着各种函数调用来调用去,会变成亿点……
关于界
在大多数题目中,我们只对多项式前若干项感兴趣,此时为所有运算设定一个次数上界,即 (mod x^n) 。不难发现有性质:
这样能避免对无穷幂级数的处理,保证复杂度。
多项式四则运算
约定 :([x^n]F(x)) 和 (F[n]) 均表示 (F(x)) 的 (n) 次项系数。
定义多项式的加减为:
多项式乘法(加法卷积)为:
使用 NTT 在 (mathcal{O}(nlog n)) 内完成计算。FFT & NTT
多项式乘法逆
如果多项式 (F(x),G(x)) 满足 (F(x)*G(x)=1) ,那么称 (F(x),G(x)) 互为乘法逆。
这里采用和卷积同样的倍增思想。假设现在已经求出了 (H(x)) 满足 (F(x)H(x)equiv 1(mod x^{lceil n/2 ceil})) ,现在需要通过 (H(x)) 求 (G(x)) .
两边平方,得(注意,平方的时候模数也要同步)
现在用 (F(x)) 消去平方项 (G(x))
移项得
这样就能直接用 NTT 递推计算了。模板题链接
如果你发现人均 500ms
,就你 1s+
,不要慌张,开个 O2
试试 /xyx
,我直接 1.22s
( o) 490ms
.
没想到老厌氧选手也有这一天 但是还是不喜 C++17
QAQ 一开准慢。
注:这里只有当前段代码,也就是针对本题删去了无用部分和封装,但事实上完整代码是……大缝合怪!
#define Clear(a,n) memset(a,0,sizeof(int)*(n))
#define Copy(a,b,n) memcpy(a,b,sizeof(int)*(n))
#define Rev(a,b) reverse(a,b)
int power( int a,int b ) { int res=1; for (;b;b>>=1,a=1ll*a*a%Mod) if (b&1) res=1ll*a*res%Mod; return res; }
inline void bmod( int &x ) { x-=Mod; x+=x>>31&Mod; }
int rev[M],rt[M],lim,lg,lasn; //反转数组,原根,上界的值和log
void Poly_Init( int n )
{
if ( n==lasn ) return; lasn=n;
for ( lg=0,lim=1; lim<=n; lg++,lim<<=1 );
for ( int i=0; i<lim; i++ ) rev[i]=(rev[i>>1]>>1) | ((i&1)<<(lg-1));
for ( int i=1,w; i<lim; i<<=1 )
{
w=power(3,(Mod-1)/(i<<1)); rt[i]=1;
for ( int j=1; j<i; j++ ) rt[i+j]=1ll*rt[i+j-1]*w%Mod;
}
}
void NTT( int *f,int opt )
{
int i,invn,len;
for ( i=0; i<lim; i++ ) if ( i>rev[i] ) swap( f[i],f[rev[i]] );
for ( i=1; i<lim; i<<=1 )
for ( int j=0; j<lim; j+=i<<1 )
for ( int k=0,t1,t2; k<i; k++ )
{
t1=f[j+k]; t2=1ll*rt[i+k]*f[i+j+k]%Mod;
bmod(f[j+k]=t1+t2); bmod(f[i+j+k]=t1+Mod-t2);
}
if ( opt ) return; invn=power(lim,Mod-2);
for ( i=0; i<lim; i++ ) f[i]=1ll*f[i]*invn%Mod;
}
void Poly_Mul( int n,int m,int *f,int *g,int *res )
{
Poly_Init(n+m); static int tf[M],tg[M];
Copy(tf,f,n); Clear(tf+n,lim-n); NTT(tf,1);
Copy(tg,g,m); Clear(tg+m,lim-m); NTT(tg,1);
for ( int i=0; i<lim; i++ ) res[i]=1ll*tf[i]*tg[i]%Mod;
Rev(res+1,res+lim); NTT(res,0);
}
void Poly_Inv( int n,int *f,int *g )
{ //G(x)equiv 2H(x)-H(x)^2F(x)(mod x^n)
static int tf[M];
if ( n==1 ) { g[0]=power(f[0],Mod-2); return; }
Poly_Inv((n+1)>>1,f,g); Poly_Init(n<<1);
Copy(tf,f,n); Clear(tf+n,lim-n); Clear(g+n,lim-n); NTT(tf,1); NTT(g,1);
for ( int i=0; i<lim; i++ ) g[i]=1ll*g[i]*(2+Mod-1ll*g[i]*tf[i]%Mod)%Mod;
Rev(g+1,g+lim); NTT(g,0); Clear(g+n,lim-n);
}
多项式除法
给定 (F(x),G(x)) ,求 (Q(x)) 使得 (F(x)=G(x)Q(x)+R(x)) ,(Q(x)) 次数为 (n-m) ,(R(x)) 次数小于 (m) .
令 (F^R(x)) 表示 (F(x)) 系数翻转之后的结果。设 (F(x)) 最高次项为 (x^{n-1}) ,那么有
于是开始推式子
由于 (Q) 的次数是 (n-m) ,所以可以进行模意义下的转化
于是就能求出 (Q^R(x)=dfrac{F^R(x)}{G^R(x)}pmod{x^{n-m+1}}) ,这个用求逆就能完成,求 (R(x)) 就直接 (F-G*Q) 即可。
void Poly_DivMod( int n,int m,int *f,int *g,int *q,int *r )
{ //F^R(x)equiv G^R(x)Q^R(x)pmod {x^{n-m+1}}
static int tf[M],tg[M],nw[M];
Copy(tf,f,n); Copy(tg,g,m); Rev(tf,tf+n); Rev(tg,tg+m);
for ( int i=n-m+1; i<(n<<1); i++ ) tg[i]=0;
Poly_Inv(n-m+1,tg,nw); Poly_Mul(n,n-m+1,tf,nw,q);
for ( int i=n-m+1; i<(n<<1); i++ ) q[i]=0;
Rev(q,q+n-m+1); Poly_Mul(n-m+1,m,q,g,r);
}
拉格朗日插值
给定 (n) 个点值,确定多项式 (f(x)) 并求出 (f(k)mod 998244353) .
拉格朗日他老人家留下的美妙遗产之一:
详细推导参见 OI-Wiki . 这样直接做就达到了 (mathcal{O}(n^2)) 比消元少了一个 (n) 呢 。在 (x) 连续的时候,可以把 (x_i) 替换成 (i) ,并令 (pre[i]=prodlimits_{j=0}^ik-j) ,(suf[i]=prodlimits_{j=i}^nk-j) ,(fac[i]=prodlimits_{j=1}^ij) ,那么有
优化到了 (mathcal{O}(n)) .
重心拉格朗日插值
回到这个式子:
令 (g=prodlimits_{i=0}^nk-x[i]) ,原式可以变为
令 (t[i]=dfrac{y[i]}{prodlimits_{i eq j}x_i-x_j}) ,有
每次修改或插入直接修改 (t[i]) 即可。注意不能处理 (k=x[i]) 的情况。模板题链接
暴力 (mathcal{O}(n^2)) 的插值 82ms
跑进第一页也是没想到啊……
//Author: RingweEH
const int N=2e3+10,Mod=998244353;
int n,k,x[N],y[N];
int power( int a,int b ) { int res=1; for (;b;b>>=1,a=1ll*a*a%Mod) if (b&1) res=1ll*a*res%Mod; return res; }
int main()
{
n=read(); k=read();
for ( int i=1; i<=n; i++ ) x[i]=read(),y[i]=read();
ll ans=0;
for ( int i=1; i<=n; i++ )
{
ll a=1,b=1;
for ( int j=1; j<=n; j++ )
{
if ( i==j ) continue;
a=a*(k-x[j]+Mod)%Mod; b=b*(x[i]-x[j]+Mod)%Mod;
}
ans=(ans+a*power(b,Mod-2)%Mod*y[i]%Mod)%Mod;
}
printf( "%lld
",ans );
return 0;
}
泰勒展开 & 牛顿迭代
原理:用一个 (n) 次多项式高度拟合一个函数,新的函数某一点 (x_0) 的 (n) 阶导数和原函数的 (n) 阶导数相同,并且都经过 ((x_0,f(x_0))) .
具体可以参考《具体数学》5.3
的 技巧2:高阶差分
关于牛顿级数和泰勒级数的阐述。
函数 (F(x)) 在 (x_0) 位置泰勒展开,令 (F^n(x)) 表示 (F(x)) 的 (n) 阶导数,有
得到这个式子的推导可以类比牛顿级数的推导,先得到 (F^n(x_0)) 为当前项的系数且在更大的时候均为 (0) (根据求导易证),然后再根据求导得到的系数部分得到下面的分子。
往往在精度满足的时候就停止,所以大部分都是有误差的。
考虑多项式求逆的过程,假设已经有 (F(x)*H(x)equiv 1(mod x^{lceil n/2 ceil})) ,函数 (G(F(x))equiv 0(mod x^n)) .
(这里显然有 (G(H(x))equiv 0pmod x^{lceil n/2 ceil}) )
将 (G(F(x))) 在 (H(x)) 泰勒展开
从第三项开始,由于 (F(x)) 和 (H(x)) 后 (leftlceildfrac{n}{2} ight ceil) 项相同,所以模意义下为 (0) .
所以有
这样就能求得 (F(x))
多项式 ln
如果不太明白为什么多项式会有自然对数这种东西,可以看 这里 (但是貌似还是需要一些高级的东西……先感性理解吧)
对 (B(x)=ln(A(x))) 两边求导
然后对 (A) 求逆,求导,卷起来,然后积分,就得到了 (B) . 题目链接
//Author: RingweEH
void Itg( int n,int *f,int *g ) //Intergral
{
for ( int i=1; i<=n; i++ ) g[i]=1ll*f[i-1]*Math::inv[i]%Mod; g[0]=0;
}
void Drv( int n,int *f,int *g ) //Derivative
{
for ( int i=0; i<n-1; i++ ) g[i]=1ll*(i+1)*f[i+1]%Mod; g[n-1]=0;
}
void Poly_ln( int n,int *f,int *g )
{
static int tf[M],tg[M];
Drv(n,f,tf); Clear(tg,n); Poly_Inv(n,f,tg);
Poly_Mul(n,n,tf,tg,tf); Itg(n,tf,g);
}
多项式 exp
已知 (A(x)) ,求 (F(x)=e^{A(x)}) ,保证 (A(0)=0) .
对两边求导得 (ln(F(x))=A(x)) . 令 (G(F(x))=ln(F(x))-A(x)) ,那么有 (G'(F(x))=dfrac{1}{F(x)}) .
由于 (A(x)) 给出,在 (x) 固定时是固定的,可以看做常量。
然后对其牛顿迭代
同样递归计算 (H(x)) 即可,边界为 (G(0)=1) . 这样的复杂度是 (mathcal{O}(nlog n)) 的。
如果 (n) 很小且模数不一定是 NTT 模数,那么可以换一种做法。
所以有
积分回去
void Poly_Exp( int n,int *f,int *g )
{ //G(x)=H(x)(1-ln(H(x))+A(x))
if ( n==1 ) { g[0]=1; return; } static int tf[M];
Poly_Exp( (n+1)>>1,f,g ); Clear(tf,n); Poly_ln(n,g,tf); tf[0]--;
for ( int i=0; i<n; i++ ) bmod( tf[i]=f[i]+Mod-tf[i] );
Poly_Mul(n,n,g,tf,g); Clear(g+n,lim-n);
}
多项式开根
给定 (A(x)) ,求 (F^2(x)equiv A(x)pmod x^n) ,保证 (a_0=1) .
设 (G(F(x))=F^2(x)-A(x)) ,那么 (G'(F(x))=2F(x)) . 对其牛顿迭代
同样递归求解即可。(G(0)=1) . 模板题链接
void Poly_Sqrt( int powe,int *f,int *g )
{//F(x)=dfrac{H(x)^2+A(x)}{2H(x)}
if ( powe==1 ) { g[0]=1; return; }
static int tf[M]; Poly_Sqrt((powe+1)>>1,f,g);
Clear(tf,powe); Poly_Inv(powe,g,tf); Poly_Mul(powe,powe,tf,f,tf);
int inv2=power(2,Mod-2); for ( int i=0; i<powe; i++ ) g[i]=1ll*(g[i]+tf[i])*inv2%Mod;
}
多项式快速幂
经典用法:(ln+exp) 把幂次转化成四则运算。
然后 (exp(ln(B(x)))) 即可。
注意:本题 (kleq 10^{10^5}) ,所以要用到:
任意质数 (p) 和任意多项式 (F(x)) , (F^p(x)equiv 1pmod{x^n}) (还得要求 (a_0=1) ),具体可以看具体数学和 这个帖子
总之就是读入的幂次要 (mod p) 啦。
void Poly_Power( int n,int k,int *f,int *g )
{
int tf[M]; Clear(tf,n); Poly_ln(n,f,tf);
for ( int i=0; i<n; i++ ) tf[i]=1ll*tf[i]*k%Mod;
Poly_Exp(n,tf,g);
}
一些技巧
差卷积
令 (B_1[n-i]=B[i]) ,那么有
所以 (C[i]=(A*B_1)_{n+i}) .
多项式平移
用 (F(x)) 得到 (F(x+c)) 的系数。设 (F(x)=sum_{i=0}^na_ix^i) .
差卷积即可。
Hints
- 如果你想对
Poly_Mul
卡常的话,可以试试小范围暴力大范围卷积。 - 不妨对手写取模加个
inline
,效果会很不错。 - 在 这里 F3 “卡常后的倍增”会得到一个
Poly_Inv
的小优化。 - 似乎
i+j+k
比i|j|k
快。 memset
应该比fill
和循环快- 指针很强
虽然我不会,如果会指令集那我只能 Orz - 一个避免大量
memset
的好方法是每次在封装函数里面开static
,而且能有效避免传指针修改外界变量的情况
习题
- 不要再忘记
Math::Init(n)
了! - 时刻记住你的卷积开闭区间情况。
力
前面那个式子可以卷积,后面那个式子反一反同样卷积,就没了。不会吧不会吧你都看过多项式除法了不会想不到 reverse 吧
礼物
设当前调整值为 (c) ,(x_i=a_i-b_i) ,先来考虑一个简化的问题:已知 (b) 序列,如何确定 (c) 使 (sum(x_i+c)^2) 最小?
整个式子可以看做是一个关于 (c) 的二次函数,那么直接枚举 (c) 就能 (mathcal{O}(m)) 求解了。
现在来考虑不固定的情况。注意到 (sum x_i) 是不会变的,那么式子依然可以拆分成两部分,一部分是最小化 (nc^2+2cX) (依然可以直接做),一部分是最大化 (sum a_ib_i) 。将 (a) 反向,构造卷积 (sum a_ib_{n-i}) ,由于 (b) 是环形,所以倍长之后取后 (n) 个中的最大值即可。复杂度是 (mathcal{O}(nlog n+m)) .
差分与前缀和
先考虑前缀和,发现其实就是和一个系数全 (1) 的多项式做一次卷积。那么 (k) 次前缀和,根据卷积的结合律,乘上去的多项式就是 (dfrac{1}{(1-x)^k}) 。
类似思考差分,一维差分要卷积的式子是 (1-x) ,那么 (k) 次就是 ((1-x)^k) ,直接多项式快速幂。
这是比较暴力的做法……所以不开 O2
就全 TLE
,否则就会 AC
。
开拓者的卓识
再次说明了“ 算贡献 ”的思想是真的很重要。
我们考虑对于一个原始序列的 (a_i) ,会对哪些 ([k,l,r]) 产生贡献。注意,这里不仅仅要选出最后的 (l,r) ,还要考虑之前的所有的 (k-1) 个区间,所以相当于求出满足 (1leq l_kleq l_{k-1}leq l_{k-2}leq cdots leq ileq cdots leq r_{k-1}leq r_k) 的 ({(l_1,r_1),(l_2,r_2),cdots ,(l_k,r_k)})集合个数,由于左右端点独立,那么可以直接插板得到
注意到组合数的下指标都是 (k-1) ,所以可以一遍递推得到,令 (displaystyle b_i=inom{k-1+i}{k-1}) ,答案就是
显然可以卷积。
学习材料
tommy0221
:多项式笔记(一)
zhouzhendong
:多项式入门教程