2011-03-13 10 views
2

에 에라토스테네스의 체를 최적화 포팅 내가 여기 파이썬에서 (빠른 타오르는) primesieve 사용 전에 어떤 시간 : 지금파이썬에서 C++

def primes2(n): 
    """ Input n>=6, Returns a list of primes, 2 <= p < n """ 
    n, correction = n-n%6+6, 2-(n%6>1) 
    sieve = [True] * (n/3) 
    for i in xrange(1,int(n**0.5)/3+1): 
     if sieve[i]: 
     k=3*i+1|1 
     sieve[  k*k/3  ::2*k] = [False] * ((n/6-k*k/6-1)/k+1) 
     sieve[k*(k-2*(i&1)+4)/3::2*k] = [False] * ((n/6-k*(k-2*(i&1)+4)/6-1)/k+1) 
    return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]] 

: Fastest way to list all primes below N

가 정확하려면,이 구현을 나는 자동으로 2, 3 등의 배수를 건너 뛰는 최적화에 대한 아이디어를 약간 이해할 수 있지만 C++에이 알고리즘을 이식 할 때는 (나는 Python과 C/C에 대한 합리적인/나쁜 이해를 잘 이해하고있다. ,하지만 락앤롤에 충분히 좋다).

내가 현재 자신을 압연 한 것은이 (isqrt()은 단순한 정수 제곱근 함수)입니다 :

template <class T> 
void primesbelow(T N, std::vector<T> &primes) { 
    T sievemax = (N-3 + (1-(N % 2)))/2; 
    T i; 
    T sievemaxroot = isqrt(sievemax) + 1; 

    boost::dynamic_bitset<> sieve(sievemax); 
    sieve.set(); 

    primes.push_back(2); 

    for (i = 0; i <= sievemaxroot; i++) { 
     if (sieve[i]) { 
      primes.push_back(2*i+3); 
      for (T j = 3*i+3; j <= sievemax; j += 2*i+3) sieve[j] = 0; // filter multiples 
     } 
    } 

    for (; i <= sievemax; i++) { 
     if (sieve[i]) primes.push_back(2*i+3); 
    } 
} 

이 구현은 괜찮은 자동으로 2 배수를 건너 뛰고,하지만 난 포트 파이썬 구현을 할 수있는 경우 나는 그것이 훨씬 빠를 수 있다고 생각한다 (50 % -30 % 정도).

우분투 10.10은 1230ms 인 Q6600에 N=100000000, g++ -O3로, 현재 실행 시간을 (이 질문은 성공적으로 응답 될 것입니다 희망에서) 결과를 비교합니다.

위의 파이썬 구현이 무엇인지 이해하거나 나를 위해 이식하는 것이 도움이 될 것입니다.

편집

내가 어려운 무엇을 발견에 대한 몇 가지 추가 정보.

수정 변수와 일반적으로 어떻게 사용되는지 기술에 문제가 있습니다. 서로 다른 Eratosthenes 최적화를 설명하는 사이트에 대한 링크 ("2, 3 및 5의 배수를 건너 뛰고 1000 개의 C 파일로 슬램을 얻는 간단한 사이트는 제외)은 최고 일 것입니다.

나는 100 % 직접 및 글자 그대로의 포트에 문제가 없을 것이라고 생각하지만, 결국 이것은 전혀 쓸모가 없을 것입니다.

편집

원래 NumPy와 버전의 코드보고 후, 실제로 매우 쉽게 구현할 수 및 일부하지 이해하기 너무 열심히 생각으로 있습니다. 이것은 내가 생각해 낸 C++ 버전입니다. 나는 2 백만 라인의 코드가 아닌 매우 효율적인 primesieve를 필요로 할 경우에 대비하여 독자들을 돕기 위해 정식 버전으로 게시하고 있습니다. 이 primesieve는 위와 같은 시스템에서 약 415ms 동안 100000000 미만의 모든 소수를 처리합니다. 3 배의 속도 향상이 예상됩니다.

#include <vector> 
#include <boost/dynamic_bitset.hpp> 

// http://vault.embedded.com/98/9802fe2.htm - integer square root 
unsigned short isqrt(unsigned long a) { 
    unsigned long rem = 0; 
    unsigned long root = 0; 

    for (short i = 0; i < 16; i++) { 
     root <<= 1; 
     rem = ((rem << 2) + (a >> 30)); 
     a <<= 2; 
     root++; 

     if (root <= rem) { 
      rem -= root; 
      root++; 
     } else root--; 

    } 

    return static_cast<unsigned short> (root >> 1); 
} 

// https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188 
// https://stackoverflow.com/questions/5293238/porting-optimized-sieve-of-eratosthenes-from-python-to-c/5293492 
template <class T> 
void primesbelow(T N, std::vector<T> &primes) { 
    T i, j, k, l, sievemax, sievemaxroot; 

    sievemax = N/3; 
    if ((N % 6) == 2) sievemax++; 

    sievemaxroot = isqrt(N)/3; 

    boost::dynamic_bitset<> sieve(sievemax); 
    sieve.set(); 

    primes.push_back(2); 
    primes.push_back(3); 

    for (i = 1; i <= sievemaxroot; i++) { 
     if (sieve[i]) { 
      k = (3*i + 1) | 1; 
      l = (4*k-2*k*(i&1))/3; 

      for (j = k*k/3; j < sievemax; j += 2*k) { 
       sieve[j] = 0; 
       sieve[j+l] = 0; 
      } 

      primes.push_back(k); 
     } 
    } 

    for (i = sievemaxroot + 1; i < sievemax; i++) { 
     if (sieve[i]) primes.push_back((3*i+1)|1); 
    } 
} 
+0

파이썬 코드는 단일 요소를 배열로 푸시하지 않습니다. 요소 블록으로 초기화 한 다음 제자리에서 수정합니다. 그것은 C++에서 훨씬 더 빠를 것입니다. 여러분이 가진 것처럼'dynamic_bitset'을 시작하고, 그것을 채우고, 하나의 루프에서 체를 작성하십시오. 아마도 Python 버전의 기술을 복사하지 않고도 훨씬 빠르게 버전을 실행할 수 있습니다. –

+0

파이썬 구현의 어느 부분을 이해하는데 문제가 있습니까? 이 정보를 질문에 추가하십시오. 또한, 사람들은 당신이 말했듯이, 당신의 이해와 많은 도움이되지 않을 것이기 때문에 정말로 당신을위한 항구로 대답해서는 안됩니다. –

+0

@Jeremiah Willcock : 정확하지 않습니다. sieving하는 동안 그것은 내부에서 수정되지만, 결국 배열 (list)'xrange (1, n)의 i에 대해 return [2,3] + [3 * i + 1 | 1]/3-correction)을 사용하면된다. 그리고 그것은 제가하는 일이기도합니다. 첫 번째 루프는 체재 중이며 이미 알고있는 소수를 푸시하고 두 번째 루프는 sievelimit의 루트 이후에 오는 소수를 푸시하는 것입니다. – orlp

답변

3

나는 최대한 많이 설명하려고 노력할 것이다. sieve 배열에는 비정상적인 색인 체계가 있습니다. 1 또는 5 mod 6과 일치하는 각 숫자에 대한 비트를 저장합니다. 따라서 2*k 숫자에 저장되고 k*6 + 5 숫자는 2*k + 1 위치에 저장됩니다. 3*i+1|1 동작은 그 반대이다 : 이는 형태 2*n의 수를 취하여 6*n + 1로 변환하고 2*n + 1 소요합니다 (+1|1 것은 05-1하고 3 변환) 6*n + 5로 변환한다.주 루프는 5 (i이 1 일 때)으로 시작하는 해당 속성의 모든 숫자를 통해 k을 반복합니다. iksieve에 해당하는 색인입니다. sieve으로 첫 번째 슬라이스를 업데이트하면 체의 모든 비트가 k*k/3 + 2*m*k (m 자연수) 형식의 인덱스로 지워집니다. 해당 인덱스의 해당 숫자는 k^2에서 시작하여 각 단계에서 6*k까지 증가합니다. 두 번째 슬라이스 업데이트는 k*(k-2*(i&1)+4)/3 (k1 mod 6k * (k+2)과 일치 함)의 번호 k * (k+4)에서 시작하고 각 단계에서 마찬가지로 숫자를 6*k 증가시킵니다.

여기에 대한 설명에서 또 다른 시도이다 : candidates 적어도 5이며, 하나 1 또는 5 모드 6 합동 모든 숫자의 집합이 될 수 있습니다. 해당 집합에서 두 요소를 곱하면 집합에 다른 요소가 생깁니다. 숫자 candidatesk에 대한 succ(k)k보다 큰 candidates의 다음 요소 (숫자 순)로 설정하십시오. 이 경우, 자체의 내부 루프 (sieve 정상 색인을 사용하여) 기본적 :

for k in candidates: 
    for l in candidates where l >= k: 
    sieve[k * l] = False 
: 때문에 소자 sieve에 저장되어있는 한계
for k in candidates: 
    for (l = k; ; l += 6) sieve[k * l] = False 
    for (l = succ(k); ; l += 6) sieve[k * l] = False 

, 즉과 동일하다 어떤 점에서 체에서 ( k 자체보다 다른) candidatesk의 배수를 모두 제거합니다

은 (중 현재 k는 이전 l로 사용하는 경우 또는 이제 k로 사용하는 경우).

+0

좋아요, 그건 기본적으로 2를 배제하기 위해 제가 체로 한 것입니다,하지만 조금 더 복잡합니다. 몇 번 읽었을 때 이것을 C++ 구현으로 바꿀 수 있는지 살펴 보겠습니다. 모든 경우 +1. – orlp

+0

@nightcracker : 이제 인덱스 -> 숫자 계산이 완료되었습니다. 코드가 까다 롭지 만, 그 이유는'sieve '에 대한 인덱싱 때문입니다. –

+0

감사합니다. – orlp

0

제쳐두고 소수를 "근사 할"수 있습니다. 대략 프라임 P. 전화 여기 몇 화학식이다 :

P = 2 * K + 1 // 나눌 수없는

P = 6 (2)에 의한 기록 * K + {1,5} // 2 나눌하지 3

P = 30 * K + 2, 3, 5

숫자의 세트의 특성을 발견하여 {1, 7, 11, 13, 17, 19, 23, 29} // divisble하지 이 공식에 따르면 P는 소수 일 수는 없지만 모든 소수는 집합 P에 있습니다. 집합 P에 속한 소수만을 테스트하면 아무 것도 놓치지 않을 것입니다.

당신은 이러한 공식을 재구성 할 수 있습니다

P = X * K + {-i, -j, -k, K, J, I}

이 당신을 위해 더 편리합니다.

Here

이 기술은 실질적으로 활용 될 수있는 정도를 나타낼 수있는, 2, 3, 5, 7

이 링크로 나눌하지 P하는 식으로이 기술을 사용하는 코드이다.

1

Howard Higgle의 반응, Howard에 피기 백킹으로 소수성을 위해 2, 3 또는 5로 나눌 수없는 모든 자연수 집합에서 숫자를 테스트 할 필요가 없습니다.배열 자체의 모든 숫자와 그 이후의 모든 숫자를 배열에 곱하면됩니다 (1을 제외하고 자기 제거). 이러한 겹치는 제품은 결정 론적 곱셈 프로세스를 확장하는 모든 점까지 배열의 모든 비원 수를 제공합니다. 따라서 배열의 첫 번째 비 프라임은 7 제곱 또는 49입니다. 두 번째, 7 번 11 또는 77 등. 전체 설명은 여기를 참조하십시오. http://www.primesdemystified.com

관련 문제