[BZOJ2004]公交线路 DP+矩阵快速幂

设计状态F[i , S]为最快的车到了第i个站台，所有车的状态为S，S必须满足二进制首位为1，且共有K个1。其实S只用记录以i截至的P个站台即可。
DP方程即为F[i, S]=Sigma{F[i-1, S’]} (S’转移到S合法)。
转移合法定义为：首先S’和S合法，然后S’中去掉一个1，右移一位与S相同，形象地理解成让某的一辆车开到最前方，因为每次都是向前走了一步，所以方案没有重复。这也就可以解释为什么只要记录P位，因为P位之后的不可能一次走到最前面。
然后这个对于一个S，它的S’都是相同的，显然可以矩阵快速幂加速。但是我们要用栈先记录所有合法状态，这样的状态最多有C(9,4)=C(9,5)=126种，保证一下复杂度。
顺便说一下，最初和最终状态goal都是高位k个1，后面(p-k)个0 。
最终复杂度：O(logN*C(p-1,k-1)^3)。
我写的代码并不优美。。。

const md=30031;
type
  jznode=array[0..140,0..140]of longint;
var
  n,k,p,i,j,l,size,top,o,goal:longint;
  st:array[0..140]of longint;
  id:array[0..1100]of longint;
  jz,g,t,r:jznode;
procedure c(var a,b:jznode;sa,sb,sc:longint);
var
  i,j,l:longint;
begin
  for i:=0 to sa do
    for j:=0 to sc do
    begin
      t[i,j]:=0;
      for l:=0 to sb do
        t[i,j]:=(t[i,j]+a[i,l]*b[l,j])mod md;
    end;
  a:=t;
end;

procedure ksm(var a:jznode;b:longint);
begin
  r:=a;
  dec(b);
  while b>0 do
  begin
    if b mod 2=1 then c(r,a,top,top,top);
    c(a,a,top,top,top);
    b:=b div 2;
  end;
  a:=r;
end;
function cal(x:longint):longint;
begin
  cal:=0;
  while x>0 do begin inc(cal); x:=x-(x and(-x)); end;
end;
begin
  readln(n,k,p);
  size:=(1<<p)-1;
  o:=(1<<(p-1));
  goal:=size-((1<<(p-k))-1);
  fillchar(jz,sizeof(jz),0);
  for i:=0 to size do
    id[i]:=-1;
  top:=0;
  for i:=o to size do
    if cal(i)=k then
    begin
      if id[i]=-1 then
      begin
        inc(top);
        st[top]:=i;
        id[i]:=top;
      end;
      for j:=1 to p do
        if ((i>>(j-1))and 1)=1 then
        begin
          l:=i;
          l:=l xor (1<<(j-1));
          l:=l>>1;
          l:=l or (1<<(p-1));
          if cal(l)=k then
          begin
            if id[l]=-1 then
            begin
              inc(top);
              st[top]:=l;
             id[l]:=top;
            end;
            jz[id[i],id[l]]:=1;
          end;
        end;
    end;
  ksm(jz,n-k);
  fillchar(g,sizeof(g),0);
  g[id[goal],1]:=1;
  c(jz,g,top,top,1);
  writeln(jz[id[goal],1]);
  readln;
end.