您的位置:首页 > 其它

hdu 5868 矩阵快速幂+burnside引理 +欧拉函数+乘法逆元

2016-10-06 23:00 357 查看
Different Circle Permutation 

方法来自官方题解:
https://async.icpc-camp.org/d/546-2016-icpc
题意:

题意:n 个点构成一个环,每个点可以染色为黑色或白色,要求任意两个相邻的点不能都是黑色,问在旋转同构意义下的染色数。

题解:

设 
f(n)f(n)
 为不考虑旋转同构时给
n 个点的环染色的方案数,显然有递推式 
f(n)
= f(n-1) + f(n - 2)f(n)=f(n−1)+f(n−2)
。这个可以靠矩阵来
O(\log
n)O(logn)
求一次 
f(n)f(n)
。(即矩阵快速幂求解)

考虑旋转则根据 burnside 引理有 
F(n)
= \frac{1}{n}\sum_{i=1}^{n}f(\gcd(i , n))F(n)=​n​​1​​∑​i=1​n​​f(gcd(i,n))


即 
F(n)
= \frac{1}{n}\sum_{d|n}{f(d) \varphi(\frac{n}{d})}F(n)=​n​​1​​∑​d∣n​​f(d)φ(​d​​n​​)
,枚举 
dd
 运算即可,复杂度为
O(\sqrt{n}\log
n)O(√​n​​​logn)


其中,“枚举d”就是枚举n的所有约数,枚举的方法。。哎,太弱了。比如求n的约数:

for(int i=1;i*i<=n;i++)
if(n%i==0){

p[pCnt++]=i;
if(i*i!=n) p[pCnt++]=n/i;
}


sqrt(n)的时间复杂度。

注意还要特判 

n=1n=1
 的情况,因为题意要求是”任意两人”,所以一个人坐在单独的一个座位上仍然是合法的。
最后,还有个很重要的地方,由于所有的求值都是模 MOD=1000000000+7的,而F(n)的公式中有个 除以n的操作,这涉及到求“除n”的乘法逆元。
又因为1000000000+7是质数,所以“除n”的乘法逆元 等于n^(MOD-2)。这个证明见我的博客:
http://blog.csdn.net/yukizzz/article/details/51105009
last and last,AC代码:

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define MOD 1000000007
typedef long long ll;

using namespace std;

struct Mat{
ll mat[2][2];
};

int p[100005];
int pCnt;

Mat operator * (Mat a,Mat b){
Mat c;
memset(c.mat,0,sizeof(c.mat));

for(int k=0;k<2;k++)
for(int i=0;i<2;i++)
for(int j=0;j<2;j++)
c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%MOD;

return c;
}

Mat operator ^ (Mat a,int k){
Mat c;

for(int i=0;i<2;i++)
for(int j=0;j<2;j++)
c.mat[i][j]=(i==j); //单位矩阵

for(;k>0;k>>=1){
if(k&1) c=c*a;
a=a*a;
}

return c;
}

//求fibo数列, O(logn)
ll f(int n){
Mat a;
a.mat[0][0]=a.mat[0][1]=a.mat[1][0]=1;
a.mat[1][1]=0;

if(n==1) return 1;
if(n==2) return 3;

Mat res;
res.mat[1][0]=1;
res.mat[0][0]=3;

res=(a^(n-2))*res;
return res.mat[0][0];
}

//求n的欧拉函数
ll eula(int n){

if(n==1) return 1;
int ret=n;

for(ll i=2;i*i<=n;i++)
if(n%i==0){
ret=ret-ret/i;
while(n%i==0) n/=i;
}

if(n!=1) {ret=ret-ret/n;}

return ret;
}

//求数字n的所有 约数。
void prime(int n){

pCnt=0;
memset(p,0,sizeof(p));

for(int i=1;i*i<=n;i++)
if(n%i==0){

p[pCnt++]=i;
if(i*i!=n) p[pCnt++]=n/i;
}

}

//考虑旋转同构的结果
ll F(int n){
ll ret=0;

prime(n);
sort(p,p+pCnt);

for(ll i=0;i<pCnt;i++)
ret=ret=(ret+( f(p[i])*eula(n/p[i]) ))%MOD;

return ret;
}

ll pow_quick(ll n ,ll k){
ll ret=1;

for(;k>0;k>>=1){

if(k&1) ret=(ret*n)%MOD;
n=(n*n)% MOD;
}
//printf("pow_quick=%lld\n",ret);
return ret;
}

int main(int argc, const char * argv[]) {

int n;

while(~scanf("%d",&n)){

//printf("fn=%lld\n",f(n));
//printf("eula= %d\n",eula(n));

if(n==1) {puts("2");continue;}

ll ans=F(n);
//printf("ans = %lld ",ans);
ans= ans*pow_quick(n,MOD-2)%MOD;
printf("%lld\n",(ans+MOD)%MOD);
}

return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: