# High throughput prime numbers

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.

Contestant Language Score
alephalpha C 6.27GiB
Neil C 4.80GiB
(example) C 49.6MiB
emanresu A Vyxal 21.8MiB
• I'm not really sure what new this adds to the site. It seems awfully similar to the challenge it was inspired by, borderline duplicate. The most significant changes such as the bans on builtins and multithreading seem to just make the challenge worse. Jan 7 at 11:41
• @WheatWizard I spent a lot of time doing the Fizz Buzz challenge. The amount of data produced by an optimized program in a second is immense (at least 10GB per second). Therefore, IO by itself is a critical part of optimization to make it not create bottlenecks. However, producing primes is a computationally very intensive task compared to Fizz Buzz. Reaching the throughput of 100MB per second by the best optimized program on best hardware will be impossible, I believe. This task is about writing an optimized algorithm for producing prime numbers. The IO overhead will be insignificant. Jan 7 at 12:00
• Either this comes down to "Do output fast" which is a challenge we already have, or "Do primes fast", which is an active area of academic research and really beyond the scope of this site. I don't think either makes for a good challenge. Jan 7 at 12:16
• @WheatWizard Beyond the scope of this site because it "is an active area of academic research"? I can come up with several counterexamples, but after all, let's see if people show interest and participate in this challenge, which is what I think this site is about. Jan 7 at 12:25
• @xiver77 The idea that “producing primes is a computationally very intensive task” is mistaken. The primes up to $n$ can be sieved in $O(n)$ time. Output will always be the slow part. For example, primesieve computes all the primes up to $10^{10}$ in about 300 milliseconds, but takes half a minute to print them out (which, by the way, it does at over 200 MB/s). Jan 8 at 1:52

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

## ~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

• You got 17.9Mib on my machine. Jan 10 at 13:02
• @xiver77 Can you test this again please, with the version of Vyxal downloadable from github? A recent update changes the behaviour of the stack, which should make this a lot faster. 2 days ago
• "stack-based" ... "which makes it incredibly difficult to optimize anything": It is kinda false. Factor is stack-based but can be JIT'ed and therefore is pretty fast among interpreted languages. 2 days ago
• It has indeed got faster (from 17.9 to 21.8MiB). Last time I think there was a release download link in github, so I downloaded the full release and ran Vyxal.py directly. This time they provided a pip installer, so I used that. The current version I used is 2.7.5. yesterday

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

• Excessive memory use isn't caused by computing primes, rather, by compiler optimisations. For me, due to compiler optimisations, seconds before outputting the binary the compiler consumes 7236M of RAM. Jan 8 at 13:06
• Maybe gcc's assembly output from gcc -S is recompilable? I'll accept that if it works. I cannot compile your program in any way since my machine only has 4GB RAM. Jan 8 at 13:41

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

• Got 4.8 Gig on my machine. Thanks for your suggestion of using parentheses in the command. Jan 9 at 3:56
• Could you explain the method you are using please. Jan 9 at 15:52
• @graffe It's just a prime sieve. The first loop computes the primes up to 100,000,000 (printing them as it goes), and then the second loop sieves the integers in batches of 100,000,000, printing the primes it finds at the end of each batch.
– Neil
Jan 9 at 16:29

# 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`.

• Wow! This is an amazing tutorial on how to speed up pari code Jan 10 at 10:04