【HDU2865】Birthday Toy-Burnside引理+数论+DP矩阵优化

测试地址：Birthday Toy
题目大意：要给一种轮状玩具着色，这种轮状玩具外围是环形的，由N(≤109)颗小珠子组成，中间有一颗大珠子，外围环上相邻的珠子之间有连边，大珠子和所有小珠子之间都有连边。每一颗珠子都要着一个颜色，颜色共有K(≤109)种，有连边的两颗珠子不能是同一种颜色，旋转后相同的着色方案视为相同，问本质不同的着色方案数有多少，对109+7取模。
做法：这一道题需要使用：Burnside引理，欧拉函数，DP+矩阵快速幂优化，乘法逆元。
这一道题大体和POJ2888类似，都限定了相邻珠子的颜色，不同的是这一题限制变得更有规律了，但是K也变得很大，我们就从POJ2888得出的结论开始推。POJ2888题解请看这里。
首先大珠子可以涂任意一种颜色，然后小珠子只能涂剩下的颜色了，为了方便，接下来我们令K=K−1。
在这一题里，矩阵M′除了对角线元素为0外，其他元素均为1，我们能不能利用这个特殊性质求出M′d的对角线元素呢？
我们可以找到一个规律：M′dij=∑k≠jM′d−1ik，用这种方法写出几个矩阵，根据观察可以推断（或者用数学归纳法证明，我懒得证了），M′d的对角线元素都相等，非对角线元素也相等，且对角线元素与非对角线元素之间的差总为1或−1，而且这两个差交替出现。所以我们设M′d的对角线元素为m′(d)，可以得到递推式：
m′(2i)=(K−1)2×m′(2(i−1))−(K−1)(K−2)
m′(2i+1)=(K−1)×(m′(2i)−1)
其中i为正整数，m′(0)=1。
其实好像还有更简单的递推式，但是我使用了更简单粗暴的分析方法，大家就将就着看吧……
那么很显然这个递推式就可以使用矩阵加速优化了，其余的关于利用欧拉函数优化计算Burnside公式的时间复杂度等内容和POJ2888相同，上文已经给了链接，这里就不赘述了。注意负数取模。
犯二的地方：N可能是一个完全平方数！！！如果是使用枚举d，然后分别计算d和N/d的答案这种方法，不加判断的话，就意味着如果枚举到了N−−√，N−−√的答案会被计算两次，我写POJ2888时没注意这个居然过了，可能是数据比较水吧（现在已经修改）。
以下是本人代码：

#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>#define ll long long#define mod 1000000007using namespace std;ll N,K,ans;struct matrix {ll s[2][2];} M[40];void exgcd(ll a,ll b,ll &x,ll &y){  ll x0=1,x1=0,y0=0,y1=1;  while(b)  {    ll tmp,q;    q=a/b;    tmp=x0,x0=x1,x1=tmp-q*x1;    tmp=y0,y0=y1,y1=tmp-q*y1;    tmp=a,a=b,b=tmp%b;  }  x=x0,y=y0;}matrix mult(matrix A,matrix B){  matrix S;  memset(S.s,0,sizeof(S.s));  for(int i=0;i<=1;i++)    for(int j=0;j<=1;j++)      for(int k=0;k<=1;k++)        S.s[i][j]=(S.s[i][j]+A.s[i][k]*B.s[k][j])%mod;  return S;}matrix power(ll x){  matrix S;  S.s[0][0]=1,S.s[0][1]=0;  S.s[1][0]=0,S.s[1][1]=1;  int i=0;  while(x)  {    if (x&1) S=mult(S,M[i]);    i++;x>>=1;  }  return S;}ll phi(ll x){  ll s=x;  for(ll i=2;i*i<=x;i++)    if (!(x%i))    {      s=s/i*(i-1);      while(!(x%i)) x/=i;    }  if (x>1) s=s/x*(x-1);  return s;}void solve(ll x){  matrix S=power(x/2);  ll tmp=S.s[0][0]+S.s[0][1];  if (x%2) tmp=((K-1)*((tmp-1)%mod))%mod;  tmp=(K*tmp)%mod;  tmp=(phi(N/x)*tmp)%mod;  ans=(ans+tmp)%mod;}int main(){  while(scanf("%lld%lld",&N,&K)!=EOF)  {    ans=0;    K--;    M[0].s[0][0]=((K-1)*(K-1))%mod;    M[0].s[0][1]=((-(K-1)*(K-2))%mod+mod)%mod;    M[0].s[1][0]=0;    M[0].s[1][1]=1;    for(int i=1;i<=35;i++) M[i]=mult(M[i-1],M[i-1]);    for(ll i=1;i*i<=N;i++)      if (!(N%i))      {        solve(i);        if (i!=N/i) solve(N/i);      }    ll x0,y0;    exgcd(N,mod,x0,y0);    x0=(x0%mod+mod)%mod;    printf("%lld\n",(((x0*ans)%mod)*(K+1))%mod);  }  return 0;}