LOJ #6183 看无可看
LOJ #6183 看无可看
优秀的数学题，Orz samjia2000
这个题显然需要将和转积来处理，这个时候就要用到特征方程的一些知识了！
其实就是这个样子：$f[x]=a\times f[x-1]+b\times f[x-2]$
那么必然可以写作：$f[x]-t\times f[x-1]=k\times (f[x-1]-t\times f[x-2])$
化简：$f[x]=(k+t)\times f[x-1]-k\times t\times f[x-2]$
即：$a=k+t,b=-kt$
消去$K$，得到：$k^2=a\times k+b$
这个东西就是这个二阶递推式的特征方程。
那么显然，你可以得到：$f[x]-k_2\times f[x-1]=k_1^{x-2}\times (f[2]-k_2\times f[1]),f[x]-k_1\times f[x-1]=k_2^{x-2}\times (f[2]-k_1\times f[1])$
那么对于上述两个方程，可以得到：
$(k_1-k_2)f[x]=k_1^{x-2}\times (f[2]-k_2\times f[1])-k_2^{x-2}\times (f[2]-k_1\times f[1])$
即：$f[x]=\frac{k_1^{x-2}\times (f[2]-k_2\times f[1])-k_2^{x-2}\times (f[2]-k_1\times f[1])}{k_1-k_2}$
设：$A=\frac{f[2]-k_2\times f[1]}{k_1-k_2},B=\frac{f[2]-k_1\times f[1]}{k_1-k_2}$
故可得：$f[x]=A \times k_1^{x-2}+B\times k_2^{x-2}$，的通项公式如上
那么根据通项公式，就可以把题目里的原始式子化为$\sum\limits_{S'\subset S ,|S'|=k}A\times k_1^{\prod\limits_{x\in S'} x}+B\times k_2^{\prod\limits_{x\in S'} x}$
那么接下来的就是裸的分治FFT了！
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <queue>
#include <iostream>
#include <bitset>
using namespace std;
#define N 100005
#define mod 99991
#define ll long long
const long double pi=acos(-1);
int a[N],n,k;
struct cp
{
long double x,y;
cp(){}
cp(long double a,long double b){x=a,y=b;}
cp operator + (const cp &a) const {return cp(x+a.x,y+a.y);}
cp operator - (const cp &a) const {return cp(x-a.x,y-a.y);}
cp operator * (const cp &a) const {return cp(x*a.x-y*a.y,x*a.y+y*a.x);}
}A[N<<2],B[N<<2];
void FFT(cp *a,int len,int flag)
{
int i,j,k,t;cp w,x,tmp;long long tt;
for(i=k=0;i<len;i++)
{
if(i>k)swap(a[i],a[k]);
for(j=len>>1;(k^=j)<j;j>>=1);
}
for(k=2;k<=len;k<<=1)
{
x=cp(cos(2*pi*flag/k),sin(2*pi*flag/k));t=k>>1;
for(i=0;i<len;i+=k)
for(w=cp(1,0),j=i;j<i+t;j++)
tmp=a[j+t]*w,a[j+t]=a[j]-tmp,a[j]=a[j]+tmp,w=w*x;
}if(flag==-1)for(i=0;i<len;i++){a[i].x/=len;tt=a[i].x+0.1;tt%=mod;a[i]=cp(tt,0);}
}
struct ploy
{
vector<cp>a;int len;
ploy(){a.clear();len=0;}
ploy(int x){a.resize(3);a[0]=cp(1,0);a[1]=cp(x,0);len=2;}
void print(){printf("%d\n",len);for(int i=0;i<len;i++)printf("%.0lf ",a[i].x);puts("");}
friend ploy operator * (const ploy &a,const ploy &b)
{
// a.print();b.print();
int len=1;while(len<a.len+b.len)len<<=1;
for(int i=0;i<a.len;i++)A[i]=a.a[i];for(int i=a.len;i<len;i++)A[i]=cp(0,0);
for(int i=0;i<b.len;i++)B[i]=b.a[i];for(int i=b.len;i<len;i++)B[i]=cp(0,0);
FFT(A,len,1);FFT(B,len,1);for(int i=0;i<len;i++)A[i]=A[i]*B[i];FFT(A,len,-1);
ploy c;c.len=min(a.len+b.len-1,k+1);c.a.resize(c.len);
for(int i=0;i<c.len;i++)c.a[i]=A[i];
// c.print();
return c;
}
}ret,ret1;
int f1,f2,c1,c2;
int q_pow(int x,int n){int ret=1;for(;n;n>>=1,x=(ll)x*x%mod)if(n&1)ret=(ll)ret*x%mod;return ret;}
int get(int x){return q_pow(3,x);}
ploy solve(int l,int r)
{
if(l==r)return ploy(get(a[l]));int m=(l+r)>>1;
return solve(l,m)*solve(m+1,r);
}
int get1(int x){return q_pow(mod-1,x);}
ploy solve1(int l,int r)
{
if(l==r)return ploy(get1(a[l]));int m=(l+r)>>1;
return solve1(l,m)*solve1(m+1,r);
}
int main()
{
// freopen("see.in","r",stdin);
// freopen("see.out","w",stdout);
scanf("%d%d",&n,&k);
for(int i=1;i<=n;i++)scanf("%d",&a[i]);
scanf("%d%d",&f1,&f2);c1=(f1+f2)*41663ll%mod;
c2=(3ll*c1-f1+mod)%mod;
c1=c1*3ll%mod;c2=((ll)c2*(mod-1))%mod;
// c1=c1*(ll)q_pow(66661ll,k-1)%mod;c2=(ll)c2*q_pow(mod-1,k-1)%mod;
ret=solve(1,n);ret1=solve1(1,n);
printf("%lld\n",(((long long)(ret.a[k].x+0.1)*c1+(long long)(ret1.a[k].x+0.1)*c2)%mod+mod)%mod);
}
