BZOJ 3610 Heoi2014 林中路径 矩阵乘法

题目大意：给定一张有向图，多次询问从S到T经过不超过K条边的所有路径的长度的平方和

首先这题一点也不麻烦

现有一带权整数集合S

我们令一个矩阵F_S表示从第i个点到第j个点，经过k条边(k∈S)的所有路径长度的平方和

令矩阵G_S表示从第i个点到第j个点，经过k条边(k∈S)的所有路径长度之和

令矩阵H_S表示从第i个点到第j个点，经过k条边(k∈S)的路径条数

定义T_S为(F_S,G_S,H_S)组成的三元组

定义T的加法为T_S1+T_S2=T_(S1+S2)

定义T的乘法为T_S1*T_S2=T_(S1*S2) 其中后面的'*'表示两个带权集合的卷积

比如T_{3}+T_{5}=T_{3,5}

T_{2,3}*T_{5,7}=T_{7,8,9,10}

现在我们来看看T的乘法如何计算

首先讨论F_(S1*S2)

比如我们现在有三个点i j k

现在我要从i经过k到达j

假设从i到达k有三条路径，长度分别为k1,k2,k3 (k1,k2,k3∈S1)

假设从k到达j有两条路径，长度分别为k4,k5 (k4,k5∈S2)

那么则有F_(S1*S2)[i][j]=(k1+k4)^2+(k1+k5)^2+(k2+k4)^2+(k2+k5)^2+(k3+k4)^2+(k3+k5)^2

=2*(k1^2+k2^2+k3^2)+3*(k4^2+k5^2)+2*(k1+k2+k3)*(k4+k5)

=F_S1[i][k]*H_S2[k][j]+H_S1[i][k]*F_S2[k][j]+2*G_S1[i][k]*G_S2[k][j]

故F_(S1*S2)=F_S1*H_S2+H_S1*F_S2+2*G_S1*G_S2

等式后面的乘号是矩阵乘法

同理可得到G_(S1*S2)=G_S1*H_S2+H_S1*G_S2 , H_(S1*S2)=H_S1*H_S2

但是要注意的一点是由于下标的运算是卷积，因此会出现重复的情况

比如我们想要用T_[0,2]和T_[0,2]做乘法得到T_[0,4]，这是不正确的，比如T_{2}会被计算3次(T_{0}*T_{2}+T_{1}*T_{1}+T_{2}*T_{0})

因此我们这么搞

对于所有的t(0<=t<=log2(K+1))，预处理出T_{2^t}和T_[0,2^t)

其中T_[0,2^t)的处理方法是T_[0,2^(t-1))*T_{2^(t-1)}+T_[0,2^(t-1))

然后将K进行二进制拆分，按照数位DP那样搞

比如K=1010010，那么我们将T_[0,1010011)拆分成下列区间:

T_[0,1010011)

=T_[0,1000000)+T_[1000000,1010000)+T_[1010000,1010010)+T_[1010010,1010011)

=T_{0}*T_[0,1000000)+T_{1000000}*T_[0,10000)+T_{1010000}*T_[0,10)+T_{1010010}*T_[0,1)

每个乘号的右侧已经预处理出来了，左侧维护一下即可

时间复杂度O(n^3logk)

其实说了这么多只有最后一段是有用的- - 代码也不是很难写- -

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 110
#define MOD 1000000007
using namespace std;
int n,m,k,q;
//f表示平方和，g表示路径长度和，h表示路径条数
long long f[35][M][M],g[35][M][M],h[35][M][M];
//f[t][i][j]表示从i到j经过[0,2^t)条边的平方和 
long long _f[35][M][M],_g[35][M][M],_h[35][M][M];
//_f[t][i][j]表示从i到j经过恰好2^t条边的平方和 
long long ans_f[M][M],ans_g[M][M],ans_h[M][M];
long long _ans_f[M][M],_ans_g[M][M],_ans_h[M][M];
void Add(   long long f1[M][M],long long g1[M][M],long long h1[M][M],
            long long f2[M][M],long long g2[M][M],long long h2[M][M])
{
    int i,j;
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
        {
            (f1[i][j]+=f2[i][j])%=MOD;
            (g1[i][j]+=g2[i][j])%=MOD;
            (h1[i][j]+=h2[i][j])%=MOD;
        }
}
void Multiplication(long long f1[M][M],long long g1[M][M],long long h1[M][M],
                    long long f2[M][M],long long g2[M][M],long long h2[M][M],
                    long long f3[M][M],long long g3[M][M],long long h3[M][M])
{
    int i,j,k;
    static long long re_f[M][M],re_g[M][M],re_h[M][M];
    memset(re_f,0,sizeof re_f);
    memset(re_g,0,sizeof re_g);
    memset(re_h,0,sizeof re_h);
    for(k=1;k<=n;k++)
        for(i=1;i<=n;i++)
            for(j=1;j<=n;j++)
            {
                (re_f[i][j]+=f1[i][k]*h2[k][j]+h1[i][k]*f2[k][j]+2*g1[i][k]*g2[k][j])%=MOD;
                (re_g[i][j]+=h1[i][k]*g2[k][j]+g1[i][k]*h2[k][j])%=MOD;
                (re_h[i][j]+=h1[i][k]*h2[k][j])%=MOD;
            }
    memcpy(f3,re_f,sizeof re_f);
    memcpy(g3,re_g,sizeof re_g);
    memcpy(h3,re_h,sizeof re_h);
}
void Binary_Decomposition(int k)
{
    static long long temp_f[M][M],temp_g[M][M],temp_h[M][M];
    int T;
    for(T=0;1<<T<=k;T++);
    for(T--;~T;T--)
        if(k&(1<<T) )
        {
            Multiplication(_ans_f,_ans_g,_ans_h,f[T],g[T],h[T],temp_f,temp_g,temp_h);
            Add(ans_f,ans_g,ans_h,temp_f,temp_g,temp_h);
            Multiplication(_ans_f,_ans_g,_ans_h,_f[T],_g[T],_h[T],_ans_f,_ans_g,_ans_h);
        }
}
int main()
{
    int T,i,x,y;
    cin>>n>>m>>k>>q;
    for(i=1;i<=n;i++)
        h[0][i][i]=1;
    for(i=1;i<=m;i++)
    {
        scanf("%d%d",&x,&y);
        _f[0][x][y]++;
        _g[0][x][y]++;
        _h[0][x][y]++;
    }
    for(T=1;1<<T<=k+1;T++)
    {
        Multiplication(_f[T-1],_g[T-1],_h[T-1],_f[T-1],_g[T-1],_h[T-1],_f[T],_g[T],_h[T]);
        Multiplication(_f[T-1],_g[T-1],_h[T-1],f[T-1],g[T-1],h[T-1],f[T],g[T],h[T]);
        Add(f[T],g[T],h[T],f[T-1],g[T-1],h[T-1]);
    }
    for(i=1;i<=n;i++)
        _ans_h[i][i]=1;
    Binary_Decomposition(k+1);
    for(i=1;i<=q;i++)
    {
        scanf("%d%d",&x,&y);
        printf("%d\n",(int)ans_f[x][y]);
    }
    return 0;
}