题解 CF506E 【Mr. Kitayuta's Gift】

shadowice1984 · 2018-10-11 09:00:33 · 题解

神仙矩阵乘法……

我矩乘今天就是要玩死你

其实这种回文串计数的题还是相当套路的

先简述一下题意

给定一个字符串S,输入一个m问你有多少个长度为|S|+m的回文串含有S这个子序列

前置芝士:矩阵快速幂

蛤？淦cfE题不会矩乘？出门左转你站膜板区,包教包会

本题题解

为了减少不必要的分情况讨论,我们以下的解法仅仅针对|S|+m是偶数的情况

首先这是一道可爱(并不)的数数题，那么既然是计数我们常用的一个套路就是dp了

这道题显然最要命的限制就是回文串了,由于我们现在仅仅考虑偶回文串的情况,所以情况还比较简单。

我们仿照一般的对字符串进行dp的套路,我们不停的一位一位的放字符来进行决策,换句话讲我们的转移就是枚举下一位放什么,这样dp有一个相当棒的好处就是我们永远不会重复数同一个串

由于我们的限制是回文串因此我们的决策是同时在字符串的头和尾放一个字符进行转移,(我们这样做一定可以dp出所有的偶回文串)

那么我们就会有一个比较naive的dp设计,我们设dp(i,j,k)表示现在字符串的前k位和后k位已经被决策完毕,当前字符串的(i,j)这段区间还没有被匹配完毕时的方案数,另设一个状态target(k)表示该字符串的前k位和后k位同时被确定,并且已经包含S这个子序列的字符串数目

那么我们状态中的"还没有匹配的区间"是什么意思呢？

我们发现一个回文串会可能会有很多个等于字符串S的子序列,所以我们的dp状态的限制维度设计的不好的话我们就会很容易将同一个串多次计数

所以我们采取这样一种策略匹配子序列从头尾两端同时添加字符,一旦可以和|S|匹配子序列就进行匹配,可以证明的是如果这个字符串含有子序列S 就肯定能被这个算法识别出一个子序列来，并且我们还可以证明的是对于同一个字符串被识别出的子序列有且仅有一种,因此我们会发现这样设计状态的话每个字符串就只会被统计一次了,所以我们的"还没有匹配的区间"指的就是这样贪心的匹配完前后k位之后

那么我们的dp转移就十分的简易了

假设我们已知dp(i,j,k)的值我们现在想要更新其他状态的值

那么我们根据S(i)==S(j)这个命题的真假来进行分类讨论进行转移

Case1:S(i)==S(j)\&\&j-i\leq1

此时我们发现这种情况下我们只需要填对S(i)这个字符我们就可以成功的完成对于S的匹配了,另外的25种填字符方案都会仅仅使步数+1

所以我们的转移方程就是

Case2:S(i)==S(j)\&\&j-i > 1

这种情况下我们发现必须填对S(i),S(j)所代表的字符才能进行转移,其他情况下仅仅只是令步数+1

Case3:S(i)!=S(j)

这种情况下我们会发现只要填对S(i)或者S(j)所代表的字符才能进行转移,其余的24种情况才会令步数+1

Case4:target

显然target状态的话可以随便转移，因此26种字符都可以填

那么对于一个偶回文串来讲我们需要做的就是提取target((|S|+m)/2)的值就可以得到总的方案数了，上面的dp显然是可以大力矩乘的,不过由于状态的总量是O(|S|^2)的这样写我们会得到一个O(|S|^{6}log(|S|+m))的严格过不去算法

那么我们想想有什么的办法可以优化这个暴力呢?我们真的利用了全部的信息吗？有没有信息可以压缩呢？

答案是显然的,上述的算法存在着严重的信息浪费我们明明刷出了一个信息量是O(|S|^{4})级别的矩阵但是我们却仅仅提取的其中一项(也就是target)的值,很明显中间的变量全部被浪费了

所以我们要想办法利用起这些中间变量来

怎么办呢？

我们观察一下我们刚才的dp转移方程,如果状态A和状态B之间存在着这样一种转移方程的话

我们就在A->B之间连C条A指向B的有向边

(注意如果两个状态仅仅是k值不同的话我们视为相同的状态,比如dp(i,j,k)和dp(i,j,k+1)视为相同状态)

那么我们发现这样建出来的图如果去掉自环的话就是一张DAG了,而我们target((|S|+m)/2)的值就是从(1,n)这个点出发,走(|S|+m)/2条边之后到达target点的方案数

很好所以我们就这样把一个dp问题转化为了图上的路径计数问题

那这个问题怎么做呢？

我们首先值得观察的一点就是图上点的自环数目总共就只有三种:24,25,26

那么这又有什么用呢？

我们观察这个图上的每一条路径,这条路径经过的边可能很多，不过经过的点的种类却很少,最多经过|S|个点(因为走过的边绝大部分都是在走自环)

那么如果我们将经过点的点集相同的路径成为"一类路径"的话我们会发现两类路径的所包含的路径条数相等当且仅当这两条路径上自环数目为24,25,26的点的个数分别相等(经过点集相同之后两个路径不同的唯一可能就是走自环方式不同了)

其中26相等的限制明显是个废话,因为整张图就一个点有26的自环它就是终点，我们统计的就是从(1,n)到终点的路径，所以终点显然会出现在每一条路径中

另一个相当有趣的性质就是如果我们已知一条路径上有a个自环数目为24的点,那么自环数目为25的点就是\lceil \frac{n-a}{2} \rceil，这个性质可以由dp的转移特性推出，因为我们每次走过一个自环为24的点仅仅消耗一个字符，不过走过自环为25却经常消耗两个字符,唯一的特例出现在i==j的时候,此时我们的自环数为25不过却消耗了一个字符,这就是我们的公式里面加上了上取整的原因

那这个性质告诉我们什么呢？图中本质不同的路径类仅仅只有O(|S|)种,换句话说两个路径类包含的路径条数不同当且仅当他们包含的自环数为24的点不同

那么我们可以设一个cnt(i,j)表示从j出发到达终点走了i个自环数为24的点的方案数，这个状态可以通过记忆化搜索的方式以O(|S|^3)的方式求出

有了这个状态我们就可以枚举走了几个24自环点,用矩阵快速幂求出这类路径的总方案数，然后读一下cnt数组找到这类路径一共有几条二者乘一下就可以了

这样我们就得到了一个O(|S|^4log(|S|+m))的过不去算法

怎么办呢？

我们发现刚才的算法之所以慢是因为你做了O(|S|)次矩阵乘法,问题来了你为什么要做矩乘呢？肯定是为了算从起点到终点的路径条数对不对啊,那么问题来了你刷出了一个O(|S|^2)的矩阵你却就提取这个矩阵的一位，剩下的中间变量又双叒叕被你扔了……算法当然会慢了

那么我们接下来该怎么办呢？

我们要抛弃传统矩乘的思路,传统矩阵乘法在统计路径条数的时候显然是只有一个起点和一个汇点的,但是我们考虑一下矩阵乘法的意义立即就会发现一个事实，假如我们最后刷出来的矩阵是C的话,C(i,j)表示的就是从i走到j的方案数

换句话说我们的矩阵乘法其实兹瓷多个源点和汇点的,只是我们一般不用而已

那么对于这道题来讲多源点和汇点的矩乘就显得十分有用了

我们这样建一张图,前|S|-1个点连成一条有向的链并且每个点有一个24的自环指向自己,之后(|S|+1)/2个点连成一条有向链并且和前面的链接起来,这些点每个点有一个25的自环,最后每个有25自环的点连出去一条边和一个有26个自环的点相连(注意我们需要开(|S|+1)/2)个有26自环的点)

你这样建一个图你就会发现首先图上任意两点的路径都是唯一的(不考虑自环),第二就是任意一个走了n个24自环点m个25自环点1个26自环的路径都可以在图上被表示出来,换句话说,这张图蕴含了所有的路径类

那么我们就可以对这个图做一遍矩乘就可以计算出所有路径类的路径条数了(这对应了图上某两个点之间的方案数)

至此我们得到了一个O(|S|^{3}log(|S|+m)的做法我们已经可以解决偶回文串的问题了

问题来了似乎还有奇回文串的问题啊？

笑容逐渐凝固

关于奇回文串的解决方案

奇回文串和偶回文串的唯一不同在于奇回文串有一个回文中心它不需要和其他的字符一样，而偶回文串没有

那其实我们的做法异常简单粗暴，我们枚举回文中心是什么字符

首先我们对于所有的字符来讲target((|S|+m)/2)一定是一个合法的状态

然后假设我们枚举的回文中心的字符是char的话那么dp(i,i,(|S|+m)/2)(其中S(i)==char)也会是一个合法的状态

而target((|S|+m)/2)的值当然可以直接按照偶回文串的解决方案求出

现在就是我们多了一些dp值需要求出

怎么办呢？

我们把dp的转移图中的边全部反向,计算一个cnt(i,j)数组表示从j点出发到(1,n)经过了i个24自环点的方案数,然后我们枚举考虑每个dp(i,i)状态到(1,n)的24自环点数目然后读取一下矩阵上的信息就可以计算了

关于卡常数

嗯~ o(*￣▽￣*)然后我们信心满满的交了这个算法

“Time limit exceeded on test 14”

就算你有6000ms时限,就算你是cf评测姬,我说你T,你还是得T

卡常数了不起啊？

sorry，卡常数真的是可以为所欲为的

以后天天卡他常数！

“Time limit exceeded,Time limit exceeded”

然后你发现你的矩阵似乎是一个上三角矩阵

这意味着我们一般的矩阵乘法

for(int i=1;i<=n;i++)
    for(int k=1;k<=n;k++)
        for(int j=1;j<=n;j++)(a[i][j]+=b[i][k]*c[k][j])%=mod;

可以换成

for(int i=1;i<=n;i++)
    for(int k=i;k<=n;k++)
        for(int j=k;j<=n;j++)(a[i][j]+=b[i][k]*c[k][j])%=mod;

常数减少至原来的1/6就可以通过这题了

上代码~

#include<cstdio>
#include<algorithm>
using namespace std;const int N=2*1e2+10;const int mod=10007;
int tr[N][N];int tp[N*N];int n;char mde[N];int m;int ctt;int trp[N<<1][N<<1];int S;int ans;
struct Grph
{
    int v[2*N*N];int x[2*N*N];int ct;int al[N*N];int dp[N][N*N];bool book[N][N*N];
    inline void add(int u,int V){v[++ct]=V;x[ct]=al[u];al[u]=ct;}
    inline int dfs(int cnt,int u)
    {
        if(cnt<0)return 0;if(book[cnt][u])return dp[cnt][u];book[cnt][u]=true;
        for(int i=al[u];i;i=x[i])(dp[cnt][u]+=dfs(cnt-tp[u],v[i]))%=mod;
        return dp[cnt][u];
    }
}g2;//转移的dag 
struct mar//矩阵类 
{
    int mp[N<<1][N<<1];
    inline int* operator [](const int x){return mp[x];}
    void operator *=(const mar& b)
    {
        for(int i=1;i<=S;i++)for(int j=1;j<=S;j++)trp[i][j]=0;
        for(int i=1;i<=S;i++)
            for(int k=i;k<=S;k++)
                for(int j=k;j<=S;j++)(trp[i][j]+=mp[i][k]*b.mp[k][j])%=mod;
        for(int i=1;i<=S;i++)for(int j=1;j<=S;j++)mp[i][j]=trp[i][j];
    }
}trs,res;int f[N];int g[N];
int main()
{
    scanf("%s",mde+1);scanf("%d",&m);for(n=1;mde[n+1]!='\0';n++);
    for(int i=1;i<=n;i++)for(int j=i;j<=n;j++)tr[i][j]=++ctt;++ctt;
    for(int i=1;i<=n;i++)
        for(int j=i;j<=n;j++)
            if(mde[i]==mde[j])//建图 
            {if(j-i<=1)g2.add(ctt,tr[i][j]);else g2.add(tr[i+1][j-1],tr[i][j]);}
            else {tp[tr[i][j]]=1;g2.add(tr[i+1][j],tr[i][j]),g2.add(tr[i][j-1],tr[i][j]);}
    for(int i=0;i<n;i++)g2.book[i][tr[1][n]]=true;g2.dp[tp[tr[1][n]]][tr[1][n]]=1;
    for(int i=0;i<n;i++)g2.dfs(i,ctt);//记忆化搜索 
    for(int i=1;i<n;i++)trs[i][i]=24;int t=(n+1)/2;//构建矩阵 
    for(int i=0;i<t;i++)trs[n+i][n+i]=25;
    for(int i=0;i<t;i++)trs[n+t+i][n+t+i]=26;
    for(int i=1;i<n+t-1;i++)trs[i][i+1]++;
    for(int i=0;i<t;i++)trs[n+i][n+t+i]++;S=n+(t<<1)-1;
    for(int i=1;i<=S;i++)res[i][i]=1;
    for(int p=(n+m)>>1;p;p>>=1,trs*=trs)if(p&1)res*=trs;
    for(int i=0,v;i<n;i++)v=(n-i+1)/2,f[i]=res[n-i][n+t+v-1],g[i]=res[n-i][n+v-1];
    if((n+m)&1)
    {
        for(int i=0;i<n;i++)(ans+=f[i]*g2.dp[i][ctt])%=mod;(ans*=26)%=mod;
        for(int i=1;i<=26;i++)
            for(int j=1;j<=n;j++)
                if(mde[j]=='a'+i-1)for(int k=0;k<n;k++)(ans+=g[k]*g2.dp[k][tr[j][j]])%=mod;
    }
    else{for(int i=0;i<n;i++)(ans+=f[i]*g2.dp[i][ctt])%=mod;}printf("%d",ans);return 0;//拜拜程序~ 
}