SDOI2013 项链

Vexoben

Mar 11, 2019

It takes 8 minutes to read this article.

题目链接： SDOI2013 项链

（刷新以获取数学公式）

题意

有一种项链，它的珠子是正三菱柱。三菱柱的侧面是正方形构成的，上面刻有数字。

给定 $n$ 和 $a$ ，求满足条件的项链种数：

1：这串项链由 $n$ 颗珠子构成的。

2：每一个珠子上面的数字 $x$ ,必须满足 $0<x≤a$ ，且珠子上面的数字的最大公约数要恰好为 $1$ 。两个珠子被认为是相同的，当且仅当他们经过旋转，或者翻转后能够变成一样的。

3：相邻的两个珠子必须不同。

4：两串项链如果能够经过旋转变成一样的，那么这两串项链就是相同的！

由于答案可能很大，所以对输出的答案 $mod \;1000000007$ 。

数据范围： $n≤10^{14}$ ， $a≤10^7$ ， $T≤10$

题解

一道模板题大综合。

先计算 $m＝$ 珠子的种类数。即计算 $\sum_{i=1}^{n}\sum_{j=i}^{n}\sum_{k=j}^{n}[gcd(i, j, k) = 1]$

记 $f(x)=\lfloor \frac{n}{x} \rfloor +2C_{\lfloor \frac{n}{x} \rfloor}^{2}+C_{\lfloor \frac{n}{x} \rfloor }^{3}$ 为 $gcd$ 是 $x$ 的倍数的方案数，直接莫比乌斯反演就可以算出 $m$ 。

本质不同项链数显然是要套Burside定理。先记 $g(x)$ 为给 $x$ 颗珠子染色，不要求本质不同，但相邻珠子不同的方案数。显然 $g(x)=(m-2)g(x-1)+(m-1)g(x-2),g(1)=0,g(2)=m(m-1)$

特征根解一下就知道 $g(x)=(m-1)^x+(-1)^{x}(m-1)$

根据Burside引理，答案就是

$\sum_{i=1}^{n} g(gcd(i,n)) = \sum_{d|n} g(d)\phi(\frac{n}{d})$

然而由于 $n$ 特别大，可能是 $mod$ 的倍数，因此要先模 $mod^2$ ，最后除以 $n$ 。

中间可能会用到 $O(1)$ 快速乘，去学习了一发：

#define LL long long
#define Ldb long double
LL Qmul(LL x, LL y) {
	return (x * y - (LL) (((Ldb) x * y + 0.5) / (Ldb) MOD) * MOD + MOD) % MOD;
}

据说是因为两个乘法都会溢出，但是溢出的减一减就刚好了。

#include<bits/stdc++.h>
#define LL long long
#define Ldb long double
using namespace std;
const int N = 1e7 + 10;
const int mod = 1e9 + 7;
const LL MOD = 1000000014000000049;
const LL inv6 = 833333345000000041;

LL a, n, m;
int tot, fcnt, pri[N], fct[10005];
short miu[N];
bitset<N> use;

void Init() {
	miu[1] = 1;
	for (int i = 2; i < N; ++i) {
		if (!use[i]) {
			pri[++tot] = i;
			miu[i] = -1;
		}
		for (int j = 1; j <= tot && pri[j] * i < N; ++j) {
			int t = pri[j] * i;
			use[t] = 1;
			miu[t] = (i % pri[j] == 0) ? 0 : -miu[i];
			if (i % pri[j] == 0) break;
		}
	}
}

LL Qmul(LL x, LL y) {
	return (x * y - (LL) (((Ldb) x * y + 0.5) / (Ldb) MOD) * MOD + MOD) % MOD;
}

LL C3(LL n) {
	if (n <= 2) return 0;
	else return Qmul(Qmul(Qmul(n, n - 1), n - 2), inv6);
}

LL calc(LL n) {
	return (C3(n) + n * n) % MOD;
}

LL phi(LL n) {
	LL ans = n;
	for (int i = 1; i <= fcnt; ++i) {
		if (n % fct[i] == 0) {
			ans /= fct[i];
			ans *= fct[i] - 1;
		}
	}
	return ans;
}

LL Qpow(LL x, LL p) {
	LL ans = 1;
	while (p) {
		if (p & 1) ans = Qmul(ans, x);
		x = Qmul(x, x);
		p >>= 1;
	}
	return ans;
}

LL Inv(LL x) {
	LL ans = 1, p = mod - 2;
	while (p) {
		if (p & 1) ans = ans * x % mod;
		x = x * x % mod;
		p >>= 1;
	}
	return ans;
}

LL f(LL n) {
	LL tmp = Qpow(m - 1, n);
	if (n & 1) tmp -= (m - 1);
	else tmp += (m - 1);
	if (tmp >= MOD) tmp -= MOD;
	if (tmp < 0) tmp += MOD;
	return tmp;
}

void Solve() {
	scanf("%lld%lld", &n, &a);
	m = 0;
	for (int i = 1; i <= a; ++i) {
		LL tmp = miu[i] == 0 ? 0 : miu[i] * calc(a / i);
		m += tmp;
		if (m > MOD) m -= MOD;
		if (m < 0) m += MOD;
	}

	if (n == 1) return (void) (printf("%lld\n", m % mod));
    
	fcnt = 0;
	LL tmp = n;
	for (int i = 1; 1LL * pri[i] * pri[i] <= tmp && i <= tot; ++i) {
		if (tmp % pri[i] != 0) contishuxuenue;
		fct[++fcnt] = pri[i];
		while (tmp % pri[i] == 0) tmp /= pri[i];
	}
	if (tmp != 1) fct[++fcnt] = tmp;

	LL ans = 0;
	for (LL i = 1; i * i <= n; ++i) {
		if (n % i != 0) continue;
		ans = (ans + Qmul(f(i), phi(n / i))) % MOD;
		if (i * i != n) ans = (ans + Qmul(f(n / i), phi(i))) % MOD;
	}
	if (n % mod != 0) {
		ans = ans % mod * Inv(n % mod) % mod;
	}
	if (n % mod == 0) {
		ans = ans / mod * Inv(n / mod) % mod;
	}
	printf("%lld\n", ans);
}

int main() {
	int T;
	scanf("%d", &T);
	Init();
	while (T--) Solve();
	return 0;
}