Project Euler 501 Eight Divisors (数论)

Posted 2020-10-24 lzyrapx

题目链接：

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;
}