[LOJ572]Misaka Network 与求和

题目

传送门 to LOJ

思路

答案为
$$\sum _{i=1}^{n}f(i{)}^k\left[2\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\phi (j)-1\right]$$

记 $h(x)=f(x{)}^k$ 则问题在于求出 $h$ 和 $h\ast \phi$ 的前缀和。

由于 $h$ 不是积性函数，所以传统的方法（如生成函数）都并不是很奏效。

用杜教筛求出 $\phi$ 的前缀和，则只需求出 $h$ 的前缀和。

本题要我们从新的角度考察 $min25$ 筛：将每个数质因数分解，只要能 把最大的质因数统一计算，复杂度就是正确的。

这当然可以做到，因为最大的质因数的贡献都是 $1$ 。

在 $min25$ 筛的 $\text{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;
}