HDU-5868 Different Circle Permutation

做这道题，需要先熟悉polya定理及欧拉函数优化，矩阵快速幂，求逆元

题目大意转化之后就是，有n个位置，可以涂黑色或白色两种颜色，黑色不能相邻

问排除旋转后相同的方案，总共有多少种涂色方案

那么和一般的涂色问题相比，多了一个黑色不能相邻的条件，只需要稍微改一下polya定理的公式

求出每种旋转的循环节个数C=gcd(n,i)后，原来的m^C改为f(C)

f(n)指不考虑旋转同构的情况下给n个位置涂色的方案数，可以推导得f(n)=f(n-1)+f(n-2)

//（由此能看出来polya定理中，m^C也相当于不考虑同构时对C个位置的涂色方案数）

f(1)=1,f(2)=3,但是当n==1时，其实有涂为黑色和白色两种情况，所以n==1时特判一下

那么总结下，n<=1e9比较大，所以用矩阵快速幂求得f(n)，再用欧拉函数优化求循环节gcd(n,i)

算出来代入修改过的polya公式中计算，中间用了扩展欧几里得求n的逆元

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
ll euler(ll n)
{
     ll res=n,a=n;
     for(int i=2;i*i<=a;i++)
	 {
         if(a%i==0)
		 {
             res=res/i*(i-1);
             while(a%i==0) a/=i;
         }
     }
     if(a>1) res=res/a*(a-1);
     return res;
}
const ll MOD=1e9+7;
const int N=2;
struct node
{
	ll a[10][10];
};
node shu,ans,mp;
//shu是输入的矩阵，ans是所求答案
node matrix(node x,node y)
{
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++){
            mp.a[i][j]=0;
            for(int p=1;p<=N;p++)
                mp.a[i][j]=(mp.a[i][j]+x.a[i][p]*y.a[p][j]+MOD)%MOD;
            //矩阵乘法 
        }
    return mp;
}
void work(ll k)
{//矩阵快速幂 
	for(int i=1;i<=N;i++)
			for(int j=1;j<=N;j++)
				ans.a[i][j]=0;
	for(int i=1;i<=N;i++) ans.a[i][i]=1;
	node t=shu;
    while(k){
        if(k&1) 
            ans=matrix(ans,t);
        k>>=1;
        t=matrix(t,t);
    }
}
ll get(ll n)
{
	if(n==1) return 1;
	if(n==2) return 3;
	shu.a[1][1]=1;
	shu.a[1][2]=1;
	shu.a[2][1]=1;
	shu.a[2][2]=0;
	work(n-2);
	return (ans.a[1][1]*3+ans.a[1][2]*1)%MOD;
}
ll e_gcd(ll a, ll b, ll& x, ll& y)
{
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    }
    ll ans = e_gcd(b, a%b, y, x);
    y -= x*(a/b);
    return ans;
}
int main()
{
	ll n;
	while(~scanf("%lld",&n))
	{
		if(n==1) //特判n==1
		{
			printf("2\n");
			continue;
		}
		ll ans=0,i;
		for(i=1;i*i<n;i++)
			if(n%i==0)
			{
				ans+=euler(i)*get(n/i)%MOD+euler(n/i)*get(i)%MOD;
			}
		if(i*i==n) ans+=euler(i)*get(i)%MOD;
		ans%=MOD;
		ll x,y;
		e_gcd(n,MOD,x,y); //扩展欧几里得求逆元
		x=(x%MOD+MOD)%MOD; //求n的逆元x
		ans=(ans*x%MOD+MOD)%MOD;
		printf("%lld\n",ans);
	}
	return 0;
}