【MOBAN】多项式算法模板和小记录
多项式乘法,多项式求逆,多项式除法,多项式开方,多项式求导和积分,多项式取模,多项式指数函数,多项式对数函数,多项式快速幂,多项式多点求值,多项式多点插值,常系数线性递推,多项式求复合逆,拆系数FFT
参考文章:https://blog.csdn.net/kscla/article/details/79356786
预备
不用long long 多件套和LJ NTTconst int maxn = 300005,mod=998244353; const int inv2 = 499122177; const int g = 3; int add(int x,int y) { x+=y;return x>=mod?x-mod:x; } int sub(int x,int y) { x-=y;return x<0?x+mod:x; } int mul(int x,int y) { return 1ll*x*y%mod; } int ksm(int a,int b) { int ans=1; for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a); return ans; } void ntt(int *a,int s,int dft) { for(int i=0,j=0;i<s;i++) { if(i<j) swap(a[i],a[j]); for(int k=(s>>1);(j^=k)<k;k>>=1); } for(int st=1;st<s;st<<=1) { int dwg = dft==1?ksm(g,(mod-1)/(st<<1)):ksm(g,mod-1-(mod-1)/(st<<1)); for(int i=0;i<s;i+=(st<<1)) { int ng = 1; for(int j=i;j<i+st;j++) { int x=a[j]; int y=mul(a[j+st],ng); ng = mul(ng,dwg); a[j]=add(x,y); a[j+st]=sub(x,y); } } } if(dft==1) return; int invs = ksm(s,mod-2); for(int i=0;i<s;i++) a[i]=mul(a[i],invs); }
多项式乘法
void MUL(int *a,int *b,int *c,int le) { for(int i=0;i<le;i++) ta[i]=b[i]; for(int j=0;j<le;j++) tb[j]=c[j]; ntt(ta,le,1); ntt(tb,le,1); for(int i=0;i<le;i++) ta[i]=mul(ta[i],tb[i]); ntt(ta,le,-1); for(int i=0;i<le;i++) a[i] = ta[i]; }
多项式求逆:
求A(x)*B(x)==1(mod x^n) 设定A(x)*G(x)==1(mod x^(n/2) 那么就有B(x) – G(x) == 0 (mod x^(n/2) ) 就有B(x)^2 + G(x)^2 – 2B(x)G(x) == 0 (mod x^(n)) 两边同乘上A,就最终可以得到 B(x) = G(x)(2-AG(x))了void ginv(int deg,int *a,int *b) { if(deg==1) { b[0]=ksm(a[0],mod-2); return; } ginv(deg>>1,a,b); for(int i=0;i<deg;i++) ta[i]=a[i]; for(int i=deg;i<(deg<<1);i++) ta[i]=0; ntt(ta,deg<<1,1); ntt(b,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i]))); ntt(b,(deg<<1),-1); for(int i=deg;i<(deg<<1);i++) b[i]=0; }
多项式除法
我们设定Ar(x) = x^n*A(1/x) 比如举个例子A(x) = x^3 + 2x^2 + 4x + 1 Ar(x) = 1+2x+4*x^3+x^3 这恰好等于A(x)的系数反转。 考虑原式我们只要将x替换成1/x,会发现一些有趣的事情。 F(1/x)x^n == Q(1/x)G(1/x)x^n + x^nR(x) Fr(x) = Qr(x) + Gr(x) + Rr(x)x^(n-m+1) 那么Ar(x) == Dr(x)Br(x) ( mod x^(n-m+1))#include<stdio.h> #include<cstdio> #include<algorithm> #include<cstring> using namespace std; const int mod = 998244353; const int g = 3; const int maxn = 700005; int mul(int x,int y) { return 1ll*x*y%mod; } int add(int x,int y) { x+=y; return x>=mod?x-mod:x;} int sub(int x,int y) { x-=y; return x<0?x+mod:x;} int ksm(int a,int b) { int ans=1; for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a); return ans; } int n,m; int tmp[maxn],ta[maxn],tb[maxn]; void ntt(int *a,int s,int dft) { for(int i=0,j=0;i<s;i++) { if(i<j) swap(a[i],a[j]); for(int k=(s>>1);(j^=k)<k;k>>=1); } for(int st=1;st<s;st<<=1) { int dwg = dft==1?ksm(g,(mod-1)/(st<<1)):ksm(g,mod-1-(mod-1)/(st<<1)); for(int i=0;i<s;i+=(st<<1)) { int ng = 1; for(int j=i;j<i+st;j++) { int x=a[j]; int y=mul(a[j+st],ng); ng = mul(ng,dwg); a[j]=add(x,y); a[j+st]=sub(x,y); } } } if(dft==1) return; int invs = ksm(s,mod-2); for(int i=0;i<s;i++) a[i]=mul(a[i],invs); } void ginv(int deg,int *a,int *b) { if(deg==1) { b[0]=ksm(a[0],mod-2); return; } ginv(deg>>1,a,b); for(int i=0;i<deg;i++) ta[i]=a[i]; for(int i=deg;i<(deg<<1);i++) ta[i]=0; ntt(ta,deg<<1,1); ntt(b,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i]))); ntt(b,(deg<<1),-1); for(int i=deg;i<(deg<<1);i++) b[i]=0; } int Q[maxn],F[maxn],G[maxn],Fr[maxn],Gr[maxn],Gri[maxn]; void MUL(int *a,int *b,int *c,int le) { for(int i=0;i<le;i++) ta[i]=b[i]; for(int j=0;j<le;j++) tb[j]=c[j]; ntt(ta,le,1); ntt(tb,le,1); for(int i=0;i<le;i++) ta[i]=mul(ta[i],tb[i]); ntt(ta,le,-1); for(int i=0;i<le;i++) a[i] = ta[i]; } int main() { scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) { scanf("%d",&F[i]); Fr[n-i]=F[i]; } for(int i=0;i<=m;i++) { scanf("%d",&G[i]); Gr[m-i]=G[i]; } int le = 1; for(;le<=(n-m);le<<=1); for(int i=n-m+1;i<=m;i++) Gr[i]=0; ginv(le,Gr,Gri); le = 1; for(;le<=(n<<1);le<<=1); memset(ta,0,sizeof ta); memset(tb,0,sizeof tb); MUL(Q,Gri,Fr,le); reverse(Q,Q+n-m+1); for(int i=n-m+1;i<le;i++) Q[i]=0; for(int i=0;i<=n-m;i++) printf("%d ",Q[i]); puts(""); le = 1; for(;le<=(n<<1);le<<=1); memset(ta,0,sizeof ta); memset(tb,0,sizeof tb); MUL(tmp,Q,G,le); for(int i=0;i<m;i++) printf("%d ",sub(F[i],tmp[i])); }
多项式取模(源自多项式除法)
我们已知除数和商和被除数,那么余就等于被除数-商×除数int Q[maxn],Fr[maxn],Gr[maxn],Gri[maxn],Yy[maxn]; void gmod(int *O,int *F,int *G,int n,int m) {//F%G==O (F n G m) for(int i=0;i<=n;i++) Fr[i] = F[n-i]; for(int i=0;i<=m;i++) Gr[i] = G[m-i]; for(int i=n-m+2;i<=m;i++) Gr[i] = 0; int le = 1; for(;le<=(n-m);le<<=1); for(int i=n-m+1;i<=m;i++) Gr[i]=0; ginv(le,Gr,Gri); le = 1; for(;le<=(n<<1);le<<=1); MUL(Q,Gri,Fr,le); reverse(Q,Q+n-m+1); for(int i=n-m+1;i<le;i++) Q[i]=0; MUL(tmp,Q,G,le); for(int i=0;i<m;i++) O[i] = sub(F[i],tmp[i]); }
多项式开方
考虑和多项式求逆用到的同样的倍增的方法 如果我已知B(x)*B(x) == A(x) mod(x^m) 求C(x)使得C(x)C(x) == A(x) mod(x^(2m)) 则(B(x)+C(x))*(B(x)-C(x))==0(mod x^(m)) 那么我们讨论B(x)-C(x)==0情况 最后又可以得到C(x)==(A(x)+B(x)^2)/(2*B(x)) (mod x^(2m)) code://代入a,b^2 = avoid ginv(int deg,int *a,int *b) { if(deg==1) { b[0]=ksm(a[0],mod-2); return; } ginv(deg>>1,a,b); for(int i=0;i<deg;i++) ta[i]=a[i]; for(int i=deg;i<(deg<<1);i++) ta[i]=0; ntt(ta,deg<<1,1); ntt(b,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i]))); ntt(b,(deg<<1),-1); for(int i=deg;i<(deg<<1);i++) b[i]=0; } void gsqrt(int deg,int *a,int *b) { if(deg==1) { b[0]=1; return; } gsqrt(deg>>1,a,b); memset(d,0,sizeof d); ginv(deg,b,d); for(int i=0;i<deg;i++) c[i]=a[i]; for(int i=deg;i<(deg<<1);i++) c[i]=0; ntt(c,deg<<1,1); ntt(b,deg<<1,1); ntt(d,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i]=mul(add(b[i],mul(c[i],d[i])),inv2); ntt(b,deg<<1,-1); for(int i=deg;i<(deg<<1);i++) b[i]=0; }
多项式求导和积分
有关微积分的参考: 某博客https://www.cnblogs.com/zwfymqz/p/9127894.html 高中数学书选修2-2 http://jiaofu.yousi.com/compontent/pdf/?url=http://jfpdf.yousi.com/160424030303585634.pdf#page=11] 若[latex] F(x) = \sum_{i=0}^{n}a_{i} * x^i [/latex] 则[latex] F'(x) = \sum_{i=1}^{n} i * a_{i} * x^{i-1} [/latex] 则[latex] \int F(x) = \sum_{i=1}^{n} \frac{a_{i}* x^{i+1} }_{i+1} [/latex] 所以我们如果把多项式看做一个函数就可以求导和积分了。int direv(int *a,int *b,int len) { //求导 for(int i=1;i<len;i++) b[i-1]=mul(a[i],i); b[len-1]=0; } void inter(int *a,int *b,int len) { //积分 for(int i=1;i<len;i++) b[i]=mul(a[i-1],ksm(i,mod-2)); b[0]=0; //这里也可以预处理逆元 }
多项式对数函数
给出 n-1 次多项式 A(x),求一个mod x^n下的多项式B(x)==ln A(x) 我们知道[latex]\ln x[/latex]的导数为[latex]\frac{1}{x}[/latex],我们设[latex]F(x) = \ln x [/latex] [latex]B(x) = F(A(x))[/latex] [latex]B(x)\prime = F(A(x))\prime * A(x)\prime[/latex] [latex]B(x) = \int \frac {A\prime}{A} [/latex] 直接代入搞就可以了void direv(int *a,int *b,int len) { //求导 for(int i=1;i<len;i++) b[i-1]=mul(a[i],i); b[len-1]=0; } void inter(int *a,int *b,int len) { //积分 for(int i=1;i<len;i++) b[i]=mul(a[i-1],ksm(i,mod-2)); b[0]=0; //这里也可以预处理逆元 } int F[maxn],G[maxn],A[maxn],B[maxn]; void getln(int deg,int *a,int *b) { direv(a,A,deg); ginv(deg,a,B); ntt(A,deg<<1,1); ntt(B,deg<<1,1); for(int i=0;i<(deg<<1);i++) A[i]=mul(A[i],B[i]); ntt(A,deg<<1,-1); inter(A,b,deg); for(int i=0;i<(deg<<1);i++) A[i]=B[i]=0; }
多项式指数指数
已知函数F,求G==e^F (mod) 由于牛顿迭代公式 xi= xi-1 - f(xi-1) / f ' (xi-1) 那么初始x0 = ea0 显然只有a0等于0时有意义。 因为lnG - F = 0 套用牛顿迭代公式得到Gi = Gi-1'(1-lnGi-1'+F) 就可以搞了。#include<stdio.h> #include<iostream> #include<cstdio> #include<algorithm> #include<cmath> using namespace std; const int mod = 998244353; const int g = 3; const int maxn = 400005; int mul(int x,int y) { return 1ll*x*y%mod; } int add(int x,int y) { x+=y; return x>=mod?x-mod:x;} int sub(int x,int y) { x-=y; return x<0?x+mod:x; } int n; int a[maxn],b[maxn],c[maxn]; int ksm(int a,int b){ int ans=1; for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a); return ans; } void ntt(int *a,int s,int dft) { for(int i=0,j=0;i<s;i++) { if(i<j) swap(a[i],a[j]); for(int k=(s>>1);(j^=k)<k;k>>=1); } for(int st=1;st<s;st<<=1) { int dwg = dft==1?ksm(g,(mod-1)/(st<<1)):ksm(g,mod-1-(mod-1)/(st<<1)); for(int i=0;i<s;i+=(st<<1)) { int ng = 1; for(int j=i;j<i+st;j++) { int x=a[j]; int y=mul(a[j+st],ng); ng = mul(ng,dwg); a[j]=add(x,y); a[j+st]=sub(x,y); } } } if(dft==1) return; int invs = ksm(s,mod-2); for(int i=0;i<s;i++) a[i]=mul(a[i],invs); } int ta[maxn]; void ginv(int deg,int *a,int *b) { if(deg==1) { b[0]=ksm(a[0],mod-2); return; } ginv(deg>>1,a,b); for(int i=0;i<deg;i++) ta[i]=a[i]; for(int i=deg;i<(deg<<1);i++) ta[i]=0; ntt(ta,deg<<1,1); ntt(b,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i]))); ntt(b,(deg<<1),-1); for(int i=deg;i<(deg<<1);i++) b[i]=0; } void direv(int *a,int *b,int len) { //求导 for(int i=1;i<len;i++) b[i-1]=mul(a[i],i); b[len-1]=0; } void inter(int *a,int *b,int len) { //积分 for(int i=1;i<len;i++) b[i]=mul(a[i-1],ksm(i,mod-2)); b[0]=0; //这里也可以预处理逆元 } int F[maxn],G[maxn],A[maxn],B[maxn]; void getln(int deg,int *a,int *b) { direv(a,A,deg); ginv(deg,a,B); ntt(A,deg<<1,1); ntt(B,deg<<1,1); for(int i=0;i<(deg<<1);i++) A[i]=mul(A[i],B[i]); ntt(A,deg<<1,-1); inter(A,b,deg); for(int i=0;i<(deg<<1);i++) A[i]=B[i]=0; } int f[maxn]; void EXP(int *a,int *b,int deg) { if(deg==1) {b[0]=1; return;} EXP(a,b,deg>>1); getln(deg,b,f); f[0]=sub(a[0]+1,f[0]); for(int i=1;i<deg;i++) f[i]=sub(a[i],f[i]); ntt(f,deg<<1,1); ntt(b,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i] = mul(b[i],f[i]); ntt(b,deg<<1,-1); for(int i=deg;i<(deg<<1);i++) b[i]=f[i]=0; } int main() { scanf("%d",&n); for(int i=0;i<n;i++) scanf("%d",&F[i]); int len = 1; for(;len<n;len<<=1); EXP(F,G,len); for(int i=0;i<n;i++) printf("%d ",G[i]); }
多项式快速幂
啥,直接二分快速幂? 不不不,这是O(nlognlogk)的,他不够优秀 我们知道F(x)^k == exp(ln(F(x)^k)) 并且ln(F(x)^k)==k*ln(F(x)) 我们对多项式求ln,然后系数乘以k之后exp还原就好 注意细节! 我们发现求ln的时候用到的求导和反积分会损失掉常数项! 所以在我们回代回exp之前记得单独处理常数项。int len = 1; for(;len<=(2*max(m,n));len<<=1); getln(len,F,OF);//求出F的k次快速幂 memset(F,0,sizeof F); for(int i=0;i<len;i++) F[i] = mul(k%mod,OF[i]); memset(OF,0,sizeof OF); OF[0] = ksm(orz,k%(mod-1));//由于求导和积分会损失掉常数项,为保证牛顿迭代正确,必须要常数项带入 EXP(F,OF,len);
常系数线性递推
求一个满足k阶齐次线性递推数列a_i的第n项。 即:[latex]a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i} [/latex] 首先矩阵卡速米的做法很显然,把这个递推矩阵写出来,然后直接爆推就可以了。然而k这么大,一次矩阵乘法就GG了。 我们如果能够构造一个数列ci使得[latex] A^{n} = \sum_{i=0}^{k} c_{i} * A_{i} [/latex]。但其实我们完全不用真的算出A^n,因为实际上我们想要得到一个递推矩阵,然后最后乘上一个列向量ST(给定的初始那k个数)。 我们把这个向量乘进去[latex] ST * A^{n} = \sum_{i=0}^{k} c_{i} * ST * A_{i} [/latex]。我们只需要他最低的那项。我们仔细观察发现这里的ST*Ai的最低项不就是STi吗?所以我们只需要构造出这样一个数列ci,我们就可以得到,答案就是[latex] ANS = \sum_{i=0}^{k-1} c_{i} * ST_{i} [/latex]。现在我们考虑怎么构造这个数列。 如果我们把一个这个递推矩阵的特征多项式手算出来:引入 特征多项式
我们已知对于一个矩阵A,若有非0行(列)向量y,和常数x,使得Ay = xy,那么就称y为A的特征向量,x为A的特征值。我们移一个项,使得|A-Ex|y=0 (E为单位矩阵),由于y行列式不为0,那么一定|A-Ex| = 0 ,同时我们称f(x) = |A-Ex|为矩阵A的特征多项式。由于某个神奇的定理(Cayley-Hamilton定理),f(A) = 0。(啥?你说为什么A是一个矩阵而不是一个常数?我们其实把这个多项式f的x替换成矩阵变量就可以了,感性理解的话|A-A*E| = 0) 所以手算出来后发现,这个多项式恰好是[latex] f(A) = (-1)^{k} * (A^k - \sum_{i=1}^{k} a_{i}*A^{k-i}) [/latex] 由于f(A)==0 , 那么 直接消掉(-1)^{k},得到[latex]f(A) = (A^k - \sum_{i=1}^{k} a_{i}*A^{k-i}) [/latex] 也就是说A^n == A^n % (f(A)) , 之后我们就将这个以矩阵为变量的多项式愉快地降低到了次数k-1。也就是说,我们只要构出了这个多项式,也就是我们要的ci数列对应构造出的多项式。多项式取模即可 luogu 线性递推// luogu-judger-enable-o2 #include<stdio.h> #include<iostream> #include<cstring> #include<algorithm> #include<cstdio> #include<cmath> using namespace std; const int maxn = 32005*6; const int mod = 998244353; const int g = 3; int add(int x,int y) { x+=y; return x>=mod?x-mod:x; } int sub(int x,int y) { x-=y; return x<0?x+mod:x; } int mul(int x,int y) { return 1ll*x*y%mod; } int ksm(int a,int b) { int ans=1; for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a); return ans; } void ntt(int *a,int s,int dft) { for(int i=0,j=0;i<s;i++) { if(i<j) swap(a[i],a[j]); for(int k=(s>>1);(j^=k)<k;k>>=1); } for(int st=1;st<s;st<<=1) { int dwg; if(dft==1) dwg=ksm(g,(mod-1)/(st<<1) ); else dwg=ksm(g,mod-1-(mod-1)/(st<<1)); for(int i=0;i<s;i+=(st<<1)) { int ng = 1; for(int j=i;j<i+st;j++) { int x = a[j]; int y = mul(a[j+st],ng); ng = mul(ng,dwg); a[j] = add(x,y); a[j+st] = sub(x,y); } } } if(dft==1) return; int invs = ksm(s,mod-2); for(int i=0;i<s;i++) a[i]=mul(a[i],invs); } int ta[maxn],tb[maxn]; void ginv(int *a,int *b,int deg) {// in>>a out<<b if(deg==1) { b[0] = ksm(a[0],mod-2); return;} ginv(a,b,deg>>1); for(int i=0;i<deg;i++) ta[i] = a[i]; for(int i=deg;i<(deg<<1);i++) ta[i]=0; ntt(b,deg<<1,1); ntt(ta,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i]))); ntt(b,deg<<1,-1); for(int i=deg;i<(deg<<1);i++) b[i]=0; } void MUL(int *a,int *b,int *c,int deg){ //in>>a>>b out<<c for(int i=0;i<deg;i++) ta[i]=a[i]; for(int i=0;i<deg;i++) tb[i]=b[i]; ntt(ta,deg,1); ntt(tb,deg,1); for(int i=0;i<deg;i++) c[i]=mul(ta[i],tb[i]); ntt(c,deg,-1); } int k; void clr(int *a) { for(int i=0;i<=2*k;i++) a[i]=0; } int ar[maxn],br[maxn],bri[maxn],Q[maxn],tmp[maxn]; void gmod(int *a,int *b,int *c,int n,int m) {//in>>a%b out<<c a n b m clr(ar); clr(br); clr(bri); clr(Q); clr(tmp); for(int i=0;i<=n;i++) ar[i]=a[n-i]; for(int i=0;i<=m;i++) br[i]=b[m-i]; for(int i=n-m+2;i<=m;i++) br[i]=0; int le=1; for(;le<=(n-m);le<<=1); ginv(br,bri,le); le = 1; for(;le<=2*n;le<<=1); MUL(bri,ar,Q,le); reverse(Q,Q+n-m+1); for(int i=n-m+1;i<le;i++) Q[i]=0; MUL(Q,b,tmp,le); for(int i=0;i<m;i++) c[i]=sub(a[i],tmp[i]); for(int i=m;i<=2*k;i++) c[i]=0; } int n; int xs[maxn],st[maxn],a[maxn]; int sg[maxn],ans[maxn]; int main() { // freopen("orz.in","r",stdin); scanf("%d%d",&n,&k); for(int i=1;i<=k;i++) scanf("%d",&xs[i]),xs[i]%=mod,xs[i]+=(xs[i]<0)?mod:0; for(int i=0;i<k;i++) { scanf("%d",&st[i]),st[i]%=mod,st[i]+=(st[i]<0)?mod:0; } for(int i=1;i<=k;i++) sg[k-i] = mod-xs[i]; sg[k]=1; int le=1; for(;le<=2*k;le<<=1); a[1]=1; ans[0]=1; while(n) { if(n&1) { MUL(ans,a,ans,le<<1); int zz=(k<<1); while(ans[--zz]==0); if(zz>=k) gmod(ans,sg,ans,zz,k); } MUL(a,a,a,le); int zz = (k<<1); while(a[--zz]==0); if(zz>=k)gmod(a,sg,a,zz,k); n>>=1; } int AS=0; for(int i=0;i<k;i++) AS=add(AS,mul(st[i],ans[i])); printf("%d\n",AS); // cout<<ans[0]<<endl; }bzoj 4161 O(k^2logn)
#include<stdio.h> #include<iostream> #include<algorithm> #include<cstdio> #include<cmath> using namespace std; const int maxn = 60006 , mod = 998244353; int k; long long n; int a[maxn],h[maxn]; int mul(int x,int y) { return 1ll*x*y%mod; } int add(int x,int y) { x+=y; return x>=mod?x-mod:x; } int b[maxn],c[maxn],tmp[maxn]; void qmul(int *x,int *y) { for(int i=0;i<k;i++) { for(int j=0;j<k;j++) tmp[i+j] = add(tmp[i+j],mul(x[i],y[j])); } for(int i=k*2-2;i>=k;i--) { for(int j=0;j<=k;j++) tmp[i-k+j] = add(tmp[i-k+j],mul(tmp[i],a[k-j])); tmp[i] = 0; } for(int i=0;i<k;i++) x[i] = tmp[i],tmp[i]=0; } int main() { scanf("%d%lld",&k,&n); n--; for(int i=0;i<k;i++) scanf("%d",&h[i]); for(int i=1;i<=k;i++) scanf("%d",&a[i]); b[0] = 1; c[1] = 1; if(n<k) { printf("%d",h[n]); return 0; } for(;n;n>>=1,qmul(c,c)) if(n&1) qmul(b,c); int ans = 0; for(int i=0;i<k;i++) ans = add(ans,mul(b[i],h[i])); printf("%d",ans); }
多项式多点求值
给定一个n次多项式,现在请你对于,求出。 我们考虑分治处理,对于 我们写成一个类似线段树的结构把插值的多项式存下来。#include<stdio.h> #include<cstdio> #include<algorithm> #include<iostream> #include<cmath> #include<vector> using namespace std; const int maxn = 64005; const int mod = 998244353; const int g = 3; int mul(int x,int y) { return 1ll*x*y%mod; } int add(int x,int y) { x+=y; return x>=mod?x-mod:x; } int sub(int x,int y) { x-=y; return x<0?x+mod:x; } int ksm(int a,int b) { int ans = 1; for(;b;b>>=1,a=mul(a,a)) if(b&1) ans = mul(ans,a); return ans; } void ntt(int *a,int deg,int dft) { for(int i=0,j=0;i<deg;i++) { if(i<j) swap(a[i],a[j]); for(int k=(deg>>1);(j^=k)<k;k>>=1); } for(int st=1;st<deg;st<<=1) { int dwg = (dft==1) ? ksm(g,(mod-1)/(st<<1)) : ksm(g,mod-1-(mod-1)/(st<<1)); for(int i=0;i<deg;i+=(st<<1)) { int ng = 1; for(int j=i;j<i+st;j++) { int x = a[j]; int y = mul(ng,a[j+st]); ng = mul(ng,dwg); a[j] = add(x,y); a[j+st] = add(x,mod-y); } } } if(dft==1) return; int invs = ksm(deg,mod-2); for(int i=0;i<deg;i++) a[i] = mul(invs,a[i]); } int F[18][maxn*4]; int n,m; int ta[maxn*4],tb[maxn*4]; void MULL(int *a,int *b,int *c,int le) { for(int i=0;i<le;i++) ta[i]=b[i]; for(int j=0;j<le;j++) tb[j]=c[j]; ntt(ta,le,1); ntt(tb,le,1); for(int i=0;i<le;i++) ta[i]=mul(ta[i],tb[i]); ntt(ta,le,-1); for(int i=0;i<le;i++) a[i] = ta[i]; } void ginv(int deg,int *a,int *b) { if(deg==1) { b[0]=ksm(a[0],mod-2); return; } ginv(deg>>1,a,b); for(int i=0;i<deg;i++) ta[i]=a[i]; for(int i=deg;i<(deg<<1);i++) ta[i]=0; ntt(ta,deg<<1,1); ntt(b,deg<<1,1); for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i]))); ntt(b,(deg<<1),-1); for(int i=deg;i<(deg<<1);i++) b[i]=0; } int Q[maxn*4],Fr[maxn*4],Gr[maxn*4],Gri[maxn*4],tmp[maxn*4]; void gmod(int *O,int *F,int *G,int n,int m) {//F%G==O (F n G m) for(int i=0;i<=n;i++) Fr[i] = F[n-i]; for(int i=0;i<=m;i++) Gr[i] = G[m-i]; for(int i=n-m+2;i<=m;i++) Gr[i] = 0; int le = 1; for(;le<=(n-m);le<<=1); for(int i=n-m+1;i<=m;i++) Gr[i]=0; ginv(le,Gr,Gri); le = 1; for(;le<=(n<<1);le<<=1); MULL(Q,Gri,Fr,le); reverse(Q,Q+n-m+1); for(int i=n-m+1;i<le;i++) Q[i]=0; MULL(tmp,Q,G,le); for(int i=0;i<m;i++) O[i] = sub(F[i],tmp[i]); for(int i=0;i<le;i++) tmp[i] = Fr[i] = Gr[i] = Gri[i] = Q[i] = 0; } vector<int>ve[maxn*2]; int tot,rt,ls[maxn*2],rs[maxn*2],X[maxn]; int tpa[maxn*4],tpb[maxn*4],tpc[maxn*4]; void maketree(int &p,int l,int r) { p = ++tot; int len = 1; for(;len<r-l+2;len<<=1); ve[p].resize(len); // len ci if(l==r) { ve[p][0] = (mod-(X[l]%mod+mod)%mod); ve[p][1] = 1; return; } int mid = (l+r)>>1; maketree(ls[p],l,mid); maketree(rs[p],mid+1,r); int ss = ve[ls[p]].size(); for(int i=0;i<ss;i++) tpa[i] = ve[ls[p]][i]; for(int i=ss;i<len;i++) tpa[i]=0; ss = ve[rs[p]].size(); for(int i=0;i<ss;i++) tpb[i] = ve[rs[p]][i]; for(int i=ss;i<len;i++) tpb[i]=0; MULL(tpc,tpa,tpb,len); for(int i=0;i<len;i++) ve[p][i] = tpc[i]; } int G[maxn*4]; int ANS[maxn]; void DC(int dep,int p,int l,int r,int cs) { int m = r-l+1; for(int i=0;i<=m;i++) G[i] = ve[p][i]; if(cs>=m) gmod(F[dep],F[dep-1],G,cs,m); else { for(int i=0;i<=m-1;i++) F[dep][i] = F[dep-1][i]; } for(int i=0;i<=m;i++) G[i] = 0; if(l==r) { ANS[l] = F[dep][0]; return; } int mid = (l+r)>>1; DC(dep+1,ls[p],l,mid,m-1); DC(dep+1,rs[p],mid+1,r,m-1); } int main() { scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) { int x; scanf("%d",&x); F[0][i] = x; } for(int i=1;i<=m;i++) scanf("%d",&X[i]); maketree(rt,1,m); for(int i=0;i<=m;i++) G[i] = ve[rt][i]; DC(1,rt,1,m,n); for(int i=1;i<=m;i++) printf("%d\n",ANS[i]); }
多项式多点插值
给定n个点(x1,y1),(x2,y2),(x3,y3)......(xn,yn)求一个n-1次多项式F满足F(xi) = yi 如果我们直接考虑拉格朗日插值就会得到一个n^2的做法 [latex] {F(x)=\sum_{i=1}^n{\frac{\prod_{{j}\neq{i}}{({x}-{x_j})}}{\prod_{{j}\neq{i}}{({x_i}-{x_j})}}{y_i}}} [/latex]多项式求复合逆
预备:拉格朗日反演
若两个多项式,都没有常数项,有A(B(x))=x,已知A,称B为A的复合逆。显然若A(B(x))=x则B(A(x)),且一定满足一次项系数互为逆元。 那么如果已知A要求B的n次项系数,有 [latex][x^n]B(x)=\frac{1}{n} [x^{n-1}] (\frac{x}{A(x)})^{n}[/latex] 同时拉格朗日反演有扩展有 [latex][x^n]A(B(x))=\frac{1}{n}[x^{n-1}]A'(x)(\frac{x}{B(x)})^n[/latex] 就这样我们就可以在nlogn的时间复杂度内完成求解多项式复合逆的某一项系数了任意模数NTT
由于某些毒瘤出题人出的模数并不棒,想毒瘤一下大家,没法NTT,FFT又精度上有点出事。多模数NTT
一般对于长度为1e5,数范围在1e9的卷积,可以想到其卷积之后范围应该是大约最大1e23(1e51e91e9)的。 因此我们可以考虑三个1e9左右的模数来NTT之后利用中国剩余定理合并。 把这三个模数背下来吧,因为他们很优秀——原根都是3. 998244353,1004535809,469762049(1<<26)*7,(1<<23)*119,(1<<21)*479
拆系数FFT
多模数NTT太慢了,那么,快乐的拆系数FFT来啦! 我们把每个数拆成aM+b的形式(M是常数),这样的话,就可以把一个A×B 转化为 (Aa + Ab) × (Ba + Bb)。然后分为三个区域AaBa , AbBa + AaBb,Aa*Bb分别对应乘M^2,乘M到M^2,不乘任何数。 一般为了方便,选择取M为(1<<15)=32768.缺点是空间消耗比较爆炸,在一些情况下精度也并不友好。 手写复数套装code:struct cd{ db x,y; }fa[maxm],fb[maxm]; cd operator+(cd aa,cd bb) { return (cd){aa.x+bb.x,aa.y+bb.y}; } cd operator-(cd aa,cd bb) { return (cd){aa.x-bb.x,aa.y-bb.y}; } cd operator*(cd aa,cd bb) { return (cd){aa.x*bb.x-aa.y*bb.y,aa.x*bb.y+aa.y*bb.x}; } cd operator/(cd aa,db bb) { return (cd){aa.x/bb,aa.y/bb}; } void fft(cd *a,int n,int dft) { for(int j=0,i=0;i<n;i++) { if(i<j) swap(a[i],a[j]); for(int k=(n>>1);((j^=k)<k);k>>=1); } for(int st=1;st<n;st<<=1) { cd dwfg=(cd){cos(dft*pi/st),sin(dft*pi/st)}; for(int i=0;i<n;i+=(st<<1)) { cd nfg=(cd){1,0}; for(int j=i;j<i+st;j++) { cd x=a[j],y=a[j+st]*nfg; nfg=nfg*dwfg; a[j]=x+y; a[j+st]=x-y; } } } if(dft==-1) for(int i=0;i<n;i++) a[i]=a[i]/n; } int fta[maxm],ftb[maxm],ftc[maxm],ftd[maxm]; int oz1[maxm],oz2[maxm],oz3[maxm],oz4[maxm]; void MUL(int *c,int *a,int *b,int deg) { for(int i=0;i<deg;i++) fa[i]=(cd){1.0*a[i],0},fb[i]=(cd){1.0*b[i],0}; fft(fa,deg,1); fft(fb,deg,1); for(int i=0;i<deg;i++) fa[i]=fa[i]*fb[i]; fft(fa,deg,-1); for(int i=0;i<deg;i++) c[i]=( (long long)(fa[i].x + 0.5) )%mod; } void anymul(int *a,int *b,int *c,int deg) { // c = a*b for(int i=0;i<deg;i++) { fta[i] = a[i]>>10; ftb[i] = a[i]&((1<<10)-1); ftc[i] = b[i]>>10; ftd[i] = b[i]&((1<<10)-1); } MUL(oz1,ftb,ftd,deg);MUL(oz2,ftb,ftc,deg); MUL(oz3,ftd,fta,deg);MUL(oz4,fta,ftc,deg); for(int i=0;i<deg;i++) c[i]=add(add(oz1[i],mul(add(oz2[i],oz3[i]),1<<10)),mul(oz4[i],1<<20)); }