1 ///
2 /// @file pmath.hpp
3 /// @brief Auxiliary math functions for primesieve.
4 ///
5 /// Copyright (C) 2018 Kim Walisch, <kim.walisch@gmail.com>
6 ///
7 /// This file is distributed under the BSD License. See the COPYING
8 /// file in the top level directory.
9 ///
10
11 #ifndef PMATH_HPP
12 #define PMATH_HPP
13
14 #include <stdint.h>
15 #include <algorithm>
16 #include <cmath>
17 #include <cstddef>
18 #include <limits>
19 #include <type_traits>
20
21 namespace {
22
23 template <typename X, typename Y>
ceilDiv(X x,Y y)24 inline X ceilDiv(X x, Y y)
25 {
26 return (X) ((x + y - 1) / y);
27 }
28
29 template <typename T>
isPow2(T x)30 constexpr bool isPow2(T x)
31 {
32 return x != 0 && (x & (x - 1)) == 0;
33 }
34
35 template <typename T>
numberOfBits()36 constexpr T numberOfBits()
37 {
38 return (T) std::numeric_limits<
39 typename std::make_unsigned<T>::type
40 >::digits;
41 }
42
43 template <typename T>
floorPow2(T x)44 inline T floorPow2(T x)
45 {
46 for (T i = 1; i < numberOfBits<T>(); i += i)
47 x |= (x >> i);
48
49 return x - (x >> 1);
50 }
51
52 template <typename T>
ilog2(T x)53 inline T ilog2(T x)
54 {
55 T log2 = 0;
56 T bits = numberOfBits<T>();
57
58 for (T i = bits / 2; i > 0; i /= 2)
59 {
60 T one = 1;
61 if (x >= (one << i))
62 {
63 x >>= i;
64 log2 += i;
65 }
66 }
67
68 return log2;
69 }
70
71 #if __cplusplus >= 201402L
72
73 /// C++14 compile time square root using binary search
74 template <typename T>
ctSqrt(T x,T lo,T hi)75 constexpr T ctSqrt(T x, T lo, T hi)
76 {
77 if (lo == hi)
78 return lo;
79
80 const T mid = (lo + hi + 1) / 2;
81
82 if (x / mid < mid)
83 return ctSqrt<T>(x, lo, mid - 1);
84 else
85 return ctSqrt(x, mid, hi);
86 }
87
88 #else
89
90 #define MID ((lo + hi + 1) / 2)
91
92 /// C++11 compile time square root using binary search
93 template <typename T>
ctSqrt(T x,T lo,T hi)94 constexpr T ctSqrt(T x, T lo, T hi)
95 {
96 return lo == hi ? lo : ((x / MID < MID)
97 ? ctSqrt<T>(x, lo, MID - 1) : ctSqrt<T>(x, MID, hi));
98 }
99
100 #endif
101
102 template <typename T>
ctSqrt(T x)103 constexpr T ctSqrt(T x)
104 {
105 return ctSqrt<T>(x, 0, x / 2 + 1);
106 }
107
108 template <typename T>
isqrt(T x)109 inline T isqrt(T x)
110 {
111 T r = (T) std::sqrt((double) x);
112
113 constexpr T maxSqrt = ctSqrt(std::numeric_limits<T>::max());
114 r = std::min(r, maxSqrt);
115
116 while (r * r > x)
117 r--;
118 while (x - r * r > r * 2)
119 r++;
120
121 return r;
122 }
123
124 /// Returns 2^64-1 if (x + y) > 2^64-1
checkedAdd(uint64_t x,uint64_t y)125 inline uint64_t checkedAdd(uint64_t x, uint64_t y)
126 {
127 if (x >= std::numeric_limits<uint64_t>::max() - y)
128 return std::numeric_limits<uint64_t>::max();
129 else
130 return x + y;
131 }
132
133 /// Returns 0 if (x - y) < 0
checkedSub(uint64_t x,uint64_t y)134 inline uint64_t checkedSub(uint64_t x, uint64_t y)
135 {
136 if (x > y)
137 return x - y;
138 else
139 return 0;
140 }
141
142 template <typename A, typename B, typename C>
inBetween(A min,B x,C max)143 inline B inBetween(A min, B x, C max)
144 {
145 using T = typename std::common_type<A, B, C>::type;
146
147 if ((T) x < (T) min)
148 return (B) min;
149 if ((T) x > (T) max)
150 return (B) max;
151
152 return x;
153 }
154
155 /// primeCountApprox(x) >= pi(x)
primeCountApprox(uint64_t start,uint64_t stop)156 inline std::size_t primeCountApprox(uint64_t start, uint64_t stop)
157 {
158 if (start > stop)
159 return 0;
160 if (stop <= 10)
161 return 4;
162
163 // pi(x) <= x / (log(x) - 1.1) + 5, for x >= 4
164 double x = (double) stop;
165 double logx = std::log(x);
166 double div = logx - 1.1;
167 double pix = (stop - start) / div + 5;
168
169 return (std::size_t) pix;
170 }
171
primeCountApprox(uint64_t stop)172 inline std::size_t primeCountApprox(uint64_t stop)
173 {
174 return primeCountApprox(0, stop);
175 }
176
177 /// Approximation of the maximum prime gap near n
178 template <typename T>
maxPrimeGap(T n)179 inline T maxPrimeGap(T n)
180 {
181 double x = (double) n;
182 x = std::max(8.0, x);
183 double logx = std::log(x);
184 double prime_gap = logx * logx;
185
186 return (T) prime_gap;
187 }
188
189 } // namespace
190
191 #endif
192