[分治FFT] LOJ#6183. 看无可看

推一推 $f_n$ 的通项

$f_n=a\times 3^n-b\times (-1)^n$

$a$ 和 $b$ 可以用 $f_0$ 和 $f_1$ 解出来

那么只要能求 $\sum_{s'\subseteq s,|s|=k} w^{\sum_{x\in s'}x}$ 也就是 $\sum_{s'\subseteq s,|s|=k} \prod_{x\in 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;
}