main.c

#include "prime.h"

typedef struct {
    uint64_t p;
    uint64_t sp;
} np_t;

static np_t fillBase(uint32_t *p, np_t n, unsigned char *f, uint64_t st,
uint64_t end, uint64_t lim) {
    #define push(a) do {\
        *p_ = st + i * 30 + a;\
        if (*p_ < L1C / 8) ++n.sp;\
        if (sq(*p_++) > lim) return (np_t){p_ - p + n.p, n.sp};\
    } while (0)
    p += n.p;
    uint32_t *p_ = p;
    unsigned i = 0;
    do {
        if (f[i] & 1 << 0) {
            push(1);
        }
        if (f[i] & 1 << 1) {
            push(7);
        }
        if (f[i] & 1 << 2) {
            push(11);
        }
        if (f[i] & 1 << 3) {
            push(13);
        }
        if (f[i] & 1 << 4) {
            push(17);
        }
        if (f[i] & 1 << 5) {
            push(19);
        }
        if (f[i] & 1 << 6) {
            push(23);
        }
        if (f[i] & 1 << 7) {
            push(29);
        }
    } while (st + ++i * 30 < end);
    return (np_t){p_ - p + n.p, n.sp};
    #undef push
}

static void sieve(unsigned char *f, uint64_t st, uint32_t *p, unsigned nsp) {
    memset(f, -1, BC);
    unsigned i = 0;
    do {
        uint64_t sti30 = st + i * 30;
        sieve_7_13(f + i, sti30);
        sieve_17(f + i, sti30);
        sieve_19(f + i, sti30);
        sieve_23(f + i, sti30);
        sieve_29(f + i, sti30);
        sieve_31(f + i, sti30);
        sieve_37(f + i, sti30);
        sieve_41(f + i, sti30);
        sieve_43(f + i, sti30);
        sieve_47(f + i, sti30);
        sieve_53(f + i, sti30);
        sieve_59(f + i, sti30);
        sieve_61(f + i, sti30);
        sieve_67(f + i, sti30);
        sieve_71(f + i, sti30);
        sieve_73(f + i, sti30);
        sieve_79(f + i, sti30);
        sieve_83(f + i, sti30);
        sieve_89(f + i, sti30);
        sieve_93_(f + i, sti30, p, nsp);
    } while ((i += L1C) < BC);
    sieve_93_(f, st, p + nsp, -1);
}

static unsigned count(unsigned char *f) {
    unsigned c = 0;
    unsigned i = 0;
    do {
        uint64_t qf;
        memcpy(&qf, f + i + 8 * 0, sizeof(qf));
        c += __builtin_popcountll(qf);
        memcpy(&qf, f + i + 8 * 1, sizeof(qf));
        c += __builtin_popcountll(qf);
        memcpy(&qf, f + i + 8 * 2, sizeof(qf));
        c += __builtin_popcountll(qf);
        memcpy(&qf, f + i + 8 * 3, sizeof(qf));
        c += __builtin_popcountll(qf);
        memcpy(&qf, f + i + 8 * 4, sizeof(qf));
        c += __builtin_popcountll(qf);
        memcpy(&qf, f + i + 8 * 5, sizeof(qf));
        c += __builtin_popcountll(qf);
        memcpy(&qf, f + i + 8 * 6, sizeof(qf));
        c += __builtin_popcountll(qf);
        memcpy(&qf, f + i + 8 * 7, sizeof(qf));
        c += __builtin_popcountll(qf);
    } while ((i += 8 * 8) < BC);
    return c;
}

static uint64_t countPrimes(uint64_t lim) {
    _Alignas(64) unsigned char f[BC + 256];
    uint32_t *p = aligned_alloc(64, sqrt(lim));
    unsigned bfsz = 256;
    np_t n = fillBase(p, (np_t){0}, countPrimes_bf, 0, bfsz * 30, -1);
    sieve(f, 0, p, n.sp);
    uint64_t c = count(f);
    unsigned r = lim % 30;
    if (!r) {
        --lim;
    } else if (1 < r && r < 7) {
        lim -= r - 1;
    } else if (7 < r && r < 11) {
        lim -= r - 7;
    } else if (13 < r && r < 17) {
        lim -= r - 13;
    } else if (r == 18) {
        --lim;
    } else if (19 < r && r < 23) {
        lim -= r - 19;
    } else if (23 < r && r < 29) {
        lim -= r - 23;
    }
    unsigned bc30 = BC * 30;
    r = lim % bc30;
    uint64_t lim_ = lim - r + bc30;
    n = fillBase(p, n, f + bfsz, bfsz * 30, bc30, lim_);
    uint64_t st = bc30;
    for (; st + bc30 < lim_; st += bc30) {
        sieve(f, st, p, n.sp);
        if (sq(st + 1) < lim_) {
            n = fillBase(p, n, f, st, st + bc30, lim_);
        }
        c += count(f);
    }
    sieve(f, st, p, n.sp);
    unsigned q = r / 30;
    memset(f + q + 1, 0, BC - q);
    f[q] &= countPrimes_rm[r % 30 / 2];
    c += count(f);
    return c + 5;
}

int main(int argc, char **argv) {
    if (argc != 2) {
        puts("must have 1 argument");
        return 1;
    }
    L1C = sysconf(_SC_LEVEL1_DCACHE_SIZE) - (1 << 12);
    BC = L1C * 8;
    uint64_t lim = strtod(argv[1], 0);
    if (lim < 100000000) {
        puts("input < 1e8");
        return 1;
    }
    printf("%"PRIu64"\n", countPrimes(lim));
    return 0;
}