LOJ #6183 看无可看

最新推荐文章于 2023-11-21 12:28:31 发布

weixin_33737774

最新推荐文章于 2023-11-21 12:28:31 发布

阅读量101

点赞数

原文链接：http://www.cnblogs.com/Winniechen/p/10623375.html

版权

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);
}

转载于:https://www.cnblogs.com/Winniechen/p/10623375.html