HS
AFO
HS的博客
P5293 [HNOI2019]白兔之舞

题目链接

洛谷

题目描述

有一张顶点数为 $(L+1)\times n$ 的有向图。这张图的每个顶点由一个二元组$(u,v)$表示$(0\le u\le L,1\le v\le n)$。 这张图不是简单图,对于任意两个顶点$(u1,v1)(u2,v2)$,如果 $u1<u2$,则从$(u1,v1)$到$(u2,v2)$一共有 $w[v1][v2]$条不同的边,如果 $u1\ge u2$ 则没有边。
白兔将在这张图上上演一支舞曲。白兔初始时位于该有向图的顶点$(0,x)$。
白兔将会跳若干步。每一步,白兔会从当前顶点沿任意一条出边跳到下一个顶点。白兔可以在任意时候停止跳舞(也可以没有跳就直接结束)。当到达第一维为 $L$ 的顶点就不得不停止,因为该顶点没有出边。
假设白兔停止时,跳了 $m$ 步,白兔会把这只舞曲给记录下来成为一个序列。序列的第 $i$ 个元素为它第 $i$ 步经过的边。
问题来了:给定正整数 $k$ 和 $y$($1\le y\le n$),对于每个 $t$($0\le t<k$),求有多少种舞曲(假设其长度为 $m$)满足$ m \mod k=t$,且白兔最后停在了坐标第二维为 $y$ 的顶点?
两支舞曲不同定义为它们的长度($m$)不同或者存在某一步它们所走的边不同。
输出的结果对 $p$ 取模。

输入格式

输入文件名为 dance.in。
第一行五个用空格隔开的整数 $n,k,L,x,y,p$。
接下来 $n$ 行,每行有 $n$ 个用空格隔开的整数,第 $i$ 行的第 $j$ 个数表示 $w[i][j]$。

输出格式

输出文件名为 dance.out。
依次输出 $k$ 行,每行一个数表示答案对 $p$ 取模的结果。

输入输出样例

输入 #1

2 2 3 1 1 998244353
2 1
1 0

输出 #1

16
18

输入 #2

3 4 100 1 3 998244353
1 1 1
1 1 1
1 1 1

输出 #2

506551216
528858062
469849094
818871911

说明/提示

【样例解释 1】
t=0:
1.路径长度为0,方案数为1。
2.路径长度为2,一共有六类路径
(0,1)->(1,1)->(2,1) 该路径有w[1][1] * w[1][1]=4条
(0,1)->(1,1)->(3,1) 该路径有w[1][1] * w[1][1]=4条
(0,1)->(2,1)->(3,1) 该路径有w[1][1] * w[1][1]=4条
(0,1)->(1,2)->(2,1) 该路径有w[1][2] * w[2][1]=1条
(0,1)->(1,2)->(3,1) 该路径有w[1][2] * w[2][1]=1条
(0,1)->(2,2)->(3,1) 该路径有w[1][2] * w[2][1]=1条
答案就为1+4+4+4+1+1+1=16
t=1:
1.路径长度为1,一共有3类路径
(0,1)->(1,1) 该路径有w[1][1]=2条
(0,1)->(2,1) 该路径有w[1][1]=2条
(0,1)->(3,1) 该路径有w[1][1]=2条
2.路径长度为3,一共有3类路径
(0,1)->(1,1)->(2,1)->(3,1) 该路径有w[1][1] * w[1][1] * w[1][1]=8条
(0,1)->(1,1)->(2,2)->(3,1) 该路径有w[1][1] * w[1][2] * w[2][1]=2条
(0,1)->(1,2)->(2,1)->(3,1) 该路径有w[1][2] * w[2][1] * w[1][1]=2条
答案就为2+2+2+8+2+2=18
【数据范围】
测试点1,2:$L<=100000$
测试点3: $n=1,w[1][1]=1$,k的最大质因子为2
测试点4: $n=1,k$ 的最大质因子为2
测试点5: $n=1,w[1][1]=1$
测试点6: $n=1$
测试点7,8: $k$的最大质因子为$2$
对于全部数据:
p为一个质数
$10^8<p<2^{30}$
$1\le n\le 3$
$1\le x\le n$
$1\le y\le n$
$0\le w[i][j]<p$
$1\le k\le 65536$
$k$为$p-1$的约数
$1\le L\le 10^8$
【编译命令】
对于c++语言: g++ -o dance dance.cpp –lm
对于c语言: gcc -o dance dance.c –lm
对于pascal语言: fpc dance.pas

分析

首先考虑$DP$,设$f_{i,j}$表示走了$i$步,到了第二维为$j$的点。

那么对于答案$g_i$($i$表示走了$i$步后走到了题目中的目标点),有
$$
g_i=\binom{L}{i}f_{i,y}
$$
考虑$f$的状态转移方程式。
$$
f_{i,j}=\sum_{k}^nf_{i-1,k}w[k][j]
$$
即,枚举上一步的第二维位置乘上路径条数。

由于$n$很小,考虑将这个转移方程写成转移矩阵,那么设这个矩阵为$S$。

那么$F_i=F_0 \times S^{i}$,其中$F_i$表示走了$i$步后的$f$构成的矩阵,$F_0$取决于读入的出发点。

那么写出答案的式子
$$
ans_t=\sum_{i=0}^L[i\bmod k==t]g_{i}
$$

$$
ans_t=\sum_{i=0}^L[(i-t) \bmod k==0]g_{i}
$$

$$
ans_t=\sum_{i=0}^L[k|(i-t)]g_{i}
$$

接着考虑单位根反演,有
$$
ans_t=\frac{1}{k}\sum_{i=0}^Lg_i\sum_{j=0}^{k-1}w_k^{(i-t)j}
$$

$$
ans_t=\frac{1}{k}\sum_{i=0}^Lg_i\sum_{j=0}^{k-1}w_k^{-tj}w_k^{ij}
$$

$$
ans_t=\frac{1}{k}\sum_{i=0}^L\sum_{j=0}^{k-1}g_iw_k^{-tj}w_k^{ij}
$$

$$
ans_t=\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^Lg_iw_k^{-tj}w_k^{ij}
$$

$$
ans_t=\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-tj}\sum_{i=0}^Lg_iw_k^{ij}
$$

因为$k|(p-1)$ ,所以在模 $p$ 意义下 $\omega_k =g^{\frac{p-1} k}$ ,其中 $g$ 为 $p$ 的原根。

带入$g$的方程
$$
ans_t=\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-tj}\sum_{i=0}^L\binom{L}{i}f_{i,y}w_k^{ij}
$$
在带入$f$的矩阵方程(注意,虽然是矩阵,但是在式子中只考虑第$y$项,即$f_{i,y}$)
$$
ans_t=\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-tj}\sum_{i=0}^L\binom{L}{i}(F_0 \times S^{i})w_k^{ij}
$$

$$
ans_t=\frac{1}{k}F_0\sum_{j=0}^{k-1}w_k^{-tj}\sum_{i=0}^L\binom{L}{i}S^{i}w_k^{ij}
$$

$$
ans_t=\frac{1}{k}F_0\sum_{j=0}^{k-1}w_k^{-tj}\sum_{i=0}^L\binom{L}{i}(Sw_k^j)^{i}
$$

放入单位矩阵,次数为$L-i$
$$
ans_t=\frac{1}{k}F_0\sum_{j=0}^{k-1}w_k^{-tj}\sum_{i=0}^L\binom{L}{i}(Sw_k^j)^{i}I^{L-i}
$$
于是凑出了二项式定理的形式
$$
ans_t=\frac{1}{k}F_0\sum_{j=0}^{k-1}w_k^{-tj}(Sw_k^j+I)^L
$$

$$
ans_t=\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-tj}F_0(Sw_k^j+I)^L
$$

$F_0(Sw_k^j+I)^L$显然可以使用矩阵快速求出,上文已经提过,虽然是矩阵,但是在式子中只考虑第$y$项,那么,设$F_0(Sw_k^j+I)^L$的第$y$项为$h_j$。
$$
ans_t=\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-tj}h_j
$$
这很像$FFT$的式子,不过这里的循环卷积长度为$k$。

那么考虑使用$\texttt{Bluestein’s Algorithm}$。

具体的,考虑分拆$w_k$的指数,设为$ij$,对于一般的情况,我们将其分拆为$\frac{(i+j)^2}{2}-\frac{i^2}{2}-\frac{j^2}{2}$,但是在本题中,分拆之后$w_k$会变成$w_{2k}$,题目中并没有保证$2k$是$p-1$的因数,而且也不一定存在模意义下的二次剩余,所以不能这样分拆。

换一个思路,将指数分拆成$\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}$。

那么式子就是
$$
ans_t=\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-\tbinom{t+j}{2}+\tbinom{t}{2}+\tbinom{j}{2}}h_j
$$

$$
ans_t=\frac{1}{k}w_k^{\tbinom{t}{2}}\sum_{j=0}^{k-1}w_k^{-\tbinom{t+j}{2}+\tbinom{j}{2}}h_j
$$

$$
ans_t=\frac{1}{k}w_k^{\tbinom{t}{2}}\sum_{j=0}^{k-1}w_k^{-\tbinom{t+j}{2}}\times w_k^{\tbinom{j}{2}}h_j
$$

记$f_i=w_k^{-\tbinom{i}{2}}$,$g_i=w_k^{\tbinom{i}{2}}h_i$。
$$
ans_t=\frac{1}{k}w_k^{\tbinom{t}{2}}\sum_{j=0}^{k-1}f_{t+j}g_j
$$
翻转一下$g$。
$$
ans_t=\frac{1}{k}w_k^{\tbinom{t}{2}}\sum_{j=0}^{k-1}f_{t+j}g^R_{k-1-j}
$$

$$
ans_t=\frac{1}{k}w_k^{\tbinom{t}{2}}\sum_{j=0}^{k-1}f_{t+j}g^R_{k-1-j}
$$

$$
ans_t=\frac{1}{k}w_k^{\tbinom{t}{2}}\sum_{j=0}^{k-1}(f \otimes g^R)_{t+k-1}
$$

所以只要求出$f$和$g$再使用$MTT$卷积一下就行了。

复杂度 $O(kn^3\log L+k\log k)$ 。

一些对我而言的坑

不要打错变量…

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=(1<<18)+10;
const double Pi=acos(-1);
template <typename T>inline T read()
{
    register T sum=0;
    register char cc=getchar();
    int sym=1;
    while(cc!='-'&&(cc>'9'||cc<'0'))cc=getchar();
    if(cc=='-')sym=-1,cc=getchar();
    sum=sum*10+cc-'0';
    cc=getchar();
    while(cc>='0'&&cc<='9')sum=sum*10+cc-'0',cc=getchar();
    return sym*sum;
}
template <typename T>inline T read(T &a)
{
    a=read<T>();
    return a;
}
template <typename T,typename... Others> inline void read(T& a, Others&... b)
{
    a=read(a);
    read(b...);
}
int n,K,L,X,Y,p,inv,omega[maxn],rev[maxn],f[maxn],g[maxn],h[maxn];
int add(int x,int y)
{
    x+=y;
    return x>=p? x-p:x;
}
void Add(int &x,int y)
{
    x+=y;
    if(x>=p)
        x-=p;
}
int mul(int x,int y)
{
    return 1ll*x*y%p;
}
int fpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)
            res=mul(res,x);
        x=mul(x,x);
        y>>=1;
    }
    return res;
}
struct Matrix
{
    int v[3][3];
    Matrix()
    {
        memset(v,0,sizeof(v));
    }
    int *operator [] (int x)
    {
        return v[x];
    }
    Matrix operator + (const Matrix &o)
    const{
        Matrix res;
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
                res[i][j]=add(v[i][j],o.v[i][j]);
        return res;
    }
    Matrix operator * (const Matrix &o)
    const{
        Matrix res;
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
            {
                res[i][j]=0;
                for(int k=0;k<3;k++)
                    Add(res[i][j],mul(v[i][k],o.v[k][j]));
            }
        return res;
    }
    Matrix operator * (const int &o)
    const{
        Matrix res;
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
                res[i][j]=mul(v[i][j],o);
        return res;
    }
}I,A;
Matrix fpow(Matrix x,int y)
{
    Matrix res=I;
    while(y)
    {
        if(y&1)
            res=res*x;
        x=x*x;
        y>>=1;
    }
    return res;
}
struct cp
{
    double r,i;
    cp(double r=0,double i=0):r(r),i(i) {}
    cp operator + (const cp &o)
    const{
        return {r+o.r,i+o.i};
    };
    cp operator - (const cp &o)
    const{
        return {r-o.r,i-o.i};
    };
    cp operator * (const cp &o)
    {
        return {r*o.r-i*o.i,r*o.i+i*o.r};
    };
    cp conj()
    {
        return {r,-i};
    }
}omg[maxn];
void DFT(cp *A,int len)
{
    for(int i=0;i<len;i++)
        if(i<rev[i])
            swap(A[i],A[rev[i]]);
    for(int i=1;i<len;i<<=1)
        for(int j=0,s=i<<1;j<len;j+=s)
            for(int k=0;k<i;k++)
            {
                cp x=A[j+k],y=A[j+i+k]*omg[len/i*k];
                A[j+k]=x+y;
                A[j+i+k]=x-y;
            }
}
void MTT(int *a,int *b,int *c,int len)
{
    static cp s1[maxn],s2[maxn],s3[maxn],s4[maxn],s5[maxn],s6[maxn];
    int n=1,bit=0;
    while(n<len)n<<=1,bit+=1;
    for(int i=1;i<n;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    for(int i=0;i<n;i++)
        omg[i]=cp(cos(Pi*i/n),sin(Pi*i/n));
    for(int i=0;i<n;i++)
        s1[i]=cp(a[i]&32767,a[i]>>15);
    for(int i=0;i<n;i++)
        s2[i]=cp(b[i]&32767,b[i]>>15);
    DFT(s1,n);
    DFT(s2,n);
    for(int i=0;i<n;i++)
    {
        int j=(n-i)&(n-1);
        cp da=(s1[i]+s1[j].conj())*cp(0.5,0);
        cp db=(s1[i]-s1[j].conj())*cp(0,-0.5);
        cp dc=(s2[i]+s2[j].conj())*cp(0.5,0);
        cp dd=(s2[i]-s2[j].conj())*cp(0,-0.5);
        s3[j]=da*dc;
        s4[j]=da*dd;
        s5[j]=db*dc;
        s6[j]=db*dd;
    }
    for(int i=0;i<n;i++)
    {
        s1[i]=s3[i]+s4[i]*cp(0,1);
        s2[i]=s5[i]+s6[i]*cp(0,1);
    }
    DFT(s1,n);
    DFT(s2,n);
    for(int i=0;i<n;i++)
    {
        long long da=(long long)(s1[i].r/n+0.5)%p;
        long long db=(long long)(s1[i].i/n+0.5)%p;
        long long dc=(long long)(s2[i].r/n+0.5)%p;
        long long dd=(long long)(s2[i].i/n+0.5)%p;
        long long val=da+((db+dc)<<15)+(dd<<30);
        c[i]=val%p;
        if(c[i]<0)
            c[i]+=p;
    }
}
int getrt()
{
    vector<int>v;
    int x=p-1;
    for(int i=2;i*i<=x;i++)
    {
        if(x%i)
            continue;
        v.push_back(i);
        while(x%i==0)
            x/=i;
    }
    for(int g=2;;g++)
    {
        int flg=true;
        for(auto i:v)
            if(fpow(g,(p-1)/i)==1)
            {
                flg=false;
                break;
            }
        if(flg)
            return g;
    }
}
int C(int x)
{
    return (1ll*x*(x-1)/2)%K;
}
int main()
{
    read(n,K,L,X,Y,p);
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            read(A[i][j]);
    for(int i=0;i<n;i++)
        I[i][i]=1;
    omega[0]=1,omega[1]=fpow(getrt(),(p-1)/K);
    for(int i=2,j=omega[1];i<=K;i++)
        omega[i]=mul(omega[i-1],j);
    for(int i=0;i<K;i++)
        g[K-1-i]=mul(fpow(A*omega[i]+I,L)[X-1][Y-1],omega[C(i)]);
    for(int i=0;i<2*K;i++)
        f[i]=omega[K-C(i)];
    MTT(f,g,h,3*K);
    inv=fpow(K,p-2);
    for(int i=0;i<K;i++)
        printf("%d\n",mul(inv,mul(omega[C(i)],h[i+K-1])));
    return 0;
}

Gluttony

文章作者

发表评论

textsms
account_circle
email

HS的博客

P5293 [HNOI2019]白兔之舞
单位根反演+DP+矩阵快速幂+MTT+Bluestein's Algorithm
扫描二维码继续阅读
2020-04-19