Project Euler 501. Eight Divisors


https://projecteuler.net/problem=501


如何一次篩法求出所有的解?拆成兩個數組,大小都小於 sqrt(n):

  • F[i]: 小於等於 i 的 PrimePi。
  • G[i]: 小於等於 n/i 的 PrimePi。

//}/* .................................................................................................................................. */

const int PMAX = 1000000 + 1;
vector<LL> P; bitset<PMAX> isC;
LL n, nn; int PP[PMAX];

bool isPrime(LL x){
    ECH(it, P){
        if (sqr(*it) > x) break;
        if (x % *it == 0) return false;
    }
    return true;
}

#define ii (i*P[j])
#include <unordered_map>
unordered_map<LL, LL> F, G;

void sieve(){
    
    nn = sqrt(n); FOR(i, 2, nn+1){
        if (!isC[i]) P.PB(i);
        for (int j=0;j<SZ(P)&&ii<nn;++j){
            isC[ii]=1; if (!(i%P[j])) break;
        }
    }
    
    REP_1(i, nn+1){
        F[i] = i-1;
        G[i] = n/i - 1;
    }
    
    FOR_1(p, 2, nn+1) if (F[p-1] != F[p]){
        LL pp = (LL) p * p;
        LL up = min(nn, n / pp);
        REP_1(q, up){
            LL d = (LL)p * q;
            if (d <= nn) G[q] -= G[d] - F[p-1];
            else G[q] -= F[n/d] - F[p-1];
        }
        DWN_1(i, nn+1, pp-1) F[i] -= F[i/p] - F[p-1];
    }
}

LL 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 pp = (LL)p*p; for(auto v: V){
            if (v < pp) break;
            dp[v] -= (dp[v/p] - dp[p-1]);
        }
    }
    return dp[n];
}


LL z = 0;

int main(){
    
#ifndef ONLINE_JUDGE
    freopen("/users/xiaodao/desktop/Exercise/in.txt", "r", stdin);
    //   freopen("out.txt", "w", stdout);
#endif
    
    
    n = 1e12; sieve(); z = 0;
    
    ECH(it, P){
        if (pow(*it, 7) > n) break;
        z += 1;
    }
    
    REP(i, P.size()){
    
        LL t = n / cub(P[i]);  if (!t) break;
        z += t <= nn ? F[t] : G[cub(P[i])];
        if (t >= P[i]) z -= 1;
        
        FOR(j, i+1, P.size()){
            LL x = P[i] * P[j];
            if (n/x > P[j]){
                z += n/x <= nn ? F[n/x] : G[x];
                z -= j+1;
            }
            else break;
        }
    }
    
    cout << z << endl;
}
// 197912312715