题目链接：
https://projecteuler.net/problem=501
题意：
\(f(n)\) be the count of numbers not exceeding \(n\) with exactly eight divisors.
就是找少于等于 \(n\)中只有8个因子的个数。
You are given \(f(100) = 10, f(1000) = 180\) and \(f(10^6) = 224427\).
Find \(f(10^{12})\)
题解:
我们知道，对于 \(n=\prod p_i^{a_i}\) ，那么，\(n\)的因子的个数有 \(\prod (a_i+1)\)个。
那么，符合题目条件的只有三种情况。
\(1.p^{7} <= n\)
\(2.p^{3} * q <= n\)
\(3.p * q * r <= n\)'
其中，\(p,q,r\)是各自不相等的质数，并且 \(p < q < r\)。
和这题套路一样。
http://codeforces.com/problemset/problem/665/F
Codefores这题代码：
#include <bits/stdc++.h>
using namespace std;
const int MAX = 2e6+5;
const int M = 7;
typedef long long ll;
vector<int> lp, primes, pi;
int phi[MAX+1][M+1], sz[M+1];
void factor_sieve()
{
lp.resize(MAX);
pi.resize(MAX);
lp[1] = 1;
pi[0] = pi[1] = 0;
for (int i = 2; i < MAX; i++) {
if (lp[i] == 0) {
lp[i] = i;
primes.emplace_back(i);
}
for (int j = 0; j < primes.size() && primes[j] <= lp[i]; j++) {
int x = i * primes[j];
if (x >= MAX) break;
lp[x] = primes[j];
}
pi[i] = primes.size();
}
}
void init()
{
factor_sieve();
for(int i = 0; i <= MAX; i++) {
phi[i][0] = i;
}
sz[0] = 1;
for(int i = 1; i <= M; i++) {
sz[i] = primes[i-1]*sz[i-1];
for(int j = 1; j <= MAX; j++) {
phi[j][i] = phi[j][i-1] - phi[j/primes[i-1]][i-1];
}
}
}
int sqrt2(long long x)
{
long long r = sqrt(x - 0.1);
while (r*r <= x) ++r;
return r - 1;
}
int cbrt3(long long x)
{
long long r = cbrt(x - 0.1);
while(r*r*r <= x) ++r;
return r - 1;
}
long long getphi(long long 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 <= primes[s-1]*primes[s-1])
{
return pi[x] - s + 1;
}
if(x <= primes[s-1]*primes[s-1]*primes[s-1] && x < MAX)
{
int sx = pi[sqrt2(x)];
long long ans = pi[x] - (sx+s-2)*(sx-s+1)/2;
for(int i = s+1; i <= sx; ++i) {
ans += pi[x/primes[i-1]];
}
return ans;
}
return getphi(x, s-1) - getphi(x/primes[s-1], s-1);
}
long long getpi(long long x)
{
if(x < MAX) return pi[x];
int cx = cbrt3(x), sx = sqrt2(x);
long long ans = getphi(x, pi[cx]) + pi[cx] - 1;
for(int i = pi[cx]+1, ed = pi[sx]; i <= ed; i++)
{
ans -= getpi(x/primes[i-1-1]) - i + 1;
}
return ans;
}
long long lehmer_pi(long long x)
{
if(x < MAX) return pi[x];
int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
int b = (int)lehmer_pi(sqrt2(x));
int c = (int)lehmer_pi(cbrt3(x));
long long sum = getphi(x, a) + (long long)(b + a - 2) * (b - a + 1) / 2;
for (int i = a + 1; i <= b; i++)
{
long long w = x / primes[i-1];
sum -= lehmer_pi(w);
if (i > c) continue;
long long lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++)
{
sum -= lehmer_pi(w / primes[j-1]) - (j - 1);
}
}
return sum;
}
long long power(long long a, long long b)
{
long long x = 1, y = a;
while(b)
{
if (b&1) x = x * y;
y = y * y;
b >>= 1;
}
return x;
}
void solve(long long n)
{
ll ans = 0;
// case 1: p^3 <= n ,p is a prime
for(int i = 0; i < (int)primes.size(); i++) {
if (power(primes[i], 3) <= n) {
//std::cout << primes[i] << '\n';
ans += 1;
}
else {
break;
}
}
// std::cout << "ans=" <<ans << '\n'<<'\n';
// case 2: p*q <= n (p < q)
// p, q is distinct primes
ll cnt = 0;
ll q = 0;
for(int i = 0; i < (int)primes.size(); i++) //p
{
ll x = (ll)primes[i]; // p
q = n / x; //q
if(q <= x)continue;
if(q == 0)continue;
//std::cout << "p=" << x << '\n';
//std::cout << "q=" << q << '\n';
cnt = lehmer_pi(q);
if (q >= primes[i]) {
cnt -= lehmer_pi(x);
}
if (cnt <= 0) break;
//std::cout << "cnt=" <<cnt << '\n';
ans += cnt;
}
// std::cout << "p*q finish!" << '\n';
std::cout << ans << '\n';
}
int main(int argc, char const *argv[])
{
ll n;
init();
scanf("%lld", &n);
solve(n);
return 0;
}
PE501代码：
利用 Lucy_Hedgehog's method. \(O(n^{3/4})\)。跑10min左右...太慢了..
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 2000010;
vector<int> mark,prime;
void init()
{
mark.resize(maxn);
mark[1] = 1;
for (int i = 2; i < maxn; i++)
{
if (mark[i] == 0)
{
mark[i] = i;
prime.emplace_back(i);
}
for (int j = 0; j < (int)prime.size() && prime[j] <= mark[i]; j++)
{
int x = i * prime[j];
if (x >= maxn) break;
mark[x] = prime[j];
}
}
}
ll check(ll v, ll n, ll ndr, ll nv) {
return v >= ndr ? (n / v - 1) : (nv - v);
}
ll primenum(ll n) // O(n^(3/4))
{
assert(n>=1);
ll r = (ll)sqrt(n);
ll ndr = n / r;
assert(r*r <= n && (r+1)*(r+1) > n);
ll nv = r + ndr - 1;
std::vector<ll> S(nv+1);
std::vector<ll> V(nv+1);
for(ll i=0;i<r;i++) {
V[i] = n / (i+1);
}
for(ll i=r;i<nv;i++) {
V[i] = V[i-1] - 1;
}
for(ll i = 0;i<nv;i++) {
S[i] = V[i] - 1; //primes number
}
for(ll p=2;p<=r;p++) {
if(S[nv-p] > S[nv-p+1]) {
ll sp = S[nv-p+1]; // sum of primes smaller than p
ll p2 = p*p;
// std::cout << "p=" << p << '\n'; // p is prime
for(ll i=0;i<nv;i++) {
if(V[i] >= p2) {
S[i] -= 1LL * (S[check(V[i] / p, n, ndr, nv)] - sp);
}
else break;
}
}
}
ll ans = S[0];
return ans;
}
ll qpower(ll a, ll b)
{
ll res = 1;
while(b)
{
if (b&1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
void solve(ll n)
{
ll ans = 0;
// case 1: p^7 <= n ,p is a prime
for(int i = 0; i < (int)prime.size(); i++) {
// for example 2^7 = 128 <=n
if (qpower(prime[i], 7) <= n) {
//std::cout << primes[i] << '\n';
ans += 1;
}
else {
break;
}
}
std::cout << "p^7 finish!" << '\n';
// case 2: p^3*q <= n (p < q)
// p, q is distinct primes
ll cnt = 0;
for(int i = 0; i < (int)prime.size(); i++) //p
{
ll x = (ll)prime[i]*prime[i]*prime[i]; // p^3
x = n / x; //q
if(x == 0)continue;
cnt = primenum(x);
if (x >= prime[i]) {
cnt -= 1;
}
if (cnt <= 0) break;
ans += cnt;
}
std::cout << "p^3*q finish!" << '\n';
//case 3: p*q*r <= n (p < q < r)
// p,q,r is distinct primes
for(int i = 0; i < (int)prime.size(); i++) // p
{
if (qpower(prime[i], 3) > n) {
break;
}
for(int j = i+1; j < (int)prime.size(); j++) // q
{
ll x = n / ((ll)prime[i]*prime[j]);
if(x <= j)continue;
if(x == 0)continue;
cnt = primenum(x); // r
cnt -= j+1;
if (cnt <= 0) break;
ans += cnt;
}
}
std::cout << "p*q*r finish!" << '\n';
std::cout << ans << '\n';
cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.\n";
}
int main(int argc, char const *argv[])
{
/* code */
init();
solve(100);
solve(1000);
solve(1000000);
solve(1e12);
return 0;
}
利用Meisell-Lehmer's method。\(O(n^{2/3})\)。跑10s..
#include <bits/stdc++.h>
using namespace std;
const int MAX = 2e6+5;
const int M = 7;
vector<int> lp, primes, pi;
int phi[MAX+1][M+1], sz[M+1];
void factor_sieve()
{
lp.resize(MAX);
pi.resize(MAX);
lp[1] = 1;
pi[0] = pi[1] = 0;
for (int i = 2; i < MAX; i++) {
if (lp[i] == 0) {
lp[i] = i;
primes.emplace_back(i);
}
for (int j = 0; j < primes.size() && primes[j] <= lp[i]; j++) {
int x = i * primes[j];
if (x >= MAX) break;
lp[x] = primes[j];
}
pi[i] = primes.size();
}
}
void init()
{
factor_sieve();
for(int i = 0; i <= MAX; i++) {
phi[i][0] = i;
}
sz[0] = 1;
for(int i = 1; i <= M; i++) {
sz[i] = primes[i-1]*sz[i-1];
for(int j = 1; j <= MAX; j++) {
phi[j][i] = phi[j][i-1] - phi[j/primes[i-1]][i-1];
}
}
}
int sqrt2(long long x)
{
long long r = sqrt(x - 0.1);
while (r*r <= x) ++r;
return r - 1;
}
int cbrt3(long long x)
{
long long r = cbrt(x - 0.1);
while(r*r*r <= x) ++r;
return r - 1;
}
long long getphi(long long 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 <= primes[s-1]*primes[s-1])
{
return pi[x] - s + 1;
}
if(x <= primes[s-1]*primes[s-1]*primes[s-1] && x < MAX)
{
int sx = pi[sqrt2(x)];
long long ans = pi[x] - (sx+s-2)*(sx-s+1)/2;
for(int i = s+1; i <= sx; ++i) {
ans += pi[x/primes[i-1]];
}
return ans;
}
return getphi(x, s-1) - getphi(x/primes[s-1], s-1);
}
long long getpi(long long x)
{
if(x < MAX) return pi[x];
int cx = cbrt3(x), sx = sqrt2(x);
long long ans = getphi(x, pi[cx]) + pi[cx] - 1;
for(int i = pi[cx]+1, ed = pi[sx]; i <= ed; i++)
{
ans -= getpi(x/primes[i-1-1]) - i + 1;
}
return ans;
}
long long lehmer_pi(long long x)
{
if(x < MAX) return pi[x];
int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
int b = (int)lehmer_pi(sqrt2(x));
int c = (int)lehmer_pi(cbrt3(x));
long long sum = getphi(x, a) + (long long)(b + a - 2) * (b - a + 1) / 2;
for (int i = a + 1; i <= b; i++)
{
long long w = x / primes[i-1];
sum -= lehmer_pi(w);
if (i > c) continue;
long long lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++)
{
sum -= lehmer_pi(w / primes[j-1]) - (j - 1);
}
}
return sum;
}
long long power(long long a, long long b)
{
long long x = 1, y = a;
while(b)
{
if (b&1) x = x * y;
y = y * y;
b >>= 1;
}
return x;
}
void solve(long long n)
{
init();
long long ans = 0, val = 0;
// case : p^7 <= n ,p is a prime
for(int i = 0; i < primes.size(); i++) {
// for example 2^7 = 128 <=n
if (power(primes[i], 7) <= n) {
//std::cout << primes[i] << '\n';
ans += 1;
}
else {
break;
}
}
// std::cout << "ans = " << ans << '\n';
// case : p^3*q <= n (assume q > p for finding unique pairs)
// p, q is distinct primes
for(int i = 0; i < primes.size(); i++) //p
{
long long x = (long long)primes[i]*primes[i]*primes[i]; // p^3
x = n / x; //q
val = lehmer_pi(x);
if (x >= primes[i]) {
val -= 1;
}
if (val <= 0) break;
ans += val;
}
//case : p*q*r <= n (assume r > q > p for finding unique pairs)
// p,q,r is distinct primes
for(int i = 0; i < primes.size(); i++) // p
{
if (power(primes[i], 3) > n) {
break;
}
for(int j = i+1; j < primes.size(); j++) // q
{
long long x = n / ((long long)primes[i]*primes[j]);
if(x < j)continue;
val = lehmer_pi(x); // r
val -= j+1; // 减去 计算 <=x 的素数个数中多余的
if (val <= 0) break;
ans += val;
}
}
std::cout << ans << '\n';
cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.\n";
}
int main(int argc, char const *argv[])
{
/* code */
solve(1e12);
return 0;
}
转载于:https://www.cnblogs.com/LzyRapx/p/8453170.html
