LOJ #6183 看无可看

LOJ #6183 看无可看

优秀的数学题，Orz samjia2000

这个题显然需要将和转积来处理，这个时候就要用到特征方程的一些知识了！

其实就是这个样子：$f[x]=a imes f[x-1]+b imes f[x-2]$

那么必然可以写作：$f[x]-t imes f[x-1]=k imes (f[x-1]-t imes f[x-2])$

化简：$f[x]=(k+t) imes f[x-1]-k imes t imes f[x-2]$

即：$a=k+t,b=-kt$

消去$K$，得到：$k^2=a imes k+b$

这个东西就是这个二阶递推式的特征方程。

那么显然，你可以得到：$f[x]-k_2 imes f[x-1]=k_1^{x-2} imes (f[2]-k_2 imes f[1]),f[x]-k_1 imes f[x-1]=k_2^{x-2} imes (f[2]-k_1 imes f[1])$

那么对于上述两个方程，可以得到：

$(k_1-k_2)f[x]=k_1^{x-2} imes (f[2]-k_2 imes f[1])-k_2^{x-2} imes (f[2]-k_1 imes f[1])$

即：$f[x]=frac{k_1^{x-2} imes (f[2]-k_2 imes f[1])-k_2^{x-2} imes (f[2]-k_1 imes f[1])}{k_1-k_2}$

设：$A=frac{f[2]-k_2 imes f[1]}{k_1-k_2},B=frac{f[2]-k_1 imes f[1]}{k_1-k_2}$

故可得：$f[x]=A imes k_1^{x-2}+B imes k_2^{x-2}$，的通项公式如上

那么根据通项公式，就可以把题目里的原始式子化为$sumlimits_{S'subset S ,|S'|=k}A imes k_1^{prodlimits_{xin S'} x}+B imes k_2^{prodlimits_{xin 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
",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
",(((long long)(ret.a[k].x+0.1)*c1+(long long)(ret1.a[k].x+0.1)*c2)%mod+mod)%mod);
}