多项式半家桶
upd on 2025-02-08 19:33:46 星期六：优化了 exp 常数，快了一倍。
upd on 2025-02-09 14:54:01 星期日：整体优化，用了 vector 和 重载运算符，预处理了单位根，但是还没写除法，别急，现在的板子是高速状态。
upd on 2025-02-09 16:20:10 星期日：更新了除法，取模，开根。
upd on 2025-02-09 16:42:10 星期日：更新了使用指南。
upd on 2025-02-09 17:33:05 星期日：修复了一个小错误。
速度详情
多项式乘法最大点 389 ms，总时间 1.13 s
多项式乘法逆最大点 58 ms，总时间 407 ms
多项式ln最大点 82 ms，总时间 432 ms
多项式exp最大点 172 ms，总时间 851 ms
多项式快速幂最大点 229 ms，总时间 3.05s
多项式除法最大点 88 ms，总时间 385ms
多项式开根最大点 168ms，总时间 2.22s
速度还是很可观的，都用 cin/cout （关闭同步流），没有刻意卡常，基本都在最优解前几页（再前边基本都是指令集老哥了，实在比不过）。
有卷积，求逆，exp，ln，除法，取模，开根，乘方。
使用指南：
对于所有传入的多项式，必须保证是个 \(n\) 项式。
如果要改模数，在 namespace 前几行有 ___P,___G,___Gi 分别是模数，原根，原根的逆元，自行更改。默认 \(10^6\)，可以更改 ___LEN，注意要保证 \(\_\_\_N\) 为 \(x+10\) 的形式。
定义一个 \(n\) 项多项式 ： poly f(n)。
- 卷积： 就是
*，返回a*b。 - 求逆： 调用
invp(f,n)，返回 \(f^{-1}\)。 - 除法： 就是
/。 - 取模： 就是
%。 - 开根： 就是
sqrtp(f,n)，返回 \(\sqrt{f}\)。 - exp： 就是
exp(f,n)。 - ln: 就是
ln(f,n)。 - +，-。就是
+,-。 - Diff 是求导，Int 是积分。
- Pow(x,k) 是求 \(x^k\)，如果不保证 \(x[0]=1\)，请调用
powp(x,to_string(k))。
调用以上所有函数都需要在主函数先调用一次 Poly_pre()，调用 ln,exp,Pow 时需要先调用一次 _invpre()。
代码在下面！！！
code
namespace PolyBasedNTT{
#define IL inline
const int __LEN = 4e6 + 10;
const int ___P = 998244353;
const int ___G = 3,___Gi = 332748118;
using poly = vector<int>;
int _t[__LEN];poly __RT[22];
int lastnnn = 0,lastlim = 0;
poly invint;int lastinvn = 0;
template<class T> T FMod(const T &x){return x >= ___P?x-___P:x;}
IL void _invpre(){
invint.resize(__LEN-9);invint[1] = 1;
for(int i = 2;i <= __LEN-10; ++i) invint[i] = 1ll*invint[___P%i]*(___P-___P/i)%___P;
}
IL int Quick_power(int a,int b,int Mod = ___P){
int res = 1;
for(;b;b >>= 1,a = 1ll*a*a%Mod)
if(b&1) res = 1ll*res*a%Mod;
return res;
}
IL int inv_int(int a){return Quick_power(a,___P-2);}
IL void Poly_pre(){
__RT[0].resize(1,1);int pp = __lg(__LEN);
for(int i = 1;i < pp; ++i){
poly &g = __RT[i];g.resize(1<<i,1);
const int t = Quick_power(___G,___P>>(i+1));
for(int j = 1;j < (1<<i); ++j) g[j] =1ll*g[j-1]*t%___P;
}
}
IL int _prepare(int n){
if(lastnnn == n) return lastlim;
lastnnn = n;int lim = 1,ct = -1;
while(lim < n) lim <<= 1,ct++;
for(int i = 0;i < lim; ++i) _t[i] = (_t[i>>1]>>1)|((i&1)<<ct);
return (lastlim = lim);
}
IL void NTT(poly &a,int n,bool type = true){
n = _prepare(n);a.resize(n);
for(int i = 1;i < n; ++i) if(i < _t[i]) swap(a[i],a[_t[i]]);
for(int mid = 1,now = 0;mid < n;mid <<= 1, ++now){
auto &G = __RT[now];
for(int Res = mid << 1,j = 0;j < n;j += Res){
for(int k = 0;k < mid; ++k){
int x = a[j+k],y = 1ll*G[k]*a[j+k+mid]%___P;
a[j+k] = FMod(x + y);
a[j+k+mid] = FMod(x - y + ___P);
}
}
}
if(!type) reverse(a.begin()+1,a.end());
}
IL void _Reverse(poly &a){reverse(a.begin(),a.end());}
IL poly operator * (poly a,poly b){
int lena = a.size(),lenb = b.size();
if(lena <= 50 && lenb <= 50){
poly res(lena+lenb-1);
for(int i = 0;i < lena; ++i) for(int j = 0;j < lenb; ++j)
res[i+j] = FMod(res[i+j] + 1ll*a[i]*b[j]%___P);
return res;
}
int len = lena+lenb,Len = _prepare(len),llen = inv_int(Len);
NTT(a,len);NTT(b,len);
for(int i = 0;i < Len; ++i) a[i] = 1ll*a[i]*b[i]%___P;
NTT(a,len,false);
for(int i = 0;i < len; ++i) a[i] = 1ll*a[i]*llen%___P;
a.resize(len-1);
return a;
}
IL poly operator + (poly a,poly b){
int lena = a.size(),lenb = b.size();
int len = max(lena,lenb);
poly res(len);
for(int i = 0;i < lena; ++i) res[i] = FMod(res[i] + a[i]);
for(int i = 0;i < lenb; ++i) res[i] = FMod(res[i] + b[i]);
return res;
}
IL poly operator - (poly a,poly b){
int lena = a.size(),lenb = b.size(),len = max(lena,lenb);
poly res(len);
for(int i = 0;i < lena; ++i) res[i] = FMod(res[i] + a[i]);
for(int i = 0;i < lenb; ++i) res[i] = FMod(res[i] - b[i] + ___P);
return res;
}
poly invp(poly x,int n){
if(n == 1) return {inv_int(x[0])};
poly y = invp(poly(x.begin(),x.begin()+((n+1)>>1)),(n+1)>>1);
int len = _prepare(n<<1);
NTT(x,n<<1);NTT(y,n<<1);
for(int i = 0;i < len; ++i) x[i] = 1ll*y[i]*FMod(2-1ll*x[i]*y[i]%___P+___P)%___P;
NTT(x,n<<1,false);int llen = inv_int(len);
for(int i = 0;i < len; ++i) x[i] = 1ll*x[i]*llen%___P;
x.resize(n);
return x;
}
IL poly operator / (poly x,poly y){
int n = x.size(),m = y.size(),len = n - m + 1;
_Reverse(x);_Reverse(y);y.resize(len);
poly res = x*invp(y,len);
res.resize(len);_Reverse(res);
return res;
}
IL poly operator % (poly x,poly y){
poly res = x - (x/y)*y;
res.resize(y.size()-1);
return res;
}
IL poly Diff(poly x){
int len = x.size() - 1;poly res(len);
for(int i = 0;i < len; ++i) res[i] = 1ll*x[i+1]*(i+1)%___P;
return res;
}
IL poly Int(poly x){
int len = x.size() + 1;poly res(len);
for(int i = 1;i < len; ++i) res[i] = 1ll*x[i-1]*invint[i]%___P;
return res;
}
poly ln(poly x,int n){
poly res = Int(invp(x,n)*Diff(x));
res.resize(n);return res;
}
poly exp(poly x,int n){
if(n == 1) return {1};
poly y = exp(poly(x.begin(),x.begin()+((n+1)>>1)),(n+1)>>1);
poly res = y*(poly{1}-ln(y,n)+x);
res.resize(n);
return res;
}
IL int BSGS(int a,int p,int n){
unordered_map<int,int> mp;
int m = ceil(sqrt(p)),tmp = 1;
for(int i=0;i<m;i++) mp[1ll*tmp*n%p] = i+1,tmp = 1ll*tmp*a%p;
int tmpp = 1;
for(int i=0;i<=m;i++){
if(mp[tmpp]) return (1ll*i*m%p-mp[tmpp]+1+p)%p;
tmpp=1ll*tmpp*tmp%p;
}
return -1;
}
poly sqrtp(poly x,int n){
if(n == 1){
int p = BSGS(3,___P,x[0]);
p = Quick_power(3,p/2);
return {min(p,___P-p)};
}
//x.resize(n);
poly y = sqrtp(poly(x.begin(),x.begin()+(n+1)/2),(n+1)/2);
poly res = (y+x*invp(y,n))*poly({(___P+1)/2});res.resize(n);
return res;
}
poly Pow(poly x,int k){
int n = x.size();
x = ln(x,n);
for(int i = 0;i < n; ++i) x[i] = 1ll*x[i]*k%___P;
return exp(x,n);
}
poly Pow(poly x,int k,int n){
x = ln(x,n);
for(int i = 0;i < n; ++i) x[i] = 1ll*x[i]*k%___P;
x = exp(x,n);
return x;
}
int __GET(string s,int Mod = ___P - 1){
int res = 0;
for(int i = 0;i < s.size(); ++i)
res = (res*10ll + s[i]-48)%Mod;
return res;
}
poly powp(poly x,string k){
int n = x.size(),pp = -1;
for(int i = 0;i < n; ++i) if(x[i]) {pp = i;break;}
for(int i = 0;i < n-pp; ++i) x[i] = x[i+pp];
int a = x[0],ia = inv_int(a);
for(auto &i:x) i = 1ll*i*ia%___P;
a = Quick_power(a,__GET(k));
if(k.size() <= 9){
int p = 0;
for(char c:k) p = p*10+c-48;
if(1ll*pp*p < n){
x.resize(n-pp*p);
x = Pow(x,p,n-pp*p);
poly res(n);
for(int i = 0;i < n-pp*p; ++i) res[i+pp*p] = 1ll*x[i]*a%___P;
return res;
}
else return poly(n);
}
if(pp) return poly(n);
x = ln(x,n) * poly{__GET(k,___P)};
x = exp(x,n)*(poly){a};
x.resize(n);
return x;
}
#undef IL
}using namespace PolyBasedNTT;
代码在上面！！！
简单说一下原理（以后就不写了）。
多项式卷积是前置知识，自己找，也可以看看我的。
求逆
即求一个多项式 \(G(x)\) 使得 \(F(x)*G(x) \equiv 1\pmod{x^n}\)。
考虑倍增，假如当前已经求出了 \(\bmod\;{x^{\left\lceil\frac{n}{2}\right\rceil}}\) 的逆元，设其为 \(G^{\prime}(x)\)。
那么有 \(F*G^{\prime}\equiv 1\pmod {x^{\left\lceil\frac{n}{2}\right\rceil}}\)。
显然 \(F*G\equiv 1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)。
那么两式联立就有 \(G-G^{\prime}\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)。
两边平方，再同乘 \(F\) 有 \(F{G^{\prime2}}-2G^{\prime}+G\equiv 0\pmod{x^n}\)。
即 \(G\equiv 2G^{\prime}-FG^{\prime2}\pmod{x^n}\)。
除法
已知 \(F(x)=\sum\limits_{i=0}^{n-1}a_ix^i\)，定义 \(F^{T}(x)=\sum\limits_{i=0}^{n-1}a_{n-1-i}x^i\)。
就是把系数翻转。
已知一个 \(n\) 项式 \(F(x)\)，\(m\) 项式 \(G(x)\)，求一个 \(n-m+1\) 项式 \(Q(x)\)，一个项数为 \(m-2\) 的多项式 \(R(x)\)（不足用 \(0\) 补齐），使得 \(F(x)=Q(x)*G(x)+R(x)\)。
考虑如何删掉余项变成多项式求逆。一般通过 \(\bmod\;x^n\) 除去 \(\ge x^n\) 的项，但是这里 \(R\) 的项数很小，无法删去。先不管这些，直接暴力推式子。
这时候，巧妙地应用刚刚定义的那个操作，就有。
然后两边同时 \(\bmod\;x^{n-m+1}\)，就有
然后求出 \(Q,R\) 即可。
平方根
前置知识：牛顿迭代，简单求导，积分。
求导只需要知道 \(y=x^a,y^{\prime}=ax^{a-1}\)，还有 \((f(x)+g(x))^{\prime}=f^{\prime}(x)+g^{\prime}(x)\)。
积分只需要知道 \(\int x^a\mathrm{d}x=\frac{x^{a+1}}{a+1}+C\)，其中 \(C\) 是常数。
设 \(G(Q(x))=Q^2(x)-F(x)\)，即求 \(G(x)\) 的零点。
考虑倍增法，假设已经求出了 \({Q_*}^2(x)\equiv F(x)\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)。
直接大力上牛顿迭代，就有 \(Q(x)=Q_*(x)-\frac{G(Q_*(x))}{G^{\prime}(Q_*(x))}\)。
拆一下就是 \(Q(x)=Q_*(x)-\frac{Q_*^2(x)-F(x)}{2Q_*(x)}\)。
直接卷积求逆跑就行。
多项式 ln
复合函数求导链式法则：\(f(g(x))=f^{\prime}(g(x))g^{\prime}(x)\)。
设 \(Q(x)=ln(x)\)，直接推式子：
就是对 \(F(x)\) 求导，求逆，卷积，求出 \(G^{\prime}\)，然后再对 \(G^{\prime}\) 积分即可。
多项式 exp
即求 \(G(x)=e^{F(x)}\)，对两边取 \(ln\) 有 \(ln(G(x))=F(x)\)。
设 \(Q(x)=ln(x)-F\)，考虑牛顿迭代+倍增。
有 \(G(x)=G_*(x)-\frac{Q(G_*(x))}{Q^{\prime}(G_*(x))}\)。
就是 \(G(x)=G_*(x)-\frac{ln(G_*(x))-F(x)}{G_*^{-1}(x)}\)。
然后就有 \(G(x)=G_*(x)(1-ln(G_*(x))-F(x))\)。
倍增求即可。
本文来自博客园，作者：CuFeO4，转载请注明原文链接：https://www.cnblogs.com/hzoi-Cu/p/18703325
浙公网安备 33010602011771号