High throughput prime numbers

Question

This challenge is inspired by the High throughput Fizz Buzz challenge.

The goal

Generate a list of prime numbers up to 10,000,000,000,000,000. The output of primes should be in decimal digits followed by a newline character '\n' in ascending order starting from the lowest prime 2. You may not skip a prime number or output a composite number.

Scoring

Your program's throughput will be measured on my Linux machine by the following command.

(timeout 1m ./your-program) | pv > /dev/null

At the timeout after 1 minute, your score will be the total size of output measured by pv.

An example

This is a simple example of a conforming program in C. It produces 49.6MiB of data in total for a minute, on my machine with 2.4GHz 4-core CPU and 4GiB RAM.

#include <stdio.h>

int main() {
    puts("2");
    for (long long i = 3; i < 10000000000000000; i += 2) {
        for (long long j = 3; j * j <= i; j += 2) {
            if (i % j == 0) {
                goto next;
            }
        }
        printf("%lld\n", i);
    next:;
    }
    return 0;
}

Rules

You should only print what's specified in the goal. You may not print garbage characters including ones that do not appear on the terminal.

The maximum size of your source code is 64Kib.

Otherwise, I'll accept any code that can be run on my Linux machine with 4 cores and AVX2 support.

---

Leaderboard

Contestant	Language	Score
xiver77	C + asm	9.71GiB
alephalpha	C + PARI/GP	6.85GiB
Neil	C	5.47GiB
Seggan	Jyxal	120MiB
(example)	C	56.2MiB
emanresu A	Vyxal	21.8MiB

Answer 1

Vyxal, ~7.3Mb

0 { 200 ( ∆Ṗ : ) W ṫ $ ⁋,

Vyxal really isn't fast.

Because Vyxal is stack-based, an operation such as 1 + (incrementing) is actually push(1), push(add(pop(),pop())), which makes it incredibly difficult to optimize anything.

Vyxal uses Sympy's isprime, which looks for small factors, then for numbers \$ n < 2^{64}\$ it runs a set of Miller-Rabin tests. In other words, no optimization there.

So, the above program optimizes by buffering the output into groups of 200, which through experimentation produces the most output.

This version

{ is loop forever. 200 ( ∆Ṗ : ) means "200 times, get the prime after the top of stack and push a copy of it". W ṫ $ gets all but the last one (which we use for the next iteration of the forever loop), then ⁋, prints that joined by newlines.

Other attempts

~10Kb

0 {›:ǎ,

This one's "forever, print the nth prime and increment n". It's kinda slow.

~165Kb

0 { ∆Ṗ …

This one's "forever, get the next prime and print it". It's an order of magnitude faster than the previous, but still slow.

~780Kb

Þp(n,

5x faster than the previous, but still not very fast, this one loops through an infinite generator of primes and prints each one.

~4.3Mb

0 { ⟨⟩ →primes 100 ( ∆Ṗ : ←primes $ J →primes) ←primes ⁋ ,

Like the main submission but uses an array instead of the stack, making it slower.

* ^{amounts approximated with online interpreter, unreliable}

Answer 2

C + assembly (nasm) 9.71GiB

I got to know a program called primesieve from a comment and was impressed to know that it can generate prime numbers so fast. Somehow the existence of the program triggered me in a weird way that I felt I should and can build a program that is faster than primesieve. I don't know why, but really somehow it felt like an easy opponent, and it wasn't. Over a month the goal of building a faster program than primesieve drained so much of my energy, but unfortunately, the best I could reach was about 0.8 times the speed of primesieve ⁽¹⁾.

Well, I could still make a faster one if I just copy the same algorithm primesieve used, and apply some of optimization techniques I gathered while trying to beat it. Maybe I will do that later in time, but for now I feel that just re-implementing the same algorithm is not a fun task. Also, one of the reasons why I didn't use the same algorithm is because that algorithm is quite complex and thus difficult to implement it efficiently. I really thought a simpler approach will be faster, but now it is apparent why primesieve took the hard route of implementing such a complex algorithm. It was obviously worth it.

Anyway the simpler route still took me more than a month of my free time and a lot of sleeplessness. The result is just a slower program than primesieve but it is still enough as a fast program here.

I didn't do any IO optimizations, but just a loop of printfs. There are many ways to make the output faster, but I don't have the energy to do it now.

---

⁽¹⁾ If you're interested this is an alternative main.c which actually counts primes like primesieve. The gap in speed is about 10~30%, which gets bigger as the sieving range gets larger, which is why my choice of algorithm was wrong and primesieve's was right.

prime.h

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <inttypes.h>
#include <math.h>
#include <unistd.h>

typedef unsigned char v16b
__attribute__((vector_size(16), aligned(16), may_alias));

extern unsigned BC;
extern unsigned L1C;
extern unsigned char PADD[];
extern unsigned char WIX[];
extern unsigned char countPrimes_bf[];
extern unsigned char countPrimes_rm[];
extern v16b sieve_7_13_M7[];
extern v16b sieve_7_13_M11[];
extern v16b sieve_7_13_M13[];

void sieve_7_13(unsigned char *, uint64_t);
void sieve_17(unsigned char *, uint64_t);
void sieve_19(unsigned char *, uint64_t);
void sieve_23(unsigned char *, uint64_t);
void sieve_29(unsigned char *, uint64_t);
void sieve_31(unsigned char *, uint64_t);
void sieve_37(unsigned char *, uint64_t);
void sieve_41(unsigned char *, uint64_t);
void sieve_43(unsigned char *, uint64_t);
void sieve_47(unsigned char *, uint64_t);
void sieve_53(unsigned char *, uint64_t);
void sieve_59(unsigned char *, uint64_t);
void sieve_61(unsigned char *, uint64_t);
void sieve_67(unsigned char *, uint64_t);
void sieve_71(unsigned char *, uint64_t);
void sieve_73(unsigned char *, uint64_t);
void sieve_79(unsigned char *, uint64_t);
void sieve_83(unsigned char *, uint64_t);
void sieve_89(unsigned char *, uint64_t);
void sieve_93_(unsigned char *, uint64_t, uint32_t *, unsigned);

static inline uint64_t sq(uint64_t a) {
    return a * a;
}

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 void print(unsigned char *f, uint64_t st) {
    #define p(a) printf("%"PRIu64"\n", st + i * 30 + a);
    uint64_t i = 0;
    do {
        if (f[i] & 1 << 0) {
            p(1);
        }
        if (f[i] & 1 << 1) {
            p(7);
        }
        if (f[i] & 1 << 2) {
            p(11);
        }
        if (f[i] & 1 << 3) {
            p(13);
        }
        if (f[i] & 1 << 4) {
            p(17);
        }
        if (f[i] & 1 << 5) {
            p(19);
        }
        if (f[i] & 1 << 6) {
            p(23);
        }
        if (f[i] & 1 << 7) {
            p(29);
        }
    } while (++i < BC);
    #undef p
}

static void printPrimes(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);
    *f &= ~1;
    print(f, 0);
    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_);
        }
        print(f, st);
    }
    sieve(f, st, p, n.sp);
    unsigned q = r / 30;
    memset(f + q + 1, 0, BC - q);
    f[q] &= countPrimes_rm[r % 30 / 2];
    print(f, st);
}

int main() {
    L1C = sysconf(_SC_LEVEL1_DCACHE_SIZE) - (1 << 12);
    BC = L1C * 8;
    puts("2\n3\n5\n7\n11\n13");
    printPrimes(1e16);
    return 0;
}

lut.s

    align 16
global PADD
PADD:
    db 1, 7, 7, 7, 7, 7, 7, 11, 11, 11, 11, 13, 13, 17, 17, 17, 17
    db 19, 19, 23, 23, 23, 23, 29, 29, 29, 29, 29, 29, 31, -1, -1
global WIX
WIX:
    db 0, -1, -1, 1, -1, 2, 3, -1, 4, 5, -1, 6, -1, -1, 7, -1
global countPrimes_bf
countPrimes_bf:
    dq 0xf93ddbb67e000000, 0x9eeda6eaf31e4fd5, 0xa559dd3bd3d30ce6
    dq 0x56a61e78bd92676a, 0x554c2ade2dade356, 0xf8a154039ff0a3d9,
    dq 0x3a13f666e944fd2e, 0x54bf11453a2b4cb8, 0x4f8cbcc8b37ac18c,
    dq 0xef17c19b71715821, 0x468c83e5081a9654, 0x87588f9265aefb72,
    dq 0xa0e3266581d892d2, 0x99eb813c26c73811, 0x4d33f3243e88518d,
    dq 0x4c58b42aa71c8b5a, 0xc383dc8219f6264e, 0x02cdcdb50238f12c,
    dq 0x307a4c570c944ab2, 0xf8246c44cbf10b43, 0x8dea735ca8950119,
    dq 0xc41e22a6502b9624, 0x9c742f3ad40648d1, 0x2e1568bf88056a07,
    dq 0x14089851b7e35560, 0x2770494d45aa5a86, 0x618302abcad593d2,
    dq 0xada9c22287ce2405, 0xb01689d1784d8c18, 0x522434c0a262c757,
    dq 0x4308218d32405aae, 0x60e119d9b6d2b634
global countPrimes_rm
countPrimes_rm:
    db 0x01, 0, 0, 0x03, 0, 0x07, 0x0f, 0, 0x1f, 0x3f, 0, 0x7f, 0, 0, 0xff, 0
global sieve_7_13_M7
sieve_7_13_M7:
    db 0xfd, 0xdf, 0xef, 0x7e, 0xf7, 0xfb, 0xbf, 0xfd, 0xdf, 0xef, 0x7e, 0xf7
    db 0xfb, 0xbf, 0xfd, 0xdf
    db 0xef, 0x7e, 0xf7, 0xfb, 0xbf, 0xfd, 0xdf, 0xef, 0x7e, 0xf7, 0xfb, 0xbf
    db 0xfd, 0xdf, 0xef, 0x7e
    db 0xf7, 0xfb, 0xbf, 0xfd, 0xdf, 0xef, 0x7e, 0xf7, 0xfb, 0xbf, 0xfd, 0xdf
    db 0xef, 0x7e, 0xf7, 0xfb
    db 0xbf, 0xfd, 0xdf, 0xef, 0x7e, 0xf7, 0xfb, 0xbf, 0xfd, 0xdf, 0xef, 0x7e
    db 0xf7, 0xfb, 0xbf, 0xfd
    db 0xdf, 0xef, 0x7e, 0xf7, 0xfb, 0xbf, 0xfd, 0xdf, 0xef, 0x7e, 0xf7, 0xfb
    db 0xbf, 0xfd, 0xdf, 0xef
    db 0x7e, 0xf7, 0xfb, 0xbf, 0xfd, 0xdf, 0xef, 0x7e, 0xf7, 0xfb, 0xbf, 0xfd
    db 0xdf, 0xef, 0x7e, 0xf7
    db 0xfb, 0xbf, 0xfd, 0xdf, 0xef, 0x7e, 0xf7, 0xfb, 0xbf, 0xfd, 0xdf, 0xef
    db 0x7e, 0xf7, 0xfb, 0xbf
global sieve_7_13_M11
sieve_7_13_M11:
    db 0xfb, 0xff, 0xef, 0xff, 0xbe, 0xff, 0x7d, 0xff, 0xf7, 0xff, 0xdf, 0xfb
    db 0xff, 0xef, 0xff, 0xbe
    db 0xff, 0x7d, 0xff, 0xf7, 0xff, 0xdf, 0xfb, 0xff, 0xef, 0xff, 0xbe, 0xff
    db 0x7d, 0xff, 0xf7, 0xff
    db 0xdf, 0xfb, 0xff, 0xef, 0xff, 0xbe, 0xff, 0x7d, 0xff, 0xf7, 0xff, 0xdf
    db 0xfb, 0xff, 0xef, 0xff
    db 0xbe, 0xff, 0x7d, 0xff, 0xf7, 0xff, 0xdf, 0xfb, 0xff, 0xef, 0xff, 0xbe
    db 0xff, 0x7d, 0xff, 0xf7
    db 0xff, 0xdf, 0xfb, 0xff, 0xef, 0xff, 0xbe, 0xff, 0x7d, 0xff, 0xf7, 0xff
    db 0xdf, 0xfb, 0xff, 0xef
    db 0xff, 0xbe, 0xff, 0x7d, 0xff, 0xf7, 0xff, 0xdf, 0xfb, 0xff, 0xef, 0xff
    db 0xbe, 0xff, 0x7d, 0xff
    db 0xf7, 0xff, 0xdf, 0xfb, 0xff, 0xef, 0xff, 0xbe, 0xff, 0x7d, 0xff, 0xf7
    db 0xff, 0xdf, 0xfb, 0xff
    db 0xef, 0xff, 0xbe, 0xff, 0x7d, 0xff, 0xf7, 0xff, 0xdf, 0xfb, 0xff, 0xef
    db 0xff, 0xbe, 0xff, 0x7d
    db 0xff, 0xf7, 0xff, 0xdf, 0xfb, 0xff, 0xef, 0xff, 0xbe, 0xff, 0x7d, 0xff
    db 0xf7, 0xff, 0xdf, 0xfb
    db 0xff, 0xef, 0xff, 0xbe, 0xff, 0x7d, 0xff, 0xf7, 0xff, 0xdf, 0xfb, 0xff
    db 0xef, 0xff, 0xbe, 0xff
    db 0x7d, 0xff, 0xf7, 0xff, 0xdf, 0xfb, 0xff, 0xef, 0xff, 0xbe, 0xff, 0x7d
    db 0xff, 0xf7, 0xff, 0xdf
global sieve_7_13_M13
sieve_7_13_M13:
    db 0xf7, 0xff, 0xff, 0xfe, 0xbf, 0xdf, 0xff, 0xfb, 0xfd, 0x7f, 0xff, 0xff
    db 0xef, 0xf7, 0xff, 0xff
    db 0xfe, 0xbf, 0xdf, 0xff, 0xfb, 0xfd, 0x7f, 0xff, 0xff, 0xef, 0xf7, 0xff
    db 0xff, 0xfe, 0xbf, 0xdf
    db 0xff, 0xfb, 0xfd, 0x7f, 0xff, 0xff, 0xef, 0xf7, 0xff, 0xff, 0xfe, 0xbf
    db 0xdf, 0xff, 0xfb, 0xfd
    db 0x7f, 0xff, 0xff, 0xef, 0xf7, 0xff, 0xff, 0xfe, 0xbf, 0xdf, 0xff, 0xfb
    db 0xfd, 0x7f, 0xff, 0xff
    db 0xef, 0xf7, 0xff, 0xff, 0xfe, 0xbf, 0xdf, 0xff, 0xfb, 0xfd, 0x7f, 0xff
    db 0xff, 0xef, 0xf7, 0xff
    db 0xff, 0xfe, 0xbf, 0xdf, 0xff, 0xfb, 0xfd, 0x7f, 0xff, 0xff, 0xef, 0xf7
    db 0xff, 0xff, 0xfe, 0xbf
    db 0xdf, 0xff, 0xfb, 0xfd, 0x7f, 0xff, 0xff, 0xef, 0xf7, 0xff, 0xff, 0xfe
    db 0xbf, 0xdf, 0xff, 0xfb
    db 0xfd, 0x7f, 0xff, 0xff, 0xef, 0xf7, 0xff, 0xff, 0xfe, 0xbf, 0xdf, 0xff
    db 0xfb, 0xfd, 0x7f, 0xff
    db 0xff, 0xef, 0xf7, 0xff, 0xff, 0xfe, 0xbf, 0xdf, 0xff, 0xfb, 0xfd, 0x7f
    db 0xff, 0xff, 0xef, 0xf7
    db 0xff, 0xff, 0xfe, 0xbf, 0xdf, 0xff, 0xfb, 0xfd, 0x7f, 0xff, 0xff, 0xef
    db 0xf7, 0xff, 0xff, 0xfe
    db 0xbf, 0xdf, 0xff, 0xfb, 0xfd, 0x7f, 0xff, 0xff, 0xef, 0xf7, 0xff, 0xff
    db 0xfe, 0xbf, 0xdf, 0xff
    db 0xfb, 0xfd, 0x7f, 0xff, 0xff, 0xef, 0xf7, 0xff, 0xff, 0xfe, 0xbf, 0xdf
    db 0xff, 0xfb, 0xfd, 0x7f
    db 0xff, 0xff, 0xef, 0xf7, 0xff, 0xff, 0xfe, 0xbf, 0xdf, 0xff, 0xfb, 0xfd
    db 0x7f, 0xff, 0xff, 0xef

section .bss
global L1C
L1C:
    resb 4
global BC
BC:
    resb 4

sieve_7_13.c

#include "prime.h"

#define am7(_0, _1, _2, _3, _4, _5, _6) do {\
    m0 = M7[_0]; m1 = M7[_1]; m2 = M7[_2]; m3 = M7[_3]; m4 = M7[_4];\
    m5 = M7[_5]; m6 = M7[_6];\
} while (0)

#define am11(_0, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10) do {\
    m0 = M11[_0]; m1 = M11[_1]; m2 = M11[_2]; m3 = M11[_3]; m4 = M11[_4];\
    m5 = M11[_5]; m6 = M11[_6]; m7 = M11[_7]; m8 = M11[_8]; m9 = M11[_9];\
    m10 = M11[_10];\
} while (0)

#define am13(_0, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12) do {\
    m0 = M13[_0]; m1 = M13[_1]; m2 = M13[_2]; m3 = M13[_3]; m4 = M13[_4];\
    m5 = M13[_5]; m6 = M13[_6]; m7 = M13[_7]; m8 = M13[_8]; m9 = M13[_9];\
    m10 = M13[_10]; m11 = M13[_11]; m12 = M13[_12];\
} while (0)

void sieve_7_13(unsigned char *f, uint64_t st) {
    static const void *L[] = {
        &&_7_0, &&_7_1, &&_7_2, &&_7_3, &&_7_4, &&_7_5, &&_7_6,
        &&_11_0, &&_11_1, &&_11_2, &&_11_3, &&_11_4, &&_11_5, &&_11_6, &&_11_7,
        &&_11_8, &&_11_9, &&_11_10,
        &&_13_0, &&_13_1, &&_13_2, &&_13_3, &&_13_4, &&_13_5, &&_13_6, &&_13_7,
        &&_13_8, &&_13_9, &&_13_10, &&_13_11, &&_13_12
    };
    uint64_t stb = st / 30;
    v16b *M7 = sieve_7_13_M7;
    v16b *M11 = sieve_7_13_M11;
    v16b *M13 = sieve_7_13_M13;
    v16b m0, m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12;
    v16b *vf;
    goto *L[stb % 7];
_7_0:
    am7(0, 1, 2, 3, 4, 5, 6);
    goto _7;
_7_1:
    am7(4, 5, 6, 0, 1, 2, 3);
    goto _7;
_7_2:
    am7(1, 2, 3, 4, 5, 6, 0);
    goto _7;
_7_3:
    am7(5, 6, 0, 1, 2, 3, 4);
    goto _7;
_7_4:
    am7(2, 3, 4, 5, 6, 0, 1);
    goto _7;
_7_5:
    am7(6, 0, 1, 2, 3, 4, 5);
    goto _7;
_7_6:
    am7(3, 4, 5, 6, 0, 1, 2);
_7:
    vf = (v16b *)f;
    do {
        vf[0] &= m0;
        vf[1] &= m1;
        vf[2] &= m2;
        vf[3] &= m3;
        vf[4] &= m4;
        vf[5] &= m5;
        vf[6] &= m6;
    } while ((unsigned char *)(vf += 7) < f + L1C);
    goto *(L + 7)[stb % 11];
_11_0:
    am11(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10);
    goto _11;
_11_1:
    am11(9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8);
    goto _11;
_11_2:
    am11(7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6);
    goto _11;
_11_3:
    am11(3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2);
    goto _11;
_11_4:
    am11(5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4);
    goto _11;
_11_5:
    am11(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0);
    goto _11;
_11_6:
    am11(10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9);
    goto _11;
_11_7:
    am11(8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7);
    goto _11;
_11_8:
    am11(6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5);
    goto _11;
_11_9:
    am11(4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3);
    goto _11;
_11_10:
    am11(2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1);
_11:
    vf = (v16b *)f;
    do {
        vf[0] &= m0;
        vf[1] &= m1;
        vf[2] &= m2;
        vf[3] &= m3;
        vf[4] &= m4;
        vf[5] &= m5;
        vf[6] &= m6;
        vf[7] &= m7;
        vf[8] &= m8;
        vf[9] &= m9;
        vf[10] &= m10;
    } while ((unsigned char *)(vf += 11) < f + L1C);
    goto *(L + 18)[stb % 13];
_13_0:
    am13(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
    goto _13;
_13_1:
    am13(9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8);
    goto _13;
_13_2:
    am13(5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4);
    goto _13;
_13_3:
    am13(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0);
    goto _13;
_13_4:
    am13(10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9);
    goto _13;
_13_5:
    am13(6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5);
    goto _13;
_13_6:
    am13(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1);
    goto _13;
_13_7:
    am13(11, 12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10);
    goto _13;
_13_8:
    am13(7, 8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6);
    goto _13;
_13_9:
    am13(3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2);
    goto _13;
_13_10:
    am13(12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11);
    goto _13;
_13_11:
    am13(8, 9, 10, 11, 12, 0, 1, 2, 3, 4, 5, 6, 7);
    goto _13;
_13_12:
    am13(4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2, 3);
_13:
    vf = (v16b *)f;
    do {
        vf[0] &= m0;
        vf[1] &= m1;
        vf[2] &= m2;
        vf[3] &= m3;
        vf[4] &= m4;
        vf[5] &= m5;
        vf[6] &= m6;
        vf[7] &= m7;
        vf[8] &= m8;
        vf[9] &= m9;
        vf[10] &= m10;
        vf[11] &= m11;
        vf[12] &= m12;
    } while ((unsigned char *)(vf += 13) < f + L1C);
}

sieve_17_89.c

#include "prime.h"

#define eb8(c0, c1, c2, c3, c4, c5, c6, c7, d1, d2, d3, d4, d5, d6, d7) do {\
    *f_ &= ~(1 << c0);\
    f_[q * 3 + d1] &= ~(1 << c1);\
    f_[q * 5 + d2] &= ~(1 << c2);\
    f_[q * 6 + d3] &= ~(1 << c3);\
    f_[q * 8 + d4] &= ~(1 << c4);\
    f_[q * 9 + d5] &= ~(1 << c5);\
    f_[q * 11 + d6] &= ~(1 << c6);\
    f_[q * 14 + d7] &= ~(1 << c7);\
} while (0)

#define eb(m, i, a) do {\
    *f_ &= ~(1 << i);\
    if ((f_ += q * m + a) >= f + c) return;\
} while (0)

#define eb0(i, a) eb(3, i, a)
#define eb1(i, a) eb(2, i, a)
#define eb2(i, a) eb(1, i, a)
#define eb3(i, a) eb(2, i, a)
#define eb4(i, a) eb(1, i, a)
#define eb5(i, a) eb(2, i, a)
#define eb6(i, a) eb(3, i, a)
#define eb7(i, a) eb(1, i, a)

static void sieve(unsigned char *f, uint64_t st, unsigned p) {
    unsigned c = L1C;
    uint64_t psq = sq(p);
    int64_t d = psq - st;
    unsigned char *f_ = f;
    unsigned w;
    if (__builtin_expect(psq > st, 0)) {
        f_ += (unsigned)d / 30;
        w = WIX[p % 30 / 2];
    } else {
        uint64_t p30 = (uint64_t)p * 30;
        uint64_t n = st / p30 * 30 + PADD[st % p30 / p];
        f_ += (unsigned)(p * n - st) / 30;
        w = WIX[n % 30 / 2];
    }
    unsigned q = p / 30 * 2;
    switch (p % 30) {
    case 1:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(0, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0);
                }
                eb0(0, 0);
            case 1:
                eb1(1, 0);
            case 2:
                eb2(2, 0);
            case 3:
                eb3(3, 0);
            case 4:
                eb4(4, 0);
            case 5:
                eb5(5, 0);
            case 6:
                eb6(6, 0);
            case 7:
                eb7(7, 1);
            }
        default:
            __builtin_unreachable();
        }
    case 7:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(1, 5, 4, 0, 7, 3, 2, 6, 1, 2, 3, 3, 4, 5, 6);
                }
                eb0(1, 1);
            case 1:
                eb1(5, 1);
            case 2:
                eb2(4, 1);
            case 3:
                eb3(0, 0);
            case 4:
                eb4(7, 1);
            case 5:
                eb5(3, 1);
            case 6:
                eb6(2, 1);
            case 7:
                eb7(6, 1);
            }
        default:
            __builtin_unreachable();
        }
    case 11:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(2, 4, 0, 6, 1, 7, 3, 5, 2, 4, 4, 6, 6, 8, 10);
                }
                eb0(2, 2);
            case 1:
                eb1(4, 2);
            case 2:
                eb2(0, 0);
            case 3:
                eb3(6, 2);
            case 4:
                eb4(1, 0);
            case 5:
                eb5(7, 2);
            case 6:
                eb6(3, 2);
            case 7:
                eb7(5, 1);
            }
        default:
            __builtin_unreachable();
        }
    case 13:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(3, 0, 6, 5, 2, 1, 7, 4, 3, 4, 5, 7, 8, 9, 12);
                }
                eb0(3, 3);
            case 1:
                eb1(0, 1);
            case 2:
                eb2(6, 1);
            case 3:
                eb3(5, 2);
            case 4:
                eb4(2, 1);
            case 5:
                eb5(1, 1);
            case 6:
                eb6(7, 3);
            case 7:
                eb7(4, 1);
            }
        default:
            __builtin_unreachable();
        }
    case 17:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(4, 7, 1, 2, 5, 6, 0, 3, 3, 6, 7, 9, 10, 13, 16);
                }
                eb0(4, 3);
            case 1:
                eb1(7, 3);
            case 2:
                eb2(1, 1);
            case 3:
                eb3(2, 2);
            case 4:
                eb4(5, 1);
            case 5:
                eb5(6, 3);
            case 6:
                eb6(0, 3);
            case 7:
                eb7(3, 1);
            }
        default:
            __builtin_unreachable();
        }
    case 19:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(5, 3, 7, 1, 6, 0, 4, 2, 4, 6, 8, 10, 12, 14, 18);
                }
                eb0(5, 4);
            case 1:
                eb1(3, 2);
            case 2:
                eb2(7, 2);
            case 3:
                eb3(1, 2);
            case 4:
                eb4(6, 2);
            case 5:
                eb5(0, 2);
            case 6:
                eb6(4, 4);
            case 7:
                eb7(2, 1);
            }
        default:
            __builtin_unreachable();
        }
    case 23:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(6, 2, 3, 7, 0, 4, 5, 1, 5, 8, 9, 13, 14, 17, 22);
                }
                eb0(6, 5);
            case 1:
                eb1(2, 3);
            case 2:
                eb2(3, 1);
            case 3:
                eb3(7, 4);
            case 4:
                eb4(0, 1);
            case 5:
                eb5(4, 3);
            case 6:
                eb6(5, 5);
            case 7:
                eb7(1, 1);
            }
        default:
            __builtin_unreachable();
        }
    case 29:
        switch (w) {
            for (;;) {
            case 0:
                for (; f_ + p < f + c; f_ += p) {
                    eb8(7, 6, 5, 4, 3, 2, 1, 0, 6, 10, 12, 16, 18, 22, 28);
                }
                eb0(7, 6);
            case 1:
                eb1(6, 4);
            case 2:
                eb2(5, 2);
            case 3:
                eb3(4, 4);
            case 4:
                eb4(3, 2);
            case 5:
                eb5(2, 4);
            case 6:
                eb6(1, 6);
            case 7:
                eb7(0, 1);
            }
        default:
            __builtin_unreachable();
        }
    default:
        __builtin_unreachable();
    }
}

void sieve_17(unsigned char *f, uint64_t st) {
    sieve(f, st, 17);
}

void sieve_19(unsigned char *f, uint64_t st) {
    sieve(f, st, 19);
}

void sieve_23(unsigned char *f, uint64_t st) {
    sieve(f, st, 23);
}

void sieve_29(unsigned char *f, uint64_t st) {
    sieve(f, st, 29);
}

void sieve_31(unsigned char *f, uint64_t st) {
    sieve(f, st, 31);
}

void sieve_37(unsigned char *f, uint64_t st) {
    sieve(f, st, 37);
}

void sieve_41(unsigned char *f, uint64_t st) {
    sieve(f, st, 41);
}

void sieve_43(unsigned char *f, uint64_t st) {
    sieve(f, st, 43);
}

void sieve_47(unsigned char *f, uint64_t st) {
    sieve(f, st, 47);
}

void sieve_53(unsigned char *f, uint64_t st) {
    sieve(f, st, 53);
}

void sieve_59(unsigned char *f, uint64_t st) {
    sieve(f, st, 59);
}

void sieve_61(unsigned char *f, uint64_t st) {
    sieve(f, st, 61);
}

void sieve_67(unsigned char *f, uint64_t st) {
    sieve(f, st, 67);
}

void sieve_71(unsigned char *f, uint64_t st) {
    sieve(f, st, 71);
}

void sieve_73(unsigned char *f, uint64_t st) {
    sieve(f, st, 73);
}

void sieve_79(unsigned char *f, uint64_t st) {
    sieve(f, st, 79);
}

void sieve_83(unsigned char *f, uint64_t st) {
    sieve(f, st, 83);
}

void sieve_89(unsigned char *f, uint64_t st) {
    sieve(f, st, 89);
}

sieve_93_.s

%define MD30 -2004318071
%define MD30Q -8608480567731124087

%define f rdi
%define st rsi
%define pp r8
%define nsp ecx
%define c r9
%define c_d r9d
%define fc r9
%define p r10
%define p_d r10d
%define f_ r11
%define q r12
%define fp r13
%define psq rbx
%define d rbp
%define d_d ebp
%define md30 rsp
%define md30_d esp
%define w rbx
%define w_d ebx

%macro eb8 16
    mov fp, f_
    add fp, p
    cmp fp, fc
    jae _%1_0_
    lea rax, [q + q * 2 + %10]
    lea rbx, [q + q * 4 + %11]
    lea rdx, [rax * 2 + %12]
    lea rbp, [q * 8 + %13]
    lea rsp, [q + q * 8 + %14]
    lea r14, [rax + q * 8 + %15]
    lea r15, [rdx + q * 8 + %16]
    jmp _%1_0_1
_%1_0_0:
    mov fp, f_
    add fp, p
    cmp fp, fc
    jae _%1_0_
_%1_0_1:
    and byte [f_], ~(1 << %2)
    and byte [f_ + rax], ~(1 << %3)
    and byte [f_ + rbx], ~(1 << %4)
    and byte [f_ + rdx], ~(1 << %5)
    and byte [f_ + rbp], ~(1 << %6)
    and byte [f_ + rsp], ~(1 << %7)
    and byte [f_ + r14], ~(1 << %8)
    and byte [f_ + r15], ~(1 << %9)
    mov f_, fp
    jmp _%1_0_0
%endmacro

%macro eb 3
    and byte [f_], ~(1 << %2)
    lea f_, [f_ + q * %1 + %3]
    cmp f_, fc
    jae next_
%endmacro

%macro eb_4 2
    and byte [f_], ~(1 << %1)
    lea rax, [q + q * 2 + %2]
    add f_, rax
    cmp f_, fc
    jae next_
%endmacro

%macro eb0 2
    eb_4 %1, %2
%endmacro

%macro eb1 2
    eb 2, %1, %2
%endmacro

%macro eb2 2
    eb 1, %1, %2
%endmacro

%macro eb3 2
    eb 2, %1, %2
%endmacro

%macro eb4 2
    eb 1, %1, %2
%endmacro

%macro eb5 2
    eb 2, %1, %2
%endmacro

%macro eb6 2
    eb_4 %1, %2
%endmacro

%macro eb7 2
    eb 1, %1, %2
%endmacro

extern L1C
extern BC
extern PADD
extern WIX
    align 64
L:
    dq _1, 0, 0, _7, 0, _11, _13, 0, _17, _19, 0, _23, 0, 0, _29, 0
    dq _1_0, _1_1, _1_2, _1_3, _1_4, _1_5, _1_6, _1_7
    dq _7_0, _7_1, _7_2, _7_3, _7_4, _7_5, _7_6, _7_7
    dq _11_0, _11_1, _11_2, _11_3, _11_4, _11_5, _11_6, _11_7
    dq _13_0, _13_1, _13_2, _13_3, _13_4, _13_5, _13_6, _13_7
    dq _17_0, _17_1, _17_2, _17_3, _17_4, _17_5, _17_6, _17_7
    dq _19_0, _19_1, _19_2, _19_3, _19_4, _19_5, _19_6, _19_7
    dq _23_0, _23_1, _23_2, _23_3, _23_4, _23_5, _23_6, _23_7
    dq _29_0, _29_1, _29_2, _29_3, _29_4, _29_5, _29_6, _29_7
    align 16
global sieve_93_
sieve_93_:
    push rbx
    push rbp
    push r12
    push r13
    push r14
    push r15
    movq xmm0, rsp
    mov r8, rdx
    mov c_d, [rel L1C]
    mov eax, [rel BC]
    cmp nsp, -1
    cmove c_d, eax
    mov p_d, [rel pp]
    mov psq, p
    imul psq, p
    mov d, psq
    sub d, st
L0:
    imul rax, c, 30
    cmp d, rax
    jge end
    mov f_, f
    mov md30_d, MD30
    cmp psq, st
    ja L00
    mov rax, st
    xor rdx, rdx
    imul rbx, p, 30
    div rbx
    imul rax, 30
    mov rbx, rax
    mov eax, edx
    shr rdx, 32
    div p_d
    movzx eax, byte [rax + PADD]
    add rax, rbx
    mov rbx, rax
    imul rax, p
    sub rax, st
    imul rax, md30
    shr rax, 36
    cmp eax, c_d
    jae next
    add f_, rax
    mov rax, MD30Q
    mul rbx
    shr rdx, 4
    imul rdx, -30
    add rdx, rbx
    shr edx, 1
    movzx w_d, byte [rdx + WIX]
    mov eax, p_d
    imul rax, md30
    shr rax, 36
    imul edx, eax, -30
    add edx, p_d
    shr edx, 1
L1:
    shl eax, 1
    mov q, rax
    add c, f
    jmp [rdx * 8 + L]
    align 16
_1:
    jmp [w * 8 + L + 16 * 8]
    align 16
_1_0:
    eb8 1, 0, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0
_1_0_:
    eb0 0, 0
_1_1:
    eb1 1, 0
_1_2:
    eb2 2, 0
_1_3:
    eb3 3, 0
_1_4:
    eb4 4, 0
_1_5:
    eb5 5, 0
_1_6:
    eb6 6, 0
_1_7:
    eb7 7, 1
    jmp _1_0
    align 16
_7:
    jmp [w * 8 + L + 24 * 8]
    align 16
_7_0:
    eb8 7, 1, 5, 4, 0, 7, 3, 2, 6, 1, 2, 1, 3, 4, 4, 3
_7_0_:
    eb0 1, 1
_7_1:
    eb1 5, 1
_7_2:
    eb2 4, 1
_7_3:
    eb3 0, 0
_7_4:
    eb4 7, 1
_7_5:
    eb5 3, 1
_7_6:
    eb6 2, 1
_7_7:
    eb7 6, 1
    jmp _7_0
    align 16
_11:
    jmp [w * 8 + L + 32 * 8]
    align 16
_11_0:
    eb8 11, 2, 4, 0, 6, 1, 7, 3, 5, 2, 4, 0, 6, 6, 6, 6
_11_0_:
    eb0 2, 2
_11_1:
    eb1 4, 2
_11_2:
    eb2 0, 0
_11_3:
    eb3 6, 2
_11_4:
    eb4 1, 0
_11_5:
    eb5 7, 2
_11_6:
    eb6 3, 2
_11_7:
    eb7 5, 1
    jmp _11_0
    align 16
_13:
    jmp [w * 8 + L + 40 * 8]
    align 16
_13_0:
    eb8 13, 3, 0, 6, 5, 2, 1, 7, 4, 3, 4, -1, 7, 8, 6, 7
_13_0_:
    eb0 3, 3
_13_1:
    eb1 0, 1
_13_2:
    eb2 6, 1
_13_3:
    eb3 5, 2
_13_4:
    eb4 2, 1
_13_5:
    eb5 1, 1
_13_6:
    eb6 7, 3
_13_7:
    eb7 4, 1
    jmp _13_0
    align 16
_17:
    jmp [w * 8 + L + 48 * 8]
    align 16
_17_0:
    eb8 17, 4, 7, 1, 2, 5, 6, 0, 3, 3, 6, 1, 9, 10, 10, 9
_17_0_:
    eb0 4, 3
_17_1:
    eb1 7, 3
_17_2:
    eb2 1, 1
_17_3:
    eb3 2, 2
_17_4:
    eb4 5, 1
_17_5:
    eb5 6, 3
_17_6:
    eb6 0, 3
_17_7:
    eb7 3, 1
    jmp _17_0
    align 16
_19:
    jmp [w * 8 + L + 56 * 8]
    align 16
_19_0:
    eb8 19, 5, 3, 7, 1, 6, 0, 4, 2, 4, 6, 0, 10, 12, 10, 10
_19_0_:
    eb0 5, 4
_19_1:
    eb1 3, 2
_19_2:
    eb2 7, 2
_19_3:
    eb3 1, 2
_19_4:
    eb4 6, 2
_19_5:
    eb5 0, 2
_19_6:
    eb6 4, 4
_19_7:
    eb7 2, 1
    jmp _19_0
    align 16
_23:
    jmp [w * 8 + L + 64 * 8]
    align 16
_23_0:
    eb8 23, 6, 2, 3, 7, 0 ,4, 5, 1, 5, 8, -1, 13, 14, 12, 13
_23_0_:
    eb0 6, 5
_23_1:
    eb1 2, 3
_23_2:
    eb2 3, 1
_23_3:
    eb3 7, 4
_23_4:
    eb4 0, 1
_23_5:
    eb5 4, 3
_23_6:
    eb6 5, 5
_23_7:
    eb7 1, 1
    jmp _23_0
    align 16
_29:
    jmp [w * 8 + L + 72 * 8]
    align 16
_29_0:
    eb8 29, 7, 6, 5, 4, 3, 2, 1, 0, 6, 10, 0, 16, 18, 16, 16
_29_0_:
    eb0 7, 6
_29_1:
    eb1 6, 4
_29_2:
    eb2 5, 2
_29_3:
    eb3 4, 4
_29_4:
    eb4 3, 2
_29_5:
    eb5 2, 4
_29_6:
    eb6 1, 6
_29_7:
    eb7 0, 1
    jmp _29_0
    align 16
next_:
    sub c, f
next:
    dec nsp
    jz end
    add pp, 4
    mov p_d, [rel pp]
    mov psq, p
    imul psq, p
    mov d, psq
    sub d, st
    jmp L0
    align 16
L00:
    mov eax, d_d
    imul rax, md30
    shr rax, 36
    add f_, rax
    mov eax, p_d
    imul rax, md30
    shr rax, 36
    imul edx, eax, -30
    add edx, p_d
    shr edx, 1
    movzx w_d, byte [rdx + WIX]
    jmp L1
    align 16
end:
    movq rsp, xmm0
    pop r15
    pop r14
    pop r13
    pop r12
    pop rbp
    pop rbx
    ret

build.sh

#!/bin/sh -x

SRC="main.c sieve_7_13.o sieve_17_89.o sieve_93_.o lut.o -lm"
C="gcc"
O="-O3 -march=native"
$C $O -c sieve_7_13.c
$C $O -c sieve_17_89.c
nasm -felf64 sieve_93_.s
nasm -felf64 lut.s
$C $O -or $SRC
rm *.o

Answer 3

Jyxal 0.4.1, 46 MiB on my computer, official score 120 MiB

{500(&›¥æߥ)W⁋₴

Read the README to see how to run this

For reference, the example program outputted 21.8 MiB on my computer.

I had a blast micro-optimizing the compiler. The old Jyxal answer was disqualified because it considered multiples of 5 prime. The new version of Jyxal is almost completely rewritten in Koltin (barring the math library) and contains a new post-compilation optimizer (powered by Proguard™). This allows Jyxal to move to Java 11. It also optimizes for loops. The needless push-pops are optimized away and the JVM operand stack is used whenever possible.

{500(&›¥æߥ)W⁋₴ # Takes no input
{               # Infinite loop
 500(      )    # Repeat 500 times
     &›         # Increment the register
       ¥æ       # And check it for primality
         ß      # If it is prime...
          ¥     # Push the number
            W   # Wrap the stack in a list
             ⁋  # Join by newlines
              ₴ # And then print the result
                # Implicitly loops back

The primality test is nothing as elaborate as Vyxal's. For numbers that are smaller than 9,223,372,036,854,775,807, the BigComplex object is converted into a long, and then the last bit is tested. If it is 0, that means the number is even and that the number is not prime. Then using an optimization I found on SO, I use a brute-force division algorithm that skips multiples of 2 and 3, lessening the number of divisors I need to check.

Compiled with

java -jar Jyxal.v.0.4.1.jar code.txt t

The t flag disables vectorisation of monads, which removes the need for MethodHandles and allows the monad to be inlined.

If you are interested, here is the bytecode decompiled into Java:

public final class Main {
    private static Object register;

    static {
        register = BigComplex.ZERO;
    }

    public static void main(String[] var0) {
        ProgramStack stack = new ProgramStack(var0);

        while(true) {
            for(int var1 = 500; var1 != 0; --var1) {
                if (RuntimeHelpers.truthValue(RuntimeMethods.isPrime(register = RuntimeMethods.increment(register)))) {
                    stack.push(register);
                }
            }

            System.out.print(RuntimeMethods.joinByNewlines(JyxalList.create(stack)));
        }
    }
}

Answer 4

C (gcc -O3)

#include <stdio.h>
#include <string.h>
#define MAX_SQRT 100000000ULL
unsigned long composites[MAX_SQRT >> 7UL] = { 1 };
unsigned long primes[10000000] = { 2 };
int main(int argc, char **argv) {
  unsigned long *last = primes;
  printf("2\n");
  for (unsigned long i = 3; i < MAX_SQRT; i += 2) {
    if (composites[i >> 7UL] & 1UL << (i >> 1UL)) continue;
    printf("%lu\n", i);
    *++last = i;
    for (unsigned long j = i * i >> 1UL; j < MAX_SQRT >> 1UL; j += i) composites[j >> 6UL] |= 1UL << j;
  }
  for (unsigned long long i = MAX_SQRT + 1; i < MAX_SQRT * MAX_SQRT; i += MAX_SQRT) {
    memset(composites, 0, sizeof(composites));
    for (unsigned long *j = primes; ++j <= last; ) {
      unsigned long k = (unsigned long)(*j - 1 - ((*j - 1 + i) % (*j * 2) >> 1UL));
      while (k < MAX_SQRT >> 1UL) {
        composites[k >> 6UL] |= 1UL << k;
        k += *j;
      }
    }
    for (unsigned long j = 0; j < MAX_SQRT >> 1UL; j++) if (!(composites[j >> 6UL] & 1UL << j)) printf("%llu\n", i + j * 2);
  }
}

Outputs 6.2GB in 1 minute on my fast PC, for a speed of 105 MB/s. Note that I timed using (timeout 1m ./a.out) | pv > /dev/null, which avoids timing out pv itself.

Answer 5

C++

#include <stdio.h>
#include <array>

// ~~ Algorithm Description ~~
// Since sqrt(max) = 100000000 = M, we should precompute primes until M,
// because all factors of N are always smaller or equal to sqrt(N).
// There are 5761455 primes smaller than M.

// After these maneuvers, assuming that `ull' is 8 bytes large, the binary
// will grow by approximately 46 megabytes.

// Notice that this solution largely makes use of the C++ compile-time features
// and the fact that the program's compile-time isn't timed nor the binary size
// isn't measured.

// BUILDING:
// g++ -O3 -std=c++20 pgen.cpp -o pgen -fconstexpr-ops-limit=1000000000000 -fconstexpr-loop-limit=1000000000

typedef unsigned long long ull;

template <std::size_t N>
constexpr std::array<ull, N> get_ptab() {
    std::array<ull, N> ptab;
    ull gen = 3;
    ptab[0] = 2; ptab[1] = 3; ptab[2] = 5;
    for(std::size_t i = 7; gen < N; i += 2) {
        bool ok = true;
        if(i % 3 == 0 || i % 5 == 0)
            continue;
        for(std::size_t j = 2; j * j <= i; j++) {
            if(i % j == 0) {
                ok = false;
                break;
            }
        }
        if(ok)
            ptab[gen++] = i;
    }
    return ptab;
}

constexpr std::array<ull, 5761455> ptab = get_ptab<5761455>();

int main(void) {
    // Step 1: print all the primes that we've hardcoded already.
    for(std::size_t i = 0; i < ptab.size(); i++)
        printf("%llu\n", ptab[i]);
    // Step 2: starting with ptab.back(), go up until max, checking if
    // N is divisible by any of the primes we've hardcoded.
    // If it is, skip it. If it isn't, print it.
    for(std::size_t i = ptab.back(); i < 10000000000000000; i+=2) {
        bool ok = true;
        for(std::size_t j = 0; j < ptab.size() && i <= ptab[j]; j++) {
            if(i % ptab[j] == 0) {
                ok = false;
                break;
            }
        }
        if(ok)
            printf("%llu\n", i);
    }
    return 0;
}

The performance is yet unknown since the code still compiles on my machine. If you manage to compile, run and time it, let me know how quickly it performs.

Answer 6

PARI/GP, 929MiB on my computer

forprime(p=2, 10^16, print(p))

The package name for PARI/GP is pari-gp on Ubuntu, Debian, Arch Linux and some other distros.

Runs as gp -qf ./file_name.gp.

---

PARI/GP + gp2c, 3.76GiB on my computer

gp2c compiles PARI/GP functions to C codes, where the generated .so file can be called in another PARI/GP script.

So there are two files:

a.gp:

a() = forprime(p=2, 10^16, print(p))

b.gp:

install("a","v","a","./a.gp.so");
a()

How to run

The package name for gp2c is pari-gp2c on Ubuntu and Debian, gp2c on Arch Linux (AUR).

First compile a.gp with gp2c:

gp2c a.gp > a.gp.c

The command to compile the generated C code is in the first comment of the generated C code. It should look like:

cc -c -o a.gp.o -g -O3 -Wall -fomit-frame-pointer -fno-strict-aliasing -fPIC -I"/usr/include/x86_64-linux-gnu" a.gp.c && cc -o a.gp.so -shared -g -O3 -Wall -fomit-frame-pointer -fno-strict-aliasing -fPIC -Wl,-shared -Wl,-z,relro a.gp.o -lc -lm -L/usr/lib/x86_64-linux-gnu -lpari

After that you can run b.gp and measure the throughput:

(timeout 1m gp -qf ./b.gp) | pv > /dev/null

---

C + PARI/GP's C library, 4.81GiB on my computer

Modified from the C code generated by gp2c.

#include <pari/pari.h>

int main()
{
    pari_init(8000000, 500000);

    forprime_t iter;
    u_forprime_init(&iter, 2, 10000000000000000);
    long p;
    while (p = u_forprime_next(&iter))
        printf("%ld\n", p);
}

If you are using Ubuntu or Debian, you need to install the package libpari-dev. For Arch Linux, pari-gp is enough.

Compiles with gcc -O3 -lpari.