[BZOJ3202][SDOI2013]项链
sol
数论多合一。
可以考虑分两步做。
1、求有多少种不同的珠子，设其数量为\(m\)。
2、用\(m\)种颜色对\(n\)个珠子的项链染色，要求相邻不同色，且旋转相同视作同一种方案，求方案数。
求有多少种不同的珠子
\(A=\sum_{i=1}^a\sum_{j=1}^a\sum_{k=1}^a[\gcd(i,j,k)=1]\)？
然而算重了。
如果\(i,j,k\)都不相等的话就多算了\(6\)次，如果有两个相等就多算了\(3\)次，如果三个都相等......那肯定是三个\(1\)啦，并没有多算。
这个问题怎么解决呢？我们再算\(B=\sum_{i=1}^a\sum_{j=1}^a[\gcd(i,j)=1]\)，那么答案就是\(\frac 16(A+3B+2)\)（把每种情况补齐到\(6\)次，再一起除以\(6\)）。
求\(A\)跟\(B\)是很简单的莫比乌斯反演，具体来说就是：
求染色方案
这里要用到\(Burnside\)定理。
这里的置换只有一种，就是旋转。设旋转\(i\)个位置，那么就会有产生\(\gcd(i,n)\)个等价类。
所以答案就是$$\frac{\sum_{i=1}^nf_{\gcd(i,n)}}{n}$$其中\(f_i\)表示\(i\)个点的环染色的方案数。也可以认为是给一个长为\(n\)的序列染色，然后保证首尾不同色的方案数。
考虑对上式做一些简化。因为\(\gcd(i,n)\)一定是\(n\)的约数，所以可以枚举某个约数\(d\)，然后答案就变成了$$\frac{\sum_{d|n}f_d\sum_{i=1}^n[\gcd(i,n)=d]}{n}=\frac{\sum_{d|n}f_d\varphi(\frac nd)}{n}$$
求\(f_n\)
考虑一下\(f_n\)的递推。
讲道理的话，前一种转移表示在末尾填一个与首尾都不相同的颜色，这样有\(m-2\)种选择；后一种转移表示先填一个和首位相同的颜色，再填一个不同色，这样有\(1\times(m-1)\)种选择。
然后就可以矩乘了。
求\(f_n\)的通项公式
有人跟我讲这题要求通项公式我就没去写矩乘。常数优秀一点说不定跑得过？
下面讲一下求通项公式
其特征方程$$x^2-(m-2)x-(m-1)=0$$
解得其特征根$$x_1=-1,x_2=m-1$$
于是设其通项公式为$$f_i=a(-1)i+b(m-1)i$$
将\(f_1=0,f_2=m(m-1)\)代入解得$$a=1-m,b=1$$
所以$$f_i=(1-m)(-1)i+(m-1)i$$
复杂度还是\(O(\log n)\)但是常数小很多。
求\(\varphi(n)\)
最好不要枚举因数然后暴力算\(\varphi(n)\)。
一种可行的方案是：先将\(n\)质因数分解，然后\(dfs\)搜出每一个因数，在搜的过程中自然可以求出它的\(\varphi\)。
这样的话总的复杂度就是\(O(a+T(\sqrt a+d(n)\log n))\)。
你以为这样就做完了吗
\(n\le10^{14}\)
这样可能会产生一个问题：\(n\)是\(P=10^9+7\)的倍数。这样算出来的分子要除以\(n\)，不好意思不存在逆元。
怎么办呢？
我们把答案对\(P^2\)取模。如果\(P|n\)，那么分子既然是\(n\)的倍数就也是\(P\)的倍数。直接除以\(P\)，接着还需要除以\(\frac nP\)，而\(\frac nP\)在模\(P\)意义下是存在逆元的（\(n\le P^2\)），于是就可以乘逆元了。
code
#include<cstdio>
#include<algorithm>
using namespace std;
#define ll long long
const ll mod = 1000000007ll;
const ll MOD = 1000000007ll*1000000007ll;
const ll Inv6 = 833333345000000041ll;
const int N = 10000005;
int T,zhi[N],pri[N/10],tot,mu[N],cnt,q[30];ll n,m,p[30],ans;
void add(ll &x,ll y){x+=y;if(x>=MOD)x-=MOD;}
ll mul(ll x,ll y){
return (x*y-(ll)(((long double)x*y+0.5)/(long double)MOD)*MOD+MOD)%MOD;
}
ll sqr(ll x){return mul(x,x);}
ll cub(ll x){return mul(mul(x,x),x);}
ll fastpow(ll a,ll b){
ll res=1;
while (b) {if (b&1) res=mul(res,a);a=mul(a,a);b>>=1;}
return res;
}
ll fpow(ll a,ll b){
ll res=1;
while (b) {if (b&1) res=res*a%mod;a=a*a%mod;b>>=1;}
return res;
}
void init(){
mu[1]=1;
for (int i=2;i<N;++i){
if (!zhi[i]) pri[++tot]=i,mu[i]=-1;
for (int j=1;i*pri[j]<N;++j){
zhi[i*pri[j]]=1;
if (i%pri[j]) mu[i*pri[j]]=-mu[i];
else break;
}
}
for (int i=1;i<N;++i) mu[i]+=mu[i-1];
}
ll cal(ll n){
ll A=0,B=0;
for (ll i=1,j;i<=n;i=j+1){
j=n/(n/i);
add(A,mul(cub(n/i),mu[j]-mu[i-1]+MOD)%MOD);
add(B,mul(sqr(n/i),mu[j]-mu[i-1]+MOD)%MOD);
}
add(A,mul(B,3));add(A,2);return mul(A,Inv6);
}
ll F(ll n){
ll res=n&1?1-m:m-1;add(res,MOD);
add(res,fastpow(m-1+MOD,n));return res;
}
void dfs(int i,ll cur,ll phi){
if (i==cnt+1) {add(ans,mul(phi,F(n/cur)));return;}
dfs(i+1,cur,phi);
cur*=p[i];phi*=p[i]-1;dfs(i+1,cur,phi);
for (int j=2;j<=q[i];++j)
cur*=p[i],phi*=p[i],dfs(i+1,cur,phi);
}
int main(){
init();
scanf("%d",&T);while (T--){
scanf("%lld%lld",&n,&m);m=cal(m);
ll x=n;cnt=0;
for (int i=1;i<=tot&&(ll)pri[i]*pri[i]<=x;++i)
if (x%pri[i]==0){
p[++cnt]=pri[i];q[cnt]=0;
while (x%pri[i]==0) x/=pri[i],++q[cnt];
}
if (x>1) p[++cnt]=x,q[cnt]=1;
ans=0;dfs(1,1,1);
if (n%mod==0){
ans/=mod;
ans=ans*fpow(n/mod,mod-2)%mod;
}else{
ans%=mod;
ans=ans*fpow(n%mod,mod-2)%mod;
}
printf("%lld\n",ans);
}
return 0;
}
