[SPOJ]DIVCNTK – Counting Divisors (general)

SPOJ

luogu Remote Judge

latex \sum\limits_{i=1}^{n}\sigma(i^k)\pmod{2^{64}},n,k\le 10^{10}

latex \sigma(i) 表示 latex i 的约数个数。

显然 latex \sigma(x) 是积性函数。
f(1)=1\\ f(p)=k+1\\ f(p^c)=kc+1

latex g(i)=i,求出 latex g(n,|P|),乘以 latex (k+1) 就是 latex \sum\limits_{i=1}^{n}[i\in P]F(i)

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

const ull Max = 233333;
ull n, m, k, sn, w[Max], g[Max];
ull cntp, p[Max], id1[Max], id2[Max];
bool notp[Max];

ull S(ull x, ull y) {
	if (x < p[y] || x <= 1) return 0;
	ull id = x <= sn ? id1[x] : id2[n / x];
	ull re = (g[id] * (k + 1) - ((y - 1) * (k + 1)));
	for (ull i = y; i <= cntp && p[i] * p[i] <= x; i++) {
		ull t = p[i], t2 = p[i] * p[i];
		for (ull j = 1; t2 <= x; j++, t = t2, t2 *= p[i]) {
			re += ((k * j + 1) * S(x / t, i + 1) + (k * (j + 1) + 1));
		}
	}
	return re;
}

void solve() {
	scanf("%lld", &n), sn = sqrt(n);
	k = 1;
	cntp = m = 0;
	for (ull i = 2; i <= sn; i++) {
		if (!notp[i]) {
			p[++cntp] = i;
		}
		for (ull j = 1; j <= cntp && i * p[j] <= sn; j++) {
			notp[i * p[j]] = 1;
			if (i % p[j] == 0) break;
		}
	}
	for (ull i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i);
		w[++m] = n / i;
		g[m] = (w[m] - 1);
		if (w[m] <= sn) id1[w[m]] = m;
		else id2[n / w[m]] = m;
	}
	for (ull i = 1; i <= cntp; i++) {
		for (ull j = 1; j <= m && p[i] * p[i] <= w[j]; j++) {
			ull d = w[j] / p[i];
			ull id = d <= sn ? id1[d] : id2[n / d];
			g[j] -= g[id] - i + 1;
		}
	}
	printf("%llu\n", S(n, 1) + 1);
}

int main() {
	int cas = 0;
	scanf("%d", &cas);
	while (cas--) solve();
	return 0;
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注