码迷，mamicode.com
首页 > 其他好文 > 详细

Project Euler 501 Eight Divisors (数论)

时间：2018-02-18 21:30:40      阅读：274      评论：0      收藏：0      [点我收藏+]

标签：resize   左右   ref   http   type   tin   air   finish   rod   

题目链接：

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 (数论)

标签：resize   左右   ref   http   type   tin   air   finish   rod   

原文地址：https://www.cnblogs.com/LzyRapx/p/8453170.html
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！