@description@
给出 N, K, 请计算下面这个式子:
[sum_{i=1}^{N}sum_{j=1}^{N}sgcd(i, j)^k
]
其中, sgcd(i, j)表示(i, j)的所有公约数中第二大的,特殊地,如果gcd(i, j) = 1, 那么sgcd(i, j) = 0。
考虑到答案太大,请输出答案对2^32取模的结果.
@solution@
记 mnp[i] 表示 i 的最小质因子。则：
[egin{aligned}
sum_{i=1}^{N}sum_{j=1}^{N}sgcd(i, j)^k =& sum_{d=1}^{N}sum_{i=1}^{N}sum_{j=1}^{N}[gcd(i, j) = d] imes (frac{d}{mnp[d]})^{k}\
=& sum_{d=1}^{N}(frac{d}{mnp[d]})^{k}sum_{i=1}^{lfloor frac{N}{d}
floor}sum_{j=1}^{lfloor frac{N}{d}
floor}[gcd(i, j) = 1]\
=& sum_{d=1}^{N}(frac{d}{mnp[d]})^{k}(sum_{i=1}^{lfloor frac{N}{d}
floor}phi(i) - 1)
end{aligned}]
后面那个 (phi) 前缀和，杜教筛基础，不讲了。前面整除分块，接下来讲讲具体怎么求。
记 (S(N) = sum_{d=1}^{N}(frac{d}{mnp[d]})^{k})，考虑怎么求 (S)。分成两部分：d 为质数，贡献为 (pi(N))（N 以内质数个数）；d 为合数。
注意到这种涉及到最小质因子的可以类比 min-25 筛的思想搞（其实就是埃氏筛法）。
我们做 min-25 筛前半部分时，有一个辅助数组 (f(N, p)) 表示 N 以内的 质数 或 最小质因子 > 第 p 个质数的数 的 K 次幂之和。
利用 (f(N, p)) 就可以得到 N 以内最小质因子 >= 第 p 个质数的数的 K 次幂之和 (g(N, p))，这一点很容易做到。
然后有：
[S(N) = pi(N) + sum_{i=1}^{pcnt}g(lfloor frac{N}{prime_i}
floor, i)
]
自然数幂和用斯特林数处理好些，因为模数不是质数，不一定有逆元。
@accepted code@
#include <cstdio>
#include <algorithm>
using namespace std;
const unsigned long long MOD = (1LL << 32);
typedef unsigned int ui;
const int MAXN = 1000000;
ui S[55][55], inv[55];
void init() {
for(int i=0;i<=50;i++) {
S[i][0] = (i == 0);
for(int j=1;j<=i;j++)
S[i][j] = S[i-1][j-1] + S[i-1][j] * ui(j);
}
inv[1] = 1;
for(int i=3;i<=51;i+=2) {
for(int j=1;j<=i;j++)
if( (j*MOD + 1) % i == 0 ) {
inv[i] = (j*MOD + 1) / i % MOD;
break;
}
}
}
ui sum(ui n, int k) {
int cnt = 0; ui ans = 0, tmp = 1; n++;
for(int i=0;i<=k;i++) {
ui del = (n - i), pw = 0;
if( del == 0 ) break;
while( del % 2 == 0 ) del /= 2, cnt++;
tmp *= del;
del = i + 1;
while( del % 2 == 0 ) del /= 2, pw++;
if( (cnt - pw) < 60 ) {
ui p = ui(1LL << (cnt - pw));
ans += S[k][i]*tmp*inv[del]*p;
}
}
return ans;
}
int prm[MAXN + 5], pcnt;
bool nprm[MAXN + 5]; ui phi[MAXN + 5];
void sieve() {
phi[1] = 1;
for(int i=2;i<=MAXN;i++) {
if( !nprm[i] ) prm[++pcnt] = i, phi[i] = i - 1;
for(int j=1;i*prm[j]<=MAXN;j++) {
nprm[i*prm[j]] = true;
if( i % prm[j] == 0 ) {
phi[i*prm[j]] = phi[i]*prm[j];
break;
}
else phi[i*prm[j]] = phi[i]*phi[prm[j]];
}
}
for(int i=1;i<=MAXN;i++)
phi[i] += phi[i - 1];
}
int n, k;
ui powk(ui b) {
ui ret = 1;
for(int i=k;i;i>>=1,b*=b)
if( i & 1 ) ret *= b;
return ret;
}
int a[MAXN + 5], id1[MAXN + 5], id2[MAXN + 5], cnt;
int id(int x) {return x <= MAXN ? id1[x] : id2[n / x];}
bool tag[MAXN + 5]; ui sp[MAXN + 5];
ui phisum(int n) {
if( n <= MAXN ) return sp[id(n)] = phi[n];
if( tag[id(n)] ) return sp[id(n)];
ui ans = sum(n, 1);
for(int i=2;i<=n;i++) {
int p = n / i, j = n / p;
ans -= phisum(p) * (sum(j, 0) - sum(i - 1, 0));
i = j;
}
tag[id(n)] = true; return sp[id(n)] = ans;
}
ui f[MAXN + 5], g[MAXN + 5], h[MAXN + 5];
void solve() {
for(int i=1;i<=cnt;i++)
f[i] = a[i] - 1, g[i] = sum(a[i], k) - 1, h[i] = 0;
ui tmp = 0;
for(int i=1;i<=pcnt;i++) {
long long p = 1LL*prm[i]*prm[i]; ui pw = powk(prm[i]);
if( p > n ) break;
for(int j=1;j<=cnt;j++) {
if( p > a[j] ) break;
h[j] += (g[id(a[j]/prm[i])] - tmp);
f[j] -= f[id(a[j]/prm[i])] - (i - 1);
g[j] -= (g[id(a[j]/prm[i])] - tmp)*pw;
}
tmp += pw;
}
for(int i=1;i<=cnt;i++) h[i] += f[i];
}
void get_id() {
for(int i=1;i<=n;i=n/(n/i)+1) {
int p = n / i;
if( p <= MAXN ) id1[p] = (++cnt);
else id2[n / p] = (++cnt);
a[cnt] = p;
}
}
int main() {
init();
sieve(), scanf("%d%d", &n, &k), get_id(), solve();
ui ans = 0;
for(int i=2;i<=n;i++) {
int p = n / i, j = n / p;
ans += (2*phisum(p) - 1)*(h[id(j)] - h[id(i - 1)]);
i = j;
}
printf("%lld
", (long long)ans);
}
@details@
妄想用线性求逆元然后发现这种方法并不适用于模数为合数的情况。当时我人就傻在那里。
然后不想打 exgcd，于是突然脑洞大开写了暴力解不定方程求逆元。