推一推 fn f n 的通项
fn=a×3n−b×(−1)n f n = a × 3 n − b × ( − 1 ) n
a a 和 可以用 f0 f 0 和 f1 f 1 解出来
那么只要能求 ∑s′⊆s,|s|=kw∑x∈s′x ∑ s ′ ⊆ s , | s | = k w ∑ x ∈ s ′ x 也就是 ∑s′⊆s,|s|=k∏x∈s′wx ∑ s ′ ⊆ s , | s | = k ∏ x ∈ s ′ w x 就行了
分治FFT
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
#define show(a) \
cerr<<#a<<" - size : "<<a.size()<<endl;\
cerr<<'\t';\
for(int i=0;i<a.size();i++) cerr<<a[i]<<' ';\
cerr<<'\n'
using namespace std;
const int N=500010,P=99991;
const double pi=acos(-1);
struct E{
double r,i;
E(double _r=0,double _i=0):r(_r),i(_i){}
friend E operator +(const E &a,const E &b){
return E(a.r+b.r,a.i+b.i);
}
friend E operator -(const E &a,const E &b){
return E(a.r-b.r,a.i-b.i);
}
friend E operator *(const E &a,const E &b){
return E(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
}
};
typedef vector<int> poly;
inline char nc(){
static char buf[100000],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}
inline void read(int &x){
char c=nc(); x=0;
for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());
}
inline int Pow(int x,int y){
int ret=1; y%=P-1;
for(;y;y>>=1,x=1LL*x*x%P) if(y&1) ret=1LL*ret*x%P;
return ret;
}
int n,k,num,a[N],f0,f1,aa,bb,rev[N];
E w[2][N];
inline void FFT(E *a,int n,int r){
for(int i=1;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1)
for(int j=0;j<n;j+=(i<<1))
for(int k=0;k<i;k++){
E x=a[j+k],y=a[j+k+i]*w[r][num/(i<<1)*k];
a[j+k]=x+y; a[j+k+i]=x-y;
}
if(!r) for(int i=0;i<n;i++) a[i].r/=n;
}
E tmp[N];
poly operator *(poly a,poly b){
//show(a); show(b);
poly ret; ret.resize(a.size()+b.size()-1);
if(ret.size()<=100){
for(int i=0;i<a.size();i++)
for(int j=0;j<b.size();j++)
ret[i+j]=(ret[i+j]+1LL*a[i]*b[j])%P;
return ret;
}
int n,L=0; for(n=1;n<=ret.size();n<<=1,L++); n<<=1;
for(int i=1;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
for(int i=0;i<a.size();i++) tmp[i]=E(a[i]*0.5,a[i]*0.5);
for(int i=a.size();i<n;i++) tmp[i]=0;
for(int i=0;i<b.size();i++) tmp[i]=tmp[i]+E(b[i]*0.5,-b[i]*0.5);
FFT(tmp,n,1);
for(int i=0;i<n;i++) tmp[i]=tmp[i]*tmp[i];
FFT(tmp,n,0);
for(int i=0;i<ret.size();i++) ret[i]=(long long)(tmp[i].r+0.5)%P;
//show(ret);
return ret;
}
poly solve(int l,int r,int b){
if(l==r){
poly ret; ret.resize(2);
ret[0]=1; ret[1]=Pow(b,a[l]);
//show(ret);
return ret;
}
int mid=l+r>>1;
poly ret=solve(l,mid,b)*solve(mid+1,r,b);
//show(ret);
return ret;
}
inline void Pre(const int &n){
num=n;
w[0][0]=w[1][0]=E(1,0);
for(int i=1;i<num;i++)
w[1][i]=E(cos(2*pi*i/num),sin(2*pi*i/num));
for(int i=1;i<num;i++)
w[0][i]=w[1][num-i];
}
int main(){
read(n); read(k);
for(int i=1;i<=n;i++) read(a[i]);
int m; for(m=1;m<=n+2;m<<=1); Pre(m<<1);
poly A=solve(1,n,3),B=solve(1,n,P-1);
//show(A); show(B);
read(f0); read(f1);
aa=24998LL*(f0+f1)%P; bb=(aa+P-f0)%P;
int ans=(1LL*aa*A[k]%P+P-1LL*bb*B[k]%P)%P;
printf("%d\n",ans);
return 0;
}