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