Sieb des Eratosthenes

[Bonita Montero]

Ich hab gestern mal jemanden im IRC erzählt, dass ich mir einen AMD
16-Kerner gegönnt habe. Der meinte darauf, dass ich dessen Performance
ja mal testen könnte indem ich das Sieb des Eratosthenes als Bash-Script
laufen lasse. Ich darauf, dass man dann doch mehr die Perfomrmance des
Interpreters teste als das eigentliche Problem zu benchmarken.
Ich mein solche Benchmarks sind schon unsinnig genug, da wollte ich der
Unsinnigkeit noch eins drauf setzen und habe das Sieb des Eratosthenes
als C++20-Programm geschrieben, und zwar wirklich so effizient wie mög-
lich. Letztlich bearbeite ich mehrfach einen Vektor aus vorab auf true
bools und lösche die Vielfachen von bisher bekannten Primzahlen und su-
che darauf die nächste Primzahl und fange wieder mit deren Vielfachen
an.
In C++ ist ein bool üblicherweise ein Byte das einfach auf Null oder
Eins gesetzt ist. Dementsprechend wäre ein vector<bool> ein Vektor aus
Bytes. In C++ lassen sich aber einzelne gernerische Klassen für Sonder-
fälle, bezogen auf bestimmte Typen, spezialisieren, und das ist auch
für den vector<bool> so, der dann jedes bool einfach als ein Bit spei-
chert. Das macht, dass der Vektor eben kompakter gespeichert wird und
das wirkt sich dann eben positiv auf die Performance aus.
Am Ende werden die Primzahlen dann alle noch in eine Datei ausgegeben.
Ich mein das könnte ich in C wie in C++ über irgendwelche ASCII-Streams
tun, aber die sind meiner Erfahung nach alle ziemlich ineffizient.
Es gibt seit C++20 die Funktion to_chars, die eine Zahl ähnlich wie bei
sprintf in einen Stück Speicher ausgeben kann. Bei der gehe ich davon
aus, dass die wirklich maximal effizient arbeitet (vor allem muss hier
ein Format-String geparst werden). Damit printe ich die Primzahlen fort-
laufend in ein char-Array, hänge ein Newline an und hänge das Ganze dann
an einen dynamisch wachsenden C++-string an.
Die C++-Strings haben wie die Vektoren die Eigenschaft, dass der allo-
kierte Speicher nicht für jedes Größenwachstum neu um-allokiert wird,
sondern nur dann wenn dessen Kapazitäts-Grenzen gebrochen werden. C++
hat für Vektoren und Strings nämlich die Eigenschaft des "amortisiert
konstanten" Overheads für jedes Einfügen oder Anhängen an einen String
oder Vektor, d.h. der Aufwand für jedes Anfügen ist konstant, egal
welche Größe der Vektor zuvor hatte. Um das umsetzen zu können muss
der Container in seiner Kapazität immer um einen konstanten Faktor
wachsen. Bei Visual C++ wächst der Container in seiner Kapazität dann
immer um 50% seiner vorherigen Größe, bei libstdc++ (g++) und libc++
(clang++) wächst er immer um 100% (is performanter da seltener reallo-
kiert wird als bei MSVC, braucht aber mehr Speicher).
Ich hatte anfänglich damit gerechnet, dass die Reallokationen den Code
maßgeblich verlangsamen. Daher habe ich die Ausgabeschleife als Funk-
tion in der Funktion, also als so genanntes Lambda, geschrieben und
diese Funktion kriegt an zwei Stellen ein unterschiedliches Funktions-
objekt übergeben, das an erster Stelle nur die Länge der Zeilen mit
den Primzahlen aufsummiert und an zweiter Stelle das dann in einen
entsprechend vor-allokierten String kopiert. Da das Lambda für jedes
dieser Übergebenen Funktionsobjekte zweimal kompiliert wird habe ich
so wirklich die maximale Performance beim Durchzählen und bei der Aus-
gabe.
Interessanterweise haut das Durchzählen bzgl. der CPU-Zeit mehr rein
als das mit logarithmisch zur Größe des Ausgabe-Strings häufigen Re-
allokieren und Umkopieren des Strings. Daher hab ich das Durchzählen
zuletzt auskommentiert.
Wenn ich alle Primzahlen im Zahlenbereich von >= 2 bis < 0x1000000000
durchrechnen lasse braucht das Ganze dann nur ca 10s auf meinem AMD
7950X. Was ich auch interessant finde ist, dass das Generieren der
Datei aus dem berechneten Array daran nur ca. 0,5s Anteil hat.


#include <iostream>
#include <vector>
#include <charconv>
#include <string_view>
#include <algorithm>
#include <fstream>
#include <cctype>

using namespace std;

int main( int argc, char **argv )
{
constexpr bool USE_CHAR = true;
using bool_t = conditional_t<!USE_CHAR, bool, char>;
size_t max = 10000;
if( argc >= 2 )
{
char const* sizeParam = argv[1];
bool hex = argv[1][0] == '0' && tolower( argv[1][1] ) == 'x';
sizeParam += 2 * hex;
if( from_chars_result fcr = from_chars( sizeParam, sizeParam + strlen(
sizeParam ), max ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;
}
vector<bool_t> sieve( max, true );
for( size_t m = 2; m < max; )
{
for( size_t i = 2 * m; i < max; i += m )
sieve[i] = false;
do
if( ++m == max )
goto finished;
while( !sieve[m] );
}
finished:
auto format = [&]( auto fn )
{
for( size_t i = 2; i != max; ++i )
if( sieve[i] )
{
char append[32];
to_chars_result tcr = to_chars( append, append + sizeof append, i );
*tcr.ptr++ = '\n';
fn( string_view( append, tcr.ptr ) );
}
};
size_t outLength = 0;
//format( [&]( string_view const &append ) { outLength +=
append.length(); } );
string print;
print.reserve( outLength );
format( [&]( string_view const &append ) { print.append( append ); } );
ofstream ofs;
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc );
ofs << print << endl;
}

Mich würde an der Stelle der Performance-Vergleich zu anderen
Implementationen interessieren. Freiwillige vor.

[Bonita Montero]

[Wendelin Uez]

> In C++ ist ein bool üblicherweise ein Byte das einfach auf Null oder
> Eins gesetzt ist.

Wäre denn eine zwei dann "weiß nicht" und eine drei "vielleicht"?
Wie schrieb schon Johann Schüttelspeer: "tubii ohr notubii".

[Bonita Montero]

Es könnten das ja auch völlig andere Zahlen sein. Du kannst
aber in C++ ein Plattform-abhängig illegales bool machen;

void fn( bool &b )
{
char c = 123;
memcpy( &b, &c, 1 );
}

(memcpy() ist die einzig legale Art in C++ Aliasing zu machen).

[Hermann Riemann]

Am 12.08.23 um 12:22 schrieb Wendelin Uez:

<https://de.wikipedia.org/wiki/Fuzzylogik>

[Bonita Montero]

Am 12.08.2023 um 14:03 schrieb Hermann Riemann:
> Am 12.08.23 um 12:22 schrieb Wendelin Uez:

> <https://de.wikipedia.org/wiki/Fuzzylogik>

Das wäre ja nicht mehr fuzzy wenn es zwei andere exakte
Werte für true und false geben würde.

[Bonita Montero]

Ich hab das gerade mal zum Euler-Sieb umgebaut was ich da gebastelt
habe. Das Euler-Sieb hält noch eine Liste der bisherigen Primzahlen
und guckt erst, ob das aktuelle Vielfache was da auf false gesetzt
werden soll nicht schon ein Vielfaches von anderen bisherigen Prim-
zahlen ist. Das kam mir anfangs noch recht aufwendig vor, aber ich
dachte mit, dass man schnell bei kleinen bisherigen Primzahlen hän-
gen bleibt da es sehr viele von deren Vielfachen gibt. Letztlich
ist der Code dann mit dieser "Optimierung" dann doch ein Vielfaches
langsamer.

#include <iostream>
#include <vector>
#include <charconv>
#include <string_view>
#include <algorithm>
#include <fstream>
#include <cctype>

#include <cstring>

#if defined(_MSC_VER)
#pragma warning(disable: 4996)
#endif

using namespace std;

int main( int argc, char **argv )
{

constexpr bool USE_BOOL = true;
constexpr size_t BUF_SIZE = 0x100000;
try
{
size_t max;
if( argc < 2 )
return EXIT_FAILURE;
char const *sizeStr = argv[1];

bool hex = argv[1][0] == '0' && tolower( argv[1][1] ) == 'x';

sizeStr += 2 * hex;
if( from_chars_result fcr = from_chars( sizeStr, sizeStr + strlen(
sizeStr ), max, !hex ? 10 : 16 ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;
using bool_t = conditional_t<USE_BOOL, bool, char>;

vector<bool_t> sieve( max, true );

vector<size_t> soFar;

for( size_t m = 2; m < max; )
{
for( size_t i = 2 * m; i < max; i += m )

{
bool multiple = false;
for( size_t sf : soFar )
if( !(i % sf) )
{
multiple = true;
break;
}
if( !multiple )
sieve[i] = false;
}
soFar.emplace_back( m );
while( ++m < max && !sieve[m] );
}
if( argc >= 3 && !*argv[2] )
return EXIT_SUCCESS;

ofstream ofs;
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc );

string buf;
buf.reserve( BUF_SIZE + 31 );

for( size_t i = 2; i != max; ++i )
if( sieve[i] )
{

char primeStr[32];
to_chars_result tcr = to_chars( primeStr, end( primeStr ), i );
*tcr.ptr = '\n';
buf.append( string_view( primeStr, tcr.ptr + 1 ) );
if( buf.size() >= BUF_SIZE )
ofs << buf,
buf.resize( 0 );
}
ofs << buf;
}
catch( bad_alloc const & )
{
cout << "out of memory" << endl;
}
catch( ios_base::failure const & )
{
cout << "I/O error" << endl;
}
}

[Bonita Montero]

Jemand aus comp.lang.c++ hat eine andere Idee für das Euler'sche Sieb
vorgeschlagen, und zwar, dass man das jeweils die nächste rechts von
der aktuellen Primzahl stehende nächsten Vielfache einer bisherigen
Primzahl in den obersten Knoten einer Heap-Datenstruktur packt.
Dann zählt man eben alle Zahlen bis zum nächsten Vielfachen durch und
"die" sind dann eben eine Primzahl. Natürlich kann es mehrere Vielfache
geben die auf die selbe Zahl fallen, dass man den Heap beim Durchlaufen
fortlaufend umorganisieren muss. D.h. wenn man irgendein Vielfaches
getroffen hat dann addiert man eben die Primzahl zu dem das ein Viel-
faches ist drauf und sortiert das im Heap neu ein. Das Umsortieren im
Heap hat dann logarithmische Laufzeit.
So wie ich das mache, dass ich zu der Primzahl eben keinen Faktor,
sondern eben das nächste Vielfache, speichere brauch ich dann also
nichtmal eine Multiplikation, geschweige denn wie im vorangegangenen
Code eine teure Division (so 20 Takte auf modernen CPUs).
Da es im Heap dann wie gesagt mehrere Vielfache unterschiedlicher
Primzahlen die gleich sind gibt muss bei diesen Überlappungen der
Heap dann mehrfach umorganisiert werden.
Ich mein das Ganze ist dann zwar super-elegant, aber eben auch nicht
so effizient wie der einfache Sieb des Eratosthenes per Bitmap. Wenn
ich alle Primzahlen von 2 bis <= 1E8 durchgehe hab ich auf meinem
Zen4-System einen Performance-Unterschied von Faktor 82.

#include <iostream>
#include <queue>
#include <vector>
#include <charconv>
#include <fstream>

using namespace std;

int main( int argc, char **argv )
{

constexpr size_t BUF_SIZE = 0x100000;

size_t max;
if( argc < 2 )
return EXIT_FAILURE;
char const *sizeStr = argv[1];

bool hex = sizeStr[0] == '0' && (sizeStr[1] == 'x' || sizeStr[1] == 'X');

sizeStr += 2 * hex;
if( from_chars_result fcr = from_chars( sizeStr, sizeStr + strlen(
sizeStr ), max, !hex ? 10 : 16 ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;

struct next_multiple { size_t p, pNext; };
auto nm_less = []( next_multiple const &left, next_multiple const
&right ) { return left.pNext > right.pNext; };
priority_queue<next_multiple, vector<next_multiple>, decltype(nm_less)>
multiples( nm_less );
multiples.emplace( 2, 4 );
auto *nextMultiple = &multiples.top();
for( size_t p = 3; p <= max; ++p )
{
if( p < nextMultiple->pNext ) [[unlikely]]
{
multiples.emplace( p, p + p );
nextMultiple = &multiples.top();
continue;
}
do
{
auto nextIncr = *nextMultiple;
multiples.pop();
nextIncr.pNext += nextIncr.p;
multiples.emplace( nextIncr );
nextMultiple = &multiples.top();
} while( p == nextMultiple->pNext );
}
vector<size_t> primes;
primes.reserve( multiples.size() );
while( !multiples.empty() )
{
auto const next = multiples.top();
primes.emplace_back( next.p );
multiples.pop();
}
sort( primes.begin(), primes.end() );

if( argc >= 3 && !*argv[2] )

return EXIT_FAILURE;

ofstream ofs;
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc );

union ndi_t { }; // non-default-initialzied
vector<ndi_t> buf( BUF_SIZE + 32 );
char
*const bufBegin = (char *)buf.data(),
*bufEnd = bufBegin,
*const bufFlush = bufBegin + BUF_SIZE;
for( size_t p : primes )
{
to_chars_result tcr = to_chars( bufEnd, bufEnd + 32, p );
bufEnd = tcr.ptr;
*bufEnd++ = '\n';
if( bufEnd >= bufFlush )
ofs << string_view( bufBegin, bufEnd ),
bufEnd = bufBegin;
}
ofs << string_view( bufBegin, bufEnd );
}

[Bonita Montero]

Ich hab mir jetzt mal gedacht, dass man die Bitmap-Lösung und die mit
dem binären Heap kombinieren könnte. Mein Grundgedanke war einfach,
dass wenn die Bitmap so klein zu machen, dass die in den Cache geht.
Die Bitmap erfasst dann jeweils nur ein Fenster des Zahlenstrahls
und in dem Heap sind alle bisherigen Primzahlen und deren nächstes
Vielfaches rechts von der nächsten Bitmap geführt, "sortiert" nach
dem nächsten Vielfachen.

#include <iostream>
#include <vector>
#include <charconv>

#include <string_view>
#include <algorithm>
#include <fstream>
#include <cctype>
#include <cstring>

#include <queue>

#if defined(_MSC_VER)
#pragma warning(disable: 26495) // uninitialized member
#endif

using namespace std;

int main( int argc, char **argv )
{

try
{

if( argc < 2 )
return EXIT_FAILURE;

auto parseParam = [&]( char const *str ) -> size_t
{
bool hex = str[0] == '0' && (str[1] == 'x' || str[1] == 'X');
str += 2 * hex;
size_t param;
if( from_chars_result fcr = from_chars( str, str + strlen( str ),
param, !hex ? 10 : 16 ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;
return param;
};
size_t
max = parseParam( argv[1] ),
bits = argc >= 3 ? parseParam( argv[2] ) : 512ull * 1024 * 1024 * 8;
union ndi_t { uint8_t c; ndi_t() {} };
vector<ndi_t> bitmap( (bits + 7) / 8 );

struct next_multiple { size_t p, pNext; };
auto nm_less = []( next_multiple const &left, next_multiple const
&right ) { return left.pNext > right.pNext; };
priority_queue<next_multiple, vector<next_multiple>, decltype(nm_less)>
multiples( nm_less );

for( size_t w = 2; w <= max; w += bits )
{
for_each( bitmap.begin(), bitmap.end(), [&]( ndi_t &n ) { n.c = 0xFF;
} );
size_t wEnd = w + bits <= max ? w + bits : max;
auto clearInterval = [&]( next_multiple &nm )
{
size_t b = nm.pNext - w;
for( ; b < wEnd - w; b += nm.p )
bitmap[b / 8].c &= ~(1u << b % 8);
nm.pNext = w + b;
};
if( multiples.size() ) [[likely]]
for( next_multiple const *pNm; (pNm = &multiples.top())->pNext < wEnd; )
{
auto nm = *pNm;
multiples.pop();
clearInterval( nm );
multiples.push( nm );
}
for( size_t b = 0, e = wEnd - w; b != e; ++b )
if( bitmap[b / 8].c >> b % 8 & 1 )
{
next_multiple nm { w + b, 2 * (w + b) };
clearInterval( nm );
multiples.push( nm );
}
}
/*while (multiples.size())
cout << multiples.top().p << endl,
multiples.pop();*/

}
catch( bad_alloc const & )
{
cout << "out of memory" << endl;
}
catch( ios_base::failure const & )
{
cout << "I/O error" << endl;
}
}

Obiges Programm nimmt als ersten Parameter den Maximalwert und als
zweiten die Anzahl der Bits in der Bitmap, ggf. mit einem 0x davor,
dass man die Angabe hexadezimal machen kann.
Auf jeden Fall ist meine Idee nicht aufgegangen und der Algorithmus
ist weiterhin langsamer als die reine Bitmap-Lösung, auch wenn die
Bitmap so groß ist, dass die in in keinen Cache passt.
Eine Ausgabe in eine Datei hab ich dann gar nicht erst programmiert
weil diese Lösung dann ja eh wegfällt.

[Bonita Montero]

Ich hab mir mal das Sieb heute morgen auf alle Kerne parallelisiert;
die meiste CPU-Zeit geht ja wirklich auf das resetten der Bits jen-
seits dem Quadrat der aktuell gefundenen Primzahl.
Der Sppedup liegt auf einem Rechner mit 16 Zen4-Kernen aber nur bei
knap Faktor 2.5. Ich kanns nicht sicher sagen, aber ich denk mal das
hängt letztlich an dem Durchsatz des Speicher-Interfaces, dass eben
irgendwann gesättigt ist.

#include <cstdlib>
#include <charconv>
#include <cstring>
#include <vector>
#include <algorithm>
#include <cmath>
#include <bit>
#include <fstream>
#include <string_view>
#include <thread>
#include <mutex>
#include <condition_variable>
#include <utility>
#include <semaphore>
#include <atomic>
#include <new>
#include <span>
#include <deque>

using namespace std;

#if defined(__cpp_lib_hardware_interference_size)
constexpr ptrdiff_t CL_SIZE = hardware_destructive_interference_size;
#else
constexpr ptrdiff_t CL_SIZE = 64;
#endif

int main( int argc, char **argv )
{

if( argc < 2 )
return EXIT_FAILURE;

char const *sizeStr = argv[1];
bool hex = sizeStr[0] == '0' && (sizeStr[1] == 'x' || sizeStr[1] == 'X');

sizeStr += 2 * hex;
size_t end;

if( from_chars_result fcr = from_chars( sizeStr, sizeStr + strlen(

sizeStr ), end, !hex ? 10 : 16 ); fcr.ec != errc() || *fcr.ptr )
return EXIT_FAILURE;
if( (ptrdiff_t)end++ < 0 )
throw bad_alloc();
unsigned hc = jthread::hardware_concurrency();
if( !hc )
return EXIT_FAILURE;
using word_t = size_t;
constexpr size_t
BITS = sizeof(word_t) * 8,
CL_BITS = CL_SIZE * 8;
union alignas(CL_SIZE) ndi_words_cl { word_t words[CL_SIZE /
sizeof(word_t)]; ndi_words_cl() {} };
size_t nCls = (end + CL_BITS - 1) / CL_BITS;
vector<ndi_words_cl> sieveCachelines( (end + CL_BITS - 1) / CL_BITS );
size_t nWords = (end + BITS - 1) / BITS;
span<word_t> sieve( &sieveCachelines[0].words[0], nWords );
auto fill = sieve.begin();
*fill++ = (word_t)0xAAAAAAAAAAAAAAACu;
if( fill != sieve.end() )
for_each( fill, sieve.end(), []( word_t &w ) { w =
(word_t)0xAAAAAAAAAAAAAAAAu; } );
size_t sqrt = (ptrdiff_t)ceil( ::sqrt( (ptrdiff_t)end ) );
mutex punchMtx;
condition_variable punchCv;
using range_t = pair<size_t, size_t>;
deque<range_t> punchQueue;
atomic_uint finishCountDown;
binary_semaphore semFinish( false );
bool terminateThreads = false;
vector<jthread> punchThreads;
punchThreads.reserve( hc - 1 );
for( size_t prime = 3; prime < sqrt; ++prime )
{
for( auto pSieve = sieve.begin() + prime / BITS; ; )
{
if( word_t w = *pSieve >> prime % BITS; w ) [[likely]]
{
prime += countr_zero( w );
break;
}
prime = prime + BITS & -(ptrdiff_t)BITS;
if( ++pSieve == sieve.end()) [[unlikely]]
break;
}
if( prime >= sqrt ) [[unlikely]]
break;
size_t
firstMultiple = prime * prime,
nBits = end - firstMultiple,
nPartitions = nBits / prime / 1'000;
nPartitions += (size_t)!nPartitions;
nPartitions = nPartitions <= hc ? nPartitions : hc;
size_t bitsPerPartition = nBits / nPartitions;
unique_lock lock( punchMtx );
unsigned nThreads = -1;
for( size_t t = nPartitions, rEnd = end, rStart; t--; rEnd = rStart )
{
rStart = firstMultiple + t * bitsPerPartition;
rStart -= rStart % prime;
size_t
rClStart = rStart & -(CL_SIZE * 8),
rFirstPrimeInCl = rClStart + (rStart - rClStart) % prime;
rStart = rFirstPrimeInCl >= firstMultiple ? rFirstPrimeInCl : rStart;
if( rStart < rEnd )
punchQueue.emplace_back( rStart, rEnd ),
punchCv.notify_one(),
++nThreads;
}
range_t firstRange( punchQueue.front() );
punchQueue.pop_front();
finishCountDown.store( nThreads, memory_order_relaxed );
lock.unlock();
auto punch = [&]( auto, size_t start, size_t end )
{
size_t multiple = start;
if( prime >= BITS ) [[likely]]
do
sieve[multiple / BITS] &= rotl( (word_t)-2, multiple % BITS );
while( (multiple += prime) < end );
else
{
auto pSieve = sieve.begin() + multiple / BITS;
do
{
word_t
word = *pSieve,
mask = rotl( (word_t)-2, multiple % BITS ),
oldMask;
do
word &= mask,
oldMask = mask,
mask = rotl( mask, prime % BITS ),
multiple += prime;
while( mask < oldMask );
*pSieve++ = word;
} while( multiple < end );
}
};
for( unsigned t = (unsigned)punchThreads.size(); t < nThreads; ++t )
punchThreads.emplace_back( [&]()
{
for( ; ; )
{
unique_lock lock( punchMtx );
punchCv.wait( lock, [&]() { return punchQueue.size() ||
terminateThreads; } );
if( terminateThreads )
return;
range_t range = punchQueue.front();
punchQueue.pop_front();
lock.unlock();
punch( [] {}, range.first, range.second );
if( finishCountDown.fetch_sub( 1, memory_order_relaxed ) == 1 )
semFinish.release( 1 );
}
} );
punch( [] {}, firstRange.first, firstRange.second );
if( nThreads )
semFinish.acquire();
prime += 0;
}
if( punchThreads.size() )
terminateThreads = true,
atomic_thread_fence( memory_order_release ),
punchCv.notify_all();
auto scan = [&]( auto fn )
{
for( size_t prime = 0; word_t w : sieve )
{
bool hasBits = w;
for( unsigned gap; w; ++prime, w >>= gap, w >>= 1 )
{
gap = countr_zero( w );
prime += gap;
if( prime >= end )
return;
fn( prime );
}
prime = prime + BITS - hasBits & -(ptrdiff_t)BITS;

}
};
if( argc >= 3 && !*argv[2] )

return EXIT_SUCCESS;

ofstream ofs;
ofs.exceptions( ofstream::failbit | ofstream::badbit );

ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::binary |
ofstream::trunc );
constexpr size_t
BUF_SIZE = 0x100000,
AHEAD = 32;
union ndi_char { char c; ndi_char() {} };
vector<ndi_char> printBuf( BUF_SIZE + AHEAD - 1 );
auto
bufBegin = printBuf.begin(),
bufEnd = bufBegin,

bufFlush = bufBegin + BUF_SIZE;

auto print = [&]() { ofs << string_view( &bufBegin->c, &to_address(
bufEnd )->c ); };
scan( [&]( size_t prime )
{
auto [ptr, ec] = to_chars( &bufEnd->c, &bufEnd[AHEAD - 1].c, prime );
if( ec != errc() ) [[unlikely]]
throw system_error( (int)ec, generic_category(), "converson failed" );
bufEnd = ptr - &bufBegin->c + bufBegin; // pointer to iterator - NOP
bufEnd++->c = '\n';
if( bufEnd >= bufFlush ) [[unlikely]]
print(), bufEnd = bufBegin;
} );
print();
}

[Bonita Montero]

Ich hab das nochmal weiter optimiert und jetzt innerhalb der Threads
noch partitioniert, dass die Partitionen in den L2-cache gehen. Dazu
wird mit CPUID die Größe des L2-Caches ermittelt.
Auf jeden Fall dauert auf meinem 7950X Zen4 16-Kerner das Berechnen
der ersten 210 Millionen Primzahlen nur 0.29 Sekunden. Auf meinem
Linux-Rechner, ein Zen2 64-Kerner gar nur 0.09 Sekunden.

#include <cstdlib>
#include <charconv>
#include <cstring>
#include <vector>
#include <algorithm>
#include <cmath>
#include <bit>
#include <fstream>
#include <string_view>
#include <thread>

#include <utility>
#include <new>
#include <span>
#if defined(_MSC_VER)
#include <intrin.h>
#elif defined(__GNUC__) || defined(__clang__)
#include <cpuid.h>
#endif

using namespace std;

#if defined(__cpp_lib_hardware_interference_size)
constexpr ptrdiff_t CL_SIZE = hardware_destructive_interference_size;
#else
constexpr ptrdiff_t CL_SIZE = 64;
#endif

size_t getL2Size();
bool smt();

int main( int argc, char **argv )
{
if( argc < 2 )
return EXIT_FAILURE;

auto parse = []( char const *str, bool hex, auto &value )
{
from_chars_result fcr = from_chars( str, str + strlen( str ), value,

!hex ? 10 : 16 );

return fcr.ec == errc() && !*fcr.ptr;

};
char const *sizeStr = argv[1];
bool hex = sizeStr[0] == '0' && (sizeStr[1] == 'x' || sizeStr[1] == 'X');
sizeStr += 2 * hex;
size_t end;

if( !parse( sizeStr, hex, end ) )

return EXIT_FAILURE;
if( (ptrdiff_t)end++ < 0 )
throw bad_alloc();
unsigned

hc = jthread::hardware_concurrency(),
nThreads;
if( argc < 4 || !parse( argv[3], false, nThreads ) )
nThreads = hc;
if( !nThreads )

return EXIT_FAILURE;
using word_t = size_t;
constexpr size_t
BITS = sizeof(word_t) * 8,
CL_BITS = CL_SIZE * 8;
union alignas(CL_SIZE) ndi_words_cl { word_t words[CL_SIZE /
sizeof(word_t)]; ndi_words_cl() {} };
size_t nCls = (end + CL_BITS - 1) / CL_BITS;
vector<ndi_words_cl> sieveCachelines( (end + CL_BITS - 1) / CL_BITS );
size_t nWords = (end + BITS - 1) / BITS;
span<word_t> sieve( &sieveCachelines[0].words[0], nWords );
auto fill = sieve.begin();
*fill++ = (word_t)0xAAAAAAAAAAAAAAACu;
if( fill != sieve.end() )
for_each( fill, sieve.end(), []( word_t &w ) { w =
(word_t)0xAAAAAAAAAAAAAAAAu; } );
size_t sqrt = (ptrdiff_t)ceil( ::sqrt( (ptrdiff_t)end ) );

auto punch = [&]( auto, size_t start, size_t end, size_t prime )
{
if( start >= end )
return;

size_t multiple = start;
if( prime >= BITS ) [[likely]]
do
sieve[multiple / BITS] &= rotl( (word_t)-2, multiple % BITS );
while( (multiple += prime) < end );
else
{
auto pSieve = sieve.begin() + multiple / BITS;
do
{
word_t
word = *pSieve,
mask = rotl( (word_t)-2, multiple % BITS ),
oldMask;
do
word &= mask,
oldMask = mask,
mask = rotl( mask, prime % BITS ),
multiple += prime;
while( mask < oldMask );
*pSieve++ = word;
} while( multiple < end );
}
};

for( size_t prime = 3; prime < sqrt; ++prime )
{
for( auto pSieve = sieve.begin() + prime / BITS; ; )
{
if( word_t w = *pSieve >> prime % BITS; w ) [[likely]]
{
prime += countr_zero( w );
break;
}
prime = prime + BITS & -(ptrdiff_t)BITS;

if( prime >= sqrt )

break;
}
if( prime >= sqrt ) [[unlikely]]
break;

punch( false_type(), prime * prime, sqrt, prime );
}
auto scan = [&]( size_t start, size_t end, auto fn )
{
for( size_t prime = start; prime < end; )
{
word_t word = sieve[prime / BITS] >> prime % BITS;
bool hasBits = word;
for( unsigned gap; word; word >>= gap, word >>= 1 )
{
gap = countr_zero( word );
if( (prime += gap) >= end ) [[unlikely]]
return;
fn( prime );
if( ++prime >= end ) [[unlikely]]
return;
}
prime = (prime + BITS - hasBits) & -(ptrdiff_t)BITS;
}
};
double bitsPerPartition = (end - sqrt) / nThreads;

using range_t = pair<size_t, size_t>;

vector<pair<size_t, size_t>> ranges;
ranges.reserve( nThreads );
for( size_t t = nThreads, rEnd = end, rStart; t--; rEnd = rStart )
{
rStart = sqrt + t * bitsPerPartition;
size_t rClStart = rStart & -(CL_SIZE * 8);
rStart = rClStart >= sqrt ? rClStart : sqrt;
if( rStart < rEnd )
ranges.emplace_back( rStart, rEnd );
}
vector<jthread> threads;
threads.reserve( ranges.size() - 1 );
static auto dbl = []( ptrdiff_t x ) { return (double)x; };
double maxPartitionBits = dbl( getL2Size() * 8 ) / 3 * (smt() &&
nThreads > hc / 2 ? 1 : 2);
for( range_t const &range : ranges )
{
auto partitionedPunch = [&]( size_t rStart, size_t rEnd )
{
double nBits = dbl( rEnd - rStart );
ptrdiff_t nPartitions = (ptrdiff_t)ceil( nBits / maxPartitionBits );
double bitsPerPartition = nBits / dbl( nPartitions );
for( ptrdiff_t p = nPartitions, pEnd = rEnd, pStart; p--; pEnd = pStart )
{
pStart = rStart + (ptrdiff_t)((double)p * bitsPerPartition);
pStart &= -(CL_SIZE * 8);
scan( 3, sqrt,

[&]( size_t prime )
{

punch( true_type(), (pStart + prime - 1) / prime * prime, pEnd,
prime );
} );
}
};
if( &range != &ranges.back() )
threads.emplace_back( partitionedPunch, range.first, range.second );
else
partitionedPunch( range.first, range.second );
}
threads.resize( 0 );

if( argc >= 3 && !*argv[2] )
return EXIT_SUCCESS;
ofstream ofs;
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::binary |
ofstream::trunc );
constexpr size_t
BUF_SIZE = 0x100000,
AHEAD = 32;
union ndi_char { char c; ndi_char() {} };
vector<ndi_char> printBuf( BUF_SIZE + AHEAD - 1 );
auto
bufBegin = printBuf.begin(),
bufEnd = bufBegin,
bufFlush = bufBegin + BUF_SIZE;
auto print = [&]() { ofs << string_view( &bufBegin->c, &to_address(
bufEnd )->c ); };

scan( 2, end,

[&]( size_t prime )
{
auto [ptr, ec] = to_chars( &bufEnd->c, &bufEnd[AHEAD - 1].c, prime );
if( ec != errc() ) [[unlikely]]
throw system_error( (int)ec, generic_category(), "converson failed" );
bufEnd = ptr - &bufBegin->c + bufBegin; // pointer to iterator - NOP
bufEnd++->c = '\n';
if( bufEnd >= bufFlush ) [[unlikely]]
print(), bufEnd = bufBegin;
} );
print();
}

unsigned cpuId( unsigned code, unsigned (&regs)[4] )
{
#if defined(_MSC_VER)
__cpuid( (int *)regs, code );
#elif defined(__GNUC__) || defined(__clang__)
__cpuid(code, regs[0], regs[1], regs[2], regs[3]);
#endif
return regs[0];
}

bool smt()
{
unsigned regs[4];
if( cpuId( 0, regs ) < 1 )
return false;
cpuId( 1, regs );
return regs[3] >> 28 & 1;
}

size_t getL2Size()
{
unsigned regs[4];
if( cpuId( 0x80000000u, regs ) < 0x80000006u )
return 512ull * 1024;
cpuId( 0x80000006u, regs );
return (regs[2] >> 16) * 1024;
}

[Bonita Montero]

Ich hab das noch weiter optimiert und speichere nur noch die ungeraden
Bits.
Grundsätzlich basiert das Ganze auf der Idee des Segmented Sieve, d.h.
wenn ich alle Primzahlen bis zur Wurzel des Maximum habe kann ich mit
wenig Aufwand jeden Bereich bis zum Quadrat berechnen. Ich berechne
halt für jeden Thread eine Partition und innerhalb eines Threads wird
das ganze noch auf Blöcke aufgeteilt die in den L2 -Cache passen.
Letzlich hängt die Performance fast 1:1 vom Durchsatz des L2-Caches
ab, was sich u.A. daran zeigt, dass auf SMT-Systemen die Performance
mit halber und ganzer logischer Kern-Zahl nahezu gleich ist.
So, daran dürfte echt nichts mehr zu verbessern sein, alle Primzahlen
die in den 32 Bit Wertebereich gehen berechnet das Ganze in 0.15 Se-
kunden auf meinem Zen4-PC 16-Kerner. Auf dem 64-kernigen Zen2-PC
unter Linux braucht das Ganze 0.05 Sekunden.

#include <cstdlib>
#include <charconv>
#include <cstring>
#include <vector>
#include <algorithm>
#include <cmath>
#include <bit>
#include <fstream>
#include <string_view>
#include <thread>
#include <utility>
#include <new>
#include <span>

#include <array>
#include <cassert>

#if defined(_MSC_VER)
#include <intrin.h>
#elif defined(__GNUC__) || defined(__clang__)
#include <cpuid.h>
#endif

#if defined(_MSC_VER)
#pragma warning(disable: 26495)

#endif

using namespace std;

#if defined(__cpp_lib_hardware_interference_size)
constexpr ptrdiff_t CL_SIZE = hardware_destructive_interference_size;
#else
constexpr ptrdiff_t CL_SIZE = 64;
#endif

size_t getL2Size();
bool smt();

int main( int argc, char **argv )
{
if( argc < 2 )
return EXIT_FAILURE;

auto parse = []( char const *str, auto &value )

{
bool hex = str[0] == '0' && (str[1] == 'x' || str[1] == 'X');
str += 2 * hex;

auto [ptr, ec] = from_chars( str, str + strlen( str ), value, !hex ?
10 : 16 );
return ec == errc() && !*ptr;
};
size_t end;
if( !parse( argv[1], end ) )
return EXIT_FAILURE;
if( end < 2 || (ptrdiff_t)end++ < 0 )

throw bad_alloc();
unsigned
hc = jthread::hardware_concurrency(),
nThreads;

if( argc < 4 || !parse( argv[3], nThreads ) )

nThreads = hc;
if( !nThreads )
return EXIT_FAILURE;
using word_t = size_t;
constexpr size_t

BITS_PER_CL = CL_SIZE * 8,
BITS = sizeof(word_t) * 8;
auto bitEnd = []( size_t end ) { return end / 2 + (end & 1 ^ 1); };
size_t nBits = bitEnd( end );

union alignas(CL_SIZE) ndi_words_cl { word_t words[CL_SIZE /
sizeof(word_t)]; ndi_words_cl() {} };

vector<ndi_words_cl> sieveCachelines( (nBits + BITS_PER_CL - 1) /
BITS_PER_CL );
span<word_t> sieve( &sieveCachelines[0].words[0], (nBits + BITS - 1) /
BITS );
fill( sieve.begin(), sieve.end(), (word_t)-1 );

size_t sqrt = (ptrdiff_t)ceil( ::sqrt( (ptrdiff_t)end ) );
auto punch = [&]( auto, size_t start, size_t end, size_t prime )
{

size_t bit = start / 2;
end = bitEnd( end );
if( bit >= end )
return;

if( prime >= BITS ) [[likely]]

do [[likely]]
sieve[bit / BITS] &= rotl( (word_t)-2, bit % BITS );
while( (bit += prime) < end );
else
{
auto pSieve = sieve.begin() + bit / BITS;
do [[likely]]
{
word_t
word = *pSieve,
mask = rotl( (word_t)-2, bit % BITS ),

oldMask;
do
word &= mask,
oldMask = mask,
mask = rotl( mask, prime % BITS ),

bit += prime;

while( mask < oldMask );
*pSieve++ = word;

} while( bit < end );
}
};
for( size_t bSqrt = bitEnd( sqrt ), bit = 1; bit < bSqrt; ++bit )
[[likely]]
{
auto pSieve = sieve.begin() + bit / BITS;
do [[likely]]
if( word_t w = *pSieve >> bit % BITS; w ) [[likely]]
{
bit += countr_zero( w );
break;
}
while( (bit = bit + BITS & -(ptrdiff_t)BITS) < bSqrt );
if( bit >= bSqrt ) [[unlikely]]
break;
size_t prime = 2 * bit + 1;

punch( false_type(), prime * prime, sqrt, prime );
}
auto scan = [&]( size_t start, size_t end, auto fn )
{

start /= 2;
end = bitEnd( end );

if( start >= end )
return;

auto pSieve = sieve.begin() + start / BITS;
for( size_t bit = start; ; )
{
word_t word = *pSieve++ >> bit % BITS;
bool hasBits = word;
for( unsigned gap; word; word >>= gap, word >>= 1 ) [[likely]]

{
gap = countr_zero( word );

if( (bit += gap) >= end ) [[unlikely]]
return;
fn( bit * 2 + 1, true_type() );
if( ++bit >= end ) [[unlikely]]
return;
}
if( bit >= end ) [[unlikely]]
break;
bit = (bit + BITS - hasBits) & -(ptrdiff_t)BITS;

}
};
static auto dbl = []( ptrdiff_t x ) { return (double)x; };

double threadTange = dbl( end - sqrt ) / nThreads;

using range_t = pair<size_t, size_t>;
vector<pair<size_t, size_t>> ranges;
ranges.reserve( nThreads );

for( size_t t = nThreads, rEnd = end, trStart; t--; rEnd = trStart )
[[likely]]
{
trStart = sqrt + (ptrdiff_t)((ptrdiff_t)t * threadTange + 0.5);
size_t trClStart = trStart & -(2 * CL_SIZE * 8);
trStart = trClStart >= sqrt ? trClStart : sqrt;
if( trStart < rEnd )
ranges.emplace_back( trStart, rEnd );
}
double maxCacheRange = dbl( getL2Size() * (8 * 2) ) / 3 * (smt() &&

nThreads > hc / 2 ? 1 : 2);

vector<jthread> threads;
threads.reserve( ranges.size() - 1 );

for( range_t const &range : ranges )
{

auto cacheAwarePunch = [&]( size_t rStart, size_t rEnd )
{
double thisThreadRange = dbl( rEnd - rStart );
ptrdiff_t nCachePartitions = (ptrdiff_t)ceil( thisThreadRange /
maxCacheRange );
double cacheRange = thisThreadRange / dbl( nCachePartitions );
for( size_t p = nCachePartitions, cacheEnd = rEnd, cacheStart; p--;
cacheEnd = cacheStart ) [[likely]]
{
cacheStart = rStart + (ptrdiff_t)((double)(ptrdiff_t)p * cacheRange
+ 0.5);
cacheStart &= -(2 * CL_SIZE * 8);
cacheStart = cacheStart >= sqrt ? cacheStart : sqrt;
scan( 3, sqrt,
[&]( size_t prime, auto )
{
size_t start = (cacheStart + prime - 1) / prime * prime;
start = start & 1 ? start : start + prime;
punch( true_type(), start, cacheEnd, prime );

} );
}
};
if( &range != &ranges.back() )

threads.emplace_back( cacheAwarePunch, range.first, range.second );
else
cacheAwarePunch( range.first, range.second );

}
threads.resize( 0 );
if( argc >= 3 && !*argv[2] )
return EXIT_SUCCESS;
ofstream ofs;
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::binary |
ofstream::trunc );
constexpr size_t
BUF_SIZE = 0x100000,
AHEAD = 32;
union ndi_char { char c; ndi_char() {} };

vector<ndi_char> rawPrintBuf( BUF_SIZE + AHEAD - 1 );
span printBuf( &rawPrintBuf.front().c, &rawPrintBuf.back().c + 1 );

auto
bufBegin = printBuf.begin(),
bufEnd = bufBegin,
bufFlush = bufBegin + BUF_SIZE;

auto print = [&]() { ofs << string_view( bufBegin, bufEnd ); };
auto printPrime = [&]( size_t prime, auto )
{
auto [ptr, ec] = to_chars( &*bufEnd, &bufEnd[AHEAD - 1], prime );

if( ec != errc() ) [[unlikely]]
throw system_error( (int)ec, generic_category(), "converson failed" );

bufEnd = ptr - &*bufBegin + bufBegin; // pointer to iterator - NOP
*bufEnd++ = '\n';

if( bufEnd >= bufFlush ) [[unlikely]]
print(), bufEnd = bufBegin;

};
printPrime( 2, false_type() );
scan( 3, end, printPrime );
print();
}

array<unsigned, 4> cpuId( unsigned code )
{
array<unsigned, 4> regs;
#if defined(_MSC_VER)
__cpuid( (int *)&regs[0], code );

#elif defined(__GNUC__) || defined(__clang__)
__cpuid(code, regs[0], regs[1], regs[2], regs[3]);
#endif

return regs;
}

bool smt()
{
if( cpuId( 0 )[0] < 1 )
return false;
return cpuId( 1 )[3] >> 28 & 1;
}

size_t getL2Size()
{
if( cpuId( 0x80000000u )[0] < 0x80000006u )
return 512ull * 1024;
return (cpuId( 0x80000006u )[2] >> 16) * 1024;
}