分析：Meisell-Lehmer算法+二分蜜汁RE，只能用这个不知其算法名的官方题解。
Code：
#include <cstdio> #include <cmath> #define MAXN 100 #define MAXM 10001 #define MAXP 40000 #define MAX 400000 #define LL long long #define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31)))) #define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31)))) #define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2)) using namespace std; namespace pcf{ LL dp[MAXN][MAXM]; unsigned int ar[(MAX >> 6) + 5] = {0}; int len = 0, primes[MAXP], counter[MAX]; void Sieve() { setbit(ar, 0), setbit(ar, 1); for(int i = 3; (i*i) < MAX; i++, i++) { if(!chkbit(ar, i)) { int k = i << 1; for(int j = (i*i); j < MAX; j += k) setbit(ar, j); } } for (int i = 1; i < MAX; i++) { counter[i] = counter[i - 1]; if(isprime(i)) primes[len++] = i, counter[i]++; } } void init() { Sieve(); for(int n = 0; n < MAXN; n++) { for(int m = 0; m < MAXM; m++) { if(!n) dp[n][m] = m; else dp[n][m] = dp[n-1][m] - dp[n-1][m/primes[n-1]]; } } } LL phi(LL m, int n) { if (n == 0) return m; if (primes[n-1] >= m) return 1; if (m < MAXM && n < MAXN) return dp[n][m]; return phi(m, n-1) - phi(m/primes[n-1], n-1); } LL Lehmer(LL m) { if (m < MAX) return counter[m]; LL w, res = 0; int i, a, s, c, x, y; s = sqrt(0.9 + m), y = c = cbrt(0.9 + m); a = counter[y], res = phi(m, a) + a - 1; for (i = a; primes[i] <= s; i++) res = res - Lehmer(m/primes[i]) + Lehmer(primes[i]) - 1; return res; } } LL solve(LL n) { int i, j, k, l; LL x, y, res = 0; for (i = 0; i < pcf::len; i++) { x = pcf::primes[i], y = n/x; if ((x*x) > n) break; res += (pcf::Lehmer(y) - pcf::Lehmer(x)); } for (i = 0; i < pcf::len; i++) { x = pcf::primes[i]; if ((x*x*x) > n) break; res++; } return res; } int main() { int n, l, r, mid, count = 0; pcf::init(); while(scanf("%d", &n) && n) { l = 1, r = 1e9; while(l <= r) { mid = (l+r)/2; if(pcf::Lehmer(mid) < n) l = mid+1; else r = mid-1; } printf("Case %d: %d\n", ++count, l); } return 0; }
51Nod的1184可通过Meisell-Lehmer算法+二分AC。
#include <cstdio> #include <cmath> #define LL long long using namespace std; const int N = 5e6 + 5; const int M = 7; const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17; int np[N]; int prime[N], pi[N]; int phi[PM+1][M+1], sz[M+1]; int getprime() { int cnt = 0; np[0] = np[1] = 1; pi[0] = pi[1] = 0; for(int i = 2; i < N; ++i) { if(!np[i]) prime[++cnt] = i; pi[i] = cnt; for(int j = 1; j <= cnt && i * prime[j] < N; ++j) { np[i * prime[j]] = 1; if(i % prime[j] == 0) break; } } return cnt; } void init() { getprime(); sz[0] = 1; for(int i = 0; i <= PM; ++i) phi[i][0] = i; for(int i = 1; i <= M; ++i) { sz[i] = prime[i] * sz[i-1]; for(int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i-1] - phi[j/prime[i]][i-1]; } } int sqrt2(LL x) { LL r = (LL)sqrt(x - 0.1); while(r*r <= x) ++r; return (int)(r-1); } int sqrt3(LL x) { LL r = (LL)cbrt(x - 0.1); //cbrt(x): x的立方根 while(r*r*r <= x) ++r; return (int)(r-1); } LL getphi(LL x, int s) { if(s == 0) return x; if(s <= M) return phi[x%sz[s]][s] + (x/sz[s]) * phi[sz[s]][s]; if(x <= prime[s]*prime[s]) return pi[x] - s + 1; if(x <= prime[s]*prime[s]*prime[s] && x < N) { int s2x = pi[sqrt2(x)]; LL ans = pi[x] - (s2x+s-2)*(s2x-s+1)/2; for(int i = s+1; i <= s2x; ++i) ans += pi[x/prime[i]]; return ans; } return getphi(x, s-1) - getphi(x/prime[s], s-1); } LL getpi(LL x) { if(x < N) return pi[x]; LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1; for(int i = pi[sqrt3(x)]+1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x/prime[i]) - i + 1; return ans; } LL lehmer_pi(LL x) { if(x < N) return pi[x]; int a = (int)lehmer_pi(sqrt2(sqrt2(x))); int b = (int)lehmer_pi(sqrt2(x)); int c = (int)lehmer_pi(sqrt3(x)); LL sum = getphi(x, a) + (LL)(b+a-2) * (b-a+1)/2; for (int i = a+1; i <= b; ++i) { LL w = x/prime[i]; sum -= lehmer_pi(w); if(i > c) continue; LL lim = lehmer_pi(sqrt2(w)); for(int j = i; j <= lim; ++j) sum -= lehmer_pi(w/prime[j]) - (j-1); } return sum; } int main() { LL n, l, r, mid; init(); while(~scanf("%lld",&n)) { l = 1, r = 1e11; while(l <= r) { mid = (l+r)/2; if(lehmer_pi(mid) < n) l = mid+1; else r = mid-1; } printf("%lld\n", l); } return 0; }
继续加油~
