【BZOJ4870】组合数问题（计数DP，快速幂）

题意：

1 ≤ n ≤ 10^9, 0 ≤ r < k ≤ 50, 2 ≤ p ≤ 2^30 − 1

思路：From http://blog.csdn.net/qq_33229466/article/details/70665582

实际上就是要我们从nk件物品里面选出若干件，使得其数量模k等于r的方案数。 
显然的dp方程f[i,j]表示前i件物品拿了若干件使得其数量模k等于j的方案数。 
那么显然有f[i,j]=f[i−1,j]+f[i−1,j−1] 
矩阵乘法优化即可。 
复杂度O(k3logn)

还有一种更棒的做法，同样是dp，但可以发现f[n∗2,i+j]+=f[n,i]∗f[n,j] 
可以理解成枚举前n个物品的选法和后n个物品的选法。 
那么直接对dp数组做快速幂即可。 
复杂度O(k2logn)

type arr=array[0..60]of int64;
var ans,c,a:arr;
    n,p,k,r,t:int64;

procedure dp(var a:arr;b:arr);
var i,j:longint;
begin
 for i:=0 to k-1 do c[i]:=0;
 for i:=0 to k-1 do
  for j:=0 to k-1 do
   c[(i+j) mod k]:=(c[(i+j) mod k]+a[i]*b[j] mod p) mod p;
 for i:=0 to k-1 do a[i]:=c[i];
end;

begin
 assign(input,'bzoj4870.in'); reset(input);
 assign(output,'bzoj4870.out'); rewrite(output);
 read(n,p,k,r);
 inc(ans[0]); inc(ans[1 mod k]);
 a:=ans; t:=n*k-1;
 while t>0 do
 begin
  if t and 1=1 then dp(ans,a);
  dp(a,a);
  t:=t>>1;
 end;
 writeln(ans[r]);

 close(input);
 close(output);
end.