I further improved my code. Now scanning for the next prime is
signifiantly more efficient because I use a rotating bit mask.

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

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

using namespace std;

int main( int argc, char **argv )
{ constexpr size_t BUF_SIZE = 0x100000;
try
{
size_t end;
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 ), end, !hex ? 10 : 16 ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;
if( size_t e = end; ++end + 7 < e )
throw bad_alloc();
using bit_vec = vector<uint8_t>;
using bit_it = bit_vec::iterator;
bit_vec sieve( (end + 7) / 8, -1 );
ofstream ofs;
union ndi_t { char c; ndi_t() {} };
using ndi_vec = vector<ndi_t>;
using ndi_vec_it = ndi_vec::iterator;
constexpr size_t AHEAD = 32;
ndi_vec printBuf;
ndi_vec_it bufBegin, bufEnd, bufFlush;
if( argc < 3 || *argv[2] )
{
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc );
printBuf.resize( BUF_SIZE + AHEAD - 1 );
bufBegin = printBuf.begin();
bufEnd = bufBegin;
bufFlush = bufBegin + BUF_SIZE;
}
auto print = [&]() { ofs << string_view( &to_address( bufBegin )->c,
&to_address( bufEnd )->c ); };
for( size_t p = 2; p != end; )
{
bit_it pSieve = sieve.begin();
for( size_t m = p; (m += p) < end; )
pSieve[m / 8] &= ~(1u << m % 8);
pSieve += ++p / 8;
for( uint8_t m = 1u << p % 8, mNew; ; ++p )
if( p == end ) [[unlikely]]
goto brk;
else if( *pSieve & m ) [[unlikely]]
break;
else
mNew = rotl( m, 1 ),
pSieve += mNew < m,
m = mNew;
if( !ofs.is_open() )
continue;
auto [ptr, ec] = to_chars( &bufEnd->c, &to_address( bufEnd + (AHEAD -
1) )->c, p );
if( (bool)ec ) [[unlikely]]
throw system_error( (int)ec, generic_category(), "converson failed" );
bufEnd = ptr - &bufBegin->c + bufBegin; // NOP
bufEnd++->c = '\n';
if( bufEnd >= bufFlush )
print(),
bufEnd = bufBegin;
}
brk:
print();
return EXIT_SUCCESS;
}
catch( bad_alloc const & )
{
cout << "out of memory" << endl;
}
catch( ios_base::failure const & )
{
cout << "I/O error" << endl;
}
catch( system_error const &se )
{
cout << se.what() << endl;
}
}

On 2023-08-28 8:23 AM, Öö Tiib wrote:
> On Monday, 28 August 2023 at 02:50:05 UTC+3, Ben Bacarisse wrote:
>> Öö Tiib <oot...@hot.ee> writes:
>>
>>> On Thursday, 17 August 2023 at 05:26:37 UTC+3, Tim Rentsch wrote:
>>>>
>>>> I suggest computing all 32 bit primes for the benchmark
>>>> tests. Shouldn't take too long.
>>>>
>>> Took 14 minutes to find those 203280221 primes single threaded
>>> on i7-8550U.
>> What do you mean by find? Do you include printing them out because that
>> will depend rather more on the OS.
>>
> I did try to invent a sieve algorithm this weekend, it worked, but bit
> sluggish, so I posted it.
> I did not measure printing or storing to file, only sieving and counting
> (to test that there are 203280221).
>
>> I wanted to try a newer algorithm, so I coded up a rather naive version
>> of Atkin & Bernstein's sieve in C. It needs some tidying and I imagine
>> a bit can be shaved off the time.
>>
>> It can generate a bitmap of the 203280221 primes and then count them in
>> about 40 seconds on an i5-8265U CPU @ 1.60GHz (single threaded). It
>> takes about 55 seconds to write them (in decimal) to a file. The result
>> is over 2GB:
>>
>> $ wc primes
>> 203280221 203280221 2178719347 primes
>>
> Thanks, so I apparently need to optimize my algorithm at least 30 times
> faster for it to be worth anything. Maybe next weekend.

The code I posted earlier, with limits changed to compute primes up to
32^2 - 1 and with `int` changed to `int64_t`, reports

28.9473339 seconds.
203280221 primes.

msinfo32 says "Processor Intel(R) Core(TM) i7-8700T CPU @ 2.40GHz, 2400
Mhz, 6 Core(s), 12 Logical Processor(s)"

so accounting for clock rate it seems a little bit slower than Ben's.

- Alf

#include <fmt/core.h>
#include "cppm-Timer.hpp"

#include <vector>
using std::vector;

#include <stdint.h> // int64_t
#include <stdlib.h> // EXIT_FAILURE
#include <stdio.h> // fopen, fclose

using Int64 = int64_t;

auto main() -> int
{ constexpr Int64 n = 4'294'967'295; // 2^32 - 1
constexpr Int64 sqrt_n = 65'535;

cppm::Timer timer;
auto is_composite = vector<bool>( n + 1 );
for( Int64 i = 4; i < n; i += 2 ) { is_composite[i] = true; }
for( Int64 i = 3; i < n; ) {
while( is_composite[i] ) { ++i; }
const Int64 prime = i;
if( prime <= sqrt_n ) {
for( Int64 j = prime*prime; j < n; j += 2*prime ) {
is_composite[j] = true;
}
}
++i;
}
timer.stop();
fmt::print( "{} seconds.\n", timer.seconds() );

Int64 count = 0;
for( Int64 i = 2; i < n; ++i ) {
if( not is_composite[i] ) {
++count;
}
}
fmt::print( "{} primes.\n", count );

auto f = fopen( "primes.alf.txt", "w" );
if( not f ) { return EXIT_FAILURE; }
for( Int64 i = 2; i < n; ++i ) {
if( not is_composite[i] ) {
fmt::print( f, "{:d}\n", i );
}
}
fclose( f );
}

On Monday, 28 August 2023 at 16:54:04 UTC+3, Alf P. Steinbach wrote:
> On 2023-08-28 8:23 AM, Öö Tiib wrote:
> > On Monday, 28 August 2023 at 02:50:05 UTC+3, Ben Bacarisse wrote:
> >> Öö Tiib <oot...@hot.ee> writes:
> >>
> >>> On Thursday, 17 August 2023 at 05:26:37 UTC+3, Tim Rentsch wrote:
> >>>>
> >>>> I suggest computing all 32 bit primes for the benchmark
> >>>> tests. Shouldn't take too long.
> >>>>
> >>> Took 14 minutes to find those 203280221 primes single threaded
> >>> on i7-8550U.
> >> What do you mean by find? Do you include printing them out because that
> >> will depend rather more on the OS.
> >>
> > I did try to invent a sieve algorithm this weekend, it worked, but bit
> > sluggish, so I posted it.
> > I did not measure printing or storing to file, only sieving and counting
> > (to test that there are 203280221).
> >
> >> I wanted to try a newer algorithm, so I coded up a rather naive version
> >> of Atkin & Bernstein's sieve in C. It needs some tidying and I imagine
> >> a bit can be shaved off the time.
> >>
> >> It can generate a bitmap of the 203280221 primes and then count them in
> >> about 40 seconds on an i5-8265U CPU @ 1.60GHz (single threaded). It
> >> takes about 55 seconds to write them (in decimal) to a file. The result
> >> is over 2GB:
> >>
> >> $ wc primes
> >> 203280221 203280221 2178719347 primes
> >>
> > Thanks, so I apparently need to optimize my algorithm at least 30 times
> > faster for it to be worth anything. Maybe next weekend.
>
> The code I posted earlier, with limits changed to compute primes up to
> 32^2 - 1 and with `int` changed to `int64_t`, reports
>
> 28.9473339 seconds.
> 203280221 primes.
>
> msinfo32 says "Processor Intel(R) Core(TM) i7-8700T CPU @ 2.40GHz, 2400
> Mhz, 6 Core(s), 12 Logical Processor(s)"
>
> so accounting for clock rate it seems a little bit slower than Ben's.
>
Quite likely Ben's "naive" is better algorithm. Pasted your algorithm to
Apple M1 Max (that is 3.2 GHz):

sieving 18.1471 seconds.
counting 1.8121 seconds.
203280221 primes.

Am 28.08.2023 um 18:03 schrieb Öö Tiib:

> Quite likely Ben's "naive" is better algorithm. Pasted your algorithm to
> Apple M1 Max (that is 3.2 GHz):
>
> sieving 18.1471 seconds.
> counting 1.8121 seconds.
> 203280221 primes.

I'd like to know how this code performs on a Apple
silicon Mac. That's while the M1 Max has a memory
bandwidh slightly below 300GB/s.

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

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

using namespace std;

int main( int argc, char **argv )
{ constexpr size_t BUF_SIZE = 0x100000;
try
{
size_t end;
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 ), end, !hex ? 10 : 16 ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;
if( !++end )
throw bad_alloc();
using word_t = uint64_t;
constexpr size_t BITS = sizeof(word_t) * 8;
size_t bEnd = (end + BITS - 1) & -(ptrdiff_t)BITS;
if( bEnd < end )
throw bad_alloc();
using bit_vec = vector<word_t>;
using bit_it = bit_vec::iterator;
bit_vec sieve( (end + BITS - 1) / BITS, -1 );
ofstream ofs;
union ndi_t { char c; ndi_t() {} };
using ndi_vec = vector<ndi_t>;
using ndi_vec_it = ndi_vec::iterator;
constexpr size_t AHEAD = 32;
ndi_vec printBuf;
ndi_vec_it bufBegin, bufEnd, bufFlush;
if( argc < 3 || *argv[2] )
{
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc );
printBuf.resize( BUF_SIZE + AHEAD - 1 );
bufBegin = printBuf.begin();
bufEnd = bufBegin;
bufFlush = bufBegin + BUF_SIZE;
}
auto print = [&]() { ofs << string_view( &to_address( bufBegin )->c,
&to_address( bufEnd )->c ); };
for( size_t p = 2; p < end; )
{
bit_it
pSieve = sieve.begin(),
pSieveEnd = sieve.end();
if( size_t m = 2 * p; m < end )
if( p >= BITS )
{
word_t mask = ~((word_t)1 << m % BITS);
do
pSieve[m / BITS] &= mask,
mask = rotl( mask, (int)p % BITS ),
m += p;
while( m < end );
}
else
{
pSieve += m / BITS;
word_t newMask = ~((word_t)1 << m % BITS);
do
{
word_t word = *pSieve, mask;
do
mask = newMask,
word &= mask,
newMask = rotl( mask, (int)p );
while( newMask < mask );
*pSieve = word;
} while( ++pSieve != pSieveEnd );
}
if( ++p == end ) [[unlikely]]
break;
word_t mask = (word_t)1 << p % BITS;
pSieve = sieve.begin() + p / BITS;
do
{
word_t word = *pSieve;
do
if( word & mask ) [[unlikely]]
goto foundNext;
else
++p,
mask = rotl( mask, 1 );
while( mask != 1 );
} while( ++pSieve != pSieveEnd );
foundNext:
if( p >= end ) [[unlikely]]
break;
if( !ofs.is_open() )
continue;
auto [ptr, ec] = to_chars( &bufEnd->c, &to_address( bufEnd + (AHEAD -
1) )->c, p );
if( (bool)ec ) [[unlikely]]
throw system_error( (int)ec, generic_category(), "converson failed" );
bufEnd = ptr - &bufBegin->c + bufBegin; // NOP
bufEnd++->c = '\n';
if( bufEnd >= bufFlush )
print(),
bufEnd = bufBegin;
}
print();
return EXIT_SUCCESS;
}
catch( bad_alloc const & )
{
cout << "out of memory" << endl;
}
catch( ios_base::failure const & )
{
cout << "I/O error" << endl;
}
catch( system_error const &se )
{
cout << se.what() << endl;
}
}

"Alf P. Steinbach" <alf.p.steinbach@gmail.com> writes:

> On 2023-08-28 8:23 AM, Öö Tiib wrote:
>> On Monday, 28 August 2023 at 02:50:05 UTC+3, Ben Bacarisse wrote:
>>> Öö Tiib <oot...@hot.ee> writes:
>>>
>>>> On Thursday, 17 August 2023 at 05:26:37 UTC+3, Tim Rentsch wrote:
>>>>>
>>>>> I suggest computing all 32 bit primes for the benchmark
>>>>> tests. Shouldn't take too long.
>>>>>
>>>> Took 14 minutes to find those 203280221 primes single threaded
>>>> on i7-8550U.
>>> What do you mean by find? Do you include printing them out because that
>>> will depend rather more on the OS.
>>>
>> I did try to invent a sieve algorithm this weekend, it worked, but bit
>> sluggish, so I posted it.
>> I did not measure printing or storing to file, only sieving and counting
>> (to test that there are 203280221).
>>
>>> I wanted to try a newer algorithm, so I coded up a rather naive version
>>> of Atkin & Bernstein's sieve in C. It needs some tidying and I imagine
>>> a bit can be shaved off the time.
>>>
>>> It can generate a bitmap of the 203280221 primes and then count them in
>>> about 40 seconds on an i5-8265U CPU @ 1.60GHz (single threaded). It
>>> takes about 55 seconds to write them (in decimal) to a file. The result
>>> is over 2GB:
>>>
>>> $ wc primes
>>> 203280221 203280221 2178719347 primes
>>>
>> Thanks, so I apparently need to optimize my algorithm at least 30 times
>> faster for it to be worth anything. Maybe next weekend.
>
> The code I posted earlier, with limits changed to compute primes up to 32^2
> - 1 and with `int` changed to `int64_t`, reports
>
> 28.9473339 seconds.
> 203280221 primes.
>
> msinfo32 says "Processor Intel(R) Core(TM) i7-8700T CPU @ 2.40GHz, 2400
> Mhz, 6 Core(s), 12 Logical Processor(s)"
>
> so accounting for clock rate it seems a little bit slower than Ben's.

On my machine, compiled with -O3, it takes just over 2 minutes. (I had
to remove the timing code and replace the fmt/core.h use with a call to
printf.)

It's not easy to compare i7 and i5 speeds from the clock rate alone.

--
Ben.

Tiib <ootiib@hot.ee> writes:

> On Thursday, 17 August 2023 at 05:26:37 UTC+3, Tim Rentsch wrote:
>
>> I suggest computing all 32 bit primes for the benchmark
>> tests. Shouldn't take too long.
>
> Took 14 minutes to find those 203280221 primes single threaded
> on i7-8550U.

I'm pleased to see that my suggestions sparked some interest.

Ben Bacarisse <ben.usenet@bsb.me.uk> writes:

> I wanted to try a newer algorithm, so I coded up a rather naive
> version of Atkin & Bernstein's sieve in C. It needs some tidying
> and I imagine a bit can be shaved off the time.

Thank you for this. Would you mind posting your code? Also do
you have a reference for the Atkin/Bernstein algorithm (or even
better, code?)?

Tim Rentsch <tr.17687@z991.linuxsc.com> writes:

> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>
>> I wanted to try a newer algorithm, so I coded up a rather naive
>> version of Atkin & Bernstein's sieve in C. It needs some tidying
>> and I imagine a bit can be shaved off the time.
>
> Thank you for this. Would you mind posting your code?

Hmm... It's a bit of a mess... but sure.

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

typedef unsigned int nat;
typedef unsigned long long int big_nat;

nat isqrt(nat n)
{ if (n < 2)
return n;
nat s = (isqrt(n >> 2) << 1) + 1;
return s - (s*s > n);
}

static inline void *bitset(nat n)
{ return calloc(1, n/8 + 1);
}

static inline void bitset_clear(void *set, nat n)
{ ((unsigned char *)set)[n/8] &= ~(1 << n%8);
}

static inline void bitset_flip(void *set, nat n)
{ ((unsigned char *)set)[n/8] ^= 1 << n%8;
}

static inline bool bitset_test(const void *set, nat n)
{ return ((const unsigned char *)set)[n/8] & (1 << n%8);
}

void *primes_up_to(nat n)
{ #define B(n) (1ull << (n))
typedef big_nat small_set;

nat sqrt_n = isqrt(n);
void *pset = bitset(n);

static small_set s1 =
B(1) | B(13) | B(17) | B(29) | B(37) | B(41) | B(49) | B(53);
for (big_nat x = 1; x <= sqrt_n; x++) {
for (nat y = 1; y <= sqrt_n; y += 2) {
big_nat m = 4*x*x + y*y;
if ((B(m % 60) & s1) && m <= n)
bitset_flip(pset, m);
}
}

static small_set s2 = B(7) | B(19) | B(31) | B(43);
for (big_nat x = 1; x <= sqrt_n; x += 2) {
for (nat y = 2; y <= sqrt_n; y += 2) {
big_nat m = 3*x*x + y*y;
if ((B(m % 60) & s2) && m <= n)
bitset_flip(pset, m);
}
}

static small_set s3 = B(11) | B(23) | B(47) | B(59);
for (big_nat x = 2; x <= sqrt_n; x++) {
for (int y = x - 1; y >= 1; y -= 2) {
big_nat m = 3*x*x - y*y;
if ((B(m % 60) & s3) && m <= n)
bitset_flip(pset, m);
}
}

static small_set s =
B(1) | B(7) | B(11) | B(13) | B(17) | B(19) | B(23) | B(29) |
B(31) | B(37) | B(41) | B(43) | B(47) | B(49) | B(53) | B(59);
for (big_nat m = 2; m <= n; m++) {
if (!(B(m % 60) & s))
continue;
big_nat mm = m*m;
if (mm > n)
break;
if (bitset_test(pset, m)) {
for (big_nat m2 = 1; m2 <= n; m2++) {
if (!(B(m2 % 60) & s))
continue;
big_nat c = mm * m2;
if (c > n)
break;
bitset_clear(pset, c);
}
}
}
bitset_flip(pset, 2);
bitset_flip(pset, 3);
bitset_flip(pset, 5);
return pset;
}

int main(int argc, char **argv)
{ int p = argc > 1 ? atoi(argv[1]) : 32;
nat n = (nat)-1 >> 32 - p;
fprintf(stderr, "p=%d, n=%u\n", p, n);

void *pset = primes_up_to(n);
nat count = 0;
while (n > 0)
if (bitset_test(pset, n--))
count += 1;
printf("%u\n", count);
free(pset);
}

I wanted to get rid of the big_nat type. I's sure it can be written
using only 32-bit types.

> Also do
> you have a reference for the Atkin/Bernstein algorithm (or even
> better, code?)?

Well, there's this: https://en.wikipedia.org/wiki/Sieve_of_Atkin

--
Ben.

I think that's my final solution. Currently I store the bits in
size_ts and I look up the next set bit with countr_zero which maps
to TZCNT on x64 and according to some small code on godbolt ARM
also has a proper instruction for that.
Resettig the bits for the composites is done by a rotating bit-mask.
If the prime is less than size_t * 8 there's a special case optimi-
zation.

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

#if defined(_MSC_VER)
#pragma warning(disable: 26495) // uninitialized member
#endif
#if defined(__llvm__)
#pragma clang diagnostic ignored "-Wdangling-else"
#endif

using namespace std;

int main( int argc, char **argv )
{ constexpr size_t BUF_SIZE = 0x100000;
try
{
size_t end;
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 ), end, !hex ? 10 : 16 ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;
if( !++end )
throw bad_alloc();
using word_t = size_t;
constexpr size_t BITS = sizeof(word_t) * 8;
size_t bEnd = (end + BITS - 1) & -(ptrdiff_t)BITS;
if( bEnd < end )
throw bad_alloc();
union ndi_word_t { word_t w; ndi_word_t( word_t w ) : w( w ) {}
};
using bit_vec = vector<ndi_word_t>;
using bit_it = bit_vec::iterator;
bit_vec sieve( (end + BITS - 1) / BITS,
(word_t)0xAAAAAAAAAAAAAAAAu );
ofstream ofs;
union ndi_t { char c; ndi_t() {} };
using ndi_vec = vector<ndi_t>;
constexpr size_t AHEAD = 32;
ndi_vec printBuf;
ndi_vec::iterator bufBegin, bufEnd, bufFlush;
if( argc < 3 || *argv[2] )
{
ofs.exceptions( ofstream::failbit | ofstream::badbit );
ofs.open( argc >= 3 ? argv[2] : "primes.txt",
ofstream::trunc );
printBuf.resize( BUF_SIZE + AHEAD - 1 );
bufBegin = printBuf.begin();
bufEnd = bufBegin;
bufFlush = bufBegin + BUF_SIZE;
}
auto print = [&]() { ofs << string_view( &to_address( bufBegin
)->c, &to_address( bufEnd )->c ); };
bit_it const pSieveEnd = sieve.end();
size_t sqr = (ptrdiff_t)ceil( sqrt( (ptrdiff_t)end ) );
for( size_t p = 2; p < end; )
{
if( ofs.is_open() )
{
auto [ptr, ec] = to_chars( &bufEnd->c, &to_address(
bufEnd + (AHEAD - 1) )->c, p );
if( (bool)ec ) [[unlikely]]
throw system_error( (int)ec, generic_category(),
"converson failed" );
bufEnd = ptr - &bufBegin->c + bufBegin; // NOP
bufEnd++->c = '\n';
if( bufEnd >= bufFlush )
ofs << string_view( &to_address( bufBegin )->c,
&to_address( bufEnd )->c ),
bufEnd = bufBegin;
}
if( ++p == end ) [[unlikely]]
break;
for( bit_it pSieve = sieve.begin() + p / BITS; ; )
{
size_t g = countr_zero( pSieve->w & (word_t)-1 << p %
BITS );
p = (p & -(ptrdiff_t)BITS) + g;
if( g != BITS )
break;
if( ++pSieve == pSieveEnd )
goto end;
}
bit_it pSieve = sieve.begin();
if( size_t m = p * p; p < sqr )
if( p >= BITS ) [[likely]]
{
word_t mask = rotl( (word_t)-2, m % BITS );
do
pSieve[m / BITS].w &= mask,
mask = rotl( mask, (int)p % BITS ),
m += p;
while( m < end );
}
else
{
pSieve += m / BITS;
word_t mask = rotl( (word_t)-2, m % BITS );
do
{
word_t punch = -1, oldMask;
do
punch &= mask,
oldMask = mask,
mask = rotl( mask, (int)p );
while( mask < oldMask );
pSieve->w &= punch;
} while( ++pSieve != pSieveEnd );
}
}
end:
if( ofs.is_open() )
print();
return EXIT_SUCCESS;
}
catch( bad_alloc const & )
{
cout << "out of memory" << endl;
}
catch( ios_base::failure const & )
{
cout << "I/O error" << endl;
}
catch( system_error const &se )
{
cout << se.what() << endl;
}
}

I implemented Alf's code according to my taste. The code has
about 40% the performance of my code if I don't output a file
with the numbers on Sklyake CPU with g++ 12 and -march=native.

#include <iostream>
#include <vector>
#include <charconv>
#include <fstream>
#include <cmath>
#include <cstring>

using namespace std;

int main( int argc, char **argv )
{ constexpr size_t BUF_SIZE = 0x100000;
size_t end;
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 ), end, !hex ? 10 : 16 ); (bool)fcr.ec || *fcr.ptr )
return EXIT_FAILURE;
if( !++end )
throw bad_alloc();
vector<bool> primes( end, true );
ofstream ofs;
string printBuf;
if( argc < 3 || *argv[2] )
ofs.exceptions( ofstream::failbit | ofstream::badbit ),
ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc ),
printBuf.reserve( BUF_SIZE + 31 );
size_t sqr = (ptrdiff_t)ceil( sqrt( (ptrdiff_t)end ) );
for( size_t p = 2; p < end; ++p )
{
while( !primes[p] )
if( ++p == end )
goto end;
if( p < sqr )
{
size_t m = p * p;
do
primes[m] = false;
while( (m += p) < end );
}
if( !ofs.is_open() )
continue;
char strPrime[32];
auto [ptr, ec] = to_chars( strPrime, std::end( strPrime ), p );
if( (bool)ec ) [[unlikely]]
throw system_error( (int)ec, generic_category(), "converson failed" );
*ptr++ = '\n';
printBuf += string_view( strPrime, ptr );
if( printBuf.size() >= BUF_SIZE )
ofs << printBuf,
printBuf.resize( 0 );
}
end:
if( ofs.is_open() )
ofs << printBuf;
}

On Wed, 30 Aug 2023 13:01:16 +0200
Bonita Montero <Bonita.Montero@gmail.com> wrote:
> if( argc < 3 || *argv[2] )
> ofs.exceptions( ofstream::failbit | ofstream::badbit ),
> ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc ),
> printBuf.reserve( BUF_SIZE + 31 );

Commas as statement seperators instead of using a bracketed block? Are you
posting this code for lols?

Am 30.08.2023 um 17:38 schrieb Muttley@dastardlyhq.com:
> On Wed, 30 Aug 2023 13:01:16 +0200
> Bonita Montero <Bonita.Montero@gmail.com> wrote:
>> if( argc < 3 || *argv[2] )
>> ofs.exceptions( ofstream::failbit | ofstream::badbit ),
>> ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc ),
>> printBuf.reserve( BUF_SIZE + 31 );

> Commas as statement seperators instead of using a bracketed block?
> Are you posting this code for lols?

I often concatenate a small number of statements with the comma-operator
to save screen space. I never had any issues with that. That's just a
matter of personal taste. If you're noticing this until now, I'm going
to assume that you haven't really looked at the previous code.
But that's just the syntax issue. My code is also superior in terms of
safety (iterator debuggin) and performance (+2/3 more performance with
g++ 12 and a Sklake-CPU).

Am 29.08.2023 um 20:34 schrieb Bonita Montero:

> size_t g = countr_zero( pSieve->w & (word_t)-1 << p % BITS );

size_t g = countr_zero( (word_t)(pSieve->w & (word_t)-1 << p % BITS) );

OMG, an integral promotion error.
If word_t is defined < unsigned the results are wrong.

Ben Bacarisse <ben.usenet@bsb.me.uk> writes:

> Tim Rentsch <tr.17687@z991.linuxsc.com> writes:
>
>> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>>
>>> I wanted to try a newer algorithm, so I coded up a rather naive
>>> version of Atkin & Bernstein's sieve in C. It needs some tidying
>>> and I imagine a bit can be shaved off the time.
>>
>> Thank you for this. Would you mind posting your code?
>
> Hmm... It's a bit of a mess... but sure. [...]

It looks quite okay. I haven't yet had time to look at it
closely but I look forward to doing so.

>> Also do
>> you have a reference for the Atkin/Bernstein algorithm (or even
>> better, code?)?
>
> Well, there's this: https://en.wikipedia.org/wiki/Sieve_of_Atkin

Excellent! And this page has a link to DJ Bernstein's fast
primes code, which is stunningly swift.

Thank you for both.

Tim Rentsch <tr.17687@z991.linuxsc.com> writes:

> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>
>> Tim Rentsch <tr.17687@z991.linuxsc.com> writes:
>>
>>> Ben Bacarisse <ben.usenet@bsb.me.uk> writes:
>>>
>>>> I wanted to try a newer algorithm, so I coded up a rather naive
>>>> version of Atkin & Bernstein's sieve in C. It needs some tidying
>>>> and I imagine a bit can be shaved off the time.
>>>
>>> Thank you for this. Would you mind posting your code?
>>
>> Hmm... It's a bit of a mess... but sure. [...]
>
> It looks quite okay. I haven't yet had time to look at it
> closely but I look forward to doing so.

About 20% of the time was spent counting the primes, so I added a
population count function to the bitset operations. And there's an
obvious halving one can do on the upper bound of the first loop. That
brought the time down to about 28 seconds.

>>> Also do
>>> you have a reference for the Atkin/Bernstein algorithm (or even
>>> better, code?)?
>>
>> Well, there's this: https://en.wikipedia.org/wiki/Sieve_of_Atkin
>
> Excellent! And this page has a link to DJ Bernstein's fast
> primes code, which is stunningly swift.

I haven't looked at it yet. It would probably have put me off!

> Thank you for both.

You're welcome.

--
Ben.

On Wed, 30 Aug 2023 18:31:32 +0200
Bonita Montero <Bonita.Montero@gmail.com> wrote:
>Am 30.08.2023 um 17:38 schrieb Muttley@dastardlyhq.com:
>> On Wed, 30 Aug 2023 13:01:16 +0200
>> Bonita Montero <Bonita.Montero@gmail.com> wrote:
>>> if( argc < 3 || *argv[2] )
>>> ofs.exceptions( ofstream::failbit | ofstream::badbit ),
>>> ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc ),
>>> printBuf.reserve( BUF_SIZE + 31 );
>
>> Commas as statement seperators instead of using a bracketed block?
>> Are you posting this code for lols?
>
>I often concatenate a small number of statements with the comma-operator
>to save screen space. I never had any issues with that. That's just a

To save space? Seriously??

>matter of personal taste. If you're noticing this until now, I'm going

Its not just personal taste, it makes the code harder to parse for anyone
who has to maintain it and could easily lead to logical errors if someone
modifies it and forgets to replace the ; with a , at the end.

>to assume that you haven't really looked at the previous code.
>But that's just the syntax issue. My code is also superior in terms of
>safety (iterator debuggin) and performance (+2/3 more performance with
>g++ 12 and a Sklake-CPU).

Yes, we know how wonderful your code is after you seem to spend multiple days
on some trivial task. Obviously no one else could ever manage what you did,
you're a coding god.

Am 31.08.2023 um 10:22 schrieb Muttley@dastardlyhq.com:

> To save space? Seriously??

It's my taste.

> Its not just personal taste, it makes the code harder to parse for anyone
> who has to maintain ...

Not true, you just feel unsecure with that for whaterver reason.

> Yes, we know how wonderful your code is after you seem to spend multiple days
> on some trivial task. ...

The task is not trivial if you consider my solution.
You won't manage to have it that fast.

Am 31.08.2023 um 10:22 schrieb Muttley@dastardlyhq.com:

> Yes, we know how wonderful your code is after you seem to spend multiple days
> on some trivial task. Obviously no one else could ever manage what you did,
> you're a coding god.

I just benchmarked my code with an upper bound of 1E3 against your code
on a Skylake CPU with g++ 12 and -march=native. My code runs about 0.9s,
your code takes 5.26s. Any questions ?

Muttley@dastardlyhq.com writes:
>On Wed, 30 Aug 2023 18:31:32 +0200
>Bonita Montero <Bonita.Montero@gmail.com> wrote:
>>Am 30.08.2023 um 17:38 schrieb Muttley@dastardlyhq.com:
>>> On Wed, 30 Aug 2023 13:01:16 +0200
>>> Bonita Montero <Bonita.Montero@gmail.com> wrote:
>>>> if( argc < 3 || *argv[2] )
>>>> ofs.exceptions( ofstream::failbit | ofstream::badbit ),
>>>> ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc ),
>>>> printBuf.reserve( BUF_SIZE + 31 );
>>
>>> Commas as statement seperators instead of using a bracketed block?
>>> Are you posting this code for lols?
>>
>>I often concatenate a small number of statements with the comma-operator
>>to save screen space. I never had any issues with that. That's just a
>
>To save space? Seriously??

cfront used to generate extensive comma laden expressions. It did
cause problems with PCC's temporary register allocation algorithms
that I had to fix back around 1989.

Am 31.08.2023 um 16:43 schrieb Scott Lurndal:

> cfront used to generate extensive comma laden expressions. It did
> cause problems with PCC's temporary register allocation algorithms
> that I had to fix back around 1989.

There's for sure no compiler that generates different code whether
you youse individual ;-separated statements or one comma-separated
statement as both are sequence points.

On Thu, 31 Aug 2023 11:52:09 +0200
Bonita Montero <Bonita.Montero@gmail.com> wrote:
>Am 31.08.2023 um 10:22 schrieb Muttley@dastardlyhq.com:
>
>> Yes, we know how wonderful your code is after you seem to spend multiple days
>
>> on some trivial task. Obviously no one else could ever manage what you did,
>> you're a coding god.
>
>I just benchmarked my code with an upper bound of 1E3 against your code
>on a Skylake CPU with g++ 12 and -march=native. My code runs about 0.9s,
>your code takes 5.26s. Any questions ?

Yes , do you have a life? My code took about 20 mins to write , you seem to
have been working on yours for over a week now.

On Thu, 31 Aug 2023 11:00:06 +0200
Bonita Montero <Bonita.Montero@gmail.com> wrote:
>Am 31.08.2023 um 10:22 schrieb Muttley@dastardlyhq.com:
>
>> To save space? Seriously??
>
>It's my taste.
>
>> Its not just personal taste, it makes the code harder to parse for anyone
>> who has to maintain ...
>
>Not true, you just feel unsecure with that for whaterver reason.

You sound like the python devs who claim there's no problem with significant
whitespace with indentation.

On Thu, 31 Aug 2023 14:43:34 GMT
scott@slp53.sl.home (Scott Lurndal) wrote:
>Muttley@dastardlyhq.com writes:
>>On Wed, 30 Aug 2023 18:31:32 +0200
>>Bonita Montero <Bonita.Montero@gmail.com> wrote:
>>>Am 30.08.2023 um 17:38 schrieb Muttley@dastardlyhq.com:
>>>> On Wed, 30 Aug 2023 13:01:16 +0200
>>>> Bonita Montero <Bonita.Montero@gmail.com> wrote:
>>>>> if( argc < 3 || *argv[2] )
>>>>> ofs.exceptions( ofstream::failbit | ofstream::badbit ),
>>>>> ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::trunc ),
>>>>> printBuf.reserve( BUF_SIZE + 31 );
>>>
>>>> Commas as statement seperators instead of using a bracketed block?
>>>> Are you posting this code for lols?
>>>
>>>I often concatenate a small number of statements with the comma-operator
>>>to save screen space. I never had any issues with that. That's just a
>>
>>To save space? Seriously??
>
>cfront used to generate extensive comma laden expressions. It did

I don't even understand why C/C++ allows comma seperation of calls in a block.
It should be restricted to a few specific locations such as inside for().

Am 31.08.2023 um 17:23 schrieb Muttley@dastardlyhq.com:

> You sound like the python devs who claim there's no problem with significant
> whitespace with indentation.

You're too compulsive for that world.

Am 31.08.2023 um 17:22 schrieb Muttley@dastardlyhq.com:
> On Thu, 31 Aug 2023 11:52:09 +0200
> Bonita Montero <Bonita.Montero@gmail.com> wrote:
>> Am 31.08.2023 um 10:22 schrieb Muttley@dastardlyhq.com:
>>
>>> Yes, we know how wonderful your code is after you seem to spend multiple days
>>
>>> on some trivial task. Obviously no one else could ever manage what you did,
>>> you're a coding god.
>>
>> I just benchmarked my code with an upper bound of 1E3 against your code
>> on a Skylake CPU with g++ 12 and -march=native. My code runs about 0.9s,
>> your code takes 5.26s. Any questions ?
>
> Yes , do you have a life? My code took about 20 mins to write , you seem to
> have been working on yours for over a week now.

You complained about my code and I explained why it was more complex.
You probably have nothing more to say about that. With the questions
you usually ask about C++, I assume that you would not be able to
write such secure and performance code.
Writing performant code is part of my job. I write software for high
frequency trading, so you have to think about which constructs result
in which code and how efficient it is.