hdu 2865

是这位大牛的思路：http://blog.csdn.net/wukonwukon/article/details/7215467

分析；首先，由于中间一个大圆与每个小圆都相连，所以大圆用去一种颜色之后，只剩下K-1种颜色。
           设K-1种颜色染N个珠子的不同方案数为M，最后就是求M×K mod 1000000007。
           方法跟pku 2888一样，但是这次矩阵的规模很大，所以不能用矩阵来存这个图形了。
           但由于此处规定相邻珠子颜色不同， 则邻接阵为对角线上元素全为0,，其余元素全为1。
           该矩阵的幂的迹，可以推导出公式 ( p - 1 ) ^ n + ( -1 ) ^ n * ( p - 1 )其中p是矩阵的阶数，也就是K-1。
    这个公式是怎么求出来的呢？？？？？？
    几乎所有的日志中都是这个公式，但是没见到有解释怎么求的这个公式，我说说我的想法：
    假设A的n-1次幂为：

           设K-1种颜色染N个珠子的不同方案数为M，最后就是求M×K mod 1000000007。
           方法跟pku 2888一样，但是这次矩阵的规模很大，所以不能用矩阵来存这个图形了。
           但由于此处规定相邻珠子颜色不同， 则邻接阵为对角线上元素全为0,，其余元素全为1。
           该矩阵的幂的迹，可以推导出公式 ( p - 1 ) ^ n + ( -1 ) ^ n * ( p - 1 )其中p是矩阵的阶数，也就是K-1。
    这个公式是怎么求出来的呢？？？？？？
    几乎所有的日志中都是这个公式，但是没见到有解释怎么求的这个公式，我说说我的想法：
    假设A的n-1次幂为：

其中x_n是对角线上的值。乘以对角线上全0，其余为1的矩阵后。
                                                                                 则1)    x_n = y_n-1*(p-1);
                                                                                    2)    y_n = x_n-1+y_n-1*(p-2);
上面两个式子可以解出来x_n = (p-2)*x_n-1 + (p-1)*x_n-2;
事实上，到这一步就能解了，利用矩阵的乘法，然后快速求出x_n的值，进而求出矩阵幂的迹。

利用特征根解出c1=(p-1)/p   ;c2=-1/p;

x_n=( ( p - 1 ) ^ n + ( -1 ) ^ n * ( p - 1 ))/p;

n个x_n相加，结果为  有n个置换的时候有多少种串珠子的方案数。

#include<iostream>
#include<cstring>
#include<queue>
#include<cstdio>
#include<cmath>
#include<algorithm>
#define ll long long
#define mod 1000000007LL
#define PP 40000
using namespace std;
int prime[PP],isprime[PP],size,p,n,k;
int sieve(int n)
{
    int p=0,m;
    memset(isprime,1,sizeof(isprime));
    m=sqrt((double)n);
    for(int i=2;i<=m;i++)
    {
      if(isprime[i])
      prime[p++]=i;
      for(int j=0;i*prime[j]<=m&&j<p;j++)
      {
          isprime[i*prime[j]]=0;
          if(i%prime[j]==0)
          break;
      }
    }
    isprime[0]=isprime[1]=0;
    return p;
}
ll phi(ll n)
{
    ll rea=n;
    for(int i=0;i<size,prime[i]*prime[i]<=n;i++)
    {
        if(n%prime[i]==0)
        rea-=rea/prime[i];
        while(n%prime[i]==0)
        {
            n/=prime[i];
        }
    }
    if(n>1)
    rea-=rea/n;
    return rea%mod;
}
ll pow(ll a,ll n)
{
    ll ans=1;
    while(n)
    {
        if(n&1)
        ans=(ans*a)%mod;
        a=(a*a)%mod;
        n>>=1;
    }
    return ans%mod;
}
ll geter(ll p,ll n)
{
    ll ans;
    ans=pow(p-1,n);
    if(n&1)
    ans=(ans+mod-(p-1))%mod;
    else
     ans=(ans+(p-1))%mod;
     return ans;
}
ll inv(ll n)
{
    return pow(n,mod-2)%mod;
}
ll polya(ll p,ll n)
{
    ll ans=0;
    for(ll i=1;i*i<=n;i++)
    {
        if(n%i==0)
     {
          ans=(ans+(geter(p,i)*phi(n/i))%mod)%mod;
      if(i*i!=n)
    ans=(ans+(geter(p,n/i)*phi(i))%mod)%mod;
     }
    }
    return ans*inv(n)%mod;
}
int main()
{
    size=sieve(mod);
    while(scanf("%I64d%I64d",&n,&k)!=EOF)
    {
        printf("%I64d\n",(k*polya(k-1,n))%mod);
    }
    return 0;
}