看无可看

题目大意

一个序列a[i]。
选择任意k个位置，记和为S。
为答案贡献f(S)。
f是类斐波那契数列，对于i>1它满足f(i)=2*f(i-1)+3*f(i-2)。

做法

这是个常系数线性齐次递推。
对于本题M(x)=x^2-2x-3。
定义(a,b)=ax+b，那么可以重载乘法（乘法在模M(x)意义下）：
(a,b)*(c,d)=(bc+ad+2ac,bd+3ac)
对于每个a[i]初始可以用快速幂求x^a[i] mod M(x)。
任选的k个位置把二元组相乘，就能够得到结果。
因此考虑如何处理出任意k个位置乘积和的结果。
经典问题考虑分治FFT。
递归写法空间要一些特技。
然后问题在于怎么对二元组FFT。
只要让单位根是一个虚二元组即可。
也就是w=A+Bi。
A和B都是二元组。
回想正常FFT的单位根，因为这个二元组有实数域（不带x那项），因此w=(0,cos(2pi/len))+(0,sin(2pi/len))i。
然后需要卡常卡精度。。
卡常可以在分治区间较小时选择暴力卷积。
卡精度就是取round时要加一个eps。
大概这样了。

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
typedef long long ll;
typedef double db;
const db pi=acos(-1),eps=0.1;
const int maxlen=400000+10,mo=99991,lim=64;
struct dong{
    db x,y;
    friend dong operator +(dong a,dong b){
        dong c;
        c.x=a.x+b.x;c.y=a.y+b.y;
        return c;
    }
    friend dong operator -(dong a,dong b){
        dong c;
        c.x=a.x-b.x;c.y=a.y-b.y;
        return c;
    }
    friend dong operator *(dong a,dong b){
        dong c;
        c.x=a.y*b.x+a.x*b.y+2*a.x*b.x;c.y=a.y*b.y+3*a.x*b.x;
        //c.x-=floor(c.x/mo)*mo;c.y-=floor(c.y/mo)*mo;
        return c;
    }
};
struct node{
    dong x,y;
    friend node operator +(node a,node b){
        node c;
        c.x=a.x+b.x;c.y=a.y+b.y;
        return c;
    }
    friend node operator -(node a,node b){
        node c;
        c.x=a.x-b.x;c.y=a.y-b.y;
        return c;
    }
    friend node operator *(node a,node b){
        node c;
        c.x=a.x*b.x-a.y*b.y;c.y=a.x*b.y+a.y*b.x;
        return c;
    }
};
node tt[maxlen],a[maxlen],b[maxlen],c[maxlen];
dong A[maxlen],f[5000+10][5000+10];
node F[maxlen];
int la[maxlen],rr[maxlen];
int v[maxlen];
int i,j,k,l,t,n,m,len,ni,ans,mx;
ll x,y;
db ce;
dong zlt,one;
node zero;
int read(){
    int x=0,f=1;
    char ch=getchar();
    while (ch<'0'||ch>'9'){
        if (ch=='-') f=-1;
        ch=getchar();
    }
    while (ch>='0'&&ch<='9'){
        x=x*10+ch-'0';
        ch=getchar();
    }
    return x*f;
}
int qsm(int x,int y){
    if (!y) return 1;
    int t=qsm(x,y/2);
    t=(ll)t*t%mo;
    if (y%2) t=(ll)t*x%mo;
    return t;
}
void quicksortmi(int y){
    if (!y){
        zlt.x=0;
        zlt.y=1;
        return;
    }
    quicksortmi(y/2);
    zlt=zlt*zlt;
    ll j,k;
    j=ll(zlt.x+eps);k=ll(zlt.y+eps);
    j%=mo;k%=mo;
    zlt.x=j;zlt.y=k;
    if (y%2){
        zlt=zlt*one;
        ll j,k;
        j=ll(zlt.x+eps);k=ll(zlt.y+eps);
        j%=mo;k%=mo;
        zlt.x=j;zlt.y=k;
    }
}
void DFT(node *a,int sig){
    int i;
    fo(i,0,len-1){
        int p=0;
        for(int j=0,tp=i;j<ce;j++,tp/=2) p=(p<<1)+(tp%2);
        tt[p]=a[i];
    }
    for(int m=2;m<=len;m*=2){
        int half=m/2;
        fo(i,0,half-1){
            node wi;
            zlt.x=0;
            zlt.y=cos(i*pi*sig/half);
            wi.x=zlt;
            zlt.y=sin(i*pi*sig/half);
            wi.y=zlt;
            for(int j=i;j<len;j+=m){
                node u=tt[j],v=tt[j+half]*wi;
                tt[j]=u+v;
                tt[j+half]=u-v;
            }
        }
    }
    if (sig==-1)
        fo(i,0,len-1){
            zlt=tt[i].x;
            zlt.x/=len;zlt.y/=len;
            tt[i].x=zlt;
            /*zlt=tt[i].x;
            ll j=round(zlt.x+0.00000001),k=round(zlt.y+0.00000001);
            j=(ll)j*ni%mo;
            k=(ll)k*ni%mo;
            zlt.x=j;zlt.y=k;
            tt[i].x=zlt;*/
        }
    fo(i,0,len-1) a[i]=tt[i];
}
void FFT(){
    ce=log(len)/log(2);
    ni=qsm(len,mo-2);
    int i;
    ll j,k;
    DFT(a,1);DFT(b,1);
    fo(i,0,len-1) a[i]=a[i]*b[i];
    DFT(a,-1);
    fo(i,0,len-1){
        zlt=a[i].x;
        j=ll(zlt.x+eps);k=ll(zlt.y+eps);
        j%=mo;k%=mo;
        zlt.x=j;zlt.y=k;
        a[i].x=zlt;
        zlt.x=zlt.y=0;
        a[i].y=zlt;
    }
}
void brute(int p){
    int l1=rr[p*2]-la[p*2]+1,l2=rr[p*2+1]-la[p*2+1]+1;
    int i,j;
    fo(i,0,len-1) c[i]=zero;
    fo(i,0,l1-1)
        fo(j,0,l2-1)
            c[i+j]=c[i+j]+a[i]*b[j];
    ll x,k;
    fo(i,0,len-1){
        zlt=c[i].x;
        x=ll(zlt.x+eps);k=ll(zlt.y+eps);
        x%=mo;k%=mo;
        zlt.x=x;zlt.y=k;
        c[i].x=zlt;
        zlt.x=zlt.y=0;
        c[i].y=zlt;
    }
}
void solve(int p,int l,int r,int L,int R){
    la[p]=L;rr[p]=R;
    if (l==r){
        F[la[p]].x=one;
        F[la[p]+1].x=A[l];
        return;
    }
    int mid=(l+r)/2,M=(L+R)/2;
    solve(p*2,l,mid,L,M);
    solve(p*2+1,mid+1,r,M+1,R);
    int i;
    len=R-L+1;
    fo(i,0,len-1) a[i]=b[i]=zero;
    fo(i,la[p*2],rr[p*2]) a[i-la[p*2]]=F[i];
    fo(i,la[p*2+1],rr[p*2+1]) b[i-la[p*2+1]]=F[i];
    if (len>lim){
        FFT();
        fo(i,L,R) F[i]=a[i-L];
    }
    else{
        brute(p);
        fo(i,L,R) F[i]=c[i-L];
    }
}
int main(){
    freopen("see.in","r",stdin);freopen("see.out","w",stdout);
    n=read();k=read();
    one.x=1;one.y=0;
    fo(i,1,n){
        t=read();
        quicksortmi(t);
        A[i]=zlt;
    }
    one.x=0;one.y=1;
    len=1;
    while (len<=n) len*=2;
    len*=2;
    solve(1,1,n,1,len);
    zlt=F[la[1]+k].x;
    x=ll(zlt.x+eps);y=ll(zlt.y+eps);
    x%=mo;y%=mo;
    j=read();t=read();
    ans=(x*t%mo+y*j%mo)%mo;
    (ans+=mo)%=mo;
    printf("%d\n",ans);
}