51NOD1965:奇怪的式子

本文探讨了一种高效的算法,用于计算特定的乘积公式,涉及到数论、组合数学以及算法优化技巧。通过引入质数因子、μ函数和σ函数,文章详细解析了如何简化和加速大数运算过程,特别关注于∏i=1nσ0(i)μ(i)和∏i=1nσ0(i)i的计算。通过预处理和数论分块等技术,实现了对大规模数据集的有效处理。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

传送门
拆开变成
∏ i = 1 n σ 0 ( i ) μ ( i ) ∏ i = 1 n σ 0 ( i ) i \prod_{i=1}^{n}\sigma_0(i)^{\mu(i)}\prod_{i=1}^{n}\sigma_0(i)^{i} i=1nσ0(i)μ(i)i=1nσ0(i)i
考虑 ∏ i = 1 n σ 0 ( i ) μ ( i ) \prod_{i=1}^{n}\sigma_0(i)^{\mu(i)} i=1nσ0(i)μ(i)
运用 μ \mu μ 的性质,设 c ( i ) c(i) c(i) 表示 i i i 的质数因子个数
那么就是
∏ i = 1 n 2 μ ( i ) c ( i ) = 2 ∑ i = 1 n μ ( i ) c ( i ) \prod_{i=1}^{n}2^{\mu(i)c(i)}=2^{\sum_{i=1}^{n}\mu(i)c(i)} i=1n2μ(i)c(i)=2i=1nμ(i)c(i)
只需要设
f ( x , j ) f(x,j) f(x,j) 表示 1 1 1 x x x 最小质因子大于等于第 j j j 个质数的合数的 μ ( i ) c ( i ) \mu(i)c(i) μ(i)c(i) 的和
g ( x , j ) g(x,j) g(x,j) 表示 1 1 1 x x x 最小质因子大于等于第 j j j 个质数的合数的 μ ( i ) \mu(i) μ(i) 的和
预处理出 c n t p ( i ) cntp(i) cntp(i) 表示 1 1 1 x x x 的质数个数就可以递推了
这一部分直接 m i n 25 min25 min25 筛即可
考虑 ∏ i = 1 n σ 0 ( i ) i \prod_{i=1}^{n}\sigma_0(i)^{i} i=1nσ0(i)i
分开考虑每个质数的贡献
s ( n , p , k ) s(n,p,k) s(n,p,k) 表示 1 1 1 n n n 中只含有 p k p^k pk 这个因子不含 p i , i ≠ k p^i,i\ne k pi,i̸=k 的数的和
P P P 为质数集合
那么
∏ i = 1 n σ 0 ( i ) i = ∏ p ∈ P ∏ k = 1 ( k + 1 ) s ( n , p , k ) \prod_{i=1}^{n}\sigma_0(i)^{i}=\prod_{p\in P}\prod_{k=1}(k+1)^{s(n,p,k)} i=1nσ0(i)i=pPk=1(k+1)s(n,p,k)
S u m ( x ) = ∑ i = 1 x i Sum(x)=\sum_{i=1}^{x}i Sum(x)=i=1xi
那么
s ( n , p , k ) = p k S u m ( ⌊ n p k ⌋ ) − p k + 1 S u m ( ⌊ n p k + 1 ⌋ ) s(n,p,k)=p^kSum(\lfloor\frac{n}{p^k}\rfloor)-p^{k+1}Sum(\lfloor\frac{n}{p^{k+1}}\rfloor) s(n,p,k)=pkSum(pkn)pk+1Sum(pk+1n)
p ≤ n p\le \sqrt{n} pn 的时候直接暴力计算
p > n p > \sqrt{n} p>n 的时候,显然 k k k 最多就是 1 1 1
那么
s ( n , p , k ) = p S u m ( ⌊ n p ⌋ ) s(n,p,k)=pSum(\lfloor\frac{n}{p}\rfloor) s(n,p,k)=pSum(pn)
m i n 25 min25 min25 筛预处理出 S ( x ) S(x) S(x) 表示 1 1 1 x x x 的质数的和,然后数论分块即可
注意常数优化QwQ

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

const ll mod(1e12 + 39);
const ll phi(1e12 + 38);
const int maxn(1e6 + 5);

inline ll Mul1(ll x, ll y) {
    register ll ret = x * y - (ll)((long double)x / mod * y + 0.5) * mod;
    return ret < 0 ? ret + mod : ret;
}

inline ll Mul2(ll x, ll y) {
    register ll ret = x * y - (ll)((long double)x / phi * y + 0.5) * phi;
    return ret < 0 ? ret + phi : ret;
}

inline ll Pow(ll x, ll y) {
	register ll ret = 1;
	for (; y; y >>= 1, x = Mul1(x, x))
		if (y & 1) ret = Mul1(ret, x);
	return ret;
}

inline void Inc1(ll &x, ll y) {
	x = x + y >= mod ? x + y - mod : x + y;
}

inline void Inc2(ll &x, ll y) {
	x = x + y >= phi ? x + y - phi : x + y;
}

inline ll Dec1(ll x, ll y) {
	return x - y < 0 ? x - y + mod : x - y;
}

inline ll Dec2(ll x, ll y) {
	return x - y < 0 ? x - y + phi : x - y;
}

int test, pr[maxn], tot, d, id1[maxn], id2[maxn], cnt;
ll n, val[maxn], cntp[maxn], f[maxn], g[maxn], s[maxn], sp[maxn];
bitset <maxn> ispr;

# define ID(x) ((x) <= d ? id1[x] : id2[n / (x)])

inline ll Sum1(ll x) {
	return (x & 1) ? Mul1(x, (x + 1) / 2) : Mul1(x / 2, x + 1);
}

inline ll Sum2(ll x) {
	return (x & 1) ? Mul2(x, (x + 1) / 2) : Mul2(x / 2, x + 1);
}

inline ll GetSum(ll x, ll v1, ll v2) {
	return Dec2(Mul2(v1, Sum2(x / v1)), Mul2(v2, Sum2(x / v2)));
}

inline void Solve() {
	while (cnt) f[cnt] = g[cnt] = 0, cnt--;
	register ll i, j, id, ans1, ans2, v, ret, lst, cur;
	for (d = sqrt(n), i = 1; i <= n; i = j + 1) {
		j = n / (n / i), val[++cnt] = n / i;
		val[cnt] <= d ? id1[val[cnt]] = cnt : id2[j] = cnt;
		cntp[cnt] = val[cnt] - 1, s[cnt] = Sum2(val[cnt]) - 1;
	}
	for (i = 1; i <= tot && pr[i] <= n / pr[i]; ++i)
		for (j = 1; j <= cnt && pr[i] <= val[j] / pr[i]; ++j) {
			id = ID(val[j] / pr[i]), cntp[j] -= cntp[id] - i + 1;
			Inc2(s[j], Mul2(pr[i], Dec2(sp[i - 1], s[id])));
		}
	for (--i; i; --i)
		for (j = 1; j <= cnt && pr[i] <= val[j] / pr[i]; ++j) {
			id = ID(val[j] / pr[i]);
			v = Dec2(cntp[id] - i, g[id]), Inc2(g[j], v);
			Inc2(v, cntp[id] - i), Inc2(f[j], Dec2(v, f[id]));
		}
	for (i = 1; i <= cnt; ++i) Inc2(f[i], phi - cntp[i]);
	ans1 = Pow(2, f[ID(n)]), ans2 = 1;
	for (i = 1; i <= tot && pr[i] <= n / pr[i]; ++i)
		for (j = 1, v = pr[i]; v <= n; v = v * pr[i], ++j) ans2 = Mul1(ans2, Pow(j + 1, GetSum(n, v, v * pr[i])));
	for (lst = sp[i - 1], ret = 0, i = pr[i]; i <= n; i = j + 1) {
		j = n / (n / i);
		Inc2(ret, Mul2(Dec2(cur = s[ID(j)], lst), Sum2(n / i))), lst = cur;
	}
	ans2 = Mul1(ans2, Pow(2, ret)), printf("%lld\n", Mul1(ans1, ans2));
}

int main() {
	register int i, j;
	ispr[1] = 1;
	for (i = 1; i < maxn; ++i) {
		if (!ispr[i]) pr[++tot] = i, sp[tot] = (sp[tot - 1] + i) % phi;
		for (j = 1; j <= tot && i * pr[j] < maxn; ++j) {
			ispr[i * pr[j]] = 1;
			if (!(i % pr[j])) break;
		}
	}
	scanf("%d", &test);
	while (test) scanf("%lld", &n), Solve(), --test;
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值