题目大意
一个序列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);
}