数学:多项式
参考文章:
博主的理解并不算深刻,有不恰当之处欢迎指出。
拉格朗日插值
通常用于解决如下问题:
- 已知一个
次函数 图像上的 个点 ,试求出这个函数 在 处的取值。
一般拉格朗日插值法
现在让我们来用一些比较 MO 的方式解决这个问题。考虑能不能用一些奇技淫巧构造这么一个函数。
考虑一个跟之前 CRT 比较类似的思路。对每个点
直接给出这个答案多项式的形式
直接实现是
upd:拉插可以
取值为连续正整数时
钦定
那么原式变为
预处理
重心拉格朗日插值法
这是什么,没用到过。
应用
求 (CF622F The Sum of the k-th Powers)
求
。
。
根据一些懂的都懂不懂的不懂的数学知识,这个玩意是一个关于
代入
ll n,k,xx[N],yy[N];
ll fac[N],fnv[N],f[N],b[N];
ll qp(ll x,ll y) { return y?(y&1?x:1)*qp(x*x%P,y>>1)%P:1; }
ll lag(ll n,ll *x,ll *y,ll k) {
fac[0]=1; for (ll i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;
fnv[n]=qp(fac[n],P-2); for (ll i=n-1;~i;i--) fnv[i]=fnv[i+1]*(i+1)%P;
ll ret=0;
f[0]=b[n+2]=1;
for (ll i=1;i<=n+1;i++) f[i]=f[i-1]*(k-x[i])%P;
for (ll i=n+1;i;i--) b[i]=b[i+1]*(k-x[i])%P;
for (ll i=1;i<=n+1;i++) (ret+=y[i]*(f[i-1]*b[i+1]%P)%P*(fnv[i-1]*fnv[n+1-i]*(n+1-i&1?-1:1)%P))%=P;
return (ret+P)%P;
}
void mian() {
scanf("%lld%lld",&n,&k);
for (ll i=1;i<=k+2;i++) xx[i]=i,yy[i]=(yy[i-1]+qp(i,k))%P;
cout<<lag(k+1,xx,yy,n);
}
求只有一个自变量的次数已知的多项式的取值(P3270 [JLOI2016] 成绩比较)
, 。
跳过一堆步骤,总之就是中间有一步要求出每个
你惊喜地发现全程只有
快速傅里叶变换(FFT)
已经不知道被这玩意劝退了多少次了。但理解之后其实不算很阴间的东西?
前置知识
多项式
我们平常使用的什么
还有一种叫点值表示法的东西。因为一个
点值表示法看上去非常反人类,但是你用系数表示法的话就反 FFT 了。
向量、三角函数、复数
每一个第一次学的人:?多项式是怎么跟复数扯上关系的?
简单来说就是实数域没有满足 FFT 需要的性质的数,所以就跑到复数域来乱搞了,因为多项式
假设你有高中数学中向量、三角函数以及此方面的相关知识。大概就是向量的定义、相加;三角函数的定义;复数可表示为
单位根
我们以复平面的原点画一个半径为
把三角函数和向量相加的技术放到复平面上,有
快速傅里叶变换(FFT)
首先我们把听了一万年的多项式乘法之类的东西暂时忘掉,不要管它们。现在我们研究的是傅里叶变换,不是多项式乘法。在 OI 里你可以把它简单地理解为就是把一个多项式的系数表示法转为点值表示法。FFT 就是快速地实现这个过程。好像在别的领域的应用没有这么简单,但我们暂时用不到。
假设现在我们要求出一个
考虑将多项式的奇偶项分开。
记两个新的多项式
那么有
如果我们要求出一个点值,我们需要带入
但是如果我们把眼光放到复数域,事情就有点意思了。
对于每一个
再代入
你发现,这时如果你算出了
那我们考虑,令我们要求的
在求的过程中,我们枚举
那么整体来看,求出
你发现其实这两部分跟现在求解的问题没有任何区别。对于前一部分,只要令
一直这样对半劈下去,没过多久就会进入
总之就是利用单位根的性质实现分治求解。非常绝妙!
快速傅里叶逆变换(IFFT)
因为点值表示法实在太过反人类,我们需要一个把多项式的点值表示法转成系数表示法的东西,即 IFFT。因为实现方式实在太过相似,通常不跟 FFT 作区分。
假设我们已知多项式
也就是我们把点值作为多项式
证明有点离谱,晚点补(
代码实现
因为涉及到复数运算,建议写一个复数类,或者直接用 STL 的 complex。
这就不可避免地给 FFT 带来了一定的常数和精度问题。
递归实现
即直接模拟上面那个东西。FFT 的过程本质上就是一个递归分治,所以比较容易按照这个思路实现。
但是递归实现的常数放到 FFT 上简直是雪上加霜,所以我们考虑把递归去掉。
迭代实现
也就是说我们需要一种方法直接知道
先跑一遍?可以,但不太优雅,也有点慢。
盗你谷第一篇题解的图。
你发现把原序列的二进制位翻转之后就得到了新序列的位置。其实可以理解,因为我们在第
直接给出递推式:设 r
是新序列每一项在原序列中的位置,那么 r[i]=(r[i>>1]>>1)|((i&1)<<(m-1))
。m
是层数。理解很简单,就是 r[i]=r[i>>1]*2+(i%2)
,用位运算技巧变成 r[i]=(r[i>>1]<<1)|(i&1)
,然后把二进制位反转了一下。
底层常数项搞定了。而对于从最底层往最上层处理的过程,回忆一下之前的式子
注意如果是 IFFT 的话最后最好四舍五入一下,FFT 的精度问题还是不可忽视的。
下面代码中,com
是手写的复数类,ini
是求出第一个大于 _n
的 mid
是此时子问题的规模的一半,w
是单位根,x
是单位根的幂,c
用于区分 FFT 和 IFFT,
struct com {
double a,b;
com(double _a=0,double _b=0) { a=_a,b=_b; }
friend com operator + (com x,com y) { return {x.a+y.a,x.b+y.b}; }
friend com operator - (com x,com y) { return {x.a-y.a,x.b-y.b}; }
friend com operator * (com x,com y) { return {x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a}; }
} a[N],b[N];
ll n=1,m=0,r[N];
void ini(ll _n) {
while (n<=_n) n<<=1,m++;
for (ll i=0;i<n;i++) r[i]=r[i>>1]>>1|(i&1)<<m-1;
}
void fft(com *a,ll c) {
for (ll i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
for (ll mid=1;mid<n;mid<<=1) {
com w(cos(pi/mid),c*sin(pi/mid));
for (int i=0;i<n;i+=mid<<1) {
com x(1,0);
for (int j=0;j<mid;j++,x=x*w) {
com y=a[i+j],z=x*a[i+mid+j];
a[i+j]=y+z,a[i+mid+j]=y-z;
}
}
}
if (c==-1) for (ll i=0;i<n;i++) a[i].x=(ll)(a[i].x/n+0.5);
}
快速数论变换(NTT/FNTT)
明显 FFT 常数和精度问题都很大。
众所周知精度的最高境界是取模,所以 NTT 来了。即在模意义下快速实现系数表示法和点值表示法之间的转换。
其实准确来说这玩意应该叫 FNTT,NTT 是数论变换,但是没人在意。
前置知识:原根
若一个数
当
快速数论变换(NTT)
原根这玩意逆天的地方在于,在模意义下,FFT 时用到的单位根的性质(
于是,若
也就是说,NTT 就是把 FFT 的单位根换成上面的原根幂。从此我们就非常轻松地学会了在取模的条件下飞快地进行系数和点值表示之间的转换!
……吗?注意这个指数。因为
那就不办了,考虑
那就没有多少个质数模数满足条件了。一个著名模数是
当然人类的智慧是无穷的,确实有一个叫作任意模数 NTT 的东西,在下面。
实现
跟 FFT 差不多。
ll n=1,m=0,a[N],r[N];
void ini(ll _n) {
while (n<=_n) n<<=1,m++;
for (ll i=0;i<n;i++) r[i]=r[i>>1]>>1|(i&1)<<m-1;
}
void ntt(ll c) {
for (ll i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
for (ll o=1;o<n;o<<=1) {
ll w=qp(c==1?3:qp(3),(P-1)/(o<<1));
for (ll j=0,x,y;j<n;j+=o<<1) for (ll k=j,v=1;k<j+o;k++,v=v*w%P)
x=a[k],y=v*a[k|o]%P,a[k]=(x+y)%P,a[k|o]=(x+P-y)%P;
}
if (c==-1) for (ll i=0,iv=qp(n);i<n;i++) a[i]=a[i]*iv%P;
}
任意模数 NTT
见下方应用。
应用
多项式相关操作
乘法:P3803 【模板】多项式乘法(FFT)
给定一个
次多项式 ,和一个 次多项式 。请求出 和 的卷积。
,输入系数 ,2s,500MB。
点值表示法简直是浑然天成地适合用来做卷积,不难知道,拿足够的点数,然后把两个多项式在这些点上的点值乘起来即是卷积后的对应点值。
所以我们先令 FFT 出的点值个数为大于
因为系数实在太小,这题把
代码实现在后面。
任意模数乘法:P4245 任意模数多项式乘法
给定两个多项式
,求 ,系数对 取模。不保证 是 NTT 模数。
, ,2s,500MB。
注意到卷积后得到的系数在不取模的结果下最大值也就
这首先要求我们记住三个 NTT 模数。一般使用
最后就是 CRT 的技术了。这是简单的。同样是单次
求逆:P4238 【模板】多项式乘法逆
给定一个
次多项式 ,求一个多项式 使得 ,系数对 取模。
,1s,125MB。
也是一个不知道前人是怎么想出来的东西。同样令
设
那么有
然后骚操作来了。两边平方得
两边同乘
那么从
注意不要在算对
:P4725 【模板】多项式对数函数(多项式 ln)
给定一个
次多项式 ,求一个多项式 使得 ,系数对 取模。
,2s,125MB。
:P4726 【模板】多项式指数函数(多项式 exp)
给定一个
次多项式 ,求一个多项式 使得 ,系数对 取模。
,2s,125MB。
快速幂:P5245 【模板】多项式快速幂
给定一个
次多项式 和 ,求一个多项式 使得 ,系数对 取模。
, ,1.5s,250MB。
多项式封装
一种来自 czk 的比较符合直觉而且兼容性还可以的封装写法,重载了括号,跟用普通数组没什么区别,常数还行。
目前只封了 NTT 和加减乘。
namespace poly {
vector<ll> r;
struct ply {
ll n; vector<ll> a;
void sz(ll x=0) {
if (x) n=x;
else while (n>1&&!a[n-1]) n--;
a.resize(n);
}
ll& operator [] (ll x) {
if (x>=n) sz(x+1);
return a[x];
}
ply(ll x=0,ll y=0) { sz(x+1),a[x]=y; }
void ntt(ll c) {
for (ll i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
for (ll o=1;o<n;o<<=1) {
ll w=qp(c==1?3:qp(3),(P-1)/(o<<1));
for (ll j=0,x,y;j<n;j+=o<<1) for (ll k=j,v=1;k<j+o;k++,v=v*w%P)
x=a[k],y=v*a[k|o]%P,a[k]=(x+y)%P,a[k|o]=(x+P-y)%P;
}
if (c==-1) for (ll i=0,iv=qp(n);i<n;i++) a[i]=a[i]*iv%P;
}
friend ply operator * (ply a,ply b) {
ply ret; ll t=1,c=0;
while (t<a.n+b.n-1) t<<=1,c++;
a.sz(t),b.sz(t),ret.sz(t),r.resize(t);
for (ll i=0;i<t;i++) r[i]=r[i>>1]>>1|(i&1)<<c-1;
a.ntt(1),b.ntt(1);
for (ll i=0;i<t;i++) ret[i]=a[i]*b[i]%P;
ret.ntt(-1);
return ret;
}
friend ply operator + (ply a,ply b) {
ply ret(max(a.n-1,b.n-1));
for (ll i=0;i<ret.n;i++) ret[i]=(a[i]+b[i])%P;
return ret;
}
ply operator - () { ply ret; for (ll i=0;i<n;i++) ret[i]=P-a[i]; return ret; }
friend ply operator - (ply a,ply b) { return a+(-b); }
} a,b,c;
} using namespace poly;
ll n,m;
void mian() {
scanf("%lld%lld",&n,&m);
for (ll i=0;i<=n;i++) scanf("%lld",&a[i]);
for (ll i=0;i<=m;i++) scanf("%lld",&b[i]);
c=a*b;
for (ll i=0;i<=n+m;i++) cout<<c[i]<<" ";
}
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探