Loj#6183. 看无可看

向右看齐 2021-12-20 10:15 160阅读 0赞

# Loj\#6183. 看无可看 #

[题目描述][Link 1]  
首先用特征根求出通项公式\\(A\_n=p\\cdot 3^n+q\\cdot(-1)^n\\)。通过给定的\\(f\_0,f\_1\\)可以解出\\(p,q\\)。

然后我们要求的就是\\(\\sum\_\{|s'|=k\}\\Pi\_\{x\\in s'\}a\_x\\)。这就是个背包。

考虑它的生成函数就是\\(\\Pi(1+a\_ix)\\)。用分治\\(FFT\\)求解。

代码：

#include<bits/stdc++.h>
    #define ll long long
    #define N 100005
    
    using namespace std;
    inline int Get() {int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}
    
    const ll mod=99991;
    int Mod(int a) {return a<mod?a:a-mod;}
    struct matrix {
        ll a[2][2];
        void Init() {memset(a,0,sizeof(a));}
    }g,f;
    matrix operator *(const matrix &x,const matrix &y) {
        static matrix tem;
        tem.Init();
        tem.a[0][0]=(x.a[0][0]*y.a[0][0]+x.a[0][1]*y.a[1][0])%mod;
        tem.a[0][1]=(x.a[0][0]*y.a[0][1]+x.a[0][1]*y.a[1][1])%mod;
        tem.a[1][0]=(x.a[1][0]*y.a[0][0]+x.a[1][1]*y.a[1][0])%mod;
        tem.a[1][1]=(x.a[1][0]*y.a[0][1]+x.a[1][1]*y.a[1][1])%mod;
        return tem;
    }
    matrix operator +(const matrix &x,const matrix &y) {
        static matrix tem;
        tem.Init();
        tem.a[0][0]=Mod(x.a[0][0]+y.a[0][0]);
        tem.a[0][1]=Mod(x.a[0][1]+y.a[0][1]);
        tem.a[1][0]=Mod(x.a[1][0]+y.a[1][0]);
        tem.a[1][1]=Mod(x.a[1][1]+y.a[1][1]);
        return tem;
    }
    
    matrix ksm(matrix g,int x) {
        static matrix ans;
        ans.Init();
        for(int i=0;i<2;i++) ans.a[i][i]=1;
        for(;x;x>>=1,g=g*g) if(x&1) ans=ans*g;
        return ans;
    }
    
    int n,k;
    int a[N];
    bool vis[N];
    matrix G[N];
    ll ans;
    struct Com {
        double a,r;
        Com() {a=0,r=0;}
        Com(double x,double y) {a=x,r=y;}
    };
    
    Com operator +(const Com &x,const Com &y) {return Com(x.a+y.a,x.r+y.r);}
    Com operator -(const Com &x,const Com &y) {return Com(x.a-y.a,x.r-y.r);}
    Com operator *(const Com &x,const Com &y) {return Com(x.a*y.a-x.r*y.r,x.a*y.r+x.r*y.a);}
    Com operator /(const Com &x,const double &y) {return Com(x.a/y,x.r/y);}
    const double pi=acos(-1);
    Com W[20][N<<2],M[20][N<<2];
    
    void FFT(Com *a,int d,int flag) {
        static int rev[N<<2];
        int n=1<<d;
        for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
        for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int s=1;s<=d;s++) {
            int len=1<<s,mid=len>>1;
            for(int i=0;i<n;i+=len) {
                for(int j=0;j<mid;j++) {
                    Com u=a[i+j],v=(flag==1)?a[i+j+mid]*W[d][j*n/len]:a[i+j+mid]*M[d][j*n/len];
                    a[i+j]=u+v;
                    a[i+j+mid]=u-v;
                }
            }
        }
        if(flag==-1) for(int i=0;i<n;i++) a[i]=a[i]/n;
    }
    
    ll ksm(ll t,ll x) {
        ll ans=1;
        for(;x;x>>=1,t=t*t%mod)
            if(x&1) ans=ans*t%mod;
        return ans;
    }
    
    ll F[N<<2];
    Com L[N<<2],R[N<<2];
    Com A[N<<2];
    Com tem[N<<2];
    ll t;
    
    void solve(int l,int r) {
        if(l==r) {
            F[l]=ksm(t,a[l]);
            return ;
        }
        int mid=l+r>>1;
        solve(l,mid),solve(mid+1,r);
        int d=ceil(log2(r-l+2));
        for(int k=0;k<1<<d;k++)
            L[k]=R[k]=Com(0,0);
        for(int i=l;i<=mid;i++) {
            L[i-l+1]=Com(F[i],0);
        }
        for(int i=mid+1;i<=r;i++) {
            R[i-mid]=Com(F[i],0);
        }
        L[0]=R[0]=Com(1,0);
        int len=1<<d;
        FFT(L,d,1);
        FFT(R,d,1);
        for(int k=0;k<1<<d;k++) {
            tem[k]=(L[k]*R[k]);
        }
        FFT(tem,d,-1);
        for(int k=l;k<=r;k++) {
            F[k]=((ll)(tem[k-l+1].a+0.5))%mod;
        }
    }
    
    int main() {
        
        n=Get(),k=Get();
        g.a[0][0]=0,g.a[0][1]=3;
        g.a[1][0]=1,g.a[1][1]=2;
        int d=ceil(log2(n*2));
        for(int i=1;i<=d;i++)
            for(int j=0;j<1<<i;j++)
                W[i][j]=Com(cos(2*j*pi/(1<<i)),sin(2*j*pi/(1<<i))),
                M[i][j]=Com(cos(-2*j*pi/(1<<i)),sin(-2*j*pi/(1<<i)));
        for(int i=1;i<=n;i++) {
            a[i]=Get();
            G[i]=ksm(g,a[i]);
        }
        f.a[0][0]=Get(),f.a[0][1]=Get();
        ll q=(f.a[0][0]+f.a[0][1])*ksm(4,mod-2)%mod,p=(3*f.a[0][0]-f.a[0][1]+mod)*ksm(4,mod-2)%mod;
    
        ll ans=0;
        t=3;
        solve(1,n);
        ans=q*F[k]%mod;
        t=mod-1;
        solve(1,n);
        (ans+=p*F[k])%=mod;
        cout<<ans;
        return 0;
    }

转载于:https://www.cnblogs.com/hchhch233/p/10458324.html

[Link 1]: https://loj.ac/problem/6183