BZOJ 1009: [HNOI2008]GT考试（DP+KMP失配+矩阵快速幂）

1009: [HNOI2008]GT考试

Time Limit: 1 Sec   Memory Limit: 162 MB
Submit: 4181   Solved: 2552
[ Submit][ Status][ Discuss]

Description

　　阿申准备报名参加GT考试，准考证号为N位数X1X2....Xn(0<=Xi<=9),他不希望准考证号上出现不吉利的数字。
他的不吉利数学A1A2...Am(0<=Ai<=9)有M位，不出现是指X1X2...Xn中没有恰好一段等于A1A2...Am. A1和X1可以为
0

Input

　　第一行输入N,M,K.接下来一行输入M位的数。 N<=10^9,M<=20,K<=1000

Output

　　阿申想知道不出现不吉利数字的号码有多少种，输出模K取余的结果.

Sample Input

4 3 100
111

Sample Output

81

好兴奋的做到这么一道题，还发现了其中个别也不知道是否能定义为错误的数据，先说说错误数据的问题，困然我两天，终于忍不住拿别人代码比较，发现矩阵快速幂出了错，相当自信的觉得自己的模板不会错，所以就开始跑数据输出，发现AC的代码会出现负数，AC代码用的数组是int，而longlong的结果跑大数完全不对，所以可能是数据有误，也可能是我理解的不对，不过我感觉出现负数都能过就有点问题了。

然后说说这道想了两天的题，以为模拟可以，但想想这是BZOJ怎么会简单模拟，然后就开始了漫长的演算。感觉数位DP应该是可以的，但是怎么难跟m个字符串联系起来，偷偷看了别人写的，说是KMP的失配函数，然后就去看了看失配数组，果然有收获，想想也是，不就是当前数位的后缀跟m个数的前缀有多少相同的吗，这个地方可以用KMP搞定了，新的问题又来了，DP怎么搞，一看10^9，应该有公式或者矩阵快速幂，然后就在推公式，用DP[i][j]表示第i位数匹配到第j位的种类数。发现DP[i+1][p]+=Dp[i][j]，P是匹配到的任意值，也就是说这一项跟前面所有DP[i]的项都有关，通俗点讲就是当前状态的可能种类只与上一状态有关，这正是可以利用矩阵快速幂的原因。next数组储存的就是当这一位数为i时，有多少种可能情况。

好吧，还是把公式打出来吧             dp[i][j]=a0*dp[i-1][0]+a1*dp[i-1][1]+a2*dp[i-1][2]+``````+a(m-1)dp[i-1][m-1]

然后就可以矩阵快速幂了，构造一个m*m的矩阵，系数a就是失配数组求得的种类，跑一遍矩阵快速幂，最后把dp[n][0]到dp[m-1]加起来的值就是结果。

代码中longlong定义的数组，还防止出现负数，如果想A到，longlong变int，然后去掉防止负数就可以了。

代码实现：

#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<queue>
#include<cstdio>
#define ll long long
#define mset(a,x) memset(a,x,sizeof(a))

using namespace std;
const double PI=acos(-1);
const int inf=0x3f3f3f3f;
const double esp=1e-6;
const int maxn=100000;
const int mod=1e9+7;
int dir[4][2]={0,1,1,0,0,-1,-1,0};
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return a/gcd(a,b)*b;}
ll inv(ll b){if(b==1)return 1; return (mod-mod/b)*inv(mod%b)%mod;}
ll fpow(ll n,ll k){ll r=1;for(;k;k>>=1){if(k&1)r=r*n%mod;n=n*n%mod;}return r;}
int next[100],n,m,q;
char map[100];
struct matrix{  
    int r,c;
    ll dp[40][40];  
    matrix(int R=0,int C=0){  
        memset(dp,0,sizeof(dp));  
        r=R;c=C;  
    }  
}a,b;

matrix operator*(matrix A,matrix B)
{  
    matrix C(A.r,B.c);
    for(int i=0;i<A.r;i++) 
	{
		for(int j=0;j<B.c;j++)  
            for(int k=0;k<A.c;k++)
			{  
                C.dp[i][j]=(C.dp[i][j]+A.dp[i][k]*B.dp[k][j])%q;
                //防止爆表出现负数 
                C.dp[i][j]%=q;  
			    C.dp[i][j]+=q;
			    C.dp[i][j]%=q;
            }
	} 
        
    return C; 
}

//得到失配函数的next数组，用来构造矩阵 
void getfail()
{
	int i=1,j=0;  
    next[1]=0;  
    while(i<=m)  
    {  
        if (j==0||map[i]==map[j]) 
		next[++i]=++j;  
        else 
		j=next[j];  
    }  
}

//矩阵快速幂 
void mat_qpow()
{
	int i,j,k;
	while(n)
	{
		if(n&1)
		a=a*b;
		
		n>>=1;
		b=b*b;
	}
	
	ll ans=0;
	for(i=0;i<m;i++)
	ans=(ans+a.dp[0][i])%q;
	
	cout<<ans<<endl;
}

int main()
{
	int i,j,k,l;
	cin>>n>>m>>q;
	scanf("%s",map+1);
	getfail();
	for(i=0;i<m;i++)                          //默认匹配到自己的数量为1
	a.dp[i][i]=1;
	
	for(i=1;i<=m;i++)                        //利用失配数组构造矩阵 
	{
		for(j=0;j<=9;j++)
		{
			int temp=i;
			while(temp&&map[temp]-'0'!=j)
			temp=next[temp];
			
			b.dp[i-1][temp]++;
		}
	}
	
	//先进行模除 
	for(i=0;i<m;i++)
	{
		for(j=0;j<m;j++)
		b.dp[i][j]%=q;
	}
	a.c=a.r=b.c=b.r=m;
	mat_qpow();
	return 0;
}
/*
100000000
11
100000000
11111111111
*/