Legendre prime counting function: Difference between revisions
Content added Content deleted
GordonBGood (talk | contribs) (→{{header|Vlang}}: correction to increase precision to extend range...) |
(→{{header|C}}: Aligned with changes to V example of which it is a translation.) |
||
Line 38: | Line 38: | ||
Using fixed width types so requires C99 or later. |
Using fixed width types so requires C99 or later. |
||
Surprisingly only half as fast as Go on the same machine (GCC -O3 Ubuntu 22.04). |
Surprisingly only half as fast as Go on the same machine (GCC -O3 Ubuntu 22.04) for a billion numbers. |
||
However, the position is reversed if we let it run through to a trillion numbers - about 100 ms for C compared to 160 ms for Go. |
|||
<syntaxhighlight lang="c">#include <stdio.h> |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <math.h> |
#include <math.h> |
||
Line 47: | Line 49: | ||
const uint8_t masks[8] = {1, 2, 4, 8, 16, 32, 64, 128}; |
const uint8_t masks[8] = {1, 2, 4, 8, 16, 32, 64, 128}; |
||
#define half(n) (((n) - 1) >> 1) |
#define half(n) ((int64_t)((n) - 1) >> 1) |
||
#define divide(nm, d) (( |
#define divide(nm, d) ((uint64_t)((double)nm / (double)d)) |
||
int64_t countPrimes(uint64_t n) { |
int64_t countPrimes(uint64_t n) { |
||
if (n < 3) return (n < 2) ? 0 : 1; |
if (n < 3) return (n < 2) ? 0 : 1; |
||
uint64_t rtlmt = (uint64_t)sqrt((double)n); |
|||
int64_t mxndx = (int64_t)((rtlmt - 1) / 2); |
|||
int arrlen = mxndx + 1; |
int arrlen = (int)(mxndx + 1); |
||
⚫ | |||
uint32_t *smalls = malloc(arrlen * 4); |
uint32_t *smalls = malloc(arrlen * 4); |
||
uint32_t *roughs = malloc(arrlen * 4); |
uint32_t *roughs = malloc(arrlen * 4); |
||
int64_t *larges = malloc(arrlen * 8); |
int64_t *larges = malloc(arrlen * 8); |
||
for (i = 0; i < arrlen; ++i) { |
for (int i = 0; i < arrlen; ++i) { |
||
smalls[i] = (uint32_t)i; |
smalls[i] = (uint32_t)i; |
||
roughs[i] = (uint32_t)(i + i + 1); |
roughs[i] = (uint32_t)(i + i + 1); |
||
larges[i] = (int64_t)((n/(uint64_t)(i + i + 1) - 1) / 2); |
larges[i] = (int64_t)((n/(uint64_t)(i + i + 1) - 1) / 2); |
||
} |
} |
||
int cullbuflen = (mxndx + 8) / 8; |
int cullbuflen = (int)((mxndx + 8) / 8); |
||
uint8_t *cullbuf = calloc(cullbuflen, 1); |
uint8_t *cullbuf = calloc(cullbuflen, 1); |
||
int64_t nbps = 0; |
|||
int rilmt = arrlen; |
int rilmt = arrlen; |
||
for (i = 1; ; ++i) { |
for (int64_t i = 1; ; ++i) { |
||
int64_t sqri = (i + i) * (i + 1); |
|||
if (sqri > mxndx) break; |
if (sqri > mxndx) break; |
||
if (cullbuf[i >> 3] & masks[i & 7]) continue; |
if (cullbuf[i >> 3] & masks[i & 7]) continue; |
||
cullbuf[i >> 3] |= masks[i & 7]; |
cullbuf[i >> 3] |= masks[i & 7]; |
||
uint64_t bp = (uint64_t)(i + i + 1); |
|||
for (int64_t c = sqri; c < (int64_t)arrlen; c += (int64_t)bp) { |
|||
int c, ori, pm; |
|||
cullbuf[c >> 3] |= masks[c & 7]; |
|||
⚫ | |||
int nri = 0; |
int nri = 0; |
||
for (ori = 0; ori < rilmt; ++ori) { |
for (int ori = 0; ori < rilmt; ++ori) { |
||
uint32_t r = roughs[ori]; |
|||
int64_t rci = (int64_t)(r >> 1); |
|||
if (cullbuf[ |
if (cullbuf[rci >> 3] & masks[rci & 7]) continue; |
||
uint64_t d = (uint64_t)r * bp; |
|||
int64_t t = (d <= rtlmt) ? larges[( |
int64_t t = (d <= rtlmt) ? larges[(int64_t)smalls[d >> 1] - nbps] : |
||
(int64_t)smalls[half(divide(n, |
(int64_t)smalls[half(divide(n, d))]; |
||
larges[nri] = larges[ori] - t + |
larges[nri] = larges[ori] - t + nbps; |
||
roughs[nri] = |
roughs[nri] = r; |
||
nri++; |
nri++; |
||
} |
} |
||
int64_t si = mxndx; |
|||
for (pm = (rtlmt/bp - 1) | 1; pm >= bp; pm -= 2) { |
for (uint64_t pm = (rtlmt/bp - 1) | 1; pm >= bp; pm -= 2) { |
||
uint32_t c = smalls[pm >> 1]; |
uint32_t c = smalls[pm >> 1]; |
||
uint64_t e = (pm * bp) >> 1; |
|||
for ( ; si >= e; --si) smalls[si] -= c - (uint32_t)nbps; |
for ( ; si >= (int64_t)e; --si) smalls[si] -= c - (uint32_t)nbps; |
||
} |
} |
||
rilmt = nri; |
rilmt = nri; |
||
Line 104: | Line 106: | ||
uint64_t p = (uint64_t)roughs[ri]; |
uint64_t p = (uint64_t)roughs[ri]; |
||
uint64_t m = n / p; |
uint64_t m = n / p; |
||
int ei = (int)smalls[half(( |
int ei = (int)smalls[half((uint64_t)m/p)] - nbps; |
||
if (ei <= ri) break; |
if (ei <= ri) break; |
||
ans -= (int64_t)((ei - ri) * (nbps + ri - 1)); |
ans -= (int64_t)((ei - ri) * (nbps + ri - 1)); |
||
Line 143: | Line 145: | ||
10^9 50847534 |
10^9 50847534 |
||
Took 0. |
Took 0.003843 seconds |
||
</pre> |
</pre> |
||