1

I'm currently trying to implement the Sieve of Eratosthenes in C using a BitSet, but I get a segmentation fault, when I try to sieve the primes up to 1,000,000 (1 million) - 100,000 (100 thousand) is still working though and I can't figure out why I get the seg-fault.

This is the code I use (I marked the line, in which the error occurs):

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>

void eSieve(uint64_t upperLimit);
int main(int argc, char *argv[]) {
  uint64_t upperLimit;

  if (argc == 2) {
    upperLimit = (uint64_t) atoll(argv[1]);
    printf("Using custom limit: %" PRIu64 "\n", upperLimit);
  } else {
    upperLimit = 1000;
    printf("Using default limit: %" PRIu64 "\n", upperLimit);
  }

  eSieve(upperLimit);

  return 0;
}

typedef uint32_t prime_t;

void eSieve(uint64_t upperLimit) {
  if (upperLimit < 2) {
    printf("FAILURE: Bad upper limit.\n");
    return;
  }

  prime_t *sieve = calloc(1, (upperLimit + sizeof(prime_t) - 1)/sizeof(prime_t));

  if (!sieve) {
    printf("FAILURE: Could not initialize sieve.\n");
    return;
  }

  sieve[0] |= 3;    // Set first and second bit (representing 0 and 1)

  uint64_t prime, number;
  for (prime = 2; prime * prime < upperLimit; ) {
    for (number = prime * prime; number < upperLimit; number += prime) {
      // Segmentation fault for prime = 2 and number = 258048
      sieve[number/sizeof(prime_t)] |= (((prime_t) 1) << (number % sizeof(prime_t)));
    }

    while ((sieve[++prime/sizeof(prime_t)] & (prime_t)1 << (prime % sizeof(prime_t))) != 0)
      ;
  }

  number = upperLimit;
  while ((sieve[--number/sizeof(prime_t)] & (((prime_t)1) << (number % sizeof(prime_t)))) != 0)
    ;

  printf("Greatest prime-number below %" PRIu64 ": %" PRIu64 "\n", 
      upperLimit, number);
}

Does anybody know why the error occurs? I'm guessing that now enough space is allocated (somehow), but I can't see how this would be possible at the moment...

Matthias
  • 175
  • 6
  • 1
    *Where* in your code do you get SIGSEGV? – Andrew Henle Jun 29 '16 at 13:35
  • @AndrewHenle I wrote a comment. – Matthias Jun 29 '16 at 13:37
  • You say you're using a *bit* array, but you seem actually to be using a *byte* array. Maybe. Your code is a bit hard to follow. If you indeed intend to use bit arrays, then it would be better to factor out your test and set operations into macros. – John Bollinger Jun 29 '16 at 13:46
  • I don't know too much about what your code is doing, but the reason it crashed is you're trying to access elements beyond the size of array sieve – Ishay Peled Jun 29 '16 at 13:47
  • @MichaelWalz It is provided minimal, complete and verifiable. I just ran his code – Ishay Peled Jun 29 '16 at 13:48
  • 1
    `prime * prime` is apparently a larger value than `upperLimit`, simple as that? Use a debugger and check. – Lundin Jun 29 '16 at 13:54
  • 1
    I'm betting that `prime * prime` overflows 64 bits. – barak manos Jun 29 '16 at 13:55
  • @Lundin As I wrote in the comment, this occurs for prime == 2. – Matthias Jun 29 '16 at 13:55
  • @barakmanos As I wrote in the comment, the error occurs for prime == 2. – Matthias Jun 29 '16 at 13:56
  • @Matthias it's an index out of bounds problem when accessing `sieve[x]`. You can easily check this by putting some `assert`s in your code. – Jabberwocky Jun 29 '16 at 13:58
  • @MichaelWalz That was my guess as well, but I was hoping for an explanation. – Matthias Jun 29 '16 at 13:59
  • @Matthias that's easy to find out yourself. Just check if the index is within bounds before each access (just for debugging purposes). – Jabberwocky Jun 29 '16 at 15:57
  • The posted code causes the compiler `gcc` to output two warning messages about 1) the `while` loop with the `++prime` operation and 2) the `while` loop with the `--number` operation. Strongly suggest you fix those two problems. – user3629249 Jun 29 '16 at 22:32
  • the posted code is allocating memory from the `heap` via the call to `calloc()`, but failing to pass the pointer returned from `calloc()` to the function: `free()`. – user3629249 Jun 29 '16 at 22:40

2 Answers2

4

You're not getting the correct bit number:

sieve[number/sizeof(prime_t)] |= (((prime_t) 1) << (number % sizeof(prime_t)));

When you do the division and mod, you need to divide/mod by the number of bits, not the number of bytes:

sieve[number/(sizeof(prime_t)*8)] |= (((prime_t) 1) << (number % (sizeof(prime_t)*8)));

And similarly:

while ((sieve[++prime/(sizeof(prime_t)*8)] & (prime_t)1 << (prime % (sizeof(prime_t)*8))) != 0)

...

while ((sieve[--number/(sizeof(prime_t)*8)] & (((prime_t)1) << (number % (sizeof(prime_t)*8)))) != 0)

EDIT:

You're also not allocating the right amount of memory. You need a number of bytes equal to the limit divided by the number of bits, plus 1 sizeof(prime_t) to round up.

prime_t *sieve = calloc(1, (upperLimit / 8) + sizeof(prime_t));

As it right now, you're allocating twice the bytes you need.

Also, if you want to defend against cases where there are more or less than 8 bits to a byte, use CHAR_BIT in the above code in place of 8. Whatever sizeof(uint64_t) evaluates to shouldn't matter, as you'll still get the proper number of bits required.

dbush
  • 205,898
  • 23
  • 218
  • 273
  • @Matthias, no, you allocated more bytes than you need by a factor of `CHAR_BIT` (which is normally 8). – John Bollinger Jun 29 '16 at 14:06
  • @JohnBollinger Thank you really much! – Matthias Jun 29 '16 at 14:06
  • @Lundin Do you know a readable line of code to do the some (not wasting 8x memory)? – Matthias Jun 29 '16 at 14:15
  • @JohnBollinger Quick question: What happens on exotic machines, where *CHAR_BIT* is not 8? What does `sizeof(uint64_t)` return? Because it then can't return an even number like 8... What is the better way to allocate (`number/(sizeof(prime_t)*8)` or `#define BITS_IN_PRIME_T 64` and then `number/BITS_IN_PRIME_T`) – Matthias Jun 29 '16 at 14:18
  • What will `sizeof(uint64_t)` return with `CHAR_BIT` != 8? – Matthias Jun 29 '16 at 14:23
  • @Matthias every type's representation consists of an integral number of units comprising `CHAR_BIT` bits each. However, these bits may include *padding* bits. If `CHAR_BIT` were, say, 12, and the system provided `uint64_t` at all, then `sizeof(uint64_t)` would be an integer greater than or equal to 6, and the representation of `uint64_t` would include some padding bits. – John Bollinger Jun 29 '16 at 14:24
  • @JohnBollinger Thank you for the clarification! I leant really much in this question :) – Matthias Jun 29 '16 at 14:26
  • @Mattias, so in fact, I'm not quite correct to say that you're off by a factor of `CHAR_BIT` instead of 8, and dbush is not quite correct to use `sizeof(prime_t)*8`. For elements of type `uint32_t`, one ought to use plain 32. But in practice, you are unlikely to run into a machine where it makes a difference. – John Bollinger Jun 29 '16 at 14:27
  • I'm still curious about the allocation: Lets consider `upperLimit` to be 64. Then we allocate `(upperLimit + 1) * sizeof(prime_t) / 8` = `65 * 8/8` = 65 bytes, which seems wrong to me (Shouldn't we want to allocate 8 or 16 bytes?). – Matthias Jun 29 '16 at 14:50
  • @Matthias Good catch. I edited the part on how much to allocate. – dbush Jun 29 '16 at 15:01
  • Your current edit won't always be a multiple of `sizeof(prime_t)` - consider `upperLimit` to be *56* your current code will resolve in 56/8 + 8 = 15 which isn't a multiple of `sizeof(prime_t)` and therefore won't return a valid pointer to a `prime_t`-array. This is what I'm thinking of - next comment (Please bear with me): – Matthias Jun 29 '16 at 15:24
  • @Matthias That's ok, it doesn't need to be a multiple. You'll get at most `sizeof(prime_t)` more bytes than you need, which is negligible. – dbush Jun 29 '16 at 15:26
  • `((upperLimit/CHAR_BIT + sizeof(prime_t) - 1)/sizeof(prime_t)) * sizeof(prime_t)`. I'm thinking this way: The final result needs to be a multiple of `sizeof(prime_t)` - therefore the trailing `* sizeof(prime_t)`. The question is how many of this `sizeof(prime_t)` sized blocks will be needed. We have `upperLimit` numbers so we'll need *at least* `upperLimit/CHAR_BIT` bytes. Since we want not the number of bytes, but the number of `sizeof(prime_t)` sized blocks, I need to divide by `sizeof(prime_t)`. But if `upperLimit/CHAR_BIT` is less than `sizeof(prime_t)` the expression will be round down. – Matthias Jun 29 '16 at 15:27
  • Therefore I add `(sizeof(prime_t) - 1)` to be sure, that always enough space is alloced. You may wonder why I end with `.../sizeof(prime_t)) * sizeof(prime_t)`, which can see pretty much like `* 1/1`, but the parentheses are important - the parentheses are used to floor. I really hope this makes sense. – Matthias Jun 29 '16 at 15:29
  • But when only 15 bytes are allocated and I access `sieve[1]` I'll access 16 bytes and therefore memory which I does not own, don't I? – Matthias Jun 29 '16 at 15:33
  • Oh, I see since we only wanted 56 bits we'll never access `sieve[1]`... – Matthias Jun 29 '16 at 15:41
  • @Matthias If you wanted to be more exact you could use `ceil(ceil((double)upperLimit / sizeof(prime_t)) / 8) * sizeof(prime_t)`, but a few extra bytes won't hurt. – dbush Jun 29 '16 at 15:44
  • @Matthias I can't be bothered to decipher this, but there is zero relation between the number of lines of source code and memory used. There is no performance what-so-ever gained from squeezing in all these operations on a single line of code. I'd recommend studying assembler and how the compiler actually translates C code into machine code. Once you got this program up and running, you could also post it at http://codereview.stackexchange.com/ for advise how to rewrite it properly. – Lundin Jun 30 '16 at 06:15
  • @Lundin I was never talking about the line of code and that it couldn't get split up. I thought you were talking about the *and*'s and bit shifts. I'm well aware of that the lines of code written don't correlate with the memory use... My code is actually split up into `blockIndex` and `blockBit`, which I think greatly supports the reader. *codereview* is a great idea though :) – Matthias Jun 30 '16 at 12:13
  • @Matthias What he means is don't try to be too clever with your code. By breaking out the increment/decrement into a separate line from the bit checking, it makes your code easier to read and maintain. It might also be a good idea to create macros for getting/setting a bit to make it more readable. – dbush Jun 30 '16 at 12:15
  • @dbush You are right. I actually had to do so otherwise this would be a case of [unsequenced modification](http://stackoverflow.com/questions/3852768/unsequenced-value-computations-a-k-a-sequence-points). – Matthias Jun 30 '16 at 12:17
2

You allocate X bytes with calloc, dividing the total by sizeof(prime_t), yet act as if you have room for X prime_t elements later on.

Edit: Or actually even, you are allocating an array of 1 element with size X.

If you want to do it the way you are using it now, you should do:

calloc(X, sizeof(prime_t)) instead.

Edit: The major other issue in your code is that you are using byte-level indexing instead of bit-level.

Note that there are sizeof(prime_t) * 8 bits in a prime_t, so in every byte you set exactly 1 bit, true. You divide by sizeof(prime_t) instead of (sizeof(prime_t) * 8) when indexing.

marcolz
  • 2,880
  • 2
  • 23
  • 28
  • Because I only need one bit per prime not a whole `prime_t`. – Matthias Jun 29 '16 at 13:52
  • You current implementation uses 1 *byte*, not one bit. – marcolz Jun 29 '16 at 13:53
  • Where does it use one *byte*? I set single bits with `sieve[number/sizeof(prime_t)] |= (((prime_t) 1) << (number % sizeof(prime_t)));` – Matthias Jun 29 '16 at 13:54
  • 2
    `sizeof` returns bytes, not bits – Karsten Koop Jun 29 '16 at 13:55
  • 1
    Ok, so in every byte you set exactly 1 bit, true. You divide by `sizeof(prime_t)` instead of `(sizeof(prime_t) * 8)` when indexing. – marcolz Jun 29 '16 at 13:55
  • See [this](http://stackoverflow.com/questions/2525310/how-to-define-and-work-with-an-array-of-bits-in-c) to verify that my implementation is a valid bit array. – Matthias Jun 29 '16 at 13:57
  • @marcolz Oh I see, thank you really much! I would be happy if you'd add this as an answer. – Matthias Jun 29 '16 at 14:00
  • @Matthias I've just verified that the bug is that you don't allocate enough memory. Overall, the algorithm is needlessly obscure, with those real ugly while loops. I suspect you have lost track of what your own program does somewhere among those messy loops. Try to split the silly one-liners into several lines. – Lundin Jun 29 '16 at 14:07