多项式板子
一、数组版本
数组版本和 poly 版本都只涵盖目录中第 \(4\sim 10\) 部分。
namespace Poly
{
int p[maxn],q[maxn],r[maxn],w[maxn];
int inum[maxn];
int qpow(int a,int k)
{
int res=1;
for(;k;a=1ll*a*a%mod,k>>=1) if(k&1) res=1ll*res*a%mod;
return res;
}
int add(int x,int y)
{
if((x+=y)>=mod) x-=mod;
return x;
}
int dec(int x,int y)
{
if((x-=y)<0) x+=mod;
return x;
}
int extend(int n)
{
return 1<<(__lg(n-1)+1);
}
void get_r(int n)
{
for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1?n>>1:0);
}
void init(int n=maxn)
{
for(int k=2,m=1;k<=n;k<<=1,m<<=1)
{
w[m]=1;
for(int i=m+1,x=qpow(3,(mod-1)/k);i<k;i++) w[i]=1ll*w[i-1]*x%mod;
}
for(int i=1;i<n;i++) inum[i]=qpow(i,mod-2);
}
void print(int *a,int n)
{
for(int i=0;i<n;i++) printf("%d%c",a[i]," \n"[i==n-1]);
}
void ntt(int *a,int n,int op)
{
for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
for(int k=2,m=1;k<=n;k<<=1,m<<=1)
for(int i=0;i<n;i+=k)
for(int j=i,*x=w+m;j<i+m;j++,x++)
{
int v=1ll*a[j+m]**x%mod;
a[j+m]=dec(a[j],v),a[j]=add(a[j],v);
}
if(op==-1)
{
reverse(a+1,a+n);
int v=qpow(n,mod-2);
for(int i=0;i<n;i++) a[i]=1ll*a[i]*v%mod;
}
}
void mul(int *a,int *b,int n,int m)
{
int len=extend(n+m-1);
for(int i=0;i<len;i++) p[i]=i<n?a[i]:0,q[i]=i<m?b[i]:0;
get_r(len),ntt(p,len,1),ntt(q,len,1);
for(int i=0;i<len;i++) a[i]=1ll*p[i]*q[i]%mod;
ntt(a,len,-1);
}
void inv(int *a,int *b,int n)
{
static int c[maxn],d[maxn];
n=extend(n),memset(b,0,4*n),b[0]=qpow(a[0],mod-2);
for(int k=2;k<=n;k<<=1)
{
for(int i=0;i<k<<1;i++) c[i]=i<k?a[i]:0,d[i]=i<k>>1?b[i]:0;
get_r(k),ntt(d,k,1);
for(int i=0;i<k;i++) d[i]=1ll*d[i]*d[i]%mod;
ntt(d,k,-1),get_r(k<<1),ntt(c,k<<1,1),ntt(d,k<<1,1);
for(int i=0;i<k<<1;i++) c[i]=1ll*c[i]*d[i]%mod;
ntt(c,k<<1,-1);
for(int i=0;i<k;i++) b[i]=(2ll*b[i]-c[i]+mod)%mod;
}
}
void diff(int *a,int *b,int n)
{
for(int i=1;i<n;i++) b[i-1]=1ll*i*a[i]%mod;
b[n-1]=0;
}
void integ(int *a,int *b,int n)
{
for(int i=1;i<n;i++) b[i]=1ll*inum[i]*a[i-1]%mod;
b[0]=0;
}
void ln(int *a,int *b,int n)
{
static int c[maxn],d[maxn];
assert(a[0]==1);
n=extend(n),inv(a,c,n),diff(a,d,n),mul(c,d,n,n),integ(c,b,n);
}
void exp(int *a,int *b,int n)
{
static int c[maxn];
assert(a[0]==0);
n=extend(n),memset(b,0,4*n),memset(c,0,4*n),b[0]=1;
for(int k=2;k<=n;k<<=1)
{
ln(b,c,k);
for(int i=0;i<k;i++) c[i]=dec(a[i],c[i]);
c[0]++,mul(b,c,k,k);
}
}
void sqrt(int *a,int *b,int n)
{
static const int inv2=(mod+1)>>1;
static int c[maxn],d[maxn];
assert(a[0]==1);
n=extend(n),memset(b,0,4*n),b[0]=1;
for(int k=2;k<=n;k<<=1)
{
memcpy(c,a,4*k),inv(b,d,k),mul(c,d,k,k);
for(int i=0;i<k;i++) b[i]=1ll*(b[i]+c[i])*inv2%mod;
}
}
/** enough for k<mod
void pow(int *a,int *b,int n,int k)
{
static int c[maxn];
memcpy(c,a,4*n),memset(b,0,4*n),b[0]=1;
for(;k;mul(c,c,n,n),k>>=1) if(k&1) mul(b,c,n,n);
}
**/
void pow(int *a,int *b,int n,string k)
{
static int c[maxn],d[maxn];
int u=0,k1=0,k2=0;
memset(b,0,4*n);
for(int i=0;i<k.size();i++) k1=(10ll*k1+k[i]-'0')%mod,k2=(10ll*k2+k[i]-'0')%(mod-1);
for(u=0;u<n&&!a[u];u++) ;
if((u&&k.size()>=5)||1ll*u*k1>=n) return ;
for(int i=u,x=qpow(a[u],mod-2);i<n;i++) c[i-u]=1ll*a[i]*x%mod;
ln(c,d,n-u);
for(int i=0;i<n-u;i++) d[i]=1ll*d[i]*k1%mod;
exp(d,c,n-u);
for(int i=u*k1,x=qpow(a[u],k2);i<n;i++) b[i]=1ll*c[i-u*k1]*x%mod;
}
void div(int *a,int *b,int *q,int *r,int n,int m)
{///len(q)=n-m+1,len(r)<=m-1
static int c[maxn],d[maxn],e[maxn];
int len=n-m+1;
for(int i=0;i<len;i++) c[i]=a[n-1-i];
for(int i=0;i<len;i++) d[i]=i<m?b[m-1-i]:0;
inv(d,e,len),mul(c,e,len,len),reverse(c,c+len),memcpy(q,c,4*len),mul(c,b,len,m);
for(int i=0;i<m;i++) r[i]=dec(a[i],c[i]);
}
}
温馨提示：
- 调用
Poly::init()函数进行初始化。 - 长度上限
maxn一般取extend(n<<1)，如 \(n=10^5\) 时，maxn=1<<18。 - 系数不能有负数，否则
add,dec失效后ntt函数会爆int。 - 传参数组长度类似
vector的size函数， \(\deg=n-1\) 。 - 多项式快速幂更推荐使用注释中的方法，双 \(\log\) 用时仅为单 \(\log\) 的三倍，但代码简洁许多。
二、vector 版本
namespace Poly
{
int r[maxn],w[maxn],inum[maxn];
int qpow(int a,int k)
{
int res=1;
for(;k;a=1ll*a*a%mod,k>>=1) if(k&1) res=1ll*res*a%mod;
return res;
}
int add(int x,int y)
{
if((x+=y)>=mod) x-=mod;
return x;
}
int dec(int x,int y)
{
if((x-=y)<0) x+=mod;
return x;
}
int extend(int n)
{
return 1<<(__lg(n-1)+1);
}
void get_r(int n)
{
for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1?n>>1:0);
}
void init(int n=maxn)
{
for(int k=2,m=1;k<=n;k<<=1,m<<=1)
{
w[m]=1;
for(int i=m+1,x=qpow(3,(mod-1)/k);i<k;i++) w[i]=1ll*w[i-1]*x%mod;
}
for(int i=1;i<n;i++) inum[i]=qpow(i,mod-2);
}
void print(poly a,int n)
{
a.resize(n);
for(int i=0;i<n;i++) printf("%d%c",a[i]," \n"[i==n-1]);
}
void ntt(poly &a,int n,int op)
{
for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
for(int k=2,m=1;k<=n;k<<=1,m<<=1)
for(int i=0;i<n;i+=k)
for(int j=i,*x=w+m;j<i+m;j++,x++)
{
int v=1ll*a[j+m]**x%mod;
a[j+m]=dec(a[j],v),a[j]=add(a[j],v);
}
if(op==-1)
{
reverse(a.begin()+1,a.begin()+n);
int v=qpow(n,mod-2);
for(int i=0;i<n;i++) a[i]=1ll*a[i]*v%mod;
}
}
poly operator*(poly a,int b)
{
for(int i=0;i<a.size();i++) a[i]=1ll*a[i]*b%mod;
return a;
}
poly operator+(poly a,poly b)
{
int n=max(a.size(),b.size());a.resize(n),b.resize(n);
for(int i=0;i<n;i++) a[i]=add(a[i],b[i]);
return a;
}
poly operator-(poly a,poly b)
{
int n=max(a.size(),b.size());a.resize(n),b.resize(n);
for(int i=0;i<n;i++) a[i]=dec(a[i],b[i]);
return a;
}
poly operator*(poly a,poly b)
{
int n=a.size(),m=b.size(),len=extend(n+m-1);
a.resize(len),b.resize(len),get_r(len),ntt(a,len,1),ntt(b,len,1);
for(int i=0;i<len;i++) a[i]=1ll*a[i]*b[i]%mod;
return ntt(a,len,-1),a.resize(n+m-1),a;
}
void operator*=(poly &a,int b)
{
a=a*b;
}
void operator+=(poly &a,poly b)
{
a=a+b;
}
void operator-=(poly &a,poly b)
{
a=a-b;
}
void operator*=(poly &a,poly b)
{
a=a*b;
}
poly inv(poly a,int n)
{
n=extend(n),a.resize(n);
poly b(n);b[0]=qpow(a[0],mod-2);
for(int k=2;k<=n;k<<=1)
{
poly c(a.begin(),a.begin()+k),d(b.begin(),b.begin()+k);
get_r(k),ntt(d,k,1);
for(int i=0;i<k;i++) d[i]=1ll*d[i]*d[i]%mod;
ntt(d,k,-1),c*=d;
for(int i=0;i<k;i++) b[i]=(2ll*b[i]-c[i]+mod)%mod;
}
return b;
}
poly diff(poly a,int n)
{
poly b(n);
for(int i=1;i<n;i++) b[i-1]=1ll*i*a[i]%mod;
return b[n-1]=0,b;
}
poly integ(poly a,int n)
{
poly b(n);
for(int i=1;i<n;i++) b[i]=1ll*inum[i]*a[i-1]%mod;
return b[0]=0,b;
}
poly ln(poly a,int n)
{
assert(a[0]==1);
n=extend(n),a.resize(n);
return integ(inv(a,n)*diff(a,n),n);
}
poly exp(poly a,int n)
{
assert(a[0]==0);
n=extend(n),a.resize(n);
poly b={1,0};
for(int k=2;k<=n;k<<=1)
{
poly c=ln(b,k);
for(int i=0;i<k;i++) c[i]=dec(a[i],c[i]);
c[0]++,b=b*c;
}
return b;
}
poly sqrt(poly a,int n)
{
static const int inv2=(mod+1)>>1;
assert(a[0]==1);
n=extend(n),a.resize(n);
poly b(n);b[0]=1;
for(int k=2;k<=n;k<<=1)
{
poly c(a.begin(),a.begin()+k);c*=inv(b,k);
for(int i=0;i<k;i++) b[i]=1ll*(b[i]+c[i])*inv2%mod;
}
return b;
}
poly operator<<(poly a,int k)
{
poly b(a.size()+k);
for(int i=0;i<a.size();i++) b[i+k]=a[i];
return b;
}
poly operator>>(poly a,int k)
{
if(a.size()<=k) return {0};
poly b(a.size()-k);
for(int i=k;i<a.size();i++) b[i-k]=a[i];
return b;
}
void operator<<=(poly &a,int k)
{
a=a<<k;
}
void operator>>=(poly &a,int k)
{
a=a>>k;
}
poly pow(poly a,int n,string k)
{
int u=0,k1=0,k2=0;
for(int i=0;i<k.size();i++) k1=(10ll*k1+k[i]-'0')%mod,k2=(10ll*k2+k[i]-'0')%(mod-1);
for(u=0;u<n&&!a[u];u++) ;
if((u&&k.size()>5)||1ll*u*k1>=n) return poly(n);
poly b(n),c(n-u);
for(int i=u,x=qpow(a[u],mod-2);i<n;i++) c[i-u]=1ll*a[i]*x%mod;
c=ln(c,n-u);
for(int i=0;i<n-u;i++) c[i]=1ll*c[i]*k1%mod;
c=exp(c,n-u);
for(int i=u*k1,x=qpow(a[u],k2);i<n;i++) b[i]=1ll*c[i-u*k1]*x%mod;
return b;
}
pair<poly,poly> div(poly a,poly b,int n,int m)
{///len(q)=n-m+1,len(r)<=m-1
a.resize(n),b.resize(m);
if(n<m) return make_pair(poly{0},a);
poly c=a,d=b;reverse(c.begin(),c.end()),c.resize(n-m+1),reverse(d.begin(),d.end());
poly q=c*inv(d,n-m+1);q.resize(n-m+1),reverse(q.begin(),q.end());
poly r=a-q*b;r.resize(m-1);
return make_pair(q,r);
}
poly operator/(poly a,poly b)
{
return div(a,b,a.size(),b.size()).first;
}
poly operator%(poly a,poly b)
{
return div(a,b,a.size(),b.size()).second;
}
void operator/=(poly &a,poly b)
{
a=a/b;
}
void operator%=(poly &a,poly b)
{
a=a%b;
}
}
using namespace Poly;
温馨提示：
-
调用
Poly::init()函数进行初始化。 -
如果不涉及到分治 \(\texttt{NTT}\) ，个人不太推荐使用
vector版本的板子。vector版本代码长度比数组版本高 \(20\%\) ，运行效率比数组版本低 \(10\%\) ，很大程度上是因为构造新的vector时间开销较大。vector版本由于涉及到重载运算符操作，所以不得不将Poly命名空间放入全局空间，但这很有可能导致变量重名。
三、快速傅里叶变换（FFT）
\(\texttt{FFT}\) 的核心思想：系数多项式转点值多项式（\(\texttt{DFT}\)），点值多项式可以直接相乘，点值多项式再转系数多项式（\(\texttt{IDFT}\)）。
暴力求点值时间复杂度依然是 \(\mathcal O(n^2)\) ，但是如果将 \(n\) 补成 \(2\) 的方幂，并且利用单位根的性质，我们可以做到 \(\mathcal O(n\log n)\) 。
记 \(w_n^k=\cos\frac{2k\pi}n+i\sin\frac{2k\pi}n\) ，目标求 \(f(w_n^0),\cdots,f(w_n^{n-1})\) 。
令：
容易发现 \(f(x)=f_1(x^2)+x\cdot f_2(x^2)\) ，对 \(0\le k\lt\frac n2\) ，分别代入 \(w_n^k\) 和 \(w_n^{k+\frac n2}\) ：
时间复杂度 \(T(n)=2\cdot T(\frac n2)+\mathcal O(n)\) ，解得 \(T(n)=\mathcal O(n\log n)\) 。
递归版本常数太大，我们希望改成迭代版本。
观察一下使用 \(a_k\) 的次序，每次我们将偶数项（二进制下末位为 \(0\)）放在前面，将奇数项（二进制下末位为 \(1\)）放在后面。
设 \(2^d=n\) ，那么 \(d-1\) 轮操作等价于将 \(a_k\) 的二进制位翻转。
对 \(\forall 2\le i\le d\) ，维护自底向上迭代 \(i\) 次后的结果，每次通过 \(f_1(w_{2^{i-1}}^k)\) 和 \(f_2(w_{2^{i-1}}^k)\) 计算 \(f(w_{2^i}^k)\) 和 \(f(w_{2^i}^{k+2^{i-1}})\) 。
上述过程也被称为蝴蝶变换。
记 \(w=w_n\) ， \(\texttt{DFT}\) 的本质是：
读者可以验证转移矩阵的逆矩阵为：
于是有了 \(\texttt{IDFT}\) 的第一种方法：用 \(w^{-1}\) 代替 \(w\) 重新做一遍 \(\texttt{FFT}\) ，最后除以 \(n\) 。
但是上面这种方法不利于卡常，究其原因是我们要同时预处理 \(w\) 和 \(w^{-1}\) 的方幂。
于是有了第二种看待 \(\texttt{DFT}\) 的视角：
根据上面的分析，可以得到：
如果将 \(g(x)=\texttt{IDFT}[f(x)]\) 视为 \(n\) 次多项式（原本是 \(n-1\) 次），则：
对比两边每一项的系数：
从而避开 \(w^{-1}\) ，完美。
四、快速数论变换（NTT）
\(\texttt{NTT}\) 使用条件： \(n\le 2^{v_2(p-1)}\) 。
单位根的性质原根全部满足，只需将 \(w_n\) 换成 \(g_n\) 即可。
重要卡常技巧：预处理所有 \(g_{2^i}^k\) 。
五、多项式求逆（inv）
给定 \(f(x)\) ，求 \(g(x)\) 满足 \(f(x)\cdot g(x)\equiv 1\pmod{x^n}\) 。
考虑倍增，假设已知 \(f(x)\cdot g_0(x)\equiv 1\pmod{x^\frac n2}\) 。
六、多项式对数函数（ln）
给定 \(f(x)\) 满足**常数项为 \(1\) **，求 \(g(x)\) 满足 \(g(x)\equiv\ln f(x)\pmod{x^n}\) 。
两边对 \(x\) 求导：
求导求逆乘起来，再积分回去，即可得到 \(g(x)\) 。
七、多项式指数函数（exp）
前置知识：牛顿迭代。
牛顿迭代用于求函数零点。
任取 \(x_0\) 作为初始值，取 \(x_1=x_0-\frac{f(x_0)}{f'(x_0)}\) 为 \(x_0\) 处的切线与 \(x\) 轴的交点，依此类推。
牛顿迭代收敛速度非常快，对于多项式，每做一次精度就会加倍。
回到原题。
给定 \(f(x)\) 满足**常数项为 \(0\) **，求 \(g(x)\) 满足 \(g(x)\equiv e^{f(x)}\pmod{x^n}\) 。
题目等价于求 \(F(g(x))=\ln g(x)-f(x)\) 的零点，这里 \(f(x)\) 为常数， \(g(x)\) 为变量。
考虑倍增，假设 \(g_0(x)\) 为 \(\bmod x^\frac n2\) 下的零点，则：
八、多项式开根（sqrt）
给定 \(f(x)\) 满足 **常数项为 \(1\) **，求 \(g(x)\) 满足 \(g^2(x)\equiv f(x)\pmod{x^n}\) 。
考虑倍增，假设 \(g_0^2(x)\equiv f(x)\pmod{x^\frac n2}\) 。
九、多项式快速幂（pow）
给定 \(f(x)\) ，求 \(g(x)\) 满足 \(g(x)\equiv f(x)^k\pmod{x^n}\) 。
方法一
指数 \(k\) 可对 \(p\) 取模，从而保证 \(k\lt p\) 。
类比普通快速幂，时间复杂度 \(\mathcal O(n\log n\log p)\) ，但常数小。
方法二
如果 \(f(x)\) 常数项为 \(1\) ，先取 \(\ln\) ，乘上 \(k\) 倍后再 \(\exp\) 即可。
如果 \(f(x)\) 常数项不为 \(1\) ，设 \(f(x)\) 的最低次项为 \(a_mx^m\) ，则：
注意如果 \(p\ge n\)，那么计算 \(a_m^k\) 时 \(k\) 对 \(\varphi(p)\) 取模，计算 \(f(x)^k\) 时 \(k\) 对 \(p\) 取模。
原因如下：
-
对于一个非 \(p\) 的倍数 \(a_m\) ，由费马小定理 \(a_m^{p-1}=1\) 。
-
对于一个常数项为 \(1\) 的多项式 \(f(x)\) ，我们证明 \(f(x)^p\equiv 1\pmod{x^n}\) 。
设 \(a_i\) 取了 \(c_i\) 项，若不存在 \(c_i=p\) ，则广义组合数 \(\binom p{c_0,\cdots,c_n}\) 是 \(p\) 的倍数。
若 \(c_i=p\) ，当 \(i\ge 1\) 时， \((a_ix^i)^p\) 次数高于 \(n\) ，可以直接舍弃。
十、多项式除法 & 取模（div）
给定 \(f(x),g(x)\) 满足 \(\deg f=n\ge m=\deg g\) ，求 \(q(x),r(x)\) 满足 \(\deg q=n-m,\deg r\le m-1\) ，且 \(f(x)=q(x)\cdot g(x)+r(x)\) 。
对于次数为 \(k\) 的多项式 \(A(x)\) ，记 \(A_R(x)=x^kA(\frac 1x)\) 。
为方便起见，记 \(\deg r=m-1\) 。
至此我们可以求出 \(q(x)\) ，最后令 \(r(x)=f(x)-q(x)\cdot g(x)\) 即可。
注意求逆是在 \(\bmod x^{n-m+1}\) 而不是 \(x^m\) 意义下进行。
本文来自博客园，作者：peiwenjun，转载请注明原文链接：https://www.cnblogs.com/peiwenjun/p/18541501
