luogu P4726 多項式指數函數(模板題FFT、多項式求逆、多項式對數函數）

手動博客搬家: 本文發表于20181127 08:39:42, 原地址https://blog.csdn.net/suncongbo/article/details/84559818

題目鏈接: https://www.luogu.org/problem/show?pid=4726

題意: 給定\(n\)次多項式\(A(x)\) 求多項式\(f(x)\)滿足\(f(x)\equiv e^{A(x)} (\mod x^n)\)

題解

這個比對數函數復雜一些。。
前鋪知識
泰勒展開
對于一個函數，我們可以用以下方式用它的高階導數進行無窮逼近:\[f(x)=f(a)+f'(a)(x-a)+f''(a)\frac{(x-a)^2}{2!}+f'''(a)\frac{(x-a)^3}{3!}+...\]
牛頓迭代
已知要求的多項式\(f\)滿足\(g(f(x))\equiv 0(\mod x^n)\) \(g\)已知，那么可以通過如下的方式倍增求解:
假設求出了原問題\(\mod x^n\)的答案\(f_0(x)\), 考慮轉移到\(f(x) \mod x^{2n}\).
根據泰勒展開公式: \[0=g(f)=g(f_0)+g'(f_0)(f-f_0)+g''(f_0)\frac{(f-f_0)^2}{2!}+...\]
但是，由于一個顯然的結論——\(f\equiv f_0(\mod x^n)\), 因此\((f-f_0)^2\equiv 0(\mod x^{2n})\), 因此公式里只需要計算前兩項，后面的項都在\(\mod x^{2n}\)意義下為\(0\)!
即\(g(f_0)+g'(f_0)(f-f_0)=0\)
\(f=f_0-\frac{g(f_0)}{g'(f_0)}\)
如此即可求解。
（太神了啊太神了啊嗚嗚嗚。。我大概永遠也搞不出這么神奇的東西吧。。）

本題題解
根據牛頓迭代的法則，令\(g(f)=\ln f-A\), 則\(f=f_0-\frac{\ln f_0-A}{\frac{1}{f_0}}=f_0(1-\ln f_0+A)\)
遞歸求解即可。
值得注意的是FFT的大小
\(A\)應該帶入\(2n\)次, \(\ln f_0\)應該帶入\(2n\)次, \(f_0\)應該帶入\(n\)次, FFT乘法因為是\(2n\)次乘以\(n\)次，所以應該取到\(4n\)個單位根。
FFT這個地方太容易出錯了！范圍大了常數大好幾倍，范圍小了就會出錯。
時間復雜度為\(T(n)=T(\frac{n}{2})+O(n\log n)\), \(T(n)=O(n\log n)\)
常數我算的應該是\(48\)倍。每次\(n\)到\(2n\)的迭代需要做一次\(2n\)的\(\ln\)和三次\(4n\)的DFT/IDFT. 因此\(18(2n\log n)+3(4n\log n)=48n\log n\).
(飛了……)

代碼

#include<cstdio> #include<cstdlib> #include<cstring> #include<algorithm> #define llong long long #define ldouble long double #define uint unsigned int #define ullong unsigned long long #define udouble unsigned double #define uldouble unsigned long double #define modinc(x) {if(x>=P) x-=P;} #define pii pair<int,int> #define piii pair<pair<int,int>,int> #define piiii pair<pair<int,int>,pair<int,int> > #define pli pair<llong,int> #define pll pair<llong,llong> #define Memset(a,x) {memset(a,x,sizeof(a));} using namespace std;const int N = 1<<19; const int LGN = 19; const int G = 3; const int P = 998244353; llong tmp1[N+3],tmp2[N+3],tmp3[N+3],tmp4[N+3]; //inv llong tmp5[N+3],tmp6[N+3],tmp7[N+3],tmp8[N+3]; //ln llong tmp9[N+3],tmp10[N+3],tmp11[N+3],tmp12[N+3]; //exp llong a[N+3],b[N+3]; int id[N+3]; int n;llong quickpow(llong x,llong y) {llong cur = x,ret = 1ll;for(int i=0; y; i++){if(y&(1ll<<i)){ret = ret*cur%P;y-=(1ll<<i);}cur = cur*cur%P;}return ret; } llong mulinv(llong x) {return quickpow(x,P-2);}void initid(int dgr) {int len = 0;for(int i=0; i<=LGN; i++) if(dgr==(1<<i)) {len = i; break;}id[0] = 0;for(int i=1; i<dgr; i++) id[i] = (id[i>>1]>>1)|((i&1)<<(len-1)); }void ntt(int dgr,int coe,llong poly[],llong ret[]) {initid(dgr);for(int i=0; i<dgr; i++) ret[i] = poly[i];for(int i=0; i<dgr; i++) if(i<id[i]) swap(ret[i],ret[id[i]]);for(int i=1; i<=(dgr>>1); i<<=1){llong tmp = quickpow(G,(P-1)/(i<<1));if(coe==-1) tmp = mulinv(tmp);for(int j=0; j<dgr; j+=(i<<1)){llong expn = 1ll;for(int k=0; k<i; k++){llong x = ret[j+k],y = ret[j+k+i]*expn%P;ret[j+k] = x+y; modinc(ret[j+k]);ret[j+i+k] = x-y+P; modinc(ret[j+i+k]);expn = (expn*tmp)%P;}}}if(coe==-1){llong tmp = mulinv(dgr);for(int j=0; j<dgr; j++) ret[j] = ret[j]*tmp%P;} }void polyinv(int dgr,llong poly[],llong ret[]) {for(int i=0; i<dgr; i++) ret[i] = 0ll;ret[0] = mulinv(poly[0]);for(int i=1; i<=(dgr>>1); i<<=1){for(int j=0; j<(i<<2); j++) tmp1[j] = j<i ? ret[j] : 0ll;for(int j=0; j<(i<<2); j++) tmp2[j] = j<(i<<1) ? poly[j] : 0ll;ntt((i<<2),1,tmp1,tmp3); ntt((i<<2),1,tmp2,tmp4);for(int j=0; j<(i<<2); j++) tmp3[j] = tmp3[j]*tmp3[j]%P*tmp4[j]%P;ntt((i<<2),-1,tmp3,tmp4);for(int j=0; j<(i<<1); j++) ret[j] = (tmp1[j]+tmp1[j]-tmp4[j]+P)%P;}for(int i=dgr; i<(dgr<<1); i++) ret[i] = 0ll; }void polyder(int dgr,llong poly[],llong ret[]) {for(int i=0; i<dgr-1; i++) ret[i] = poly[i+1]*(i+1)%P; }void polyint(int dgr,llong poly[],llong ret[]) {for(int i=1; i<dgr; i++) ret[i] = poly[i-1]*mulinv(i)%P; }void polyln(int dgr,llong poly[],llong ret[]) {for(int i=0; i<dgr; i++) ret[i] = 0ll;polyder(dgr,poly,tmp5);polyinv(dgr,poly,tmp6);ntt((dgr<<1),1,tmp5,tmp7); ntt((dgr<<1),1,tmp6,tmp8);for(int i=0; i<(dgr<<1); i++) tmp7[i] = tmp7[i]*tmp8[i]%P;ntt((dgr<<1),-1,tmp7,tmp8);polyint(dgr,tmp8,ret); }void polyexp(int dgr,llong poly[],llong ret[]) {for(int i=0; i<dgr; i++) ret[i] = i==0;for(int i=1; i<=(dgr>>1); i<<=1){for(int j=0; j<(i<<2); j++) tmp9[j] = j>=(i<<1) ? 0ll : ((j==0)+poly[j])%P;for(int j=0; j<(i<<2); j++) tmp10[j] = j>=i ? 0ll : ret[j];polyln((i<<1),tmp10,tmp11);for(int j=0; j<(i<<1); j++) {tmp9[j] = tmp9[j]-tmp11[j]+P; modinc(tmp9[j]);}ntt((i<<2),1,tmp10,tmp12); ntt((i<<2),1,tmp9,tmp11);for(int j=0; j<(i<<2); j++) tmp12[j] = tmp12[j]*tmp11[j]%P;ntt((i<<2),-1,tmp12,tmp11);for(int j=0; j<(i<<1); j++) ret[j] = tmp11[j];} }int main() {scanf("%d",&n); int dgr = 1; while(dgr<=n) dgr<<=1;for(int i=0; i<n; i++) scanf("%lld",&a[i]);polyexp(dgr,a,b);for(int i=0; i<n; i++) printf("%lld ",b[i]);return 0; } 發表于 2019-01-23 20:22 suncongbo 閱讀(...) 評論(...) 編輯 收藏 刷新評論刷新頁面返回頂部

總結

以上是生活随笔為你收集整理的luogu P4726 多项式指数函数(模板题FFT、多项式求逆、多项式对数函数）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。