1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:

int nKiebartLUT[8]               = { 1, 7, 11, 13, 17, 19, 23, 29 };
unsigned char bKiebartInvLUT[15] = { 1, 0, 0, 2, 0, 4, 8, 0, 16, 32, 0, 64, 0, 0, 128 };

unsigned char *eratosthenes (int nSize, int *pnFlagArraySize)
{
    if (nSize < 2)
    {
        return NULL;
    }
    else
    {
        int nFlagArraySize = (nSize + 29) / 30;     // Feldgroesse geeignet aufrunden
                                                    // Variable ?Odd repraesentiert die ungerade natuerliche Zahl (2*?Odd+1):
        int nOdd = 15 * nFlagArraySize;             // n = obere Grenze (ist sicher >nSize)
        int pOdd = 3;                               // Start p = 7
        int qOdd = 24;                              // Start q = p^2 = 49
        unsigned char *bSieve = (unsigned char *)malloc(nFlagArraySize);
        memset(bSieve, 0xFF, nFlagArraySize);       // Grundinitialisierung: Alle weder durch 2,3,5 teilbaren Zahlen sind zunaechst unter Primzahlverdacht
        bSieve[0] &= ~1;                            // 1 ist keine Primzahl
        while (qOdd < nOdd)
        {
            if ((bSieve[pOdd / 15] & bKiebartInvLUT[pOdd % 15]) > 0)
            {
                int p = 2 * pOdd + 1;               // naechste Primzahl im Sieb
                int rOdd = qOdd;                    // r = p^2
                while (rOdd < nOdd)                 // solange r < n)
                {
                    unsigned char bMask = bKiebartInvLUT[rOdd % 15];
                    if (bMask != 0)
                    {
                        // hierher zum Streichen wird nur verzweigt,
                        // wenn r weder durch 2,3,5 teilbar ist
                        bSieve[rOdd / 15] &= ~bMask;
                    }
                    rOdd += p;                      // r = r + 2p
                }
            }
            pOdd++;                                 // p = p+2
            qOdd += 4 * pOdd;                       // q = q + 4p , entspricht q = p^2
        }
        *pnFlagArraySize = nFlagArraySize;
        return bSieve;
    }
}

void printSieve (unsigned char *bSieve, int nFlagArraySize)
{
    printf("2,3,5");
    for (int nIndex = 0; nIndex < nFlagArraySize; nIndex++)
    {
        unsigned char bMask = bSieve[nIndex];
        for (int nBitPos = 0; nBitPos < 8; nBitPos++)
        {
            if (bMask & 1)
            {
                printf(",%d", 30 * nIndex + nKiebartLUT[nBitPos]);
            }
            bMask >>= 1;
        }
    }
    printf("\n");
}