【51Nod1847】奇怪的数学题

【题目链接】

• 点击打开链接

【思路要点】

• 令 $f(i)=(\frac{i}{Min(i)}{)}^k(i>1)$ ，即 $f(i)$ 表示 $i$ 次大的因子的 $k$ 次方，特别规定 $f(1)=0$ 。

• 那么原式即为 ${\sum }_{i=1}^N{\sum }_{j=1}^{N}f(gcd(i,j))$ 。

• 考虑枚举 $gcd(i,j)$ ，那么原式

  $={\sum }_{g=1}^{N}f(g){\sum }_{i=1}^{\lfloor \frac{N}{g}\rfloor }{\sum }_{j=1}^{\lfloor \frac{N}{g}\rfloor }ϵ(gcd(i,j))$

  $={\sum }_{g=1}^{N}f(g){\sum }_{i=1}^{\lfloor \frac{N}{g}\rfloor }{\sum }_{j=1}^{\lfloor \frac{N}{g}\rfloor }{\sum }_{d/i,d/j}\mu (d)$

  $={\sum }_{g=1}^{N}f(g){\sum }_{d=1}^{\lfloor \frac{N}{g}\rfloor }\mu (d)\lfloor \frac{N}{gd}{\rfloor }^2$

  $={\sum }_{gd=1}^N(f\ast \mu )(gd)\lfloor \frac{N}{gd}{\rfloor }^2$

• 对 $\lfloor \frac{N}{gd}\rfloor$ 整数分块，我们剩余的工作是计算函数 $(f\ast \mu )$ 在每一个 $\lfloor \frac{N}{gd}\rfloor$ 处的前缀和。

• 令 $g(N)={\sum }_{i=1}^N(f\ast \mu )(i)$ ，注意到 $(f\ast \mu )\ast 1=f\ast (\mu \ast 1)=f$ ，因此

  ${\sum }_{i=1}^{N}f(i)={\sum }_{i=1}^N(f\ast \mu \ast 1)(i)={\sum }_{i=1}^{N}g(\lfloor \frac{N}{i}\rfloor )$

  $g(N)={\sum }_{i=1}^{N}f(i)-{\sum }_{i=2}^{N}g(\lfloor \frac{N}{i}\rfloor )$

• 上面的式子显然可以杜教筛，剩余的唯一问题在于计算函数 $f$ 在每一个 $\lfloor \frac{N}{i}\rfloor$ 处的前缀和。

• 分开考虑质数和合数对 ${\sum }_{i=1}^{N}f(i)$ 的影响，有

  ${\sum }_{i=1}^{N}f(i)=cnt(N)+{\sum }_{pisaprime,p\le \sqrt{N}}sum(\lfloor \frac{N}{p}\rfloor ,p)$

  其中 $cnt(N)$ 表示 $N$ 以内质数的个数，

  $sum(N,p)$ 表示 $N$ 以内最小质因子大于等于 $p$ 的合数的 $k$ 次方和。

  用 $Min25$ 筛解决即可。

• 还有一个小问题，由于模数非质数，求解自然数的 $k$ 次方和是需要用第二类斯特林数转化一下。

  ${\sum }_{i=1}^N{i}^k={\sum }_{i=1}^N{\sum }_{j=1}^{k}S(k,j)\ast i^{\underline{j}}$

  $={\sum }_{j=1}^{k}S(k,j){\sum }_{i=1}^N{i}^{\underline{j}}$

  $={\sum }_{j=1}^{k}S(k,j)\ast j!{\sum }_{i=1}^N\left(\genfrac{}{}{0ex}{}{i}{j}\right)$

  $={\sum }_{j=1}^{k}S(k,j)\ast j!\left(\genfrac{}{}{0ex}{}{N+1}{j+1}\right)$

  $={\sum }_{j=1}^{k}S(k,j)\frac{(N+1{)}^{\underline{j+1}}}{j+1}$

• 时间复杂度 $O(N^{\frac{3}{4}}+K^2\sqrt{N})$ 。

【代码】

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 100005;
const int MAXK = 55;
template <typename T> void chkmax(T &x, T y) {x = max(x, y); }
template <typename T> void chkmin(T &x, T y) {x = min(x, y); } 
template <typename T> void read(T &x) {
	x = 0; int f = 1;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
	for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
	x *= f;
}
template <typename T> void write(T x) {
	if (x < 0) x = -x, putchar('-');
	if (x > 9) write(x / 10);
	putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {
	write(x);
	puts("");
}
unsigned n, k, tot, limit, pre[MAXN], prime[MAXN], f[MAXN], primek[MAXN];
unsigned m, val[MAXN], home1[MAXN], home2[MAXN], cnt[MAXN], sum[MAXN];
unsigned pref[MAXN], s[MAXK][MAXK], ans[MAXN];
unsigned power(unsigned x, unsigned y) {
	if (y == 0) return 1;
	unsigned tmp = power(x, y / 2);
	if (y % 2 == 0) return tmp * tmp;
	else return tmp * tmp * x;
}
unsigned calc(unsigned n) {
	unsigned ans = 0;
	for (unsigned i = 1; i <= k && i <= n; i++) {
		unsigned now = 1;
		for (unsigned j = n + 1; j >= n - i + 1; j--)
			if (j % (i + 1) == 0) now *= j / (i + 1);
			else now *= j;
		ans += now * s[k][i];
	}
	return ans;
}
void init(unsigned n) {
	s[0][0] = 1;
	for (unsigned i = 1; i < MAXK; i++)
	for (unsigned j = 1; j <= i; j++)
		s[i][j] = s[i - 1][j - 1] + s[i - 1][j] * j;
	for (unsigned i = 2; i <= n; i++) {
		if (f[i] == 0) {
			prime[++tot] = f[i] = i;
			primek[tot] = power(i, k);
			pre[tot] = pre[tot - 1] + primek[tot];
		}
		for (unsigned j = 1; j <= tot && prime[j] <= f[i]; j++) {
			unsigned tmp = prime[j] * i;
			if (tmp > n) break;
			f[tmp] = prime[j];
		}
	}
}
int main() {
	read(n), read(k);
	limit = sqrt(n);
	init(limit);
	for (unsigned i = 1, nxt; i <= n; i = nxt) {
		unsigned tmp = n / i;
		val[++m] = tmp;
		cnt[m] = tmp - 1;
		sum[m] = calc(tmp) - 1;
		nxt = n / tmp + 1;
		if (tmp <= limit) home1[tmp] = m;
		else home2[n / tmp] = m;
	}
	for (unsigned i = 1; i <= tot; i++) {
		for (unsigned j = 1; prime[i] * prime[i] <= val[j]; j++) {
			unsigned tmp = val[j] / prime[i], pos;
			if (tmp <= limit) pos = home1[tmp];
			else pos = home2[n / tmp];
			cnt[j] -= cnt[pos] - (i - 1);
			pref[j] += sum[pos] - pre[i - 1];
			sum[j] -= primek[i] * (sum[pos] - pre[i - 1]);
		}
	}
	for (unsigned i = 1; i <= m; i++)
		pref[i] += cnt[i];
	unsigned finalans = 0;
	for (unsigned i = m; i >= 1; i--) {
		ans[i] = pref[i];
		for (unsigned j = 2, nxt; j <= val[i]; j = nxt) {
			unsigned tmp = val[i] / j, pos;
			if (tmp <= limit) pos = home1[tmp];
			else pos = home2[n / tmp];
			nxt = val[i] / tmp + 1;
			ans[i] -= ans[pos] * (nxt - j);
		}
		finalans += (ans[i] - ans[i + 1]) * (n / val[i]) * (n / val[i]);
	}
	writeln(finalans);
	return 0;
}