【SDOI2013】项链 题解

【SDOI2013】项链

Problem

众所周知。

Solution

将原问题分为两个问题求解。

Part 1

首先求珍珠的种类数。
设$f_i$表示满足$gcd=i$的本质不同珍珠个数,
$g_i$表示满足$gcd$为$i$的倍数的本质不同珍珠个数
则$f_1$就是答案
由定义可得$$g(i)=\sum _{i|d}f(d)$$
$mobius$反演得到$$f(i)=\sum _{i|d}\mu (\frac{d}{i})g(d)$$
所以$$f(1)=\sum _{i=1}^a\mu (i)g(i)$$
下面计算$g(i)$
$i$的倍数有$\lfloor \frac{a}{i}\rfloor$个
考虑旋转和翻折产生的六个置换
由$Burnside$引理有$$g(i)=\frac{\lfloor \frac{a}{i}{\rfloor }^3+3\lfloor \frac{a}{i}{\rfloor }^2+2\lfloor \frac{a}{i}\rfloor }{6}$$
数论分块即可

Part 2

知道了珍珠的种类数后求本质不同的项链数
先设珍珠的种类数为$c$
不考虑本质不同的条件,先设$h(i)$为相邻不同，首尾不同项链数
初值$$h(1)=0,h(2)=c(c-1)$$
考虑转移 讨论放入一个珍珠，左右两边是否相同
若左右相同,则贡献为$(c-1)h(i-2)$
若左右不同,则贡献为$(c-2)h(i-1)$
综上$$h(i)=(c-2)h(i-1)+(c-1)h(i-2)$$
注意到这是一个二阶常系数递推关系,特征方程法求得通解为$$h(i)=(c-1{)}^i+(-1{)}^i(c-1)$$
考虑本质不同的方案数
设旋转i位,$x_j$为第j位的种类
则$$x_j=x_{(j+i)modn}$$
将其看作第$j$个点向第$(j+i)modn$个点连一条边
则整个项链成为许多个环,环内个点选的种类相同
环的长度为$\frac{in}{lcm(i,)}=\frac{n}{gcd(i,n)}$
每个环互不相交,因此只需考虑$\frac{n}{\frac{n}{gcd(i,n)}}=gcd(i,n)$个珍珠的$h()$值
由Burnside引理
$$ans=\frac{1}{n}\sum _{i=1}^{n}h(gcd(i,n))$$
枚举$gcd(i,n)$,上式成为
$$\begin{array}{rl}ans& =\frac{1}{n}\sum _{d|n}\sum _{i=1}^n[gcd(i,n)==d]h(d)\\ & =\frac{1}{n}\sum _{d|n}\sum _{i=1}^n[gcd(\frac{n}{d},\frac{i}{d})==1]h(d)\\ & =\frac{1}{n}\sum _{d|n}\phi (\frac{n}{d})h(d)\end{array}$$
大概的思路到这里就结束了.但是有些细节需要注意.比如$\phi (\frac{n}{d})$无法直接求,需要分解质因数.模数可能会整除$n$
丑码贴在下面

#include <cstdio>
#include <iostream>
#include <cstring>
using namespace std;
typedef long long LL;
const LL mod = 1e9 + 7;
const LL MOD = mod * mod;
LL c, n, cnt, ans, a, y[10000][2];
int vis[10000010], prime[10000000], cnt2, T, mu[10000010], sum[10000010];
LL read() {
	LL res, flag = 1; char ch = getchar();
	for(; ch > '9' || ch < '0'; ch = getchar())
		if(ch == '-') flag = -1;
	for(res = 0; ch >= '0' && ch <= '9'; res = res * 10 + LL(ch - '0'), ch = getchar());
	return res * flag;
}
LL mul(LL x, LL y, LL mo) {
	LL z = (long double) x * y / mo + 0.5; 
	z = x * y - z * mo;
	return (z % mo + mo) % mo;
}
LL Qpower(LL x, LL pnt, LL mo){
	LL res = 1;
	x = x % mo;
	while(pnt) {
		if(pnt & 1) res = mul(res, x, mo);
		x = mul(x, x, mo);
		pnt >>= 1;
	}
	return res;
}
LL F(LL x) {return (Qpower(c - 1, x, MOD) + (MOD + (x & 1 ? -1 : 1) * (c - 1)) % MOD) % MOD;}
void dfs(int k, LL p, LL x) {
	if(k > cnt2) {
		ans = (ans + mul(p, F(n / x), MOD)) % MOD;
		return ;
	}
	dfs(k + 1, p, x);
	for(int i = 1; i < y[k][1]; ++i) {
		p /= y[k][0];
		x /= y[k][0];
		dfs(k + 1, p, x);
	}
	dfs(k + 1, p / (y[k][0] - 1), x / y[k][0]);
}
int main() {
	T = read();
	mu[1] = 1, sum[1] = 1;
	for(int i = 2; i <= 1e7; ++i) {
		if(!vis[i]) prime[++cnt] = i, mu[i] = -1;
		for(int j = 1; j <= cnt && i * prime[j] <= 1e7; ++j) {
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) break;
			mu[i * prime[j]] = - mu[i];
		}
		sum[i] = sum[i - 1] + mu[i];
	}
	while(T--) {
		ans = c = cnt2 = 0;
		memset(y, 0, sizeof(y));
		n = read(), a = read();
		for(LL l = 1, r; l <= a; l = r + 1) {
			LL m = a / l;
			r = a / m;
			c = (c + mul((sum[r] - sum[l - 1] + MOD) % MOD, mul(((mul(m, mul(m, m, MOD), MOD) + 3 * mul(m, m, MOD) % MOD) % MOD + m * 2 % MOD) % MOD, 833333345000000041ll, MOD), MOD)) % MOD;
		}
		LL tmp = n, p = n;
		for(int i = 1; i <= cnt && prime[i] <= tmp; ++i) {
			if(tmp % prime[i] == 0) {
				y[++cnt2][0] = prime[i];
				p = p / prime[i] * (prime[i] - 1);
				while(tmp % prime[i] == 0) {
					y[cnt2][1] ++;
					tmp /= prime[i];
				}
			}
		}
		if(tmp > 1) {
			p = p / tmp * (tmp - 1);
			y[++cnt2][0] = tmp;
			y[cnt2][1] = 1;
		}
		dfs(1, p, n);
		if(n % mod == 0) ans = ans / mod * Qpower(n / mod, mod - 2, mod) % mod;
		else ans = ans % mod * Qpower(n, mod - 2, mod) % mod;
		printf("%lld\n",ans);
	}
	return 0;
}