【Luogu5293】[HNOI2019] 白兔之舞

最新推荐文章于 2019-09-03 22:05:16 发布

weixin_30611509

最新推荐文章于 2019-09-03 22:05:16 发布

阅读量83

点赞数

原文链接：http://www.cnblogs.com/NeosKnight/p/10679151.html

版权

题目链接

题目描述

略

Sol

考场上暴力 \(O(L)\) 50分真良心。

简单的推一下式子，对于一个 t 来说，答案就是：
\[\sum_{i=0}^{L} [k|(i-t)] {L\choose i}F(i)\]

就是对于所有 mod k 的结果是 t 的 i 的后面那一坨东西的和。

\(F(i)\) 表示走了 \(i\) 步从纵坐标 x 到纵坐标为 y 的方案数，这个东西显然非常好递推。

所以 \(O(L)\) 的暴力做法就直接递推组合数就行了。

然后是正解。

\([k|(i-t)]\)这玩意可以套用单位根的性质: \([k|n]=\frac{1}{k}\sum_{i=0}^kw_k^{in}\)
正确性结合单位根性质和等比数列求和就可以得到。

所以式子就变成了:
\[\frac{1}{k}\sum_{i=0}^{L} {L\choose i}F(i)\sum_{j=0}^k w_k^{j(i-t)}\]

显而易见的把里面拆了然后交换求和顺序:
\[\frac{1}{k}\sum_{j=0}^k w_k^{-jt}\sum_{i=0}^{L}w_k^{ji} {L\choose i}F(i)\]

后面的东西是个套路，具体见 bzoj3328 PYXFIB

我们把 \(F(i)\) 写乘一个矩阵幂次的形式,因为是可以矩阵乘法递推的。
假设转移矩阵是 \(T\)
\[\frac{1}{k}\sum_{j=0}^k w_k^{-jt}\sum_{i=0}^{L}w_k^{ji} {L\choose i}T^i\]
真正的 \(F(i)\) 就是 \(T^i[x][y]\)

后面就是一个组合数+幂次项，是一个的二项式定理的形式，合并一下就是：

\[\frac{1}{k}\sum_{j=0}^k w_k^{-jt}(Tw_k^j+I)^L\]

\(I\) 是单位矩阵。

后面的东西对于一个 \(j\) 而言显然是固定的，我们设 \(G(w_k^j)=(Tw_k^j+I)^L\)

所以要求的是:
\[\frac{1}{k}\sum_{j=0}^k (w_k^{-t})^j G(w_k^j)\]

这玩意看上去好眼熟啊。不就是个 \(IDFT\) 吗?
我们 \(IDFT\) 就是对于求出点值后的多项式代入单位根的倒数,最后结果乘上 \(\frac{1}{n}\)。
所以就是这么回事了。就是让我们对求完后的 \(G(x)\) 给 \(IDFT\) 一下回去就是所有 \(t\) 的答案了

所以问题变成求解任意长度的 \(DFT\) ，这里我们就没有办法用性能优化那题的方法 \(DFT\) 了

然后就可以使用 \(Bluestein’s\; Algorithm\) 了。(具体可参见我的性能优化的题解)

注意矩阵里面的运算要卡常，先用 long long 存下答案再取模，不然常数大不开 O2 过不了。~~(出这种毒瘤题也就算了怎么还卡常啊QAQ)~~

code:

#include<bits/stdc++.h>
using namespace std;
#define Set(a,b) memset(a,b,sizeof(a))
template<class T>inline void init(T&x){
    x=0;char ch=getchar();bool t=0;
    for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=1;
    for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
    if(t) x=-x;return;
}typedef long long ll;
int mod;const int MAXN=65536;
template<class T>inline void Inc(T&x,int y){x+=y;if(x>=mod) x-=mod;}
template<class T>inline void Dec(T&x,int y){x-=y;if(x < 0 ) x+=mod;}
template<class T>inline int fpow(int x,T k){int ret=1;for(;k;k>>=1,x=(ll)x*x%mod)if(k&1) ret=(ll)ret*x%mod;return ret;}
inline int Sum(int x,int y){x+=y;if(x>=mod) x-=mod;return x;}
inline int Dif(int x,int y){x-=y;if(x < 0 ) x+=mod;return x;}
int n,k,L,x,y,w[3][3];
int F[MAXN],g,wn[MAXN],iwn[MAXN];
inline void Getroot(){
    int phi=mod-1;int x=phi;
    static int pri[50],cnt=0;
    for(int i=2;i*i<=x;++i) {
        if(x%i==0) {pri[++cnt]=i,x/=i;while(x%i==0) x/=i;}
    }
    for(g=2;;++g){bool fl=1;
        for(int i=1;i<=cnt;++i) if(fpow(g,phi/pri[i])==1) {fl=0;break;}
        if(fl)return;
    }
}
namespace Work1{
    struct matrix{
        ll a[3][3];matrix(){Set(a,0);}
        inline ll* operator [](int x){return a[x];}
        inline matrix operator +(matrix b){matrix c;for(int i=0;i<n;++i) for(int j=0;j<n;++j) c[i][j]=Sum(a[i][j],b[i][j]);return c;}
        inline matrix operator *(matrix b){
            matrix c;
            for(int i=0;i<n;++i)
                for(int j=0;j<n;++j){
                    for(int k=0;k<n;++k) c[i][j]+=a[i][k]*b[k][j];
                    c[i][j]%=mod;
                }
            return c;
        }
    }T,I;
    inline matrix matrixfpow(matrix x,int k){matrix ret=I;for(;k;k>>=1,x=x*x)if(k&1)ret=ret*x;return ret;}
    inline void work(){
        for(int i=0;i<n;++i) I[i][i]=1;
        for(int i=0;i<k;++i) {T=matrix();
            for(int u=0;u<n;++u)for(int v=0;v<n;++v) T[u][v]=(ll)w[u][v]*wn[i]%mod;
            T=T+I;T=matrixfpow(T,L);F[i]=T[x][y];
        }
        return;
    }
}
namespace Work2{
    const int MAXM=MAXN<<2;
    typedef double db;int rader[MAXM];
    inline int Init(int n){
        int len=1,up=-1;while(len<=n) len<<=1,++up;
        for(int i=1;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
        return len;
    }
    namespace MTT{
        const db PI=acos(-1);
        struct Complex{
            db x,y;Complex(db _x=0.0,db _y=0.0){x=_x,y=_y;}
            inline Complex operator +(const Complex B){return Complex(x+B.x,y+B.y);}
            inline Complex operator -(const Complex B){return Complex(x-B.x,y-B.y);}
            inline Complex operator *(const Complex B){return Complex(x*B.x-y*B.y,x*B.y+y*B.x);}
        }w[MAXM];
        inline void Calc(int n){for(int i=1;i<n;i<<=1) for(int j=0;j<i;++j) w[n/i*j]=Complex(cos(PI/i*j),sin(PI/i*j));return;}
        inline void FFT(Complex*A,int n,int f){
            for(int i=0;i<n;++i) if(rader[i]>i) swap(A[rader[i]],A[i]);
            for(int i=1;i<n;i<<=1)
                for(int j=0,p=i<<1;j<n;j+=p)
                    for(int k=0;k<i;++k){
                        Complex X=A[j|k],Y=A[j|k|i]* ((~f)? w[n/i*k]:Complex(w[n/i*k].x,-w[n/i*k].y));
                        A[j|k]=X+Y,A[j|k|i]=X-Y;
                    }
            if(!~f) for(int i=0;i<n;++i) A[i].x/=(db)n;return;
        }
        inline void Mul(int*A,int*B,int*C,int len){
            static Complex A1[MAXM],A2[MAXM],B1[MAXM],B2[MAXM];
            int MO=sqrt(mod);
            for(int i=0;i<len;++i) A1[i]=Complex(A[i]/MO,0.0),B1[i]=Complex(A[i]%MO,0.0),A2[i]=Complex(B[i]/MO,0.0),B2[i]=Complex(B[i]%MO,0.0);
            FFT(A1,len,1),FFT(A2,len,1),FFT(B1,len,1),FFT(B2,len,1);
            for(int i=0;i<len;++i) {Complex X;
                X=A1[i]*A2[i],A2[i]=A2[i]*B1[i];
                B1[i]=B1[i]*B2[i];B2[i]=B2[i]*A1[i];
                A1[i]=X,A2[i]=A2[i]+B2[i];
            }int MOD=MO*MO%mod;
            FFT(A1,len,-1),FFT(B1,len,-1),FFT(A2,len,-1);
            for(int i=0;i<len;++i) {
                int X=(ll)(A1[i].x+0.5)%mod,Y=(ll)(B1[i].x+0.5)%mod,Z=(ll)(A2[i].x+0.5)%mod;
                int ans=(ll)MOD*X%mod;Inc(ans,(ll)MO*Z%mod);Inc(ans,Y);
                C[i]=ans;
            }return;
        }
    }using MTT::Calc;
    inline int Co(int x){return (ll)x*(x-1)/2%k;}
    inline void DFT(int*A,int n,int len){
        int m=2*n-1;static int F[MAXM],G[MAXM];
        for(int i=0;i<n;++i) F[i]=(ll)A[i]*wn[Co(i)]%mod;for(int i=n;i<len;++i) F[i]=0;
        for(int i=0;i<m;++i) G[i]=iwn[Co(i)];for(int i=m;i<len;++i) G[i]=0;
        reverse(F,F+n);MTT::Mul(F,G,F,len);
        for(int k=0,i=n-1;i<m;++i,++k) A[k]=(ll)F[i]*wn[Co(k)]%mod;
        for(int i=0,inv=fpow(n,mod-2);i<n;++i) A[i]=(ll)A[i]*inv%mod;
        return;
    }
    void work(){
        int len=Init(3*k-3);Calc(len);DFT(F,k,len);
        for(int i=0;i<k;++i) printf("%d\n",F[i]);
        return;
    }
}
int main()
{
    freopen("dance.in","r",stdin);
    freopen("dance.out","w",stdout);
    init(n),init(k),init(L),init(x),init(y),init(mod);--x,--y;Getroot();
    wn[0]=iwn[0]=1,wn[1]=fpow(g,(mod-1)/k),iwn[1]=fpow(wn[1],mod-2);
    for(int i=2;i<k;++i) wn[i]=(ll)wn[i-1]*wn[1]%mod,iwn[i]=(ll)iwn[i-1]*iwn[1]%mod;
    for(int i=0;i<n;++i) for(int j=0;j<n;++j) init(w[i][j]);
    Work1::work();Work2::work();
    return 0;
}

转载于:https://www.cnblogs.com/NeosKnight/p/10679151.html

weixin_30611509

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
【Luogu5293】[HNOI2019] 白兔之舞

题目链接题目描述略Sol考场上暴力 \(O(L)\) 50分真良心。简单的推一下式子，对于一个 t 来说，答案就是：\[\sum_{i=0}^{L} [k|(i-t)] {L\choose i}F(i)\]就是对于所有 mod k 的结果是 t 的 i 的后面那一坨东西的和。\(F(i)\) 表示走了 \(i\) 步从纵坐标 x 到纵坐标为 y 的方案数，这个东西显然非常好递推。...
复制链接

扫一扫