1 ///
2 /// @file  generate.cpp
3 ///
4 /// Copyright (C) 2021 Kim Walisch, <kim.walisch@gmail.com>
5 ///
6 /// This file is distributed under the BSD License. See the COPYING
7 /// file in the top level directory.
8 ///
9 
10 #include <generate.hpp>
11 #include <isqrt.hpp>
12 
13 #include <stdint.h>
14 #include <limits>
15 #include <vector>
16 
17 using std::vector;
18 using std::numeric_limits;
19 
20 namespace primecount {
21 
22 /// Generate a vector with the prime counts <= max
23 /// using the sieve of Eratosthenes
24 ///
25 vector<int32_t> generate_pi(int64_t max)
26 {
27   int64_t sqrt = isqrt(max);
28   int64_t size = max + 1;
29   vector<char> sieve(size, 1);
30 
31   for (int64_t i = 2; i <= sqrt; i++)
32     if (sieve[i])
33       for (int64_t j = i * i; j < size; j += i)
34         sieve[j] = 0;
35 
36   vector<int32_t> pi(size, 0);
37   int32_t pix = 0;
38 
S2(T x,int64_t y,int64_t z,int64_t c,T s2_approx,int threads,bool is_print)39   for (int64_t i = 2; i < size; i++)
40   {
41     pix += sieve[i];
42     pi[i] = pix;
43   }
44 
45   return pi;
46 }
47 
48 /// Generate a vector with Möbius function values.
49 /// This implementation is based on code by Rick Sladkey:
50 /// https://mathoverflow.net/q/99545
51 ///
52 vector<int32_t> generate_moebius(int64_t max)
53 {
54   int64_t sqrt = isqrt(max);
55   int64_t size = max + 1;
56   vector<int32_t> mu(size, 1);
57 
58   for (int64_t i = 2; i <= sqrt; i++)
59   {
60     if (mu[i] == 1)
61     {
62       for (int64_t j = i; j < size; j += i)
63         mu[j] *= (int32_t) -i;
64       for (int64_t j = i * i; j < size; j += i * i)
pi_deleglise_rivat_64(int64_t x,int threads,bool is_print)65         mu[j] = 0;
66     }
67   }
68 
69   for (int64_t i = 2; i < size; i++)
70   {
71     if (mu[i] == i)
72       mu[i] = 1;
73     else if (mu[i] == -i)
74       mu[i] = -1;
75     else if (mu[i] < 0)
76       mu[i] = 1;
77     else if (mu[i] > 0)
78       mu[i] = -1;
79   }
80 
81   return mu;
82 }
83 
84 /// Generate a vector with the least prime factors
85 /// of the integers <= max.
86 /// @Examples: lfp(2) = 2, lpf(15) = 3
87 ///
88 vector<int32_t> generate_lpf(int64_t max)
89 {
90   int64_t sqrt = isqrt(max);
91   int64_t size = max + 1;
92   vector<int32_t> lpf(size, 1);
93 
94   // By convention lfp(1) = +Infinity. Note that lpf(n) is
95   // named pmin(n) in Tomás Oliveira e Silva's paper:
96   // "Computing π(x): the combinatorial method".
97   // The reason why lfp(1) is defined to be +Infinity is
98   // that phi(x / 1, c) contributes to the ordinary leaves
99   // (S1) in the Lagarias-Miller-Odlyzko and
100   // Deleglise-Rivat prime counting algorithms. And
101   // lfp(1) = +Infinity allows to simplify that algorithm.
102   if (lpf.size() > 1)
103     lpf[1] = numeric_limits<int32_t>::max();
pi_deleglise_rivat_128(int128_t x,int threads,bool is_print)104 
105   for (int64_t i = 2; i <= sqrt; i++)
106     if (lpf[i] == 1)
107       for (int64_t j = i * i; j < size; j += i)
108         if (lpf[j] == 1)
109           lpf[j] = (int32_t) i;
110 
111   for (int64_t i = 2; i < size; i++)
112     if (lpf[i] == 1)
113       lpf[i] = (int32_t) i;
114 
115   return lpf;
116 }
117 
118 /// Generate a vector with the largest prime factors
119 /// of the integers <= max.
120 /// @Examples: mfp(2) = 2, mpf(15) = 5
121 ///
122 vector<int32_t> generate_mpf(int64_t max)
123 {
124   int64_t size = max + 1;
125   vector<int32_t> mpf(size, 1);
126 
127   for (int64_t i = 2; i <= max; i++)
128     if (mpf[i] == 1)
129       for (int64_t j = i; j < size; j += i)
130         mpf[j] = (int32_t) i;
131 
132   return mpf;
133 }
134 
135 } // namespace
136