试试以二级标题为主的格式。
多项式相关
注：本篇博客不包含(FFT)基础姿势。如果您想要阅读本篇博客，请确保自己对(FFT,NTT)有基本的认识并且能够独立写出代码。
多项式是什么？
左转数学七年级上册课本。
多项式的两种表示法
系数表示法和点值表示法，过于基础不多解释。
多项式的四则运算及扩展
多项式加减法
同类项直接进行加减即可，过于简单不多解释。
多项式乘法
朴素算法：(O(n^2))。
分治乘法：没写过，不常用，不知道。
快速傅里叶变换（(FFT)）：(O(nlogn))。
快速数论变换（(NTT)）：(O(nlogn))。
（一般认为(NTT)的常数比(FFT)大，两者的适用情况略有不同。）
多项式求逆
给定一个多项式(F(x))，请求出一个多项式(G(x))，满足(F(x) imes G(x) equiv 1 (mod x^n))。系数对(998244353)取模。
推倒过程
显然(n=1)时，直接算乘法逆元即可。
假设我们已经求出了(H(x))满足(F(x) imes H(x) equiv 1 (mod x^{lceil frac{n}{2} ceil}))，我们现在想求出(G(x))满足(F(x) imes G(x) equiv 1 (mod x^n))。
显然有(F(x) imes G(x) equiv 1 (mod x^{lceil frac{n}{2} ceil}))成立，所以我们可以得到(G(x)-H(x) equiv 0 (mod x^{lceil frac{n}{2} ceil}))。
同余号两边同时平方，并且我们发现模数也可以平方，得：
[(G(x)-H(x))^2 equiv 0 (mod x^n)
]
[G(x)^2-2 imes G(x) imes H(x)+H(x)^2 equiv 0 (mod x^n)
]
同余号两边同时乘(F(x))，得：
[G(x)-2 imes H(x)+H(x)^2 imes F(x) equiv 0 (mod x^n)
]
[G(x) equiv 2 imes H(x)-H(x)^2 imes F(x) (mod x^n)
]
从上面的推倒过程我们可以发现，想要求出(G(x))，要先求出(H(x))。求出(H(x))相当于是这个问题的一个规模缩小的子问题，于是我们可以考虑递归求解。
时间复杂度：(O(nlogn))。
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,a) for(int i=head[(a)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
inline int read(){
int x=0;char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
return x;
}
const int MAXN=100005;
const int MAXLEN=270005;
const LL MOD=998244353;
const LL G=3;
const LL INVG=332748118;
int n,m;
int a[MAXLEN];
int len,rev[MAXLEN];
LL invn,A[MAXLEN],B[MAXLEN];
inline LL qpow(LL x,LL y){
LL ret=1,tt=x%MOD;
while(y){
if(y&1) ret=ret*tt%MOD;
tt=tt*tt%MOD;
y>>=1;
}
return ret;
}
inline void ntt(LL *c,int dft){
rin(i,0,n-1) if(i<rev[i])
std::swap(c[i],c[rev[i]]);
for(int mid=1;mid<n;mid<<=1){
int r=(mid<<1);
LL u=qpow(dft>0?G:INVG,(MOD-1)/r);
for(int l=0;l<n;l+=r){
LL v=1;
for(int i=0;i<mid;i++,v=v*u%MOD){
LL x=c[l+i],y=c[l+mid+i]*v%MOD;
c[l+i]=x+y;
if(c[l+i]>=MOD) c[l+i]-=MOD;
c[l+mid+i]=x-y;
if(c[l+mid+i]<0) c[l+mid+i]+=MOD;
}
}
}
if(dft<0) rin(i,0,n-1)
c[i]=c[i]*invn%MOD;
}
void solve(int Floor){
if(Floor==1){
B[0]=qpow(a[0],MOD-2);
return;
}
solve((Floor+1)>>1);
m=Floor-1+((((Floor+1)>>1)-1)<<1);
len=0;
for(n=1;n<=m;n<<=1) len++;
rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
rin(i,0,n-1) A[i]=(i<Floor?a[i]:0);
invn=qpow(n,MOD-2);
ntt(A,1);
ntt(B,1);
rin(i,0,n-1) B[i]=((2*B[i]-A[i]*B[i]%MOD*B[i])%MOD+MOD)%MOD;
ntt(B,-1);
rin(i,Floor,n-1) B[i]=0;
}
int main(){
n=read();
n--;
rin(i,0,n){
a[i]=read();
}
int nn=n;
solve(n+1);
rin(i,0,nn){
printf("%lld ",B[i]);
}
printf("
");
return 0;
}
多项式除法/多项式取模
给定一个(n)次多项式(F(x))和一个(m)次多项式(G(x))，请求出多项式(Q(x),R(x))，满足以下条件：
-
(Q(x))次数为(n−m)，(R(x))次数小于(m)。
-
(F(x)=Q(x) imes G(x)+R(x))。
所有的运算在模(998244353)意义下进行。
推倒过程
首先我们要知道对于任意一个(n)次多项式(f(x))，将其系数翻转后得到的多项式为(f_R(x)=x^nf(frac{1}{x}))。
回到题目，将等式中的所有多项式的自变量改为(frac{1}{x})，等号左右同时乘(x)，得：
[x^nF(frac{1}{x})=x^{n-m}Q(frac{1}{x}) imes x^mG(frac{1}{x})+x^{n-m+1} imes x^{m-1}R(frac{1}{x})
]
[F_R(x)=Q_R(x) imes G_R(x)+x^{n-m+1} imes R_R(x)
]
[F_R(x) equiv Q_R(x) imes G_R(x) (mod x^{n-m+1})
]
[Q_R(x) equiv F_R(x) imes G_R^{-1}(x) (mod x^{n-m+1})
]
利用多项式求逆可以求出(Q_R(x))，然后根据：
[R(x)=F(x)-Q(x) imes G(x)
]
可以求出(R(x))。
时间复杂度：(O(nlogn))。
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,a) for(int i=head[(a)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
inline int read(){
int x=0;char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
return x;
}
const int MAXN=100005;
const int MAXLEN=270005;
const LL MOD=998244353;
const LL ROOT=3;
const LL INVROOT=332748118;
int n,m;
int len,rev[MAXLEN];
int F[MAXN],G[MAXN];
LL Q[MAXN],R[MAXN];
LL invn;
LL A[MAXLEN],B[MAXLEN];
inline LL qpow(LL x,LL y){
LL ret=1,tt=x%MOD;
while(y){
if(y&1) ret=ret*tt%MOD;
tt=tt*tt%MOD;
y>>=1;
}
return ret;
}
inline void ntt(LL *c,int dft){
rin(i,0,n-1) if(i<rev[i])
std::swap(c[i],c[rev[i]]);
for(int mid=1;mid<n;mid<<=1){
int r=(mid<<1);
LL u=qpow(dft>0?ROOT:INVROOT,(MOD-1)/r);
for(int l=0;l<n;l+=r){
LL v=1;
for(int i=0;i<mid;i++,v=v*u%MOD){
LL x=c[l+i],y=c[l+mid+i]*v%MOD;
c[l+i]=x+y;
if(c[l+i]>=MOD) c[l+i]-=MOD;
c[l+mid+i]=x-y;
if(c[l+mid+i]<0) c[l+mid+i]+=MOD;
}
}
}
if(dft<0) rin(i,0,n-1)
c[i]=c[i]*invn%MOD;
}
void getinv(int Floor){
if(Floor==1){
B[0]=qpow(G[0],MOD-2);
return;
}
getinv((Floor+1)>>1);
m=Floor-1+((((Floor+1)>>1)-1)<<1);
len=0;
for(n=1;n<=m;n<<=1) len++;
rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
rin(i,0,n-1) A[i]=(i<Floor?G[i]:0);
invn=qpow(n,MOD-2);
ntt(A,1);
ntt(B,1);
rin(i,0,n-1) B[i]=((2*B[i]-A[i]*B[i]%MOD*B[i])%MOD+MOD)%MOD;
ntt(B,-1);
rin(i,Floor,n-1) B[i]=0;
}
int main(){
n=read(),m=read();
rin(i,0,n) F[i]=read();
rin(i,0,m) G[i]=read();
std::reverse(F,F+n+1);
std::reverse(G,G+m+1);
int nn=n,mm=m;
getinv(n-m+1);
//求Q
n=nn,m=mm;
memset(A,0,sizeof A);
rin(i,0,n) A[i]=F[i];
m=n+n-m;
len=0;
for(n=1;n<=m;n<<=1) len++;
rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
invn=qpow(n,MOD-2);
ntt(A,1);
ntt(B,1);
rin(i,0,n-1) A[i]=A[i]*B[i]%MOD;
ntt(A,-1);
n=nn,m=mm;
rin(i,0,n-m) Q[i]=A[i];
//求R
std::reverse(F,F+n+1);
std::reverse(G,G+m+1);
std::reverse(Q,Q+n-m+1);
memset(A,0,sizeof A);
memset(B,0,sizeof B);
rin(i,0,n-m) A[i]=Q[i];
rin(i,0,m) B[i]=G[i];
m=n;
len=0;
for(n=1;n<=m;n<<=1) len++;
rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
invn=qpow(n,MOD-2);
ntt(A,1);
ntt(B,1);
rin(i,0,n-1) A[i]=A[i]*B[i]%MOD;
ntt(A,-1);
rin(i,0,m-1) R[i]=(F[i]-A[i]+MOD)%MOD;
n=nn,m=mm;
rin(i,0,n-m) printf("%lld ",Q[i]);
printf("
");
rin(i,0,m-1) printf("%lld ",R[i]);
printf("
");
return 0;
}
多项式开方
关于这一部分博主还不能很透彻地理解，所以先放出简要推倒过程。在给出多项式的常数项不为(0)或(1)时，求开方后多项式的常数项时需要用到二次剩余。但是由于博主目前仍未学习二次剩余（太菜了），所以只能给出给出多项式的常数项为(0)或(1)时的代码，不久之后会完善这一部分。
对于一个多项式(A(x))，如果存在一个多项式(B(x))，满足(B(x)^2 equiv A(x) (mod x^n))，则称(B(x))为(A(x))在(mod x^n)下的平方根。
所有的运算在模(998244353)意义下进行。
推倒过程
(B(x))的常数项可以由(A(x))的常数项二次剩余求出。特别的，当(A(x))的常数项为(0)或(1)时，(B(x))的常数项与(A(x))相等。
类似多项式求逆，假设我们已经求出了(C(x)^2 equiv A(x) (mod x^{lceil frac{n}{2} ceil}))，现在我们想要求(B(x)^2 equiv A(x) (mod x^n))。
显然(B(x)^2 equiv A(x) (mod x^{lceil frac{n}{2} ceil}))。
所以有(B(x)^2-C(x)^2 equiv 0 (mod x^{lceil frac{n}{2} ceil}))。
两边同时平方，得：
[(B(x)^2-C(x)^2)^2 equiv 0 (mod x^n)
]
[(B(x)^2+C(x)^2)^2 equiv 4B(x)^2C(x)^2 (mod x^n)
]
[B(x)^2+C(x)^2 equiv 2B(x)C(x) (mod x^n)
]
[A(x)+C(x)^2 equiv 2B(x)C(x) (mod x^n)
]
[B(x) equiv (A(x)+C(x)^2) imes (2C(x))^{-1} (mod x^n)
]
然后多项式求逆即可。
时间复杂度：(O(nlogn))。
Update on 2018/12/13：一段更靠谱的证明：
[(B(x)^2-C(x)^2)^2 equiv 0 (mod x^n)
]
[(B(x)^2+C(x)^2)^2 equiv 4B(x)^2C(x)^2 (mod x^n)
]
[frac{(B(x)^2+C(x)^2)^2}{4 imes C(x)^2} equiv B(x)^2 (mod x^n)
]
[B(x) equiv frac{B(x)^2+C(x)^2}{2 imes C(x)} (mod x^n)
]
[B(x) equiv frac{A(x)+C(x)^2}{2 imes C(x)} (mod x^n)
]
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,a) for(int i=head[(a)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
inline int read(){
int x=0;char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
return x;
}
const int MAXN=100005;
const int MAXLEN=270005;
const LL MOD=998244353;
const LL ROOT=3;
const LL INVROOT=332748118;
const LL INV2=499122177;
int n,m;
int len,rev[MAXLEN];
int a[MAXN];
LL invn,A[MAXLEN],B[MAXLEN],C[MAXLEN];
inline LL qpow(LL x,LL y){
LL ret=1,tt=x%MOD;
while(y){
if(y&1) ret=ret*tt%MOD;
tt=tt*tt%MOD;
y>>=1;
}
return ret;
}
inline void ntt(LL *c,int dft){
rin(i,0,n-1) if(i<rev[i])
std::swap(c[i],c[rev[i]]);
for(int mid=1;mid<n;mid<<=1){
int r=(mid<<1);
LL u=qpow(dft>0?ROOT:INVROOT,(MOD-1)/r);
for(int l=0;l<n;l+=r){
LL v=1;
for(int i=0;i<mid;i++,v=v*u%MOD){
LL x=c[l+i],y=c[l+mid+i]*v%MOD;
c[l+i]=x+y;
if(c[l+i]>=MOD) c[l+i]-=MOD;
c[l+mid+i]=x-y;
if(c[l+mid+i]<0) c[l+mid+i]+=MOD;
}
}
}
if(dft<0) rin(i,0,n-1)
c[i]=c[i]*invn%MOD;
}
void getinv(int Floor){
if(Floor==1){
C[0]=qpow(B[0],MOD-2);
return;
}
getinv((Floor+1)>>1);
m=Floor-1+((((Floor+1)>>1)-1)<<1);
len=0;
for(n=1;n<=m;n<<=1) len++;
invn=qpow(n,MOD-2);
rin(i,0,n-1) A[i]=(i<Floor?B[i]:0);
rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
ntt(A,1);
ntt(C,1);
rin(i,0,n-1) C[i]=(((2*C[i]-A[i]*C[i]%MOD*C[i])%MOD)+MOD)%MOD;
ntt(C,-1);
rin(i,Floor,n-1) C[i]=0;
}
void solve(int Floor){
if(Floor==1){
B[0]=a[0];
return;
}
solve((Floor+1)>>1);
getinv(Floor);
m=Floor-1+((((Floor+1)>>1)-1)<<1);
len=0;
for(n=1;n<=m;n<<=1) len++;
invn=qpow(n,MOD-2);
rin(i,0,n-1) A[i]=(i<Floor?a[i]:0);
rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
// ntt(A,1);
// ntt(B,1);
// ntt(C,1);
// rin(i,0,n-1) B[i]=(A[i]+B[i]*B[i])%MOD*C[i]%MOD*INV2%MOD;
// ntt(B,-1);
ntt(A,1);
ntt(C,1);
rin(i,0,n-1) A[i]=A[i]*C[i]%MOD;
ntt(A,-1);
rin(i,0,Floor-1) B[i]=(A[i]+B[i])*INV2%MOD;
rin(i,Floor,n-1) B[i]=0;
rin(i,0,n-1) C[i]=0;
}
int main(){
n=read();
n--;
rin(i,0,n) a[i]=read();
int nn=n;
solve(n+1);
rin(i,0,nn) printf("%lld ",B[i]);
printf("
");
return 0;
}
拉格朗日插值
由初中知识可知，(n)个点((x_i,y_i))可以唯一地确定一个多项式。
现在，给定(n)个点，请你确定这个多项式，并将(k)代入求值，无需输出多项式每项的系数。
求出的值对(998244353)取模。
推倒过程
令(f(x)=sum_{i=1}^{n}y_il_i(x))，(l_i(x)=prod_{j=1,j eq i}^{n}frac{x-x_j}{x_i-x_j})。
容易发现(l_i(x_i)=1)，(l_i(x_j)=0 (j eq i))。
显然(f(x))过所有(n)个点，(f(x))为所求多项式，求出(f(k))即可。
时间复杂度：(O(n^2))。
Update on 2019/2/14：如果点值连续或等差，时间复杂度可以优化到(O(n))。
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,a) for(int i=head[(a)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
inline int read(){
int x=0;char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
return x;
}
const int MAXN=100005;
const LL MOD=998244353;
int n;
LL k,X[MAXN],Y[MAXN];
inline LL qpow(LL x,LL y){
LL ret=1,tt=x%MOD;
while(y){
if(y&1) ret=ret*tt%MOD;
tt=tt*tt%MOD;
y>>=1;
}
return ret;
}
int main(){
n=read(),k=read();
rin(i,1,n){
X[i]=read();
Y[i]=read();
}
LL ans=0;
rin(i,1,n){
LL son=1,mot=1;
rin(j,1,n){
if(i==j) continue;
son=son*(k-X[j]+MOD)%MOD;
mot=mot*(X[i]-X[j]+MOD)%MOD;
}
ans=(ans+Y[i]*son%MOD*qpow(mot,MOD-2))%MOD;
}
printf("%lld
",ans);
return 0;
}
所有给出点的横坐标是连续的话，还可以做到(O(n))求解。但是博主太菜了，你们懂的。
任意模数(NTT)
三模数(NTT)
由于两个多项式卷积后每项系数的大小是(O(na_i^2))级别的，所以我们只需选取两到三个大小合适的可以进行(NTT)的素数分别(NTT)后(CRT)合并答案即可，可能需要用到快速乘防止爆(long long)。
需要做(9)遍(NTT)，常数巨大。
拆系数(FFT)
我们可以把给出的多项式系数(a_i)写成(a_i=p_i imes 32768+q_i)。同理(b_i=r_i imes 32768+s_i)。
那么就有：
[a(x) imes b(x)=(p(x) imes 32768+q(x)) imes (r(x) imes 32768+s(x))
]
整理，得：
[a(x) imes b(x)=p(x) imes r(x) imes 1073741824+(p(x) imes s(x)+q(x) imes r(x)) imes 32768+q(x) imes s(x)
]
这样只需要做(4)遍(DFT)和(3)遍(IDFT)了。
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,a) for(int i=head[(a)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
inline int read(){
int x=0;char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
return x;
}
const int MAXN=100005;
const int MAXLEN=270005;
const long double pi=std::acos(-1);
int n,m;LL MOD;
int len,rev[MAXLEN];
int a[MAXN],b[MAXN];
struct Complex{
long double real,imag;
inline friend Complex operator + (Complex x,Complex y){
return (Complex){x.real+y.real,x.imag+y.imag};
}
inline friend Complex operator - (Complex x,Complex y){
return (Complex){x.real-y.real,x.imag-y.imag};
}
inline friend Complex operator * (Complex x,Complex y){
return (Complex){x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real};
}
};
Complex A[MAXLEN],B[MAXLEN],C[MAXLEN],D[MAXLEN],E[MAXLEN];
inline void fft(Complex *c,int dft){
rin(i,0,n-1) if(i<rev[i])
std::swap(c[i],c[rev[i]]);
for(int mid=1;mid<n;mid<<=1){
int r=(mid<<1);
Complex u=(Complex){std::cos(pi/mid),dft*std::sin(pi/mid)};
for(int l=0;l<n;l+=r){
Complex v=(Complex){1,0};
for(int i=0;i<mid;i++,v=v*u){
Complex x=c[l+i],y=c[l+mid+i]*v;
c[l+i]=x+y;
c[l+mid+i]=x-y;
}
}
}
if(dft<0) rin(i,0,n-1)
c[i].real/=n;
}
int main(){
n=read(),m=read(),MOD=read();
rin(i,0,n){
a[i]=read();
A[i].real=(a[i]>>15);
B[i].real=(a[i]&32767);
}
rin(i,0,m){
b[i]=read();
C[i].real=(b[i]>>15);
D[i].real=(b[i]&32767);
}
int nn=n,mm=m;
m+=n;
for(n=1;n<=m;n<<=1) len++;
rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
fft(A,1);
fft(B,1);
fft(C,1);
fft(D,1);
rin(i,0,n-1){
E[i]=A[i]*D[i]+B[i]*C[i];
A[i]=A[i]*C[i];
B[i]=B[i]*D[i];
}
fft(A,-1);
fft(B,-1);
fft(E,-1);
n=nn,m=mm;
rin(i,0,n+m) printf("%lld ",((((LL)(A[i].real+0.5)%MOD)<<30)%MOD+(((LL)(E[i].real+0.5)%MOD)<<15)%MOD+(LL)(B[i].real+0.5))%MOD);
printf("
");
return 0;
}
拆系数(FFT)（(4)遍(DFT)）
利用(FFT)的虚部空间可以实现一遍(DFT)对两个多项式进行求值，一遍(IDFT)对两组点值进行插值。
只需做(2)遍(DFT)和(2)遍(IDFT)即可完成，代码另外放出。
之后要补充的东西
牛顿迭代法
多项式多点求值
多项式指数，对数函数
*多项式快速插值
update on 2019.3.23：填了一些坑，看这里。