題目鏈接：點擊查看

題目大意：給出一個長度為 $n$ 的數列 $a$，現在需要求$$\prod _{b_1<b_2<?<b_k}(a_{b_1}+a_{b_2}+?+a_{b_k})$$

需要滿足：

$1\le b_i\le n$

$b_1<b_2<?<b_k$

$k\ge 1$

簡單題意：需要求 $a$ 數列的所有子序列之和之積

題目分析：一個非常重要的轉換就是，設 $c[k]$ 為子序列之和為 $k$ 的子序列的方案數，那么答案顯然就是 $\prod k^{c[k]}$，又因為題目保證了 $\sum a_i\le 1e5$，也就是 $k$ 的上限不超過 $1e5$，所以我們現在的問題是如何求出數組 $c[k]$ 。

一個很簡單的背包思想就是 $dp[i][j]$ 代表前 $i$ 個數字所選子序列之和為 $j$ 的方案數，$dp[i][j]=dp[i?1][j]+dp[i?1][j?a[i]]$，轉移是 $O(n^2)$ 的，最終 $c[i]=dp[n][i]$。

考慮寫成生成函數的形式，即 $F(x)=(1+x^{a_1})(1+x^{a_2})?(1+x^{a_n})$，那么 $c[i]$ 就是 $F(x)$ 中 $x^i$ 的系數了，暴力求解的話仍然是 $O(n^2)$ 的

考慮分治FFT優(yōu)化，這樣復雜度就下降為 $O(nlo{g}^{2}n)$ 了

需要注意的是，我們最終的答案是需要求解 $\prod k^{c[k]}(\mathrm{m}\mathrm{o}\mathrm{d}p)$，而這里的 $c[k]$ 在計算過程中也是需要取模的。因為本題的 $p=998244353$ 是質數，根據費馬小定理降冪，得到 $\prod k^{c[k]}(\mathrm{m}\mathrm{o}\mathrm{d}p)=\prod k^{c[k](\mathrm{m}\mathrm{o}\mathrm{d}p?1)}(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 推得 $c$ 數組是需要對 $p?1$ 取模的，但這并不是一般的模數，所以不能用 $NTT$，只能用 $MTT$ 去實現

代碼也是參考了 $std$ 的結構體封裝的動態(tài)分配數組，實現的太巧妙了

2021.8.13 UPDATE

更新了一發(fā) $vector$ 合并的版本，跑的莫名快

代碼：

// Problem: Yiwen with Formula // Contest: Virtual Judge - HDU // URL: https://vjudge.net/problem/HDU-7054 // Memory Limit: 524 MB // Time Limit: 16000 ms // // Powered by CP Editor (https://cpeditor.org)// #pragma GCC optimize(2) // #pragma GCC optimize("Ofast","inline","-ffast-math") // #pragma GCC target("avx,sse2,sse3,sse4,mmx") #include<iostream> #include<cstdio> #include<string> #include<ctime> #include<cmath> #include<cstring> #include<algorithm> #include<stack> #include<climits> #include<queue> #include<map> #include<set> #include<sstream> #include<cassert> #include<bitset> #define lowbit(x) x&-x using namespace std; typedef long long LL; typedef unsigned long long ull; template<typename T> inline void read(T &x) {T f=1;x=0;char ch=getchar();while(0==isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(0!=isdigit(ch)) x=(x<<1)+(x<<3)+ch-'0',ch=getchar();x*=f; } template<typename T> inline void write(T x) {if(x<0){x=~(x-1);putchar('-');}if(x>9)write(x/10);putchar(x%10+'0'); } const int inf=0x3f3f3f3f; const int N=1e5+100; const int mod=998244353; const int BASE=1<<15; const long double Pi=acos(-1.0); struct complex {long double x,y;complex (long double xx=0,long double yy=0){x=xx,y=yy;} }P1[N<<2],P2[N<<2],Q[N<<2]; complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);} complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);} complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);} int limit,r[N<<2]; LL q_pow(LL a,LL b) {LL ans=1;while(b) {if(b&1) {ans=ans*a%mod;}a=a*a%mod;b>>=1;}return ans; } void FFT(complex *A,int type) {for(int i=0;i<limit;i++) if(i<r[i]) swap(A[i],A[r[i]]); for(int mid=1;mid<limit;mid<<=1) {complex Wn( cos(Pi/mid) , type*sin(Pi/mid) );for(int R=mid<<1,j=0;j<limit;j+=R) {complex w(1,0);for(int k=0;k<mid;k++,w=w*Wn) {complex x=A[j+k],y=w*A[j+mid+k];A[j+k]=x+y;A[j+mid+k]=x-y;}}} } void init(int n) {limit=1;while(limit<=n) limit<<=1;for(int i=1;i<limit;i++) r[i]=r[i>>1]>>1|((i&1)?limit>>1:0); } int MTT(int *a,int *b,int n,int m,int *res,int mod) {init(n+m);for(int i=0;i<n;i++) {P1[i]={a[i]/BASE,a[i]%BASE};P2[i]={a[i]/BASE,-a[i]%BASE};}for(int i=n;i<limit;i++) {P1[i]={0,0},P2[i]={0,0};}for(int i=0;i<m;i++) {Q[i]={b[i]/BASE,b[i]%BASE};}for(int i=m;i<limit;i++) {Q[i]={0,0};}FFT(P1,1),FFT(P2,1),FFT(Q,1);for(int i=0;i<limit;i++) {Q[i].x/=limit,Q[i].y/=limit;P1[i]=P1[i]*Q[i],P2[i]=P2[i]*Q[i];}FFT(P1,-1),FFT(P2,-1);for(int i=0;i<n+m-1;i++) {long long a1b1,a1b2,a2b1,a2b2;a1b1=(long long)floor((P1[i].x+P2[i].x)/2+0.5)%mod;a1b2=(long long)floor((P1[i].y+P2[i].y)/2+0.5)%mod;a2b1=(long long)floor((P1[i].y-P2[i].y)/2+0.5)%mod;a2b2=(long long)floor((P2[i].x-P1[i].x)/2+0.5)%mod;res[i]=((a1b1*BASE+(a1b2+a2b1))*BASE+a2b2)%mod;res[i]=(res[i]+mod)%mod;}return n+m-1; } int a[N<<2],tot; struct Node {int *p,len;void init(int len) {this->len=len;p=a+tot;for(int i=0;i<len;i++) {p[i]=0;}p[0]++,p[len-1]++;tot+=len;}void mul(const Node& t) {len=MTT(p,t.p,len,t.len,p,mod-1);} }; Node solve(int l,int r) {Node ans;if(l==r) {int x;scanf("%d",&x);ans.init(x+1);} else {int mid=(l+r)>>1;ans=solve(l,mid);ans.mul(solve(mid+1,r));}return ans; } int main() { #ifndef ONLINE_JUDGE // freopen("data.in.txt","r",stdin); // freopen("data.out.txt","w",stdout); #endif // ios::sync_with_stdio(false);int w;cin>>w;while(w--) {tot=0;int n;scanf("%d",&n);Node res=solve(1,n);if(res.p[0]>1) {puts("0");} else {LL ans=1;for(int i=1;i<res.len;i++) {ans=ans*q_pow(i,res.p[i])%mod;}printf("%lld\n",ans);}}return 0; } // Problem: Yiwen with Formula // Contest: Virtual Judge - HDU // URL: https://vjudge.net/problem/HDU-7054 // Memory Limit: 524 MB // Time Limit: 16000 ms // // Powered by CP Editor (https://cpeditor.org)// #pragma GCC optimize(2) // #pragma GCC optimize("Ofast","inline","-ffast-math") // #pragma GCC target("avx,sse2,sse3,sse4,mmx") #include<iostream> #include<cstdio> #include<string> #include<ctime> #include<cmath> #include<cstring> #include<algorithm> #include<stack> #include<climits> #include<queue> #include<map> #include<set> #include<sstream> #include<cassert> #include<bitset> #define lowbit(x) x&-x using namespace std; typedef long long LL; typedef unsigned long long ull; template<typename T> inline void read(T &x) {T f=1;x=0;char ch=getchar();while(0==isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(0!=isdigit(ch)) x=(x<<1)+(x<<3)+ch-'0',ch=getchar();x*=f; } template<typename T> inline void write(T x) {if(x<0){x=~(x-1);putchar('-');}if(x>9)write(x/10);putchar(x%10+'0'); } const int inf=0x3f3f3f3f; const int N=1e5+100; const int mod=998244353; const int BASE=1<<15; const long double Pi=acos(-1.0); struct complex {long double x,y;complex (long double xx=0,long double yy=0){x=xx,y=yy;} }P1[N<<2],P2[N<<2],Q[N<<2]; complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);} complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);} complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);} int limit,r[N<<2]; LL q_pow(LL a,LL b) {LL ans=1;while(b) {if(b&1) {ans=ans*a%mod;}a=a*a%mod;b>>=1;}return ans; } void FFT(complex *A,int type) {for(int i=0;i<limit;i++) if(i<r[i]) swap(A[i],A[r[i]]); for(int mid=1;mid<limit;mid<<=1) {complex Wn( cos(Pi/mid) , type*sin(Pi/mid) );for(int R=mid<<1,j=0;j<limit;j+=R) {complex w(1,0);for(int k=0;k<mid;k++,w=w*Wn) {complex x=A[j+k],y=w*A[j+mid+k];A[j+k]=x+y;A[j+mid+k]=x-y;}}} } void init(int n) {limit=1;while(limit<=n) limit<<=1;for(int i=1;i<limit;i++) r[i]=r[i>>1]>>1|((i&1)?limit>>1:0); } vector<int>MTT(const vector<int>&a,const vector<int>&b) {int n=a.size(),m=b.size();vector<int>res(n+m-1);int mod=::mod-1;init(n+m);for(int i=0;i<n;i++) {P1[i]={a[i]/BASE,a[i]%BASE};P2[i]={a[i]/BASE,-a[i]%BASE};}for(int i=n;i<limit;i++) {P1[i]={0,0},P2[i]={0,0};}for(int i=0;i<m;i++) {Q[i]={b[i]/BASE,b[i]%BASE};}for(int i=m;i<limit;i++) {Q[i]={0,0};}FFT(P1,1),FFT(P2,1),FFT(Q,1);for(int i=0;i<limit;i++) {Q[i].x/=limit,Q[i].y/=limit;P1[i]=P1[i]*Q[i],P2[i]=P2[i]*Q[i];}FFT(P1,-1),FFT(P2,-1);for(int i=0;i<n+m-1;i++) {long long a1b1,a1b2,a2b1,a2b2;a1b1=(long long)floor((P1[i].x+P2[i].x)/2+0.5)%mod;a1b2=(long long)floor((P1[i].y+P2[i].y)/2+0.5)%mod;a2b1=(long long)floor((P1[i].y-P2[i].y)/2+0.5)%mod;a2b2=(long long)floor((P2[i].x-P1[i].x)/2+0.5)%mod;res[i]=((a1b1*BASE+(a1b2+a2b1))*BASE+a2b2)%mod;res[i]=(res[i]+mod)%mod;}return res; } vector<int>solve(int l,int r) {if(l==r) {int x;scanf("%d",&x);vector<int>ans(x+1);ans[0]++,ans[x]++;return ans;}int mid=(l+r)>>1;return MTT(solve(l,mid),solve(mid+1,r)); } int main() { #ifndef ONLINE_JUDGE // freopen("data.in.txt","r",stdin); // freopen("data.out.txt","w",stdout); #endif // ios::sync_with_stdio(false);int w;cin>>w;while(w--) {int n;scanf("%d",&n);vector<int>res=solve(1,n);if(res[0]>1) {puts("0");} else {LL ans=1;for(int i=1;i<(int)res.size();i++) {ans=ans*q_pow(i,res[i])%mod;}printf("%lld\n",ans);}}return 0; }

總結

以上是生活随笔為你收集整理的2021HDU多校7 - 7054 Yiwen with Formula(分治MTT优化dp)的全部內容，希望文章能夠幫你解決所遇到的問題。

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