ps:소수의_개수_구하기
소수의 개수 구하기
- x이하의 소수의 개수를 π(x)라고 표현하고, 이것을 Prime-counting function이라고 부른다
- π(x)이 x/ln(x)로 근사된다는 것은 잘 알려진 내용이다. 바로 그 유명하고도 중요한 소수 정리이다.
- 하지만 근사값이 아니라 특정 n에 대해서 정확한 π(n)의 값을 계산해야 할 필요가 있다. 여기에서는 정확한 값을 구하는 알고리즘들을 정리한다.
알고리즘
- 정확한 소수 개수를 카운팅하는 알고리즘에 대해서는 통합적이고 완전하게 정리된 자료를 찾기가 어려웠다.
- 한국어로 적힌 가장 자세하고 친절한 자료는 https://rkm0959.tistory.com/189 에 있다. 여기에서는 Lucy-Hedgehog Algorithm의 O(n^(3/4)) 구현법과 O(n^(2/3)) 구현법. Meissel-Lehmer Algorithm 의 그리고 O(n^(2/3)) 구현법을 다룬다.
- 하지만, 여기에서 Lucy-Hedgehog 알고리즘이라고 불리는 이 방법은, 학계에서 발표된 방법이 아니라 PS 커뮤니티 (프로젝트 오일러)쪽에서 처음 전래된 방법이라서, 학술적인 자료에서는 언급을 찾을 수가 없다.
- 그렇다. Lucy-Hedgehog은 Lucy와 Hedgehog 이라는 두 사람의 이름이 아니라, Lucy_Hedgehog이라는 유저의 닉네임이다. Lucy Hedgehog은 실제로는 소닉에 등장하는 캐릭터 이름이고..
- 다른 외국 블로그 중에서는 Lucy's algorithm이라고 언급하는 곳도 있다.
- 게다가 인터넷에 돌아다니는 구현법도 매우 다양하고, 시간복잡도 분석도 다양하고 (rkm의 자료에서는 시간복잡도를 설명할때 로그텀을 떼고 지수승만 고려해서 표기했다. 자료들에 따라서는 로그텀까지 고려해서 분석한 경우도 있다.), 그래서 자료들간에 매칭이 잘 안되는 경우가 많다.
- 우선 학술적인 관점에서 알고리즘을 정리해보자
- 위키피디아에 항목이 있는 알고리즘은 Meissel–Lehmer_algorithm 이다. 당연히 출처가 되는 논문들도 모두 정리되어있다.
- 소수 카운팅 세계신기록을 갖고있는 프로그램의 개발자인 kimwalisch가 요약한 소수 카운팅 알고리즘의 역사는 아래처럼 된다고 한다.
- 다른 디스커션들에서 찾은 요약은 https://github.com/Bodigrim/arithmoi/issues/60 이다. https://github.com/vitaly-t/prime-lib/discussions/6 도 참고.
- 1798 : Legendre ⇒ O(n^(3/4))
- 1870 : Meissel ⇒ O(n/log^3(n))
- 1870 : Lehmer ⇒ O(n/log^4(n))
- 1985 : LMO ⇒ O(n^(2/3))
- 1996 : M. Deleglise, J. Rivat ⇒ O(n^(2/3)/log^2(n))
- 2001 : Xavier Gourdon ⇒ O(n^(2/3)/log^2(n)) 의 상수 최적화
- 2015 : Douglas Staple ⇒ O(n^(2/3)/log^2(n)) 의 상수 최적화
- 사실 관련 서적이나 논문을 찾아보면 가장 정확한 히스토리를 알 수 있겠지만.. 그냥 이 정도의 흘러다니는 내용을 갖고 정리하면, 르장드르의 공식을 시작으로 해서, Miessel과 Lehmer은 그 공식을 좀 더 변형한 공식을 제안했고, 이후의 LMO 는 Lehmer의 공식을 알고리즘적으로 효율적으로 구현하는 방식을 제안한 것 같다. 그리고 Deleglise-Rivat 과 Gordon은 LMO의 큰 틀에서 부분방식을 조금 개선한 정도. 그래서 Meissel-Lehmer 알고리즘이라 하면 Lehmer의 공식을 LMO의 최적화를 적용해서 구현하는 방식을 말하는듯 하다.
- 그리고 Lucy_hedgehog의 알고리즘은 Lehmer의 공식이 아닌, 르장드르의 공식에 LMO의 구현방식을 결합시킨 알고리즘인것 같다.
- 그래서 학술적으로는 새롭거나 큰 가치가 있지는 않을것 같다.
- 그러면 우리가 사용해야 할 알고리즘은 무엇인가.
- 다른 알고리즘들과 마찬가지로, 이론적으로 가장 빠른 알고리즘과 PS에서 가장 효율적인 알고리즘은 차이가 있다.
- 상수항을 무시하는 점근적 복잡도와 PS의 인풋 범위에서의 실제 수행 속도의 차이, 구현 난이도와 코드 길이에 대한 고려, 메모리 사용량에 대한 고려 등등의 추가적인 요소들이 영향을 미친다.
- 보통 PS에서 주어지는 최대 x는 10^12 정도이다. 이때, x^(3/4) 와 x^(2/3)은 10배 차이밖에 안난다. 이 정도면 앞의 분석에서 종종 무시되었던 log팩터 하나가 붙어도 비슷해지고, 구현의 효율성에 따라서 얼마든지 실제 성능은 뒤집힐수 있는 차이이다.
- 그러면 '이해와 구현이 가장 간단한 알고리즘'과 'PS범위에서 가장 빠른 구현'을 알아두면 충분할 것 같다.
- 여기에서 '이해와 구현이 가장 간단한 알고리즘'은 Lucy_Hedgehog 알고리즘이다.
Lucy_Hedgehog algorithm
- 학계에서 쓰이는 이름은 아니지만, 그냥 PS 커뮤니티에 알려진대로 Lucy_Hedgehog algorithm 이라고 부르겠다.
- 어차피 PS커뮤니티에서는 그냥 유저의 닉네임을 갖고 이름을 붙이는 경우도 흔하다 (키타마사법, xudyh's sieve, Arpa's trick 등등). 그게 사실은 원래 존재하던 알고리즘을 재발견하거나 재정리한 것이었던 경우도 있고.
- 알고리즘 설명은 https://rkm0959.tistory.com/189 에 상세하고 친절하게 나와있으므로 그쪽을 참고.
- 이 방법의 장점은, 이해가 간단하고 소수의 개수뿐 아니라 다른 값들을 구하는 용도로도 쉽게 변형해서 사용할 수 있다는 것이다.
- 좀 더 정확히 표현하면, 어떤 함수 $f$가 completely multiplicative function 이고, $ \sum_{i \leq x} f(i) $ 를 쉽게 구할수 있다면, 이 알고리즘을 변형해서 $ \sum_{p \leq x, p is prime} f(p) $ 을 쉽게 구할수 있다.
- 대표적으로 x이하 소수의 합, x이하 소수의 제곱의 합 등등을 구할수 있다.
- dp 테이블의 초기값을 $ S_f(v, 1) = f(2) + f(3) + \ldots + f(v) $ 으로 해주고, 점화식을 $ S_f(v, p) = S_f(v, p-1) - f(p)\left[S_f(v/p, p-1) - S_f(p-1, p-1)\right] $ 형태로만 바꿔주면 된다.
구현
기본 구현
- https://rkm0959.tistory.com/189의 처음 방법을 dict로 dp table을 사용해서 그대로 구현하면 이런 형태가 된다. 굉장히 깔끔해서 코드를 이해하고 변형하기에도 쉽다.
def count_prime_v0(n): sq = math.isqrt(n) beg = sq - 1 if n // sq == sq else sq vals = [n // x for x in range(1, sq + 1)] + list(range(beg, 0, -1)) dp = {v: v - 1 for v in vals} for m in numtheory.prime_list(sq): for v in vals: if m * m > v: break dp[v] -= dp[v // m] - dp[m - 1] return dp[n]- 로컬에서 n=10^10 을 계산하는데에 1.998s, n=10^11을 계산하는데에 12.233초가 걸렸다.
- 이 코드를, dict 대신에 두개의 list에 dp값을 저장하도록 바꾸고, 약간의 구현 최적화를 거치면 다음처럼 된다.
- 그리고 굳이 sqrt(x)범위에 대한 sieve를 미리 돌려서 소수들을 구해놓지 않더라도, 알고리즘 내에서 m값의 소수여부를 알 수 있기 때문에 sieve가 필요가 없다는 것을 알게되었다. (속도에는 별 차이 없다)
def count_prime_v1(n): sqrt = math.isqrt(n) a = [0] + list(range(sqrt)) # a[i] = dp[i] b = [0] + [n // x - 1 for x in range(1, sqrt + 1)] # b[i] = dp[n//i] for m in range(2, sqrt + 1): pi = a[m - 1] if a[m] > pi: # m is prime m_sq = m * m for j in range(1, sqrt + 1): v = n // j if v < m_sq: break x = v // m b[j] -= (a[x] if x <= sqrt else b[n // x]) - pi for v in range(sqrt, 0, -1): if v < m_sq: break a[v] -= a[v // m] - pi return b[1]- 1010 ⇒ 1.520s, 1011 ⇒ 7.921s
- 이정도만 구현해도 pypy로 제출할 경우 웬만큼은 통과된다. 펜윅트리를 이용한 최적화도 굳이 필요없다.
최적화
- 실제로 최대한 빠른 구현체를 만드는 것이 목표이다.
- 우선 현재 있는 구현체들을 찾아보자.
- rkm의 블로그에서 추천하는 것은 codeforce의 Four Divisor 문제에 제출된 소스코드들을 찾아보는 것이다. 하지만 현 시점에서는 Library Checker에서 찾아보는게 더 간단하다
- 찾아보면 가장 빠른 구현체는 (cpp, python, pypy3를 막론하고) 는 이와 유사한 형태의 코드이다
- 파이썬 버전은 https://judge.yosupo.jp/submission/34113. 이보다 더 빠른 파이썬 코드들은 numpy를 이용하는 버전이다
- 코드를 열어봤을때 smalls/larges/roughs의 세개의 배열을 선언해서 쓰는 코드는, 전부 이 방식의 변형이라고 추측하고 있다
- 문제는 이 코드가 어떤 원리로 돌아가는 코드인지에 대해 설명된 곳이 없었고, 출처가 어디인지도 찾을수 없었다. 펜윅트리 같은 것을 안쓴것은 확실한데.. rkm 블로그에서 마지막에 언급하는 '매우 빠른 최적화 방법'이 여기에 해당되는 것인지 다른 방법인지도 알 수 없었고, 이게 lucy_hedgehog 알고리즘인지, Miessel-Lehmer 알고리즘인지 조차도 명확하게 알수 없었다.
- https://judge.yosupo.jp/submission/141200 에 있는 커멘트에 따르면 시간복잡도는 O(N^(3/4) / logN) 이라고 한다
- 결국 코드를 뜯어보고 알고리즘을 분석하기로 했다..
- Library Checker에서 python 1등인 abUma의 코드와 pypy 1등인 fuyuru의 코드를 가져와서 분석했다. 둘다 기본적으로는 동일한 방식이고 구현 디테일만 다르다. 아마 같은 cpp 코드를 각자 방식대로 파이썬으로 번역한 느낌이다.
- prime_count_abUma : 0.367 / 1.728
- prime_count_fuyuru : 0.347 / 1.62
- 기본적으로는 lucy_hedgehog 알고리즘에 몇가지 최적화를 추가한 형태이다.
- DP 테이블은 2개의 배열에 저장된다. 1,2,.., sqrt(n)의 dp값을 저장하는 배열이 small, n/1, n/2, …, n/sqrt(n) 의 dp값을 저장하는 배열이 large이다
- 3가지 최적화가 들어간다.
- large값 갱신 최적화
- large 값을 갱신할때, 이후에 더이상 사용되지 않을 값들은 갱신하지 않는 것이다.
- 더이상 사용되지 않을 값들은
- 이것이 rough에 저장되는 값이다. 변수명이 rough인것은 이런 수를 부르는 rough number 라는 용어가 있기 때문이다
- 구현된 방식은 한번 바깥 루프를 돌때마다, 에라토스테네스 체의 한 스텝을 진행하면서 수들을 제거하고 이것들을 skip 배열에 기록한다. skip 배열을 기반으로, 제거되지 않은 수들을 모아서 rough 배열을 새로 갱신한다.
- 사실 cpp라면 이 방식이 효율적이긴 할텐데, 코드가 많이 복잡해지고 파이썬에서 이 구현을 그대로 옮겨 쓰기는 좀 번거롭다. 파이썬에서 이 아이디어를 간단하게 구현하는 방법은, large 배열을 그냥 다시 dict로 만들고, rough 배열을 갱신하는 것은 dict에서 해당하는 key들을 제거하는 방식으로 하는 것이다.
- dict는 처음에 key값들을 정렬된 순서대로 넣어줬다면, key들을 삭제해도 여전히 정렬된 순서로 이터레이션이 가능하다.
- 아래처럼 간단하게 구현된다.
def count_prime_v2(n): sqrt = math.isqrt(n) a = [0] + list(range(sqrt)) # a[i] = dp[i] b = {x: n // x - 1 for x in range(1, sqrt + 1)} # b[i] = dp[n//i] for m in range(2, sqrt + 1): pi = a[m - 1] if a[m] > pi: # m is prime m_sq = m * m for j in b: v = n // j if v < m_sq: break x = v // m b[j] -= (a[x] if x <= sqrt else b[n // x]) - pi for x in range(m, sqrt, m): if x in b: del b[x] for v in range(sqrt, 0, -1): if v < m_sq: break a[v] -= a[v // m] - pi return b[1]- 1010 ⇒ 0.626s, 1011 ⇒ 3.675s
- small값 갱신 최적화
- 간단한 아이디어이다
- %%i에 대해서 루프를 돌면서 a[i] -= a[i m] - pi 를 수행하게 되는데, 이때 빼줘야 하는 a[i m] - pi 값은 [j*m, j*(m+1) - 1) 범위의 모든 i에 대해서 동일하다. 그래서, 빼줘야 하는 값을 한번씩만 계산하고, 그 값을 빼줘야 하는 구간에 루프를 돌면서 업데이트를 하는 식이다. 결과적으로는 이중루프로 구현되지만, 전체 반복횟수는 똑같고, 나눗셈 연산을 비롯한 중복된 계산들을 제거시킬수 있다.
- 2-wheel factorization
- 2를 제외한 소수는 모두 홀수이므로, 2에 대해서 먼저 전처리 후에, 루프를 돌때에는 홀수에 대해서만 도는 것. 에라토스테네스 체를 비롯해서 소수관련 알고리즘을 구현할때 흔히 사용되는 방법이다.
- 이것들을 모두 추가한 코드는
def count_prime_o3(n): sqrt = math.isqrt(n) a = [0] + list(range(sqrt)) # a[i] = dp[i] for v in range(sqrt, 3, -1): a[v] -= a[v // 2] b = {j: (n // j) - (n // (j + j)) for j in range(1, sqrt + 1, 2)} for m in range(3, sqrt + 1, 2): pi = a[m - 1] if a[m] > pi: # m is prime m_sq = m * m for j in b: v = n // j if v < m_sq: break b[j] -= (a[v // m] if m * j > sqrt else b[m * j]) - pi for x in range(m, sqrt, m): if x in b: del b[x] for j in range((sqrt - 1) // m, m - 1, -1): c = a[j] - pi e = min((j + 1) * m, sqrt) for i in range(j * m, e): a[i] -= c return b[1]
- 여기에 추가적으로 적용할수 있는 최적화가 있다.
- large 배열에서 key를 제거할때, n / (m*m)을 넘어가는 값들은 굳이 제거하지 않는것이다. 어차피 key값을 이터레이션 할때 그 범위를 넘어가면 종료하기 때문에, 굳이 제거하지 않아도 상관이 없다.
- 이것까지 포함시키면 이제 1.810 으로 진짜로 실제로 존재하는 가장 빠른 구현체들과 근접한 성능이 나온다.
- teflib에서는 이 외에 사소한 최적화들을 포함시켜서 적용하였다
- large 배열을, key는 계속 dict에 저장하되, 값을 따로 list에 저장하게 해서, 원소의 억세스 타임을 조금 더 줄였다.
- 루프에서 탈출 조건을 체크하는 부분을 조금 최적화해서 나눗셈 연산을 줄였다.
- 등등
구현체들..
3 - 둘다 기본적으로는 동일한 방식. 구현 디테일만 다르다. c++의 최적 코드를 파이썬으로 번역한 느낌
prime_count_v0_0: TLE - 알고리즘 그대로 옮김. 토글링도 안함
count_prime_v0: 2.160 / 12.182s 토글링 적용. dict로 dp 테이블을 저장
count_prime_v1: 1.520 / 7.921s 리스트로 dp 테이블 저장. 에라체 전처리 제거.
count_prime_lucy: 1.263s / 6.664s V1을 기술적으로 최적화.
count_prime_o1(n) 0.626s / 3.675s v1에서 불필요한 b배열 계산 제거. b배열은 dict로 구현
count_prime_o2 ??? / 2.941s o1에서 2-wheel 적용
count_prime_o3 ??? / 2.489s o2에서 small 배열 계산 최적화
count_prime_o4 ??? / 1.828s o3에서 del 범위 줄임
ps/소수의_개수_구하기.txt · 마지막으로 수정됨: 2024/01/11 01:17 저자 teferi
토론