Min25筛学习 + 【51nod1847】奇怪的数学题（Min_25筛+杜教筛）
Min25筛小结——alpc_qleonardo的博客
讲的非常清楚，不过其中大\(S(n,j)\)表示的应该是从\(1\)累加到\(n\)的\(F(i)\)，\(i\)要么是质数，要么最小质因子大于等于\(j\)。这样才满足那个递推式。
然后来看一道巧妙的例题。
题意
给出 \(N,K\) ,请计算下面这个式子:
\(\sum_{i=1}^N\sum_{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}n\sum_{j=1}nf(gcd(i,j))k\=\sum_{d=1}nf(d)k\sum_{i=1}\sum_{j=1}^{\lfloor \frac nd\rfloor}[(i,j)==1]\=\sum_{d=1}nf(d)k(2\sum_{i =1}^{\lfloor \frac nd\rfloor}\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\)次方之和。中间的转移有这么一步：
答案就是\(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\n", Ans);
}
