多米诺覆盖

描述 Description
	JSOI2005夏令营的一天，JYY和lhc一起玩多米诺游戏。  lhc提出一个令JYY抓头的问题：有一个m 行n 列的矩形方格棋盘，1<=m,n<=10^9，用1*2 的骨牌(可横放或竖放)完全覆盖，骨牌不能重叠，有多少种不同的覆盖的方法。JYY暗想，这题范围这么大，公式又那么难找，这怎么办呢？  于是他要求lhc做出让步，lhc不得不做出让步，现在1<=n<=10^9，1<=m<=5。于是JYY高兴地在5分钟内解决了这个问题，不幸的是他把这个程序弄丢了，但他从不把做过的东西再做一遍，于是他把任务交给了你。  请你帮助JYY编程解决:)  lhc不想看到长串的高精度数，你只需要求出覆盖方法总数 mod p 的值即可。

输入格式 Input Format
	一行，三个整数数n，m，p，1<=n<=10^9，1<=m<=5，1<=p<=10000。

输出格式 Output Format
	一个整数，总数 mod p 的结果。

样例输入 Sample Input

7 2 10

样例输出 Sample Output

1

来源 Source
	JSOI06夏令营测试题增强版

原来的夏令营测试题的n的范围貌似只有1000，只要for一遍DP就行了。这题n这么大，加个矩阵快速幂就可以了。

我的思想是这样的：把每一行的状态表示成一个二进制数，因为m只有5，所以状态数最多只有64种。每一种状态中，0表示没有被覆盖，1表示被覆盖。考虑当前行在放置骨牌后会转移到下一行的什么样的状态，这样就可以得到一个矩阵。因为要覆盖满，所以把第0行的状态当做都被覆盖，第n+1行的状态当做都没有被覆盖，求出矩阵的n+1次幂，输出最后状态为零的情况数即可。

 

AC CODE

 

program fu_1380;
var t,tt,tmp:Array[0..64,0..64] of longint;
    d,a,b:array[0..64] of longint;
    i,j,n,m,p,k:longint;
//============================================================================
procedure ins;
var k,j:longint;
begin
  k:=0; for j:=1 to m do k:=k+b[j]*d[j];
  inc(T[i,k]); t[i,k]:=t[i,k] mod p;
end;
//============================================================================
procedure dfs(now,step:longint);
begin
  if now=m+1 then ins else
  if a[now]=1 then dfs(now+1,1) else    //当前格被覆盖，则下一行的对应格通过转移就一定不会被覆盖。
  begin
    b[now]:=1; dfs(now+1,1); b[now]:=0;    //竖放，占两行。
    if a[now+1]=0 then dfs(now+2,2);    //有连续两个空格就须考虑平放在一行的情况。则下一行这两格都不会被覆盖。
  end;
end;
//============================================================================
procedure quick_power(x:longint);    //矩阵快速幂。
var i,j,k:longint;
begin
  if x=1 then
  begin
    for i:=0 to d[m+1]-1 do
      for j:=0 to d[m+1]-1 do tt[i,j]:=t[i,j];
    exit;
  end; quick_power(x div 2);
  for i:=0 to d[m+1]-1 do
    for j:=0 to d[m+1]-1 do
      for k:=0 to d[m+1]-1 do
        tmp[i,j]:=(tmp[i,j]+tt[i,k]*tt[k,j]) mod p;
  for i:=0 to d[m+1]-1 do
    for j:=0 to d[m+1]-1 do
    begin
      tt[i,j]:=tmp[i,j]; tmp[i,j]:=0;
    end;
  if x mod 2=1 then
  begin
    for i:=0 to d[m+1]-1 do
      for j:=0 to d[m+1]-1 do
        for k:=0 to d[m+1]-1 do
          tmp[i,j]:=(tmp[i,j]+tt[i,k]*t[k,j]) mod p;
    for i:=0 to d[m+1]-1 do
      for j:=0 to d[m+1]-1 do
      begin
        tt[i,j]:=tmp[i,j]; tmp[i,j]:=0;
      end;
  end;
end;
//============================================================================
begin
  readln(n,m,p); d[1]:=1;
  for i:=2 to m+1 do d[i]:=(d[i-1]*2) mod p;  //预处理2的k次幂。
  a[1]:=-1; a[m+1]:=-maxlongint;
  for i:=0 to d[m+1]-1 do
  begin
    inc(a[1]); k:=1;

    while a[k]>1 do
      begin a[k]:=0; inc(a[k+1]); inc(k); end;      //枚举当前行的状态。
    for k:=1 to m do b[i]:=0;
    dfs(1,1);    //深搜枚举可以转移到的状态。
  end; quick_power(n+1);    //矩阵快速幂。
  writeln(tt[0,d[m+1]-1]);    //本来是用第一行的状态乘以转移矩阵的n+1次幂的。得到最后一行的状态，再输出状态为0的情况数，但是观察发现就是tt[0,d[m+1]-1]...
end.