Description
环上均匀分布n个点,每个点可以坐人和不坐人,相邻两点不能同时坐人,旋转后相等视为同一种方案,问有多少种方案
Input
多组用例,每组用例输入一个整数n,以文件尾结束输入
(不超过50组用例,1<=n<=1e9)
Output
对于每组用例,输出合法方案数
Sample Input
4
7
10
Sample Output
3
5
15
Solution
问题为从环上的n个点中选一些点染色,相邻两点不能同时染色
先不考虑旋转同构,用g[i]表示在不考虑旋转同构的前提下给i个点的环染色的合法方案数,先考虑一条链上的情况,用dp[i][0]表示第i个点不染色的方案数,dp[i][1]表示第i个点染色的方案数,那么由以下转移方程:
dp[1][0]=dp[1][1]=1
dp[i][0]=dp[i-1][0]+dp[i-1][1],dp[i][1]=dp[i-1][0],i=2,3,…,n
显然可以发现dp[i][0]是斐波那契数列的第i+1项,考虑第n个点是否染色,如果不染色那么可以将这个链连成一个环,方案数dp[n][0];如果染色,那么第1个点必然不染,而从第2个点开始到第n个点还是可以看做一条链上的染色,方案数dp[n-1][1]=dp[n-2][0]
故有g[n]=f[n-1]+f[n+1],f[i]表示斐波那契数列第i项,进一步由f[i]=f[i-1]+f[i-2]知g[i]=g[i-1]+g[i-2],其中g[1]=1,g[2]=3
之后考虑旋转,由Burnside引理知合法方案数
枚举d,矩阵快速幂O(logn)求g,求一个数的欧拉函数也是O(logn),故总时间复杂度
注意n=1时特判答案是2,因为一个人的时候第1个人和第n个人可以同时染色
Code
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
#define mod 1000000007ll
typedef long long ll;
struct Mat
{
ll mat[3][3];//矩阵
int row,col;//矩阵行列数
};
Mat mod_mul(Mat a,Mat b,int p)//矩阵乘法
{
Mat ans;
ans.row=a.row;
ans.col=b.col;
memset(ans.mat,0,sizeof(ans.mat));
for(int i=0;i<ans.row;i++)
for(int k=0;k<a.col;k++)
if(a.mat[i][k])
for(int j=0;j<ans.col;j++)
{
ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j]%p;
ans.mat[i][j]%=p;
}
return ans;
}
Mat mod_pow(Mat a,int k,int p)//矩阵快速幂
{
Mat ans;
ans.row=a.row;
ans.col=a.col;
for(int i=0;i<a.row;i++)
for(int j=0;j<a.col;j++)
ans.mat[i][j]=(i==j);
while(k)
{
if(k&1)ans=mod_mul(ans,a,p);
a=mod_mul(a,a,p);
k>>=1;
}
return ans;
}
#define maxn 55555
int euler[maxn],prime[maxn],res;
void Get_euler(int n)//求n以内所有数的欧拉函数值,顺便得到n以内所有素数共res个
{
memset(euler,0,sizeof(euler));
euler[1]=1;
res=0;
for(int i=2;i<=n;i++)
{
if(!euler[i])euler[i]=i-1,prime[res++]=i;
for(int j=0;j<res&&prime[j]*i<=n;j++)
{
if(i%prime[j]) euler[prime[j]*i]=euler[i]*(prime[j]-1);
else
{
euler[prime[j]*i]=euler[i]*prime[j];
break;
}
}
}
}
ll f[maxn];
void init()
{
Get_euler(50000);
f[1]=1,f[2]=3;
for(int i=3;i<=50000;i++)f[i]=(f[i-1]+f[i-2])%mod;
}
ll get_f(int n)
{
if(n<=50000)return f[n];
Mat A;
A.row=A.col=2;
A.mat[0][0]=1,A.mat[0][1]=1,
A.mat[1][0]=1,A.mat[1][1]=0;
A=mod_pow(A,n-2,mod);
return 3ll*A.mat[0][0]+A.mat[0][1];
}
ll get_euler(int n)
{
if(n<=50000)return euler[n];
ll ans=n;
for(int i=2;i*i<=n;i++)
if(n%i==0)
{
ans=ans/i*(i-1);
while(n%i==0)n/=i;
}
if(n>1)ans=ans/n*(n-1);
return ans;
}
ll Mod_pow(ll a,ll b)
{
a%=mod;
ll ans=1ll;
while(b)
{
if(b&1)ans=ans*a%mod;
a=a*a%mod;
b>>=1;
}
return ans;
}
int main()
{
init();
int n;
while(~scanf("%d",&n))
{
if(n==1)
{
printf("2\n");
continue;
}
ll ans=0;
for(int i=1;i*i<=n;i++)
if(n%i==0)
{
ll temp=(get_f(i))%mod*get_euler(n/i)%mod;
ans=(ans+temp)%mod;
if(i*i!=n)
{
int j=n/i;
temp=get_f(j)*get_euler(n/j)%mod;
ans=(ans+temp)%mod;
}
}
printf("%I64d\n",ans*Mod_pow(n,mod-2)%mod);
}
return 0;
}