题目187：半质数

欧拉计划

Semiprimes

A composite is a number containing at least two prime factors. For example, 15 = 3 × 5; 9 = 3 × 3; 12 = 2 × 2 × 3.

There are ten composites below thirty containing precisely two, not necessarily distinct, prime factors: 4, 6, 9, 10, 14, 15, 21, 22, 25, 26.

How many composite integers, n < 10⁸, have precisely two, not necessarily distinct, prime factors?

题目：

一个合数是至少两个质数因子的乘积，比如 15 = 3 × 5； 9 = 3 × 3； 12 = 2 × 2 × 3。

30 以内的数中有 10 个合数，它们有且仅有两个质数因子，这两个因子可能相同：4, 6, 9, 10, 14, 15, 21, 22, 25, 26。

请问，对于 n < 10⁸，有多少个数字只有两个质数因子？这两个质数因子可能相同。

jerryxjr1220

1. import numpy as np
   
2. from numba import jit
   
3. @jit
   
4. def solve(n):
   
5.         sieve=np.ones(n+1,dtype=bool)
   
6.         for i in range(2, int((n+1)**0.5+1)):
   
7.                 if sieve[i]:
   
8.                         sieve[i*i::i]=False
   
9.         plist = np.nonzero(sieve)[0][2:]
   
10.         count = 0
    
11.         for i in range(len(plist)):
    
12.                 j = i
    
13.                 while plist[i]*plist[j] < n:
    
14.                         count += 1
    
15.                         j += 1
    
16.         return count
    
17. print(solve(10**8))

复制代码

17427258
[Finished in 4.5s]

永恒的蓝色梦想

C++ 322ms

1. #include<iostream>
   
2. #include<vector>
   
3. #include<cstring>
   
4. using namespace std;
   
5. 
   
6. 
   
7. 
   
8. template<size_t size, typename T>
   
9. void euler_sieve(bool (&is_prime)[size], vector<T>& primes) {
   
10.     T i, k;
    
11.     memset(is_prime, -1, size);
    
12.     is_prime[0] = is_prime[1] = 0;
    
13. 
    
14.     for (T i = 2; i < size; i++) {
    
15.         if (is_prime[i]) {
    
16.             primes.push_back(i);
    
17.         }
    
18. 
    
19.         for (T j : primes) {
    
20.             k = i * j;
    
21. 
    
22.             if (k >= size) {
    
23.                 break;
    
24.             }
    
25. 
    
26.             is_prime[k] = false;
    
27. 
    
28.             if (!(i % j)) {
    
29.                 break;
    
30.             }
    
31.         }
    
32.     }
    
33. }
    
34. 
    
35. 
    
36. 
    
37. int main() {
    
38.     ios::sync_with_stdio(false);
    
39.    
    
40.     constexpr unsigned int UPPER_BOUND = 100000000;
    
41.     static bool is_prime[UPPER_BOUND / 2];
    
42.     vector<unsigned int> primes;
    
43.     unsigned int i, j, result = 0;
    
44. 
    
45.     euler_sieve(is_prime, primes);
    
46. 
    
47. 
    
48.     for (i = 0; primes[i] * primes[i] < UPPER_BOUND; result += ++i);
    
49. 
    
50.     for (j = i; i < primes.size(); i++) {
    
51.         while (primes[i] * primes[j] >= UPPER_BOUND) {
    
52.             j--;
    
53.         }
    
54. 
    
55.         result += j;
    
56.         result++;
    
57.     }
    
58. 
    
59. 
    
60.     cout << result << endl;
    
61.     return 0;
    
62. }

复制代码

guosl

优化来了。

1. /*
   
2. 答案：17427258
   
3. 耗时：1.41251
   
4.       0.643071 (8核)
   
5. */
   
6. #include <iostream>
   
7. #include <vector>
   
8. #include <cmath>
   
9. #include <omp.h>
   
10. using namespace std;
    
11. 
    
12. const int N = 100000000;
    
13. const int M = int(sqrt(double(N)));
    
14. 
    
15. char cp[N];//cp[i]==1表示i为素数
    
16. bool pm[N];//pm[i]==true表示i为两个素数的乘积
    
17. 
    
18. int main(void)
    
19. {
    
20.   double t = omp_get_wtime();
    
21.   memset(cp, 1, sizeof(cp));
    
22.   int nCount = 0;
    
23.   bool bContinue = true;
    
24. #pragma omp parallel firstprivate(bContinue) shared(cp,pm) reduction(+:nCount)
    
25.   {
    
26.     for (int i = 2; i <= M; ++i)
    
27.     {
    
28. #pragma omp single copyprivate(bContinue)
    
29.       {
    
30.         bContinue = true;
    
31.         if (cp[i] == 0)
    
32.           bContinue = false;
    
33.       }
    
34.       if (!bContinue)
    
35.         continue;
    
36. #pragma omp for
    
37.       for (int j = i * i; j < N; j += i)
    
38.         cp[j] = 0;
    
39.     }
    
40. #pragma omp barrier
    
41. #pragma omp for schedule(dynamic)
    
42.     for (int i = 2; i <= M; ++i)
    
43.     {
    
44.       if (cp[i] == 0)
    
45.         continue;
    
46.       for (int j = i; i*j < N; ++j)
    
47.       {
    
48.         if (cp[j] == 1)
    
49.           pm[i*j] = true;
    
50.       }
    
51.     }
    
52. #pragma omp for
    
53.     for (int i = 4; i < N; ++i)
    
54.     {
    
55.       if (pm[i])
    
56.         ++nCount;
    
57.     }
    
58.   }
    
59.   t = omp_get_wtime() - t;
    
60.   cout << nCount << endl << t << endl;
    
61.   return 0;
    
62. }
    

复制代码

代号-K

这题本质还是素数的求解，只是在求解的过程中加入了一些判断。

代号-K

永恒的蓝色梦想 发表于 2020-5-17 20:09
谢谢指导

有个优化过的帖子，不知道在哪里了，这个帖子不是优化的，重新写一个吧。