Luogu3307:[SDOI2013]项链

传送门
求每个珠子的方案数
即有序的求三元组 ( x , y , z ) , x , y , z ≤ a (x,y,z),x,y,z\le a (x,y,z),x,y,za 满足 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=1aia2μ(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=1aia3μ(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=1nf(gcd(i,n))=n1dnf(d)φ(dn)
求f
容斥不好做
朴素想法是枚举开始和结尾的颜色,显然也不好做
考虑增量算 f ( x ) f(x) f(x)
首先可以断开 x − 1 x-1 x1 的链
x − 1 x-1 x1 的首尾不同,贡献为 ( m − 2 ) f ( x − 1 ) (m-2)f(x-1) (m2)f(x1)
x − 1 x-1 x1 的首尾相同,贡献为 ( m − 1 ) f ( x − 2 ) (m-1)f(x-2) (m1)f(x2)
那么 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)=(m2)f(x1)+(m1)f(x2)
本题应该是默认 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)=(m2)F(x)x+(m1)F(x2)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(m1)xm11+xm1
f ( x ) = ( m − 1 ) x − ( − 1 ) x − 1 ( m − 1 ) f(x)=(m-1)^x-(-1)^{x-1}(m-1) f(x)=(m1)x(1)x1(m1)
最后
注意到 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;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值