jzoj 3580. 【NOI2014模拟】矩阵染色

Description

Input

输入共一行，包含三个正整数n、m、p。

Output

输出共一行，即所求答案。

Sample Input

2 5 100000007

Sample Output

266880

【样例解释】

当m=5 时，f(1)=20，f(2)=260。所以答案为20*2^5+260*4^5=266880。

Data Constraint

这题是真的烦……

分析：
｛x,x｝
｛y,z｝
对于第一列，方案数为m*(m-1){x不等于y}，第二列首先保证第二列的两个不同，方案为m*(m-1){同理}。

显然第一行相同方案有m-1种，显然不能x=z，也就是说当z的位置选了x，那么上面的x就不能选,因为已经被排除，第二行相同同理，同时加上两个都相同的一种，容斥一下，即
｛x,x｝
｛y,y｝

方案为

m*(m-1)*(m*(m-1)-2*(m-1)+1)=m*(m-1)*(m^2-3m+3)

使用第二列推第三列，同理，乘上一个m*(m-1)*(m^2-3m+3)
有

f[i]=m*(m-1)*(m^2-3m+3)^i-1 

答案为

sum{m*(m-1)*(m^2-3m+3)^(i-1)*2^m*i^m} {1<=i<=n}

暴力求解可以得到20%分数。

假设t=m^2-3m+3,
把m*(m-1)*2^m提出来,在给里面乘上t,外面乘1/t,得

[m*(m-1)*2^m/t]*[sum(t^i*i^m)] {1<=i<=n}

前面的为常数，设为a，a=m*(m-1)*2^m/t

关键求sum(t^i*i^m) {1<=i<=n}

---

尝试使用递推……
我们设s[n,m]=sum(t^i*i^m) {1<=i<=n}。

s[n+1,m]=sum(t^i*i^m) ｛1<=i<=n+1｝

把当i=1的状态提出来,然后把i范围从[1,n+1]缩小到[2,n+1],再把区间前移为[1,n],有

s[n+1,m]=t+sum(t^(i+1)*(i+1)^m) {1<=i<=n}

我们知道(i+1)^m拆开后，各项的系数是杨辉三角的某一行，应该都知道,把t^(i+1)提一个t出来，有

通过另一种方式，我们可以得到

观察最后一个等式，我们设x=S[n,m],对于右边,当j=m时,也为x,显然是关于x的一元一次方程，可用O(m^2)解决。

复杂度：O(m^2)
说了这么多，只能解决另外的60%。

我们发现还有20%，n和m特别大，不过p很小。根据欧拉定理,可以知道 t^i*i^m余数会循环，循环节长度不超过p^2,可以过。

｛这道题其实是3个程序结合｝

代码：

var
 n,m,p,t,a,x,pow,sum,y,g:int64;
 i,j,k:longint;
 s:array [0..1001] of int64;
 re,num:array [-1..5000001] of int64;
 c:array [-1..1001,-1..1001] of int64;

function ksm(x,y:int64):int64;
var e:int64;
 begin
  if y=1 then exit(x);
  e:=ksm(x,y shr 1);
  e:=e*e mod p;
  if y and 1=1 then e:=e*x mod p;
  exit(e);
 end;

begin
 readln(n,m,p);
 t:=(m*m-3*m+3) mod p;
 a:=ksm(2,m)*(m*(m-1) mod p) mod p;
 y:=ksm(t,p-2);
 a:=(a*y) mod p;
 if n<=100000 then
  begin
   pow:=1;
   for i:=1 to n do
    begin
     pow:=(pow*t) mod p;
     sum:=(sum+pow*ksm(i,m)) mod p;
    end;
   writeln(a*sum mod p);
   exit;
  end;
 if m<=1000 then
  begin
   x:=ksm(t,n+1);
   g:=ksm(t-1,p-2);
   if t<>1 then
     s[0]:=((x+p-t) mod p*g) mod p
   else s[0]:=n;
   pow:=1;
   c[0,0]:=1;
   for i:=1 to m do
    begin
     pow:=(pow*(n+1)) mod p;
     sum:=0;
     for j:=0 to i-1 do
      begin
       c[i,j]:=(c[i-1,j-1]+c[i-1,j]) mod p;
       sum:=(sum+c[i,j]*s[j]) mod p;
      end;
     c[i,i]:=1;
     s[i]:=(((x*pow) mod p+p-(t+t*sum) mod p) mod p*g) mod p;
    end;
   writeln((a*s[m]) mod p);
   exit;
  end;
  pow:=1; j:=0;
   for i:=1 to n do
    begin
     pow:=(pow*t) mod p;
     y:=pow*ksm(i,m) mod p;
     re[i]:=y mod p;
     num[i]:=(num[i-1]+re[i]) mod p;
    if (re[i]=re[j+1]) and (i<>j+1) then inc(j)
     else
      begin
       j:=0;
       k:=i;
      end;
     if j=k then
      begin
       x:=n div k;
       sum:=((num[k]*x) mod p+num[n mod k]) mod p;
       writeln(a*sum mod p);
       exit;
      end;
    end;
end.