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

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

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll long long
#define mod 1000000007
using 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;
}