Legendre prime counting function: Difference between revisions

→‎{{header|C}}: Aligned with changes to V example of which it is a translation.
(→‎{{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:
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) 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>
#include <math.h>
Line 47 ⟶ 49:
const uint8_t masks[8] = {1, 2, 4, 8, 16, 32, 64, 128};
 
#define half(n) ((int64_t)((n) - 1) >> 1)
 
#define divide(nm, d) ((intuint64_t)((double)nm / (double)d))
 
int64_t countPrimes(uint64_t n) {
if (n < 3) return (n < 2) ? 0 : 1;
intuint64_t rtlmt = (intuint64_t)sqrt((double)n);
intint64_t mxndx = (int64_t)((rtlmt - 1) / 2);
int arrlen = (int)(mxndx + 1);
int i;
uint32_t *smalls = malloc(arrlen * 4);
uint32_t *roughs = malloc(arrlen * 4);
int64_t *larges = malloc(arrlen * 8);
for (int i = 0; i < arrlen; ++i) {
smalls[i] = (uint32_t)i;
roughs[i] = (uint32_t)(i + i + 1);
larges[i] = (int64_t)((n/(uint64_t)(i + i + 1) - 1) / 2);
}
int cullbuflen = (int)((mxndx + 8) / 8);
uint8_t *cullbuf = calloc(cullbuflen, 1);
intint64_t nbps = 0;
int rilmt = arrlen;
for (int64_t i = 1; ; ++i) {
intint64_t sqri = (i + i) * (i + 1);
if (sqri > mxndx) break;
if (cullbuf[i >> 3] & masks[i & 7]) continue;
cullbuf[i >> 3] |= masks[i & 7];
intuint64_t bp = (uint64_t)(i + i + 1);
for (int64_t c = sqri; c < (int64_t)arrlen; c += (int64_t)bp) {
int c, ori, pm;
for (c = sqri; c < arrlen; c += bp) cullbuf[c >> 3] |= masks[c & 7];
int i; }
int nri = 0;
for (int ori = 0; ori < rilmt; ++ori) {
intuint32_t qr = (int)roughs[ori];
intint64_t pirci = q(int64_t)(r >> 1);
if (cullbuf[pirci >> 3] & masks[pirci & 7]) continue;
intuint64_t d = q(uint64_t)r * bp;
int64_t t = (d <= rtlmt) ? larges[(intint64_t)smalls[d >> 1] - nbps] :
(int64_t)smalls[half(divide(n, (uint64_t)d))];
larges[nri] = larges[ori] - t + (int64_t)nbps;
roughs[nri] = (uint32_t)qr;
nri++;
}
intint64_t si = mxndx;
for (uint64_t pm = (rtlmt/bp - 1) | 1; pm >= bp; pm -= 2) {
uint32_t c = smalls[pm >> 1];
intuint64_t e = (pm * bp) >> 1;
for ( ; si >= (int64_t)e; --si) smalls[si] -= c - (uint32_t)nbps;
}
rilmt = nri;
Line 104 ⟶ 106:
uint64_t p = (uint64_t)roughs[ri];
uint64_t m = n / p;
int ei = (int)smalls[half((intuint64_t)m/p)] - nbps;
if (ei <= ri) break;
ans -= (int64_t)((ei - ri) * (nbps + ri - 1));
Line 143 ⟶ 145:
10^9 50847534
 
Took 0.003942003843 seconds
</pre>
 
9,476

edits