Project Euler 501 Eight Divisors (數論)
阿新 • • 發佈:2018-02-18
sum scan 就是 ont tinc distinct class 套路 pair
題目鏈接:
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(106) = 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;
}
Project Euler 501 Eight Divisors (數論)