某岛

… : "…アッカリ~ン . .. . " .. .
July 1, 2015

Project Euler 521. Smallest prime factor


https://projecteuler.net/problem=521


拆解成两个子问题,

  • 小于等于 n 的素数之和。
  • 小于等于 n 的合数的最小素因子之和。

回忆 PE 10,我们可以用 $$O(n^{3/4})$$ 的复杂度去统计小于等于 n 的素数的 k 次方和。
那么前者是 k=1,后者可以用 k=0 稍微修改得到。(因此每次筛 p 过程中,dp[n] 的变化,就是所有最小素因子等于 p 的合数的个数。)


Int f(LL n){
    LL r = sqrt(n); vector<LL> V; REP_1(i, r+1) V.PB(n/i);
    int up = int(*V.rbegin()); DWN_1(i, up-1, 1) V.PB(i);
    
    Int z = 0;
    
    unordered_map<LL, LL> dp;
    
    REP(i, V.size()){
        LL x = V[i];
        dp[x] = x-1;
    }
    
    
    FOR_1(p, 2, r+1) if (dp[p] > dp[p-1]){
        LL prev = dp[n];
        LL pp = (LL)p*p; for(auto v: V){
            if (v < pp) break;
            dp[v] -= (dp[v/p] - dp[p-1]);
        }
        z += Int(p) * Int(prev - dp[n]);
    }
    return z;
}



Int g(LL n){
    LL r = sqrt(n); vector<LL> V; REP_1(i, r+1) V.PB(n/i);
    int up = int(*V.rbegin()); DWN_1(i, up-1, 1) V.PB(i);
    
    unordered_map<LL, Int> dp;

    REP(i, V.size()){
        LL x = V[i];
        if (x&1) dp[x] = Int(x) * Int((x+1)/2) - 1;
        else dp[x] = Int(x/2) * Int(x+1) - 1;
    }
    
    FOR_1(p, 2, r+1) if (dp[p] != dp[p-1]) {
        LL pp = (LL)p*p; for(auto v: V){
            if (v < pp) break;
            dp[v] -= Int(p)*(dp[v/p] - dp[p-1]);
        }
    }
    return dp[n];
}



int main(){
    
#ifndef ONLINE_JUDGE
    freopen("/users/xiaodao/desktop/Exercise/in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
#endif
    
    LL n = 1e12;
    cout <<  f(n)  + g(n) << endl;
    
    //44389811
}