сегментированные премьер сито

Question

Наткнулся этого эффективного сегментированного простого решета в Интернете, пожалуйста, помогите мне понять, работая особенно использование следующего вектора

Как конкретный выбор размера сегмента влияет на производительность?

const int L1D_CACHE_SIZE = 32768; 
void segmented_sieve(int64_t limit, int segment_size = L1D_CACHE_SIZE) 
{ 
    int sqrt = (int) std::sqrt((double) limit); 
    int64_t count = (limit < 2) ? 0 : 1; 
    int64_t s = 2; 
    int64_t n = 3; 

    // vector used for sieving 
    std::vector<char> sieve(segment_size); 

    // generate small primes <= sqrt 
    std::vector<char> is_prime(sqrt + 1, 1); 
    for (int i = 2; i * i <= sqrt; i++) 
     if (is_prime[i]) 
      for (int j = i * i; j <= sqrt; j += i) 
       is_prime[j] = 0; 

    std::vector<int> primes; 
    std::vector<int> next; 

    for (int64_t low = 0; low <= limit; low += segment_size) 
    { 
     std::fill(sieve.begin(), sieve.end(), 1); 

     // current segment = interval [low, high] 
     int64_t high = std::min(low + segment_size - 1, limit); 

     // store small primes needed to cross off multiples 
     for (; s * s <= high; s++) 
     { 
      if (is_prime[s]) 
      { 
       primes.push_back((int) s); 
       next.push_back((int)(s * s - low)); 
      } 
     } 
     // sieve the current segment 
     for (std::size_t i = 1; i < primes.size(); i++) 
     { 
      int j = next[i]; 
      for (int k = primes[i] * 2; j < segment_size; j += k) 
       sieve[j] = 0; 
      next[i] = j - segment_size; 
     } 

     for (; n <= high; n += 2) 
      if (sieve[n - low]) // n is a prime 
       count++; 
    } 

    std::cout << count << " primes found." << std::endl; 
} 

Answer 1

Вот более краткая формулировка того же алгоритма, которая должна сделать принцип более понятным (часть full, runnable .cpp for segment size timings @ pastebin). Это инициализирует упакованное (только для шансов) сито вместо подсчета простых чисел, но задействованные принципы одинаковы. Загрузите и запустите .cpp, чтобы увидеть влияние размеров сегмента. В принципе, оптимальный размер должен быть вокруг размера кеша L1 для вашего процессора. Слишком мало, и накладные расходы из-за увеличения количества раундов начинают доминировать; слишком большой, и вы будете наказаны медленными таймингами кэшей L2 и L3. См. Также How does segmentation improve the running time of Sieve of Eratosthenes?.

void initialise_packed_sieve_4G (void *data, unsigned segment_bytes = 1 << 15, unsigned end_bit = 1u << 31) 
{ 
    typedef std::vector<prime_and_offset_t>::iterator prime_iter_t; 
    std::vector<prime_and_offset_t> small_factors; 

    initialise_odd_primes_and_offsets_64K(small_factors); 

    unsigned segment_bits = segment_bytes * CHAR_BIT; 
    unsigned partial_bits = end_bit % segment_bits; 
    unsigned segments  = end_bit/segment_bits + (partial_bits != 0); 

    unsigned char *segment = static_cast<unsigned char *>(data); 
    unsigned bytes = segment_bytes; 

    for (; segments--; segment += segment_bytes) 
    { 
     if (segments == 0 && partial_bits) 
     { 
     segment_bits = partial_bits; 
     bytes = (partial_bits + CHAR_BIT - 1)/CHAR_BIT; 
     } 

     std::memset(segment, 0, bytes); 

     for (prime_iter_t p = small_factors.begin(); p != small_factors.end(); ++p) 
     { 
     unsigned n = p->prime; 
     unsigned i = p->next_offset; 

     for (; i < segment_bits; i += n) 
     { 
      set_bit(segment, i); 
     } 

      p->next_offset = i - segment_bits; 
     } 
    } 
} 

Если смещения не запомнились от сегмента к сегменту, то они должны быть пересчитаны каждый раз, используя, по меньшей мере, одно деление и умножение один за пересчитывается индекса, а также условными или серьезный битную фокусы. При просеивании полного диапазона номеров 2^32 (8192 сегмента по 32 Кбайт каждый), это не менее 53 583 872 медленных делений и столько же несколько более быстрых умножений; это примерно одна секунда, добавленная к времени, необходимому для инициализации полного сита 2^32 (2^31 бит для только для коэффициентов Eratosthenes).

Вот некоторый реальный код от одного из моих старых сит, который использует эту «преобразовывающую» математику:

for (index_t k = 1; k <= max_factor_bit; ++k) 
{ 
    if (bitmap_t::traits::bt(bm.bm, k)) continue; 

    index_t n = (k << 1) + 1;  // == index_for_value(value_for_index(k) * 2) == n 
    index_t i = square(n) >> 1; // == index_for_value(square(n)) 

    if (i < offset) 
    { 
     i += ((offset - i)/n) * n; 
    } 

    for (; i <= new_max_bit; i += n) 
    { 
     bitmap_t::traits::bts(bm.bm, i); 
    } 
} 

Это занимает около 5,5 секунд для полного ситя (VC++); сначала показан код, который занимает всего 4,5 секунды с тем же компилятором или 3,5 секунды с использованием gcc 4.8.1 (MinGW64).

Вот GCC тайминги:

sieve bits = 2147483648 (equiv. number = 4294967295) 

segment size 4096 (2^12) bytes ... 4.091 s 1001.2 M/s 
segment size 8192 (2^13) bytes ... 3.723 s 1100.2 M/s 
segment size 16384 (2^14) bytes ... 3.534 s 1159.0 M/s 
segment size 32768 (2^15) bytes ... 3.418 s 1198.4 M/s 
segment size 65536 (2^16) bytes ... 3.894 s 1051.9 M/s 
segment size 131072 (2^17) bytes ... 4.265 s 960.4 M/s 
segment size 262144 (2^18) bytes ... 4.453 s 919.8 M/s 
segment size 524288 (2^19) bytes ... 5.002 s 818.9 M/s 
segment size 1048576 (2^20) bytes ... 5.176 s 791.3 M/s 
segment size 2097152 (2^21) bytes ... 5.135 s 797.7 M/s 
segment size 4194304 (2^22) bytes ... 5.251 s 780.0 M/s 
segment size 8388608 (2^23) bytes ... 7.412 s 552.6 M/s 

digest { 203280221, 0C903F86, 5B253F12, 774A3204 } 

Примечания: дополнительная второй может быть выбрито с этого времени с помощью трюка под названием «presieving», т.е. подрывного предварительно вычисленный шаблон в битовые карты вместо обнуления его в начале. Это приводит к отключению gcc до 2.1 с для полного сита с this hacked copy of the earlier .cpp. Этот трюк работает очень хорошо вместе с сегментированным просеиванием в кусках размером с кеш.

Answer 2

Я не являюсь специалистом в этом, но моя кишка говорит мне так:

1. граница поиска решето стол

   вписываться в L1 КЭШ процессора , чтобы получить все выгоды от повышение производительности современных аппаратных архитектур

2. next

   если вы хотите segmentate решета , то вы должны помнить последний индекс за каждый штрих уже просеивают, например:

   • просеивают штрихи: 2,3,5

   • размера сегмента: 8

     |0 1 2 3 4 5 6 7|0 1 2 3 4 5 6 7|0 1 2 3 4 5 6 7| // segments 
     ----------------------------------------------- 
     2|- x x x x x x x x x x x 
     3|-  x  x  x  x  x  x  x 
     5|-   x   x   x   x  
     ----------------------------------------------- 
     |    ^   ^   ^
               // next value offset for each prime 
     

   поэтому при заполнении следующего сегмента вы продолжаете плавно ...