答案为
∑
i
=
1
n
f
(
i
)
k
[
2
∑
j
=
1
⌊
n
i
⌋
φ
(
j
)
−
1
]
\sum_{i=1}^{n}f(i)^k\left[2\sum_{j=1}^{\lfloor{n\over i}\rfloor}\varphi(j) -1\right]
i=1∑nf(i)k⎣
⎡2j=1∑⌊in⌋φ(j)−1⎦
⎤
记 h ( x ) = f ( x ) k h(x)=f(x)^k h(x)=f(x)k 则问题在于求出 h h h 和 h ∗ φ h\ast\varphi h∗φ 的前缀和。
由于 h h h 不是积性函数,所以传统的方法(如生成函数)都并不是很奏效。
用杜教筛求出 φ \varphi φ 的前缀和,则只需求出 h h h 的前缀和。
本题要我们从新的角度考察 m i n 25 \tt min25 min25 筛:将每个数质因数分解,只要能 把最大的质因数统一计算,复杂度就是正确的。
这当然可以做到,因为最大的质因数的贡献都是 1 1 1 。
在 m i n 25 \tt min25 min25 筛的 step two \text{step two} step two 中,分别返回 “已经确定次大质因数” 和 “单个大质数” 两种情况的结果,容易转移。
哦,对了,这玩意儿是不是该叫洲阁筛了?
人尽皆知 { ⌊ n x ⌋ : x ∈ [ 1 , n ] } = { 1 , 2 , … , m , ⌊ n m ⌋ , ⌊ n m − 1 ⌋ , … , ⌊ n 1 ⌋ } \{\lfloor{n\over x}\rfloor:x\in[1,n]\}=\{1,2,\dots,m,\lfloor{n\over m}\rfloor,\lfloor{n\over m-1}\rfloor,\dots,\lfloor{n\over 1}\rfloor\} {⌊xn⌋:x∈[1,n]}={1,2,…,m,⌊mn⌋,⌊m−1n⌋,…,⌊1n⌋} 。
#include <cstdio>
#include <algorithm> // Almighty XJX yyds!!
#include <cstring> // Who can tell me why I'm so weak!
#include <cctype> // rainybunny root of the evil.
using llong = long long;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
# define rep0(i,a,b) for(int i=(a); i!=(b); ++i)
inline int readint(){
int a = 0, c = getchar(), f = 1;
for(; !isdigit(c); c=getchar()) if(c == '-') f = -f;
for(; isdigit(c); c=getchar()) a = a*10+(c^48);
return a*f;
}
const int MAXN = 2000000000, PRELEN = 1587401;
int primes[PRELEN], primes_size;
unsigned phi[PRELEN+1];
bool is_prime[PRELEN+1];
void sieve(int n = PRELEN){
memset(is_prime+2, true, n-1);
phi[1] = 1; // as the definition
for(int i=2,&len=primes_size; i<=n; ++i){
if(is_prime[i]) primes[++len] = i, phi[i] = i-1;
for(int j=1; j<=len; ++j){
const llong to = llong(i)*primes[j];
if(to > n) break;
is_prime[to] = false;
if(!(i%primes[j])){
phi[to] = phi[i]*primes[j]; break;
}
else phi[to] = phi[i]*(primes[j]-1);
}
}
}
const int SQRTN = 44721;
double inv[SQRTN+1]; /** Barret Reduction */
unsigned coe[SQRTN]; /** pre-computed power of primes */
int w[SQRTN+1]; /** only for big ones */
unsigned cntp[2][SQRTN+1]; /** 0 for big, 1 for small */
unsigned sumf[2][SQRTN+1];
int min25(int n){
int m = 1; /** floor( sqrt( @p n ) ) */
while(llong(m+1)*(m+1) <= n) ++ m;
rep(i,1,m) w[i] = int(n*inv[i]), cntp[0][i] = w[i]-1;
rep(i,1,m) cntp[1][i] = i-1; // small
for(int j=1; j<=primes_size; ++j){
const llong bound = llong(primes[j])*primes[j];
if(bound > n) break; // non-sense
rep(i,1,m){ // from big to small
if(w[i] < bound) break;
int todiv = i*primes[j]; unsigned* from;
if(todiv <= m) from = cntp[0]+todiv;
else from = cntp[1]+int(w[i]*inv[primes[j]]);
cntp[0][i] -= (*from)-unsigned(j-1);
}
drep(i,m,bound) cntp[1][i] -= cntp[1]
[int(i*inv[primes[j]])]-unsigned(j-1);
}
drep(j,primes_size,1){
const llong bound = llong(primes[j])*primes[j];
if(bound > n) continue; // not essential
rep(i,1,m){ // from big to small
if(w[i] < bound) break;
int todiv = i, tow = int(w[i]*inv[primes[j]]);
for(todiv*=primes[j]; todiv<=m; todiv*=primes[j]){
tow = int(tow*inv[primes[j]]); // keep moving
sumf[0][i] += sumf[0][todiv]+(cntp[0][todiv]-j+1)*coe[j];
}
for(; tow>=primes[j]; tow=int(tow*inv[primes[j]]))
sumf[0][i] += sumf[1][tow]+(cntp[1][tow]-j+1)*coe[j];
}
drep(i,m,bound) for(int tow=i; true; ){
tow = int(tow*inv[primes[j]]);
if(tow < primes[j]) break; // a bit too long
sumf[1][i] += sumf[1][tow]+(cntp[1][tow]-j+1)*coe[j];
}
}
return m;
}
unsigned sumphi[SQRTN];
unsigned qkpow(unsigned b, int q){
unsigned a = 1;
for(; q; q>>=1,b*=b) if(q&1) a *= b;
return a;
}
int main(){
int n = readint(), k = readint();
sieve(); rep(i,1,PRELEN) phi[i] += phi[i-1];
rep(i,1,SQRTN) inv[i] = double(1+1e-13)/i;
const int bound = n/(PRELEN+1);
for(int i=bound,m=1; i; --i){
const int pos = int(n*inv[i]);
while(llong(m+1)*(m+1) <= pos) ++ m;
sumphi[i] = unsigned(llong(pos+1)*pos>>1)+m*phi[m];
rep(j,1,m) sumphi[i] -= (phi[j]-phi[j-1])*unsigned(pos*inv[j]);
int j = 2, todiv = i<<1;
for(; todiv<=bound; ++j,todiv+=i) sumphi[i] -= sumphi[todiv];
for(; j<=m; ++j) sumphi[i] -= phi[int(pos*inv[j])];
}
for(int i=1; primes[i]<=SQRTN; ++i)
coe[i] = qkpow(primes[i],k);
const int m = min25(n);
rep(o,0,1) rep(i,1,m) sumf[o][i] += cntp[o][i];
unsigned ans = -phi[m]*sumf[1][m];
rep(i,1,bound) ans += sumphi[i]*(sumf[1][i]-sumf[1][i-1]);
rep(i,bound+1,m) ans += (sumf[1][i]-sumf[1][i-1])*phi[int(n*inv[i])];
rep(i,1,m) ans += (phi[i]-phi[i-1])*sumf[0][i];
ans = (ans<<1)-sumf[0][1];
printf("%u\n",ans);
return 0;
}