传送门
求每个珠子的方案数
即有序的求三元组
(
x
,
y
,
z
)
,
x
,
y
,
z
≤
a
(x,y,z),x,y,z\le a
(x,y,z),x,y,z≤a 满足
g
c
d
(
x
,
y
,
z
)
=
1
gcd(x,y,z)=1
gcd(x,y,z)=1
设
G
i
G_i
Gi 表示
i
i
i 个小于等于
a
a
a 的有序数字,满足
g
c
d
=
1
gcd=1
gcd=1 的方案数
容斥得到要求的
1
6
(
G
3
+
2
G
2
+
3
G
1
)
\frac{1}{6}(G_3+2G_2+3G_1)
61(G3+2G2+3G1)
然后
G
1
=
1
G_1=1
G1=1
运用简单莫比乌斯反演得到
G
2
=
∑
i
=
1
a
⌊
a
i
⌋
2
μ
(
i
)
G_2=\sum_{i=1}^{a}\lfloor\frac{a}{i}\rfloor^2\mu(i)
G2=i=1∑a⌊ia⌋2μ(i)
G
3
=
∑
i
=
1
a
⌊
a
i
⌋
3
μ
(
i
)
G_3=\sum_{i=1}^{a}\lfloor\frac{a}{i}\rfloor^3\mu(i)
G3=i=1∑a⌊ia⌋3μ(i)
求项链条数
运用
P
o
l
y
a
Polya
Polya 定理
设
f
(
x
)
f(x)
f(x) 表示
x
x
x 的点的环,选择上面求出的
m
m
m 种颜色,同色不相邻的方案数
那么要求的就是
1
n
∑
i
=
1
n
f
(
g
c
d
(
i
,
n
)
)
=
1
n
∑
d
∣
n
f
(
d
)
φ
(
n
d
)
\frac{1}{n}\sum_{i=1}^{n}f(gcd(i,n))=\frac{1}{n}\sum_{d|n}f(d)\varphi(\frac{n}{d})
n1i=1∑nf(gcd(i,n))=n1d∣n∑f(d)φ(dn)
求f
容斥不好做
朴素想法是枚举开始和结尾的颜色,显然也不好做
考虑增量算
f
(
x
)
f(x)
f(x)
首先可以断开
x
−
1
x-1
x−1 的链
x
−
1
x-1
x−1 的首尾不同,贡献为
(
m
−
2
)
f
(
x
−
1
)
(m-2)f(x-1)
(m−2)f(x−1)
x
−
1
x-1
x−1 的首尾相同,贡献为
(
m
−
1
)
f
(
x
−
2
)
(m-1)f(x-2)
(m−1)f(x−2)
那么
f
(
x
)
=
(
m
−
2
)
f
(
x
−
1
)
+
(
m
−
1
)
f
(
x
−
2
)
f(x)=(m-2)f(x-1)+(m-1)f(x-2)
f(x)=(m−2)f(x−1)+(m−1)f(x−2)
本题应该是默认
f
(
1
)
=
0
f(1)=0
f(1)=0,不然过不了样例
直接构造生成函数
F
(
x
)
=
∑
i
=
0
x
f
(
i
)
x
i
F(x)=\sum_{i=0}^{x}f(i)x^i
F(x)=∑i=0xf(i)xi
那么
F
(
x
)
=
(
m
−
2
)
F
(
x
)
x
+
(
m
−
1
)
F
(
x
−
2
)
x
2
+
f
(
2
)
F(x)=(m-2)F(x)x+(m-1)F(x-2)x^2+f(2)
F(x)=(m−2)F(x)x+(m−1)F(x−2)x2+f(2)
所以
F
(
x
)
=
m
−
1
1
−
(
m
−
1
)
x
−
m
−
1
1
+
x
F(x)=\frac{m-1}{1-(m-1)x}-\frac{m-1}{1+x}
F(x)=1−(m−1)xm−1−1+xm−1
f
(
x
)
=
(
m
−
1
)
x
−
(
−
1
)
x
−
1
(
m
−
1
)
f(x)=(m-1)^x-(-1)^{x-1}(m-1)
f(x)=(m−1)x−(−1)x−1(m−1)
最后
注意到
n
n
n 可能是
1
0
9
+
7
10^9+7
109+7 的倍数
可以考虑对
(
1
0
9
+
7
)
2
(10^9+7)^2
(109+7)2 取模
如果是倍数,初始答案算出来和
n
n
n 一起除去
m
o
d
mod
mod 再求逆元即可
# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn(1e7 + 5);
const int mod(1e9 + 7);
const ll dmod((ll)mod * mod);
const ll inv6(833333345000000041ll);
inline void Inc(ll &x, ll y) {
x = x + y >= dmod ? x + y - dmod : x + y;
}
inline ll Mul(ll x, ll y) {
return (x * y - (ll)(((long double)x * y + 0.5) / (long double)dmod) * dmod + dmod) % dmod;
}
inline ll Pow1(ll x, ll y) {
register ll ret = 1;
for (x %= mod, y %= mod - 1; y; y >>= 1, x = x * x % mod)
if (y & 1) ret = ret * x % mod;
return ret;
}
inline ll Pow2(ll x, ll y) {
register ll ret = 1;
for (; y; y >>= 1, x = Mul(x, x)) if (y & 1) ret = Mul(ret, x);
return ret;
}
int test, pr[maxn / 10], tot, mu[maxn], cnt;
bitset <maxn> ispr;
ll n, a, ret, ans, d[maxn], ct[maxn];
inline ll Calc(ll x) {
register ll v = Pow2(ret, x);
(x & 1) ? Inc(v, dmod - ret) : Inc(v, ret);
return v;
}
void Dfs(int x, ll v, ll phi) {
if (x > cnt) {
Inc(ans, Mul(phi, Calc(n / v)));
return;
}
register int i;
Dfs(x + 1, v, phi), v = v * d[x], phi = phi * (d[x] - 1), Dfs(x + 1, v, phi);
for (i = 2; i <= ct[x]; ++i) v *= d[x], phi *= d[x], Dfs(x + 1, v, phi);
}
inline void Solve() {
register ll i, j, x;
scanf("%lld%lld", &n, &a), ans = 0, ret = 2;
for (i = 1; i <= a; i = j + 1) {
j = a / (a / i);
Inc(ret, Mul(Mul(Mul((a / i) + 3, a / i), a / i), (mu[j] - mu[i - 1] + dmod) % dmod));
}
ret = Mul(ret, inv6), Inc(ret, dmod - 1);
for (i = 1, cnt = 0, x = n; i <= tot && pr[i] <= x / pr[i]; ++i)
if (x % pr[i] == 0) {
d[++cnt] = pr[i], ct[cnt] = 0;
while (x % pr[i] == 0) x /= pr[i], ++ct[cnt];
}
if (x > 1) d[++cnt] = x, ct[cnt] = 1;
Dfs(1, 1, 1);
if (n % mod) ans = Mul(ans, Pow2(n, dmod - 2)), ans %= mod;
else ans = (ans / mod) % mod * Pow1(n / mod, mod - 2) % mod;
printf("%lld\n", ans);
}
int main() {
register int i, j;
ispr[1] = 1, mu[1] = 1;
for (i = 2; i < maxn; ++i) {
if (!ispr[i]) pr[++tot] = i, mu[i] = -1;
for (j = 1; j <= tot && i * pr[j] < maxn; ++j) {
ispr[i * pr[j]] = 1;
if (i % pr[j]) mu[i * pr[j]] = -mu[i];
else {
mu[i * pr[j]] = 0;
break;
}
}
}
for (i = 2; i < maxn; ++i) mu[i] += mu[i - 1];
scanf("%d", &test);
while (test) Solve(), --test;
return 0;
}