Min25筛学习 + 【51nod1847】奇怪的数学题（Min_25筛+杜教筛）

Min25筛小结——alpc_qleonardo的博客
讲的非常清楚，不过其中大(S(n,j))表示的应该是从(1)累加到(n)的(F(i))，(i)要么是质数，要么最小质因子大于等于(j)。这样才满足那个递推式。

然后来看一道巧妙的例题。

---

题意

给出 (N,K) ,请计算下面这个式子:
(sum_{i=1}^Nsum_{j=1}^Nsgcd(i,j)^k)
其中,(sgcd(i, j))表示((i, j))的所有公约数中第二大的,特殊地,如果(gcd(i, j) = 1), 那么(sgcd(i, j) = 0)。
考虑到答案太大,请输出答案对(2^{32})取模的结果.
(1≤N≤10^9,1≤K≤50)

题解

显然有(sgcd(i,j)=frac{gcd(i,j)}{minp(gcd(i,j))})，(minp(x))表示(x)的最小质因子。那么设(f(x)=frac x{minp(x)})，则$$Ans=sum_{i=1}^nsum_{j=1}^nf(gcd(i,j))^k=sum_{d=1}^nf(d)^ksum_{i=1}^{lfloor frac nd floor}sum_{j=1}^{lfloor frac nd floor}[(i,j)==1]=sum_{d=1}^nf(d)^k(2sum_{i =1}^{lfloor frac nd floor}varphi(i)-1)$$

后面的(varphi)直接整除分块+杜教筛，然后问题就是如何求前面的前缀和，及(sum_{i=1}^nf(i)^k)。

考虑(min25)筛的过程，把数分为质数，合数和(1)来计算。

(f(1)=0)，无需统计。
质数：(f(p)=1)，只需要统计质数个数就行了。
合数比较麻烦，见下：

设(g(n,j))表示在([1,n])中的所有的质数或者最小质因数大于(P_j)的合数的(k)次方之和。中间的转移有这么一步：

[g(n,j)=g(n,j-1)-P_j^kcdotleft[g(lfloorfrac n{P_j} floor,j-1)-sum_{i=1}^{j-1}P_i^k ight]$$考虑后面减去的部分，代表了所有最小质因子等于$P_j$的合数的$k$次方之和。那么中括号里的值就代表这些数的$k$次方之和除以$P_j^k$，这不就是这些数的$f(i)^k$之和吗？那么我们就通过计算$g$可以统计出所有合数的$f(i)^k$之和，就是$g(n,|P|)$ 质数的个数统计类比于$g$，转移式子为： $$f(n,j)=f(n,j-1)-left[f(lfloorfrac n{P_j} floor,j-1)-(j-1) ight]]

答案就是(f(n,|P|))

最终答案就加起来就行了。

求幂和要用第二类斯特林数。

详见代码。

CODE

#include <bits/stdc++.h>
using namespace std;
#define int unsigned
const int N = 100005;
int m, k, sqr, cnt, p[N], phi[N], pw[N], sp[N];
int s2[55][55], id[N], id2[N], ans[2][N], w[N], g[N];
bool vis[N];

inline int qpow(int a, int b) {
	int re = 1;
	while(b) {
		if(b&1) re *= a;
		a *= a; b >>= 1;
	}
	return re;
}

void init(int N) {
	phi[1] = 1;
	for(int i = 2; i <= N; ++i) {
		if(!vis[i]) p[++cnt] = i, sp[cnt] = sp[cnt-1] + (pw[cnt]=qpow(i, k)), phi[i] = i-1;
		for(int j = 1, k; j <= cnt && i * p[j] <= N; ++j) {
			vis[k = i * p[j]] = 1;
			if(i % p[j] == 0) {
				phi[k] = phi[i] * p[j];
				break;
			}
			phi[k] = phi[i] * (p[j]-1);
		}
	}
	for(int i = 1; i <= N; ++i) phi[i] += phi[i-1];
	s2[0][0] = 1;
	for(int i = 1; i <= k; ++i) {
		s2[i][1] = 1;
		for(int j = 2; j <= i; ++j)
			s2[i][j] = s2[i-1][j] * j + s2[i-1][j-1];
	}
}

map<int,int>mp;
int Phi(int n) { //杜教筛求phi和
	if(n <= sqr) return phi[n];
	if(mp.count(n)) return mp[n];
	int re = n*(n+1)>>1;
	for(int i = 2, j; i <= n; i = j+1)
		j = n/(n/i), re -= (j-i+1) * Phi(n/i);
	return mp[n] = re;
}

int calc(int n) { //第二类斯特林数求幂和
	int re = 0, tmp;
	for(int i = 1; i <= k && i <= n; ++i) {
		tmp = s2[k][i];
		for(int j = n-i+1; j <= n+1; ++j)
			j%(i+1) ? tmp *= j : tmp *= j/(i+1);
		re += tmp;
	}
	return re;
}

signed main () {
	int n;
	scanf("%u%u", &n, &k); init(sqr=sqrt(n));
	for(int i = 1, j; i <= n; i = j+1) {
		j = n/(n/i), w[++m] = n/i; //给每个n/i下去整的值标号
		w[m] <= sqr ? id[w[m]] = m : id2[n/w[m]] = m; //大于根号的无法直接开数组存，用n除一下再存
		g[m] = calc(w[m])-1, ans[0][m] = w[m]-1; //g初值就是2~w[m]所有的k次方之和，-1是为了除去1
												//ans[0]就是f,初值就是2~w[m]的数的个数
	}
	for(int j = 1, ID; j <= cnt; ++j)
		for(int i = 1; 1ll*p[j]*p[j] <= w[i]; ++i) { //w[i]是降序的
			ID = w[i]/p[j] <= sqr ? id[w[i]/p[j]] : id2[n/(w[i]/p[j])];
			g[i] -= pw[j] * (g[ID] - sp[j-1]); //g转移
			ans[0][i] -= ans[0][ID] - (j-1); //f转移
			ans[1][i] += g[ID] - sp[j-1]; //合数答案累计
		}
	int Ans = 0;
	for(int i = 2, j, lst=0, now, ID; i <= n; lst = now, i = j+1) {
		j = n/(n/i), ID = j <= sqr ? id[j] : id2[n/j];
		now = ans[0][ID] + ans[1][ID];
		Ans += ((Phi(n/i)<<1)-1) * (now - lst);
	}
	printf("%u
", Ans);
}