# Improved Sieve of Eratosthenes

How can I get to a perfect coding for this algorithm?
I made the mathematical theorem, which is a development of the Sieve of Eratosthenes algorithm. I might use some help in coding.

You can find the description of the code at my paper:
Development of Sieve of Eratosthenes and Sieve of Sundaram's proof

For more understanding you can check this paper: SEQUENCE ELIMINATION FUNCTION AND THE FORMULAS OF PRIME NUMBERS

For the next development see Next level Improved Sieve of Eratosthenes

#include <iostream>
#include <math.h>
#include <cstdlib>
using namespace std;
int D2SOE(int n1_m) {
int n1 = 0, g = 0, z = 0, p = 0,f1=0,f3=0;
cout << "\n2 3 ";
bool* array = new bool[n1_m];
// Initialising the D2SOE array with false values
for (int i = 0; i < n1_m; i++)
array[i] = false;
// The main elimination theorem
for (n1=1;n1<=ceil((sqrt(2*floor((3*n1_m+1)/2.0)+1)-2)/3.0);n1++)
{ if (array[n1] != 0)
continue;
z = ((3 * n1 + 1) / 2.0);
p = ((2 * z) + 1);
f1 = (ceil(((7*p)-2) / 3.0) - ceil((p - 2) / 3.0));
f3 = (ceil(((7*p)-2) / 3.0) - ceil(((5 * p) - 2) / 3.0));
for (g = ceil(((4*((z*z) + z)) - 1) / 3.0); g < n1_m;g+=f1)
{ array[g] = true; }
if ((p +1) % 3 == 0) {
for (g = ceil(((4*((z*z) + z)) - 1) / 3.0) + f3; g < n1_m;g+=f1)
{ array[g] = true; } }
else {
for (g = ceil(((4*((z*z) + z)) - 1) / 3.0) - f3; g<n1_m;g+=f1)
{ array[g] = true; }}}
// printing for loop
for (int n1 = 1; n1 < n1_m; n1++)
if (!array[n1]) {
z = (((3 * n1) + 1) / 2.0);
cout << (2 * z) + 1 << " "; }
return 0; }
// driver program to test above
int main() {
int n1_m, N = 0;
cout << "\n Enter limit : ";
cin >> N;
n1_m = ceil((N - 2) / 3.0);
D2SOE(n1_m);
return 0; }

• Comments are not for extended discussion; this conversation has been moved to chat.
– Mast
Aug 25, 2022 at 13:10
• In practice, yours is a base {2,3} wheel sieve so we have p=r+6*j with r=1 or r=-1 and j>0 for find the prime numbers greater than 3 and then use boolean array of size 2*n_limit/6. Here gist.github.com/user140242/9d142187b89d66a105ccc87ea292bd59 you will find a similar description of the sieve you used and here gist.github.com/user140242/ed2f0b8b93e0257ffd8e5f9bafdd20f7 a segmented example Sep 8, 2022 at 10:10
• You can write the algorithm more compactly by simplifying all the equations. If you look at the following algorithm, it uses almost the same procedure but uses two arrays of size N/6 instead of one of N/3. You can find a version here and a more understandable Python version in the first part here Sep 8, 2022 at 12:59

• Separation of concerns

We dont't sieve primes for the sake of sieving primes. We want primes because we want to do something with them, not just print. Let D2SOE return the array it computed, and

  int main()
int n1_m, N = 0;
cout << "\n Enter limit : ";
cin >> N;
n1_m = ceil((N - 2) / 3.0);
bool array = D2SOE(n1_m);
print_primes(array, n1_m);
return 0;
}

• Overall impression: Unreadable.

• Please, don't }}}. Indent your code properly.
• ceil(((4*((z*z) + z)) - 1) / 3.0) seems very important, as it is repeated 3 times. Figure out a good name for it, and compute it once.
• f1 and f3 deserve better names too.
• After sieving, I'd expect true for primes.
• Correctness Didn't check it. However,

• A floating point math in an elementary number-theoretical problem is totally out of place. It may bite you hard when n1_m1 grows large enough.
• Comparing a boolean array[n1] to 0, while technically correct, gives a yucky taste.
• I don't see why to #include <cstdlib>.
• using namespace std; is always wrong.
• Don't new. Use std::vector.
• I think unreadibility is the most serious issue. Also, the comments do not help. Aug 24, 2022 at 14:49
• I think the }}} stuff is probably an attempt to make every line follow the correct indenting, however it's not a good fit for C++ syntax. Aug 24, 2022 at 16:33
• The }}} reminds me of usual Lisp/Scheme style. Aug 24, 2022 at 16:48
• I respectfully disagree with "using namespace std; is always wrong". It depends a lot on the size and domain of a program. Small-ish, self-contained programs, like this one, benefit from it. Scattering the program with std:: makes code less readable and more redundant. Good compilers warn if and when you accidentally shadow names from std. Aug 26, 2022 at 11:23
• @Peter-ReinstateMonica but that's a distinction lost on novices. Instead, using namespace std; just becomes a standard bit of boilerplate that is always at the top of every program they wrote when they were learning to code, and eventually finds its way to the top of the header file for a library, and now some other poor soul has to deal with those shadowing compiler warnings. "It's always wrong" is a better lesson at this stage. Aug 26, 2022 at 21:38

Sorry, this code is utterly unintelligible, and I wouldn't accept it in any code-base I maintain.

• Why are we using <math.h> rather than <cmath>?
• Avoid using namespace at global scope.
• All the identifiers (including the function name D2SOE()) are completely meaningless. They may mean something if you've read the paper linked from the question, but that's not even mentioned in the code comments, so how is anyone meant to know that?
• The code layout is unlike any I've ever seen, which makes it harder to read. Use an accepted whitespace style to more effectively convey your meaning to other coders.
• Why does the function accept a signed value as argument? What does it mean to pass a negative value? It looks to me that it will likely throw a std::bad_alloc if you do that.
• And why does it return an integer? What's the significance of the return value?
• Why are we allocating a raw array for storage, rather than a C++ container? Where does it get deallocated? Yes, you have a memory leak there. N.B. beware of std::vector<bool> - it doesn't behave as a standard container, and you may consider std::vector<char> a good alternative for storing booleans; alternatively, make sure you know how to use std::vector<bool> safely.
• Why do we have a loop to clear the array, rather than using std::fill()? We don't need this anyway, if we value-initialise the array by adding () to the new expression (or by moving to a real container).
• There's lots of implicit conversions from double to int that trigger compiler warnings. Avoid writing code that generates warnings, because they all need to be examined to give any confidence, and screeds of warnings can cause readers to skim-read, and potentially overlook serious problems.
• For positive integers that can be exactly represented as double, static_cast<int>(ceil((n-2) / 3.0)) is the same as n / 3. There's no need to use floating-point there. Similarly, other ceiling divisions can be simply converted to unsigned-integer (floor) divisions.
• Writing to standard output stream isn't very useful for a function you might want to make use of in a real program - better to return the primes instead.
• We're completely missing any checking that we could read N from std::cin. If that fails for any reason (stream closed, or user entered non-numeric text), we probably don't want to proceed.
• Always include your unit-tests when posting code for review.

Some thoughts following 2 other good reviews.

Avoid floating point math for an integer problem

Using size_t for array sizing.

Consider only using unsigned types for this task.

Code gets more range and less UB potential. Risks: code needs care to avoid making a "negative".

Big bool?

As bool might be more than 1 byte and bool[n1_m] could be a very large array, consider unsigned char[n1_m].

array in bool* array does not convey much info. Further: array is not an array here. It is a pointer.

Rather than use a name that only describes the type, consider a name the IDs what the object represents, like sieve for array.

Improve formatting

Use an auto-formatter to reform into a more common formatting style.

Avoid naked acronyms

In code, detail what D2SOE implies.

Long output

stdout can have environmental limitations on line length. Perhaps print a '\n' once in a while. Even better, do not print, but let the caller display the result list.

Minor corner case

When N < 3, code still prints 3.

Efficiency ideas:

For n > 30, there are at most 8 primes for each step of 30 (note 2*3*5=30). I found flagging the primes [30-59], [60-89], [90-119], ... into a byte made for efficient packing and the masking to set/get not too bad.

• You say bool might be more than 1 byte, but char suffers from the same problem and can be more than 8 bits wide. Aug 26, 2022 at 21:32
• @Nayuki a char is always one byte, by definition. One byte can be more than 8 bits, indeed. But there is no type in C++ that is smaller than a char. Aug 28, 2022 at 1:54
• @Nayuki A char is always 1 "byte" be that an 8-bit byte or more.. A bool may be a narrow as a char or multiple bytes, like an int. It is implementation defined. Aug 28, 2022 at 2:52

## 2.0 float or not 2 float

Don't use floats when you don't need to. Floating point operations are slower than integer operations and avoiding them can give you some worthwhile performance gains. For example the second snippet runs ~40% faster

Looking at the main loop, for all but the one square root (done only once) all the math can be integer math.

for (n1=1;n1<=ceil((sqrt(2*floor((3*n1_m+1)/2.0)+1)-2)/3.0);n1++)
{ if (array[n1] != 0)
continue;
z = ((3 * n1 + 1) / 2.0);
p = ((2 * z) + 1);
f1 = (ceil(((7*p)-2) / 3.0) - ceil((p - 2) / 3.0));
f3 = (ceil(((7*p)-2) / 3.0) - ceil(((5 * p) - 2) / 3.0));
for (g = ceil(((4*((z*z) + z)) - 1) / 3.0); g < n1_m;g+=f1)
{ array[g] = true; }
if ((p +1) % 3 == 0) {
for (g = ceil(((4*((z*z) + z)) - 1) / 3.0) + f3; g < n1_m;g+=f1)
{ array[g] = true; } }
else {
for (g = ceil(((4*((z*z) + z)) - 1) / 3.0) - f3; g<n1_m;g+=f1)
{ array[g] = true; }}}


you use float almost in every line. You only need a float for the sqrt.

Using unsigned ints you can use the following

result.push_back(2);
result.push_back(3);
const unsigned n = sqrt(((3 * n1_m + 1) & UINT_BOTTOM_BIT_MASK) + 1) / 3;
for (unsigned n1 = 1; n1 <= n; n1++) {
if (primeRoots[n1]) {
const unsigned z  = (3 * n1 + 1) / 2;
const unsigned p  = 2 * z + 1;
const unsigned f1 = 7 * p / 3 - p / 3;
const unsigned f3 = 7 * p / 3 - 5 * p / 3;
const unsigned start  = (4 * (z * z + z) + 1) / 3;
const unsigned start1 = (p + 1) % 3 ? start - f3 : start + f3;
for (unsigned g = start; g < n1_m; g += f1) { primeRoots[g] = 0; }
for (unsigned g = start1; g < n1_m; g += f1) { primeRoots[g] = 0; }
}
}
for (unsigned n1 = 1; n1 < n1_m; n1++) {
if (primeRoots[n1]) { result.push_back((((3 * n1) + 1) & UINT_BOTTOM_BIT_MASK) + 1); }
}


Note that the & UINT_BOTTOM_BIT_MASK remove the lowest bit and is the same as integer math / 2) * 2. The mask is defined as constexpr unsigned UINT_BOTTOM_BIT_MASK{UINT_MAX - 1};

## Additional integer optimisations.

The values of f1 and f3 can be simplified from

const unsigned f1 = 7 * p / 3 - p / 3;     // (7p-p)/3 = (6p)/3 = 2p

const unsigned f3 = 7 * p / 3 - 5 * p / 3; // (7p-5p)/3 = (2p)/3


to

const unsigned f1 = 2 * p;
const unsigned f3 = (2 * p + 1) / 3; // the + 1 to adjust flooring of integer


or remove the /3 with a /2 (shift >> 1)

const unsigned f1 = 2 * p;
const unsigned f3 = (p + n1 + 1) / 2;


Note that multiplying or dividing by power of 2 will use shift operators when working with integers.

• You also don't need to use sqrt(); you just need to find a value such that if you square it, it is equal to or larger than 2*floor((3*n1_m+1)/2.0)+1. Since n1_m only increases by one each iteration of the loop, then either the previous square root value is still valid or you need to increase that one by 1 or 2 at most. Aug 24, 2022 at 17:41
• Perhaps revised code's 2 intense loops for (unsigned g = start; g < n1_m; g += f1) { primeRoots[g] = 0; } for (unsigned g = start1; g < n1_m; g += f1) { primeRoots[g] = 0; } could be merged into 1: int d = start1-start; for (unsigned g = start; g < n1_m; g += f1) { primeRoots[g] = 0; primeRoots[g+d] = 0;} or the like. May need some edge case proofing. Aug 25, 2022 at 8:23
• @G.Sliepen Sorry my bad, the answer is targeting performance, I did not make that clear. The sqrt is insignificant in regards to n (n1_m = n / 3). Moving additional logic into the loop to avoid the sqrt call would quickly consume any benefit. Aug 25, 2022 at 20:52
• I'm not saying you should put that logic inside the inner loop, you only need it in the outer loop. But either way you are right that it becomes insignificant for large values of n1_m. On the other hand, OP asked for the "perfect coding"? Aug 25, 2022 at 21:23
• I used floating points because it was giving me wrong results, I don't know why! Aug 25, 2022 at 21:59

You can't ask for "perfect". You can ask for "lets me identify many primes, in a very short time".

First, you need to have different implementations and measure.

Write an inline function "bit_is_set" and one "bit_is_cleared". Then try different implementations. Implement the sieve as a bit array stored in an array of 8 bit unsigned int, or a bit array stored in an array of 32 or 64 bit unsigned int, or an array of unsigned char, or an array of bool as you did. Using these inline functions that is an implementation detail that can be hidden completely or almost completely.

Decide for which numbers you store a boolean. Not for all. Maybe for odd ones, or for ones not divisible by 2 or 3, or 2 or 3 or 5, or even 2 or 3 or 5 or 7. Only 48 out of 210 consecutive integers are not divisible by 2, 3, 5, or 7. 8 out of 30 are not divisible by 2, 3 or 5. The effect is that you may need less storage which is faster, and that you set fewer bits which is faster.

Try to write functions that return true if a number is prime, or return the next prime above a number, or the last prime before a number. That makes it a lot more useful.

Try to arrange things in a way that you can handle a sieve that doesn't fit into memory. For example, I'd like to examine all primes up to 10^18 for some project, and there's no way this fits into RAM.

When you run the sieve, take care of cache sizes. L1 cache is faster than L2 cache which is faster than L3 cache which is faster than RAM. So you will do one pass where you only remove primes in a range of say 25 KB which is faster because everything is done in L1 cache. Then in a range of 200 KB which fits into L2 cache and so on. Your code will run a lot faster.

That's just a few tricks to make it faster.

• The standard way to scale up sieves is to use segmented sieves. However even storing all the primes up to 10^18 is a monumental task and I don't know of an algorithm to actually generate all of them given the Prime Number Theorem gives a very good approximation of how many would have to be generated (there are sub-linear ways to count primes much faster than generating them)
– qwr
Aug 26, 2022 at 20:40

Perfect coding of a given algorithm sounds almost achievable compared to proving an algorithm optimal for a given problem. (Which is why rate of growth is an accepted criterion.)

Problems start with defining quality measures for source code, finally a single figure of merit.
I expect readability to win with (maintenance?) programmers, followed by the wider maintainability.
Usability/documentation is important with "application" programmers (can I be sure it solves the problem, do I use it right?), resource consumption, stability.

So, if development two of the Sieve of Eratosthenes was to be used as a building block of something larger (say, factorisation), it dearly needed documentation.
Documentation and source code do not get separated when the former is embedded in the latter:
Have a look at tools like doxygen.

I take readability to be about don't make me think:

• stick to established practices
- code layout and formatting
• use suggestive names
• separate concerns
• comment what's not "obvious" (from the source code)
• where in doubt, code the way you think about solution&problem
(expect this to change over time)

There is a lot of ceil((<expression>-2) / 3.0):
give such a name (index()?), make it an inline function.
To keep details of keeping compound marks out of identify primes, define&use bool marked_compound(int i) & void mark_compound(int i).
There are three "mark loops" where one should do.
(I failed to find suggestive names/figure out what things were meant for more often than not.)

    int delta = index(7*p) - index(5*p);
if (p % 3 != 2)
delta = -delta;
for (int g = ceil(((4 * ((z * z) + z)) - 1) / 3.0) ;
g < index_limit ; g += f1) {
mark_compound(g);
mark_compound(g + delta);
}


(ceil(((4 * ((z * z) + z)) - 1) / 3.0) might be index(4 * ((z * z) + z)) + 1))

Is your sieve actually better?

I modified your code to at least count the number of primes instead of outputting and it appears correct for N=10^8, 10^9 (see table). There are a lot of floating point calculations going on, compared to the standard Sieve of Eratosthenes, and both still fit the whole array into memory for N <= 2*10^9. Here I present the standard sieve, with the common optimization being just skipping even numbers instead of optimizing 2 and 3. But what I did optimize is use half the array space, so if it is possible to apply this optimization to your version, you definitely should to fit more values into cache.

Disclaimer: I only tested for round values like powers of 10. I haven't thoroughly checked edge cases for N.

int sieve(int n)
{
std::vector<bool> nums(n/2+1, 0);
for (int i=3; i*i <= n; i+=2)
{
if (nums[i/2] == 0)
{
for (int j=i*i; j<n; j+=2*i)
nums[j/2] = 1;
}
}

int s = 1;  // include 2 as prime
for (int i=3; i<n; i+= 2)
{
if (nums[i/2] == 0) ++s;
}
return s;
}


My simple benchmarks (Ubuntu 20.04, g++ 9.4 with -O3, Intel i7-7700HQ) were computed using time:

• N=10^8: d2soe 0.32s, skip2: 0.23s
• N=10^9: d2soe 3.8s, skip2: 3.8s

Testing N=10^10 and beyond will require rewriting the sieve to be segmented.

I was able to improve your N=10^9 from 3.8s to 3.1s just be replacing the C-style bool array with vector<bool> (which implementations will usually pack 8 bits per byte instead of just 1 with a bool array) so that's definitely an easy optimization you should make.

• I imagine (no proof) that this shorter code benefits from better CPU pipelining and branch prediction compared to the original. Is there a way to actually measure that? Aug 25, 2022 at 5:37
• I think so but I couldn't tell you how. Also there is almost no math computation being done here, pretty much addition and memory access. As a general principle, compilers can optimize simple code quite well.
– qwr
Aug 25, 2022 at 5:57
• Somehow for (int j=i*i; j<n; j+=2*i) nums[j/2] = 1; looks like it can be converted to the form something like for (int j_div_2=i*i/2; j_div_2<n_div_2; j_div_2+=i) nums[j] = 1;. Thus simplifying the innermost loop. Aug 25, 2022 at 8:30
• You're also using a packed btimap (std::vector<bool>), not 1 byte per sieve entry with bool* array. That saves on memory bandwidth in the early stages with small strides (more locality between crossing off), as you have 512 entries per cache line instead of 64. (Larger primes will still stride by more than 1 whole cache lines. In some size ranges it doesn't matter too much whether that's every 5 cache lines or every 40, although HW prefetch probably works better with smaller strides. And for very large strides you start getting more TLB misses with 8x the footprint.) Aug 25, 2022 at 11:28
• @chux-ReinstateMonica you're probably right but I when I threw this together I didn't want to complicate anything. Maybe the compiler already figured this out!
– qwr
Aug 25, 2022 at 15:57