[bzoj3202] [SDOI2013]项链

Description

项链是人体的装饰品之一，是最早出现的首饰。项链除了具有装饰功能之外，有些项 链还具有特殊显示作用，如天主教徒的十字架
链和佛教徒的念珠。 从古至今人们为了美化人体本身，也美 化环境，制造了各种不同风格，不同特点、不同式样的项链，满足了不同肤色、不同民族、不同审美观的人的审美需要。就材料而论，首
饰市场上的项链有黄金、白银、珠宝等几种。珍珠项链为珍珠制成的饰品，即将珍珠 钻孔后用线串在一起，佩戴于项间。天然珍珠项链具有一定的护养作用。

最近，铭铭迷恋上了一种项链。与其他珍珠项链基本上相同，不过这种项链的珠子却 与众不同，是正三菱柱的泰山石雕刻而成的。三菱柱的侧面是正方形构成的，上面刻有数字。 能够让铭铭满意的项链必须满足下面的条件：
1：这串项链由n颗珠子构成的。
2：每一个珠子上面的数字x,必须满足0<x<=a，且珠子上面的数字的最大公约数要恰 好为1。两个珠子被认为是相同的，当且仅当他们经过旋转，或者翻转后能够变成一样的。 3：相邻的两个珠子必须不同。
4：两串项链如果能够经过旋转变成一样的，那么这两串项链就是相同的！ 铭铭很好奇如果给定n和a，能够找到多少不同串项链。由于答案可能很大，所以对输 出的答案mod 1000000007。

Input

数据由多组数据构成：
第一行给定一个T<=10，代表由T组数据。
接下来T行，每行两个数n和a。

Output

对于每组数据输出有多少不同的串。

Sample Input

1 
2 2

Sample Output

3 

Solution

丧心病狂的数论二合一题。。

可以分成两部分来考虑，珠子的种类数和组成环的方案数。

---

珠子的种类数

设\(s_3\)为最多三元环的种类数，即：

\[s_3=\sum_{i=1}^a\sum_{j=1}^a\sum_{k=1}^a[\gcd(i,j,k)=1] \]

\(s_2,s_1\)依次类推，即：

\[s_2=\sum_{i=1}^a\sum_{j=1}^a[\gcd(i,j)=1],s_1=1 \]

可以容斥得到：

总种类数\(=1+(3s_2-3)/3+(s_3-(3s_2-3)-1)/6=(s_3+3s_2+2)/6.\)

然后莫比乌斯反演求一下\(s\)就行了，式子：

\[s_3=\sum_{d=1}^a\mu(d)\lfloor\frac{a}{d}\rfloor^3,s_2=\sum_{d=1}^a\mu(d)\lfloor\frac{a}{d}\rfloor^2 \]

---

环的方案数

设先前求出来的方案数为\(m\)。

看到环旋转前后属于同一种方案，可以想到\(polya\)定理：

\[ans=\sum_{i=1}^{n}f(\gcd(n,i))=\sum_{d|n}\varphi(\frac{n}{d})f(d)=\sum_{d|n}\varphi(d)f(\frac{n}{d}) \]

由于有限制：相邻的颜色不同，设\(f(n)\)为\(n\)个点的环满足条件的涂色方案，则有：

\[f(n)=(m-2)f(n-1)+(m-1)f(n-2) \]

这个式子可以理解为：枚举一个合法方案的断点，两边颜色一定不同，这时候要么当前方案是由\((n-1)\)个点组成的，只有\((m-2)\)种颜色可选，要么由\((n-2)\)个点组成的，加一个与一边相同的点，然后只有\((m-1)\)中颜色选。

然后发现这是一个二阶齐次线性方程，可以搬出特征方程的套路，设这是一个等比数列，可得特征方程：

\[x^2=(m-2)x+m-1 \]

解得：

\[x_1=m-1,x_2=-1 \]

设当前数列通项公式为：

\[f(n)=\alpha x_1^n+\beta x_2^n \]

由于：

\[f(1)=0,f(2)=m^2-m \]

列出方程：

\[\begin{cases} \alpha(m-1)-\beta=0\\ \alpha(m-1)^2+\beta=m^2-m \end{cases} \]

解得\(\alpha=1,\beta=m-1\)，所以\(f\)的通项公式为：

\[f(n)=(m-1)^n+(-1)^n(m-1) \]

答案是：

\[\sum_{d|n}\varphi(\frac{n}{d})f(d) \]

所以直接一路暴力就好了。

记得\(long\,\,long\)会乘爆，用爆\(long\,\,long\)小技巧就好了，具体看代码\(mul\)函数。

由于\(n\)有可能是\(1e9+7\)的倍数，所以特判下，如果是就模\((1e9+7)^2\)，最后除掉\(1e9+7\)即可。

然后细节还有挺多的，，注意下写法，不然调试很恶心。

#include<bits/stdc++.h>
using namespace std;

#ifdef ONLINE_JUDGE
#define getchar() ((p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin)),p1==p2)?EOF:*p1++)
#endif

namespace fast_IO {
	char buf[1<<21],*p1=buf,*p2=buf;

	template <typename T> inline void read(T &x) {
		x=0;T f=1;char ch=getchar();
		for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
		for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
	}
	template <typename T,typename... Args> inline void read(T& x,Args& ...args) {
		read(x),read(args...);
	}

	char buf2[1<<21],a[80];int p,p3=-1;

	inline void flush() {fwrite(buf2,1,p3+1,stdout),p3=-1;}
	template <typename T> inline void write(T x) {
		if(p3>(1<<20)) flush();
		if(x<0) buf2[++p3]='-',x=-x;
		do {a[++p]=x%10+48;} while(x/=10);
		do {buf2[++p3]=a[p];} while(--p);
		buf2[++p3]='\n';
	}
	template <typename T,typename... Args> inline void write(T x,Args ...args) {
		write(x),write(args...);
	}
}

using fast_IO :: read;
using fast_IO :: write;
using fast_IO :: flush;

#define sqr(x) mul(x,x)
#define cul(x) mul(sqr(x),x)

#define ll long long 
#define lf long double
#define ull unsigned long long

const int maxn = 1e7+10;
const int MOD = 1e9+7;
const int MAXN = 1e4+10;

ll mod;

int pri[maxn],vis[maxn],mu[maxn],tot;

void sieve(int N) {
	mu[1]=1;
	for(int i=2;i<=N;i++) {
		if(!vis[i]) pri[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pri[j]<=N;j++) {
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
			mu[i*pri[j]]=-mu[i];
		}
	}
	for(int i=1;i<=N;i++) mu[i]+=mu[i-1];
}

ll mul(ll a,ll b) {
	lf res=(lf)a*b/mod;
	ll t=a*b,t2=(t-(ull)res*mod+mod)%mod;return t2;
}

ll qpow(ll a,ll x) {
	ll res=1;
	for(;x;x>>=1,a=mul(a,a)) if(x&1) res=mul(res,a);
	return res;
}

ll a[MAXN],p[MAXN],cnt;

ll m,ans,inv6,nown;

void mobius(int n) {
	m=0;int T=1;
	while(T<=n) {
		int pre=T;T=n/(n/T);
		m=(m+1ll*mul((mu[T]-mu[pre-1]),cul(n/T)))%mod;
		m=(m+1ll*mul((mu[T]-mu[pre-1]),sqr(n/T))*3ll%mod)%mod;T++;
	}m=(m+2)%mod;m=mul(m,inv6)%mod;
}

void prepare(ll n) {
	cnt=0;
	for(int i=2;1ll*i*i<=n;i++) {
		if(n%i) continue;
		a[++cnt]=i;p[cnt]=0;
		while(n%i==0) n/=i,p[cnt]++;
	}
	if(n!=1) a[++cnt]=n,p[cnt]=1;
}

ll f(ll n) {
	ll Ans=qpow(m-1,n);
	if(n&1) Ans=(Ans+mod-m+1)%mod;
	else Ans=(Ans+m-1)%mod;
	return Ans;
}

void polya(int now,ll d,ll phi) {
	if(now==cnt+1) return ans=(ans+mul(phi,f(nown/d)))%mod,void();
	polya(now+1,d,phi);
	d*=a[now],phi*=a[now]-1;
	polya(now+1,d,phi);
	for(int i=2;i<=p[now];i++)
		d*=a[now],phi*=a[now],polya(now+1,d,phi);
}

ll inn[20],ina[20],mx,inv[10];

void pre_inv() {
	inv[1]=1;
	for(int i=2;i<=6;i++) inv[i]=mul(mod-mod/i,inv[mod%i]);
	inv6=inv[6];
}

int main() {
	int t;read(t);
	for(int i=1;i<=t;i++) read(inn[i],ina[i]),mx=max(mx,ina[i]);
	sieve(mx);
	for(int i=1;i<=t;i++) {
		nown=inn[i];
		if(inn[i]%MOD) mod=MOD;
		else mod=1ll*MOD*MOD;
		pre_inv();
		mobius(ina[i]);
		prepare(inn[i]);ans=0;
		polya(1,1,1);
		if(inn[i]%MOD) ans=mul(ans,qpow(inn[i],mod-2));
		else ans=mul(ans/MOD,qpow(inn[i]/MOD,MOD-2))%MOD;
		write(ans);
	}
	flush();
	return 0;
}