多项式基础学习笔记(1)
写这玩意儿的感想就是：一开始常数比你小了一点的，到后来随着各种函数调用来调用去，会变成亿点……
关于界
在大多数题目中，我们只对多项式前若干项感兴趣，此时为所有运算设定一个次数上界，即 \(\bmod x^n\) 。不难发现有性质：
这样能避免对无穷幂级数的处理，保证复杂度。
多项式四则运算
约定 ：\([x^n]F(x)\) 和 \(F[n]\) 均表示 \(F(x)\) 的 \(n\) 次项系数。
定义多项式的加减为：
多项式乘法（加法卷积）为：
使用 NTT 在 \(\mathcal{O}(n\log n)\) 内完成计算。FFT & NTT
多项式乘法逆
如果多项式 \(F(x),G(x)\) 满足 \(F(x)*G(x)=1\) ，那么称 \(F(x),G(x)\) 互为乘法逆。
这里采用和卷积同样的倍增思想。假设现在已经求出了 \(H(x)\) 满足 \(F(x)H(x)\equiv 1(\bmod x^{\lceil n/2\rceil})\) ，现在需要通过 \(H(x)\) 求 \(G(x)\) .
两边平方，得（注意，平方的时候模数也要同步）
现在用 \(F(x)\) 消去平方项 \(G(x)\)
移项得
这样就能直接用 NTT 递推计算了。模板题链接
如果你发现人均 500ms ，就你 1s+ ，不要慌张，开个 O2 试试 /xyx ，我直接 1.22s \(\to\) 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)(\bmod 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)\bmod 998244353\) .
拉格朗日他老人家留下的美妙遗产之一：
详细推导参见 OI-Wiki . 这样直接做就达到了 \(\mathcal{O}(n^2)\) 比消元少了一个 \(n\) 呢 。在 \(x\) 连续的时候，可以把 \(x_i\) 替换成 \(i\) ，并令 \(pre[i]=\prod\limits_{j=0}^ik-j\) ，\(suf[i]=\prod\limits_{j=i}^nk-j\) ，\(fac[i]=\prod\limits_{j=1}^ij\) ，那么有
优化到了 \(\mathcal{O}(n)\) .
重心拉格朗日插值
回到这个式子：
令 \(g=\prod\limits_{i=0}^nk-x[i]\) ，原式可以变为
令 \(t[i]=\dfrac{y[i]}{\prod\limits_{i\neq 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\n",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(\bmod x^{\lceil n/2\rceil})\) ，函数 \(G(F(x))\equiv 0(\bmod x^n)\) .
（这里显然有 \(G(H(x))\equiv 0\pmod x^{\lceil n/2\rceil}\) ）
将 \(G(F(x))\) 在 \(H(x)\) 泰勒展开
从第三项开始，由于 \(F(x)\) 和 \(H(x)\) 后 \(\left\lceil\dfrac{n}{2}\right\rceil\) 项相同，所以模意义下为 \(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}(n\log 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)))\) 即可。
注意：本题 \(k\leq 10^{10^5}\) ，所以要用到：
任意质数 \(p\) 和任意多项式 \(F(x)\) ， \(F^p(x)\equiv 1\pmod{x^n}\) （还得要求 \(a_0=1\) ），具体可以看具体数学和 这个帖子
总之就是读入的幂次要 \(\bmod 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}(n\log 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\) 个区间，所以相当于求出满足 \(1\leq l_k\leq l_{k-1}\leq l_{k-2}\leq \cdots \leq i\leq \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=\binom{k-1+i}{k-1}\) ，答案就是
显然可以卷积。
学习材料
tommy0221 ：多项式笔记（一）
zhouzhendong ：多项式入门教程
