Решето Эратосфена для малых лимитов


Я написала решето Эратосфена, я думаю, что это достаточно быстро для мелких лимитов. Любые предложения, как улучшить это? Особенно для небольших пределах, скажем пределов < 50,000,000. (Также закодированные параллельной сито (в C#), которая находит простые числа < 2^32 в 1,6 секунды, но это относительно медленный для мелких лимитов).

Как это работает

Давайте найдем простые <= 75. Странно композитов отмечены в беззнаковый инт массив "х", каждый бит представляет собой нечетное число. Первый (шестнадцатеричное) число в X:

х[0] = 0x9B4B3491

  • Двоичное: 1001 1011 0100 1011 0011 1001 0100 0001
  • Наименьший бит комплект: 1 Не премьер. Следующий бит не установлен: 3-простое, и т. д.

х[1] = 0xFFFFFFE5

  • Бинарные: 1111 1111 1111 1111 1111 1111 1110 0101
  • Простые числа: 67, 71, 73

Считать простые числа, 0-биты в Х, с __popcnt(unsigned int ...), (населения, число бит). __popcnt(x[0]) дает число бит, число нечетных композитов. __popcnt(~x[0]) дает число нулевых битов, то число нечетных простых чисел.

~х[0] = 0x64B4CB6E

  • Двоичная: 0110 0100 1011 0100 0110 1100 1011 1110

~х[1] = 0x0000001A

  • Бинарные: 0000 0000 0000 0000 0000 0000 0001 1010

~х[0] имеет 17 бит, ~х[1] и 3, 20 нечетных простых чисел. Размер массива Prime "п" до 21. Это 20 нечетных простых чисел и даже премьер-1 (2).

// Put the primes in "p", p[0] = 2
// Initialize: v = 1
// Count trailing zeros (tz), shift xi right by tz + 1:
//       xi = ~x[0] = ....101101101110  tz = 1, p[1] = v += tz << 1 =  3
//                                                 xi >>= tz + 1, v += 2
//                          1011011011  tz = 0, p[2] = v += tz << 1 =  5
//                                                 xi >>= tz + 1, v += 2
//                           101101101  tz = 0, p[3] = v += tz << 1 =  7
//                                                 xi >>= tz + 1, v += 2
//                            10110110  tz = 1, p[4] = v += tz << 1 = 11
//                                                 xi >>= tz + 1, v += 2
//                              101101  tz = 0, p[5] = v += tz << 1 = 13
//                                                 xi >>= tz + 1, v += 2
//                               10110  tz = 1, p[6] = v += tz << 1 = 17
//                                                 xi >>= tz + 1, v += 2
// 
// Results (Intel i7-4790 CPU @3.6 GHz from 2015):
// primes <= 100   109 ns
//        <= 200   187 ns
//        <= 400   312 ns
//        <= 800   578 ns
//        < 2^31   7.4 seconds
//        < 2^32  15.8 seconds

#include "stdafx.h"
#include <vector>
#include <iostream>
#include <ctime>
#include <Windows.h>
using namespace std;
typedef unsigned int u32;
vector<u32> x; vector<u32> p;

void buildX(u32 m)  // mark odd composites
{
    // make it safe:  1 << (int)(a >> 1)  =>    1 << (int)(a >> 1 & 31) ?
    //                1 << (int)d         =>    1 << (int)(d & 31)      ?
    //              ~0u << (int)m         =>  ~0u << (int)(m & 31)      ?

    m -= m / ~0u; m += m & 1; m >>= 1;
    x.resize((m >> 5) + 1); x[0] = 1;
    for (u32 a = 3, b = 4, c = 4, d; b < m; a += 2, b += c += 4)
        if ((x[a >> 6] & 1 << (int)(a >> 1)) == 0)
            for (d = b; d < m; d += a) x[d >> 5] |= 1 << (int)d;
    x[m >> 5] |= ~0u << (int)m;
}

void countPrimes()
{
    u32 c = 1; int i = x.size() - 1;
    while (i >= 0) c += __popcnt(~x[i--]);
    p.resize(c);
}

void primes(u32 m)  // primes <= m
{
    if (m > 1)
    {
        buildX(m); countPrimes(); p[0] = 2;
        u32 u = 1, v = 1, xi; DWORD tz;
        for (int i = 0, j = p.size(), n = 1;;)
        {
            xi = ~x[i++];
            while (xi)
            {
                _BitScanForward(&tz, xi); xi >>= tz; xi >>= 1;
                p[n++] = v += tz << 1; v += 2;
            }
            if (n >= j) break;
            v = u += 64;
        }
    }
}

int main()
{
    for (u32 m = 25; m <= 6400; m <<= 1)  // for (u32 m = 0; m < 9; m++)
    {
        primes(m);
        if (m > 1)
        {
            u32 maxP = p[p.size() - 1];
            cout << "largest prime <= " << m << " : " << maxP << "  ";
        }
        clock_t clock0 = clock();
        for (int i = 1000000; i; i--)
        {
            x.resize(0); p.resize(0); primes(m);
        }
        cout << clock() - clock0 << " ns: ";
        cout << p.size() << " primes" << endl << endl;
    }
    x.resize(0); p.resize(0); u32 m = ~0u;
    clock_t clock0 = clock(); primes(m);
    cout << (clock() - clock0) * 1e-3 << " s: ";      //  15.709 s
    cout << p.size() << " primes <= " << m << endl;

    x.resize(0); p.resize(0);
    clock0 = clock(); primes(m);
    cout << (clock() - clock0) * 1e-3 << " s: ";      //  15.522 s

    // 15.709 - 15.522 = 0.187 => memory allocation makes it ~ 1.2 % slower?

    cout << p.size() << " primes <= " << m << endl;
    for (int j = p.size(), i = j - 5; i < j;) cout << p[i++] << endl;
    getchar();
}


348
4
задан 7 марта 2018 в 08:03 Источник Поделиться
Комментарии
2 ответа

В целом код немного трудно читать с нескольких инструкций по строкам исходного кода, несколько переменных (различных типов) на одну строку исходного кода, а не описательные имена переменных.

Вы используете много жестко константы, при условии, что вы работаете с 32-битными целыми числами. Использование именованных констант для этих бы было понятнее, где некоторые из ваших номеров, а также позволяет для более легкого расширения до 64-разрядной компиляции.

Выражение 1 << (int)d используется в buildX результаты в неопределенное поведение, когда d больше, чем число битов, используемых в левый операнд (32 в данном случае). Множество аппаратных реализаций сдвиг влево уменьшить счетчик сдвига по модулю число битов в операнде перемещается, так это, скорее всего, делать то, что вы ожидаете, но это не гарантируется.

Поскольку векторы уже есть все памяти они должны выделяться, когда вы входите в свой 1 млн. итерации цикла сроки, время, предназначены исключительно для вашего алгоритма и не учитывают время, необходимое для выделения памяти как векторы расширения. Сроки для вашего максимального делу, так как он запускается только один раз, включает время, затраченное растет векторов. (Это просто 2 ассигнования, но разница в таймингах.)

3
ответ дан 8 марта 2018 в 01:03 Источник Поделиться

Быстрее "buildX(u32 м)". Нечетные, кратные 3 написано
три 32-битных слова на "х", начиная с координат x[1],
простое решение для первого слова: х[0] = 0x9b4b3491;
Аналогичным образом следующий премьер(ы) 5,7.. может быть отмечен,
но с возвратами deminishing.

/*                   old-ns-new               old-s-new
primes < 25 47 47 p < 2^31 7.4 7.1
< 50 62 62 < 2^32 15.8 15.2
< 100 109 94
< 200 187 156
< 400 312 250
< 800 578 436
< 1600 1092 843
< 3200 2152 1700
< 6400 4274 3385 */

void buildX(u32 m) // mark odd composites
{
m -= m / ~0u; m += m & 1; m >>= 1;
u32 a = 1, b = 2, c = 3, d = m >> 5;
x.resize(d + 1); x[0] = 0x9b4b3491;
for (; a <= d; a += 3) x[a] = 0x24924924;
for (; b <= d; b += 3) x[b] = 0x49249249;
for (; c <= d; c += 3) x[c] = 0x92492492;
for (a = 5, b = 12, c = 8; b < m; a += 2, b += c += 4)
if ((x[a >> 6] & 1 << (a >> 1)) == 0)
for (d = b; d < m; d += a) x[d >> 5] |= 1 << d;
x[m >> 5] |= ~0u << m;
}

Последняя версия: https://pastebin.com/JMdTxbeJ
Праймы < 2^32 за 4,2 секунды.

0
ответ дан 11 марта 2018 в 04:03 Источник Поделиться