Сегментированное сито эратосфена?
Достаточно просто сделать простое сито:
for (int i=2; i<=N; i++){
if (sieve[i]==0){
cout << i << " is prime" << endl;
for (int j = i; j<=N; j+=i){
sieve[j]=1;
}
}
cout << i << " has " << sieve[i] << " distinct prime factors\n";
}
Но что делать, когда N очень велико, и я не могу держать этот массив в памяти? Я искал подходы с сегментированными ситами, и они, кажется, связаны с поиском простых чисел вплоть до sqrt(N), но я не понимаю, как это работает. Что если N очень большое (скажем, 10^18)?
6 ответов
Основная идея сегментированного сита состоит в том, чтобы выбрать простые числа просеивания, меньшие, чем квадратный корень из n, выбрать достаточно большой размер сегмента, который, тем не менее, помещается в память, а затем просеять каждый из сегментов по очереди, начиная с наименьшего. В первом сегменте вычисляется наименьшее кратное каждого простого сита, находящегося в сегменте, затем кратные значения простого сита отмечаются как составные обычным образом; когда все простые числа просеивания были использованы, оставшиеся немаркированные числа в сегменте являются простыми. Затем, для следующего сегмента, для каждого простого сита вы уже знаете первое кратное в текущем сегменте (это было кратное число, которое закончило просеивание для этого простого в предыдущем сегменте), так что вы просеиваете каждое простое сито, и так далее пока ты не закончил.
Размер n не имеет значения, за исключением того, что большее n займет больше времени для просеивания, чем меньшее n; Размер, который имеет значение, - это размер сегмента, который должен быть настолько большим, насколько это удобно (скажем, размер основного кэша памяти на компьютере).
Вы можете увидеть простую реализацию сегментированного сита здесь. Обратите внимание, что сегментированное сито будет намного быстрее, чем сито с приоритетной очередью О'Нила, упомянутое в другом ответе; если вам интересно, здесь есть реализация.
РЕДАКТИРОВАТЬ: Я написал это для другой цели, но я покажу это здесь, потому что это может быть полезно:
Хотя сито Эратосфена очень быстрое, оно требует O(n) пространства. Это может быть уменьшено до O(sqrt(n)) для простых чисел просеивания плюс O(1) для растровой матрицы путем выполнения просеивания в последовательных сегментах. В первом сегменте вычисляется наименьшее кратное каждого простого сита, которое находится внутри сегмента, затем кратные значения простого сита просчитываются составным образом обычным образом; когда все простые числа просеивания были использованы, оставшиеся немаркированные числа в сегменте являются простыми. Затем для следующего сегмента наименьшее кратное каждого простого сита является кратным, которое закончило просеивание в предыдущем сегменте, и поэтому просеивание продолжается до завершения.
Рассмотрим пример сита от 100 до 200 в сегментах по 20. Пять простых сит: 3, 5, 7, 11 и 13. В первом сегменте с 100 по 120 битрейн имеет десять слотов, причем слот 0 соответствует 101 слот k соответствует 100+2k+1, а слот 9 соответствует 119. Наименьшее кратное 3 в сегменте равно 105, что соответствует слоту 2; слоты 2+3=5 и 5+3=8 также кратны 3. Наименьшее кратное 5 равно 105 в слоте 2, а слот 2+5=7 также кратно 5. Наименьшее кратное 7 равно 105 в слоте 2 и слоте 2+7=9 также кратно 7. И так далее.
Функция primesRange принимает аргументы lo, hi и delta; lo и hi должны быть четными, с lo При вызове primesRange(100, 200, 10) простые числа просеивания ps равны [3, 5, 7, 11, 13]; qs изначально [2, 2, 2, 10, 8] соответствует наименьшим коэффициентам 105, 105, 105, 121 и 117 и сбрасывается для второго сегмента на [1, 2, 6, 0, 11], соответствующий наименьшему кратны 123, 125, 133, 121 и 143. Вы можете увидеть эту программу в действии на http://ideone.com/iHYr1f. И в дополнение к ссылкам, показанным выше, если вы заинтересованы в программировании с простыми числами, я скромно рекомендую это эссе в моем блоге.function primesRange(lo, hi, delta)
function qInit(p)
return (-1/2 * (lo + p + 1)) % p
function qReset(p, q)
return (q - delta) % p
sieve := makeArray(0..delta-1)
ps := tail(primes(sqrt(hi)))
qs := map(qInit, ps)
while lo < hi
for i from 0 to delta-1
sieve[i] := True
for p,q in ps,qs
for i from q to delta step p
sieve[i] := False
qs := map(qReset, ps, qs)
for i,t from 0,lo+1 to delta-1,hi step 1,2
if sieve[i]
output t
lo := lo + 2 * delta
Просто мы делаем сегментацию с помощью сита, которое у нас есть. Основная идея заключается в том, что нам нужно найти простые числа от 85 до 100. Мы должны применять традиционное сито, но так, как описано ниже:
Итак, мы берем первое простое число 2, делим начальное число на 2(85/2) и, округляя до меньшего числа, получаем p=42, теперь снова умножаем на 2 и получаем p=84, с этого момента начинаем добавлять 2 до последнего числа. Итак, мы сделали, что мы удалили все факторы 2(86,88,90,92,94,96,98,100) в диапазоне.
Мы берем следующее простое число 3, делим начальное число на 3(85/3) и, округляя до меньшего числа, получаем p=28, теперь снова умножаем на 3 и получаем p=84, с этого момента начинаем прибавлять 3 до последнее число. Итак, мы сделали, что мы удалили все факторы 3(87,90,93,96,99) в диапазоне.
Возьмите следующее простое число =5 и т. Д................... Продолжайте выполнять вышеуказанные шаги. Вы можете получить простые числа (2,3,5,7,...) с помощью традиционного сита до sqrt(n). И затем используйте его для сегментированного сита.
Существует версия Sieve, основанная на приоритетных очередях, которая дает столько простых чисел, сколько вы запрашиваете, а не все из них вплоть до верхней границы. Это обсуждается в классической статье "Подлинное сито эратосфена", и поиск "очереди с приоритетом ситарафосфена" приводит к появлению множества реализаций на разных языках программирования.
Если кто-то хотел бы увидеть реализацию C++, вот мое:
void sito_delta( int delta, std::vector<int> &res)
{
std::unique_ptr<int[]> results(new int[delta+1]);
for(int i = 0; i <= delta; ++i)
results[i] = 1;
int pierw = sqrt(delta);
for (int j = 2; j <= pierw; ++j)
{
if(results[j])
{
for (int k = 2*j; k <= delta; k+=j)
{
results[k]=0;
}
}
}
for (int m = 2; m <= delta; ++m)
if (results[m])
{
res.push_back(m);
std::cout<<","<<m;
}
};
void sito_segment(int n,std::vector<int> &fiPri)
{
int delta = sqrt(n);
if (delta>10)
{
sito_segment(delta,fiPri);
// COmpute using fiPri as primes
// n=n,prime = fiPri;
std::vector<int> prime=fiPri;
int offset = delta;
int low = offset;
int high = offset * 2;
while (low < n)
{
if (high >=n ) high = n;
int mark[offset+1];
for (int s=0;s<=offset;++s)
mark[s]=1;
for(int j = 0; j< prime.size(); ++j)
{
int lowMinimum = (low/prime[j]) * prime[j];
if(lowMinimum < low)
lowMinimum += prime[j];
for(int k = lowMinimum; k<=high;k+=prime[j])
mark[k-low]=0;
}
for(int i = low; i <= high; i++)
if(mark[i-low])
{
fiPri.push_back(i);
std::cout<<","<<i;
}
low=low+offset;
high=high+offset;
}
}
else
{
std::vector<int> prime;
sito_delta(delta, prime);
//
fiPri = prime;
//
int offset = delta;
int low = offset;
int high = offset * 2;
// Process segments one by one
while (low < n)
{
if (high >= n) high = n;
int mark[offset+1];
for (int s = 0; s <= offset; ++s)
mark[s] = 1;
for (int j = 0; j < prime.size(); ++j)
{
// find the minimum number in [low..high] that is
// multiple of prime[i] (divisible by prime[j])
int lowMinimum = (low/prime[j]) * prime[j];
if(lowMinimum < low)
lowMinimum += prime[j];
//Mark multiples of prime[i] in [low..high]
for (int k = lowMinimum; k <= high; k+=prime[j])
mark[k-low] = 0;
}
for (int i = low; i <= high; i++)
if(mark[i-low])
{
fiPri.push_back(i);
std::cout<<","<<i;
}
low = low + offset;
high = high + offset;
}
}
};
int main()
{
std::vector<int> fiPri;
sito_segment(1013,fiPri);
}
Основываясь на ответе Swapnil Kumar, я сделал следующий алгоритм на C. Он был построен с помощью mingw32-make.exe.
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
int main()
{
const int MAX_PRIME_NUMBERS = 5000000;//The number of prime numbers we are looking for
long long *prime_numbers = malloc(sizeof(long long) * MAX_PRIME_NUMBERS);
prime_numbers[0] = 2;
prime_numbers[1] = 3;
prime_numbers[2] = 5;
prime_numbers[3] = 7;
prime_numbers[4] = 11;
prime_numbers[5] = 13;
prime_numbers[6] = 17;
prime_numbers[7] = 19;
prime_numbers[8] = 23;
prime_numbers[9] = 29;
const int BUFFER_POSSIBLE_PRIMES = 29 * 29;//Because the greatest prime number we have is 29 in the 10th position so I started with a block of 841 numbers
int qt_calculated_primes = 10;//10 because we initialized the array with the ten first primes
int possible_primes[BUFFER_POSSIBLE_PRIMES];//Will store the booleans to check valid primes
long long iteration = 0;//Used as multiplier to the range of the buffer possible_primes
int i;//Simple counter for loops
while(qt_calculated_primes < MAX_PRIME_NUMBERS)
{
for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
possible_primes[i] = 1;//set the number as prime
int biggest_possible_prime = sqrt((iteration + 1) * BUFFER_POSSIBLE_PRIMES);
int k = 0;
long long prime = prime_numbers[k];//First prime to be used in the check
while (prime <= biggest_possible_prime)//We don't need to check primes bigger than the square root
{
for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
if ((iteration * BUFFER_POSSIBLE_PRIMES + i) % prime == 0)
possible_primes[i] = 0;
if (++k == qt_calculated_primes)
break;
prime = prime_numbers[k];
}
for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
if (possible_primes[i])
{
if ((qt_calculated_primes < MAX_PRIME_NUMBERS) && ((iteration * BUFFER_POSSIBLE_PRIMES + i) != 1))
{
prime_numbers[qt_calculated_primes] = iteration * BUFFER_POSSIBLE_PRIMES + i;
printf("%d\n", prime_numbers[qt_calculated_primes]);
qt_calculated_primes++;
} else if (!(qt_calculated_primes < MAX_PRIME_NUMBERS))
break;
}
iteration++;
}
return 0;
}
Он устанавливает максимальное количество простых чисел, которые будут найдены, затем массив инициализируется известными простыми числами, такими как 2, 3, 5...29. Таким образом, мы создаем буфер, который будет хранить сегменты возможных простых чисел, этот буфер не может быть больше, чем мощность наибольшего начального простого числа, которое в данном случае равно 29.
Я уверен, что для повышения производительности можно выполнить множество оптимизаций, например распараллелить процесс анализа сегментов и пропустить числа, кратные 2, 3 и 5, но это служит примером низкого потребления памяти.
Число является простым, если ни одно из меньших простых чисел не делится на него. Поскольку мы перебираем простые числа по порядку, мы уже пометили все числа, которые делятся хотя бы на одно из простых чисел, делящимися. Следовательно, если мы достигнем ячейки, и она не отмечена, то она не делится ни на какое меньшее простое число и, следовательно, должна быть простой.
Запомните эти пункты:
// Generating all prime number up to R
// creating an array of size (R-L-1) set all elements to be true: prime && false: composite
#include<bits/stdc++.h>
using namespace std;
#define MAX 100001
vector<int>* sieve(){
bool isPrime[MAX];
for(int i=0;i<MAX;i++){
isPrime[i]=true;
}
for(int i=2;i*i<MAX;i++){
if(isPrime[i]){
for(int j=i*i;j<MAX;j+=i){
isPrime[j]=false;
}
}
}
vector<int>* primes = new vector<int>();
primes->push_back(2);
for(int i=3;i<MAX;i+=2){
if(isPrime[i]){
primes->push_back(i);
}
}
return primes;
}
void printPrimes(long long l, long long r, vector<int>*&primes){
bool isprimes[r-l+1];
for(int i=0;i<=r-l;i++){
isprimes[i]=true;
}
for(int i=0;primes->at(i)*(long long)primes->at(i)<=r;i++){
int currPrimes=primes->at(i);
//just smaller or equal value to l
long long base =(l/(currPrimes))*(currPrimes);
if(base<l){
base=base+currPrimes;
}
//mark all multiplies within L to R as false
for(long long j=base;j<=r;j+=currPrimes){
isprimes[j-l]=false;
}
//there may be a case where base is itself a prime number
if(base==currPrimes){
isprimes[base-l]= true;
}
}
for(int i=0;i<=r-l;i++){
if(isprimes[i]==true){
cout<<i+l<<endl;
}
}
}
int main(){
vector<int>* primes=sieve();
int t;
cin>>t;
while(t--){
long long l,r;
cin>>l>>r;
printPrimes(l,r,primes);
}
return 0;
}