多项式 ln exp
单独写篇博客...
用的vector写板子，比数组慢一倍左右，安全性较高
多项式ln
[ln (F(x))=G(x)
]
[G'(x) = frac {dln (F(x))} {dF(x)} frac{d(F(X))}{dx}=frac{F'(x)}{F(x)}
]
[G(x)=intfrac{F'(x)}{F(x)} dx
]
#include<bits/stdc++.h>
using namespace std;
int read(){
int x=0,pos=1;char ch=getchar();for(;!isdigit(ch);ch=getchar()) if(ch=='-') pos=0;
for(;isdigit(ch);ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';return pos?x:-x;
}
#define FOR(i,a,b) for(int i=a;i<=b;++i)
const int N = 200201;
#define ROF(i,a,b) for(int i=b;i>=a;--i)
int n,m;
int a[N*2];
typedef vector<int> poly;poly b[N];
const int mod = 998244353;
const int LSZ=20,SZ=1<<LSZ;
int ksm(int a,int b){
int res=1;
while(b){
if(b&1) res=1ll*res*a%mod;
a=1ll*a*a%mod,b>>=1;
}return res;
}
int tr[SZ],omg[SZ][2];
void init(int n){
int lim=0;while((1<<lim)<n) lim++;
FOR(i,0,n-1) tr[i]=tr[i>>1]>>1|((i&1)<<(lim-1));
omg[0][0]=omg[0][1]=1;
omg[1][0]=ksm(3,(mod-1)/n);
omg[1][1]=ksm(omg[1][0],mod-2);
FOR(i,2,n) omg[i][0]=1ll*omg[i-1][0]*omg[1][0]%mod,omg[i][1]=1ll*omg[i-1][1]*omg[1][1]%mod;
}
void dft(poly &f,int n,int opt){
//assert((int)f.size()<=n);
f.resize(n,0);
for(int i=0;i<n;++i) if(i<tr[i]) swap(f[i],f[tr[i]]);
for(int l=2;l<=n;l<<=1){ int m=l/2;
for(int g=0;g!=n;g+=l){
for(int i=0;i<m;i++){
int t=1ll*omg[n/l*i][opt]*(f[g+i+m])%mod;
f[g+i+m]=(f[g+i]-t+mod)%mod;f[g+i]=(f[g+i]+t)%mod;
}
}
}
if(opt) for(int i=0,iv=ksm(n,mod-2);i<n;i++) f[i]=1ll*f[i]*iv%mod;
}
void pop(poly &a){
if(!a.size()) return;
while(a[a.size()-1]==0) a.pop_back();
}
void dao(poly& a){
for(int i=1;i<a.size();i++) a[i-1]=1ll*a[i]*i%mod;a.pop_back();
}
void jifen(poly &a){
for(int i=a.size()-1;i>=0;--i) a[i]=1ll*a[i-1]*ksm(i,mod-2)%mod;a[0]=0;
}
poly operator+(const poly &a,const poly &b){
poly res(a);
res.resize(max(a.size(),b.size()),0);
for(int i=0;i<b.size();++i) res[i]=(res[i]+b[i])%mod;
return res;pop(res);
}
poly operator*(const poly &a,const poly &b){
poly A=a,B=b;//pop(A),pop(B);
int len=1;while(len<A.size()+B.size()) len*=2;init(len);
poly C(len,0);
dft(A,len,0),dft(B,len,0);
FOR(i,0,len-1) C[i]=1ll*A[i]*B[i]%mod;
dft(C,len,1);pop(C);return C;
}
void inv(const poly &a,poly &b,int n){
if(n==1) return b[0]=ksm(a[0],mod-2),void();
inv(a,b,(n+1)/2);
int len=1;while(len<n*2) len*=2;init(len);
if(b.size()<len) b.resize(len,0);FOR(i,(n+1)/2,len-1) b[i]=0;
poly c=a;if(c.size()<len) c.resize(len,0);FOR(i,n,len-1) c[i]=0;
dft(c,len,0);dft(b,len,0);
FOR(i,0,len-1) b[i]=(2ll-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
dft(b,len,1);
}
void getln(const poly &a,poly &ans,int n){
poly b;b.resize(n,0);
inv(a,b,n);pop(b);
poly c=a;dao(c);
ans=b*c;
jifen(ans);pop(ans);
}
int main(){
//freopen("P4725_8.in","r",stdin);
//freopen("8.txt","w",stdout);
n=read();poly a;poly ans;
a.resize(n,0);ans.resize(n,0);
for(int i=0;i<n;i++){
a[i]=read();
}
getln(a,ans,n);
for(int i=0;i<n;i++){
printf("%lld ",ans[i]);
}
return 0;
}
多项式exp
[G(x)=e^{F(x)}
]
[ln(G(x))=F(x)
]
[H(G(x))=ln(G(x))-F(x)=0
]
由牛顿迭代，
[G(x)=G_0(x)-frac{ln(G_0(x))-F(x)}{G_0(x)^{-1}}
]
[=G_0(x)(1-ln(G_0(x))+F(x))
]
边界，(G(x)[x^0]=e^{F(x)[X^0]})
#include<bits/stdc++.h>
using namespace std;
int read(){
int x=0,pos=1;char ch=getchar();for(;!isdigit(ch);ch=getchar()) if(ch=='-') pos=0;
for(;isdigit(ch);ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';return pos?x:-x;
}
#define FOR(i,a,b) for(int i=a;i<=b;++i)
const int N = 200201;
#define ROF(i,a,b) for(int i=b;i>=a;--i)
int n,m;
int a[N*2];
typedef vector<int> poly;poly b[N];
const int mod = 998244353;
const int LSZ=20,SZ=1<<LSZ;
int ksm(int a,int b){
int res=1;
while(b){
if(b&1) res=1ll*res*a%mod;
a=1ll*a*a%mod,b>>=1;
}return res;
}
int tr[SZ],omg[SZ][2];
void init(int n){
int lim=0;while((1<<lim)<n) lim++;
FOR(i,0,n-1) tr[i]=tr[i>>1]>>1|((i&1)<<(lim-1));
omg[0][0]=omg[0][1]=1;
omg[1][0]=ksm(3,(mod-1)/n);
omg[1][1]=ksm(omg[1][0],mod-2);
FOR(i,2,n) omg[i][0]=1ll*omg[i-1][0]*omg[1][0]%mod,omg[i][1]=1ll*omg[i-1][1]*omg[1][1]%mod;
}
void dft(poly &f,int n,int opt){
//assert((int)f.size()<=n);
f.resize(n,0);
for(int i=0;i<n;++i) if(i<tr[i]) swap(f[i],f[tr[i]]);
for(int l=2;l<=n;l<<=1){ int m=l/2;
for(int g=0;g!=n;g+=l){
for(int i=0;i<m;i++){
int t=1ll*omg[n/l*i][opt]*(f[g+i+m])%mod;
f[g+i+m]=(f[g+i]-t+mod)%mod;f[g+i]=(f[g+i]+t)%mod;
}
}
}
if(opt) for(int i=0,iv=ksm(n,mod-2);i<n;i++) f[i]=1ll*f[i]*iv%mod;
}
void pop(poly &a){
if(!a.size()) return;
while(a[a.size()-1]==0) a.pop_back();
}
void dao(poly& a){
for(int i=1;i<a.size();i++) a[i-1]=1ll*a[i]*i%mod;a.pop_back();
}
void jifen(poly &a){
if(!a.size()) return; for(int i=a.size()-1;i>=0;--i) a[i]=1ll*a[i-1]*ksm(i,mod-2)%mod;a[0]=0;
}
poly operator+(const poly &a,const poly &b){
poly res(a);
res.resize(max(a.size(),b.size()),0);
for(int i=0;i<b.size();++i) res[i]=(res[i]+b[i])%mod;
return res;pop(res);
}
poly operator*(const poly &a,const poly &b){
poly A=a,B=b;//pop(A),pop(B);
int len=1;while(len<A.size()+B.size()) len*=2;init(len);
poly C(len,0);
dft(A,len,0),dft(B,len,0);
FOR(i,0,len-1) C[i]=1ll*A[i]*B[i]%mod;
dft(C,len,1);pop(C);return C;
}
void inv(const poly &a,poly &b,int n){
if(n==1) return b[0]=ksm(a[0],mod-2),void();
inv(a,b,(n+1)/2);
int len=1;while(len<n*2) len*=2;init(len);
if(b.size()<len) b.resize(len,0);FOR(i,(n+1)/2,len-1) b[i]=0;
poly c=a;if(c.size()<len) c.resize(len,0);FOR(i,n,len-1) c[i]=0;
dft(c,len,0);dft(b,len,0);
FOR(i,0,len-1) b[i]=(2ll-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
dft(b,len,1);
}
void getln(const poly &a,poly &ans,int n){
poly b;b.resize(n,0);
inv(a,b,n);pop(b);
poly c=a;dao(c);
ans=b*c;
jifen(ans);
}
void getexp(const poly &a,poly &ans,int n){
if(n==1) return ans[0]=1,void();
getexp(a,ans,(n+1)/2);
ans.resize(n,0);
poly tmp;getln(ans,tmp,n);tmp.resize(n,0);
tmp[0]=(1ll+a[0]-tmp[0]+mod)%mod;
FOR(i,1,n-1) tmp[i]=(0ll+a[i]-tmp[i]+mod)%mod;
ans=ans*tmp;
ans.resize(n,0);
}
int main(){
n=read();poly a;poly ans;
a.resize(n,0);ans.resize(n,0);
for(int i=0;i<n;i++){
a[i]=read();
}
getexp(a,ans,n);
for(int i=0;i<n;i++){
printf("%d ",ans[i]);
}
return 0;
}
练习
以后有其他的这类题也会放在这里。
集训队作业2013 城市规划
考虑ln和exp的本质。
[G(F(x))=e^{F(x)} =sum_{ige 0} frac{F^i(x)}{i!}
]
也就是说，G(F(x))是i个F(x)的EGF乘起来的和，也就是把G的方案分成了i个F的方案。
本题中，G表示无向图个数的EGF，F表示无向连通图的EGF，显然满足上式。G是很好求的，F直接ln即可。
那么要把i个F组合成G，相当于反过来，直接exp。这两个技巧可用于许多计数题。
#include<bits/stdc++.h>
using namespace std;
int read(){
int x=0,pos=1;char ch=getchar();for(;!isdigit(ch);ch=getchar()) if(ch=='-') pos=0;
for(;isdigit(ch);ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';return pos?x:-x;
}
#define FOR(i,a,b) for(int i=a;i<=b;++i)
const int N = 200201;
#define ROF(i,a,b) for(int i=a;i>=b;--i)
int n,m;
typedef vector<int> poly;
const int mod = 1004535809;
const int LSZ=19,SZ=1<<LSZ;
int ksm(int a,int b){
int res=1;
while(b){
if(b&1) res=1ll*res*a%mod;
a=1ll*a*a%mod,b>>=1;
}return res;
}
int tr[SZ],omg[SZ][2];
void init(int n){
int lim=0;while((1<<lim)<n) lim++;
FOR(i,0,n-1) tr[i]=tr[i>>1]>>1|((i&1)<<(lim-1));
omg[0][0]=omg[0][1]=1;
omg[1][0]=ksm(3,(mod-1)/n);
omg[1][1]=ksm(omg[1][0],mod-2);
FOR(i,2,n) omg[i][0]=1ll*omg[i-1][0]*omg[1][0]%mod,omg[i][1]=1ll*omg[i-1][1]*omg[1][1]%mod;
}
void dft(poly &f,int n,int opt){
//assert((int)f.size()<=n);
f.resize(n,0);
for(int i=0;i<n;++i) if(i<tr[i]) swap(f[i],f[tr[i]]);
for(int l=2;l<=n;l<<=1){ int m=l/2;
for(int g=0;g!=n;g+=l){
for(int i=0;i<m;i++){
int t=1ll*omg[n/l*i][opt]*(f[g+i+m])%mod;
f[g+i+m]=(f[g+i]-t+mod)%mod;f[g+i]=(f[g+i]+t)%mod;
}
}
}
if(opt) for(int i=0,iv=ksm(n,mod-2);i<n;i++) f[i]=1ll*f[i]*iv%mod;
}
void dao(poly& a){
for(int i=1;i<a.size();i++) a[i-1]=1ll*a[i]*i%mod;a.pop_back();
}
void jifen(poly &a){
if(!a.size()) return; for(int i=a.size()-1;i>=0;--i) a[i]=1ll*a[i-1]*ksm(i,mod-2)%mod;a[0]=0;
}
poly operator+(const poly &a,const poly &b){
poly res(a);
res.resize(max(a.size(),b.size()),0);
for(int i=0;i<b.size();++i) res[i]=(res[i]+b[i])%mod;
return res;
}
poly operator*(const poly &a,const poly &b){
poly A=a,B=b;//pop(A),pop(B);
int len=1;while(len<A.size()+B.size()) len*=2;init(len);
poly C(len,0);
dft(A,len,0),dft(B,len,0);
FOR(i,0,len-1) C[i]=1ll*A[i]*B[i]%mod;
dft(C,len,1);return C;
}
void inv(const poly &a,poly &b,int n){
if(n==1) return b[0]=ksm(a[0],mod-2),void();
inv(a,b,(n+1)/2);
int len=1;while(len<n*2) len*=2;init(len);
if(b.size()<len) b.resize(len,0);FOR(i,(n+1)/2,len-1) b[i]=0;
poly c=a;if(c.size()<len) c.resize(len,0);FOR(i,n,len-1) c[i]=0;
dft(c,len,0);dft(b,len,0);
FOR(i,0,len-1) b[i]=(2ll-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
dft(b,len,1);
}
void getln(const poly &a,poly &ans,int n){
poly b;b.resize(n,0);
inv(a,b,n);
poly c=a;dao(c);
ans=b*c;
jifen(ans);
}
void getexp(const poly &a,poly &ans,int n){
if(n==1) return ans[0]=1,void();
getexp(a,ans,(n+1)/2);
ans.resize(n,0);
poly tmp;getln(ans,tmp,n);tmp.resize(n,0);
tmp[0]=(1ll+a[0]-tmp[0]+mod)%mod;
FOR(i,1,n-1) tmp[i]=(0ll+a[i]-tmp[i]+mod)%mod;
ans=ans*tmp;
ans.resize(n,0);
}
int f[N],g[N];
void init1(){
f[0]=f[1]=1;
FOR(i,2,n) f[i]=1ll*f[i-1]*i%mod;
g[n]=ksm(f[n],mod-2);ROF(i,n-1,0) g[i]=1ll*g[i+1]*(i+1)%mod;
}
int main(){
n=read();poly a;a.resize(n+1,0);init1();
int iv2=ksm(2,mod-2);
a[0]=a[1]=1;
for(int i=2;i<=n;i++){
a[i]=1ll*ksm(2,1ll*i*(i-1)/2%(mod-1))*g[i]%mod;
}
poly ans;getln(a,ans,n+1);
printf("%lld",1ll*ans[n]*f[n]%mod);
return 0;
}