Legendre prime counting function: Difference between revisions

=={{header|C}}==

{{trans|Vlang}}

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).

<syntaxhighlight lang="c">#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#include <stdint.h>

#include <time.h>

const uint8_t masks[8] = {1, 2, 4, 8, 16, 32, 64, 128};

#define half(n) (((n) - 1) >> 1)

#define divide(nm, d) ((int)((double)nm / (double)d))

int64_t countPrimes(uint64_t n) {

if (n < 3) return (n < 2) ? 0 : 1;

int rtlmt = (int)sqrt((double)n);

int mxndx = (rtlmt - 1) / 2;

int arrlen = mxndx + 1;

int i;

uint32_t *smalls = malloc(arrlen * 4);

uint32_t *roughs = malloc(arrlen * 4);

int64_t *larges = malloc(arrlen * 8);

for (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 = (mxndx + 8) / 8;

uint8_t *cullbuf = calloc(cullbuflen, 1);

int nbps = 0;

int rilmt = arrlen;

for (i = 1; ; ++i) {

int sqri = (i + i) * (i + 1);

if (sqri > mxndx) break;

if (cullbuf[i >> 3] & masks[i & 7]) continue;

cullbuf[i >> 3] |= masks[i & 7];

int bp = i + i + 1;

int c, ori, pm;

for (c = sqri; c < arrlen; c += bp) cullbuf[c >> 3] |= masks[c & 7];

int nri = 0;

for (ori = 0; ori < rilmt; ++ori) {

int q = (int)roughs[ori];

int pi = q >> 1;

if (cullbuf[pi >> 3] & masks[pi & 7]) continue;

int d = q * bp;

int64_t t = (d <= rtlmt) ? larges[(int)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)q;

nri++;

}

int si = mxndx;

for (pm = (rtlmt/bp - 1) | 1; pm >= bp; pm -= 2) {

uint32_t c = smalls[pm >> 1];

int e = (pm * bp) >> 1;

for ( ; si >= e; --si) smalls[si] -= c - (uint32_t)nbps;

}

rilmt = nri;

nbps++;

}

int64_t ans = larges[0] + (int64_t)((rilmt + 2*(nbps - 1)) * (rilmt - 1) / 2);

int ri, sri;

for (ri = 1; ri < rilmt; ++ri) ans -= larges[ri];

for (ri = 1; ; ++ri) {

uint64_t p = (uint64_t)roughs[ri];

uint64_t m = n / p;

int ei = (int)smalls[half((int)m/p)] - nbps;

if (ei <= ri) break;

ans -= (int64_t)((ei - ri) * (nbps + ri - 1));

for (sri = ri + 1; sri < ei + 1; ++sri) {

ans += (int64_t)smalls[half(divide(m, (uint64_t)roughs[sri]))];

}

}

free(smalls);

free(roughs);

free(larges);

free(cullbuf);

return ans + 1;

}

int main() {

uint64_t n;

int i;

clock_t start = clock();

for (i = 0, n = 1; i < 10; ++i, n *= 10) {

printf("10^%d %ld\n", i, countPrimes(n));

}

clock_t end = clock();

printf("\nTook %f seconds\n", (double) (end - start) / CLOCKS_PER_SEC);

return 0;

}</syntaxhighlight>

{{out}}

<pre>

10^0 0

10^1 4

10^2 25

10^3 168

10^4 1229

10^5 9592

10^6 78498

10^7 664579

10^8 5761455

10^9 50847534

Took 0.003942 seconds

</pre>