【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 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;
}