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