Implementation of Lucy-Hedgehog algorithm

우리는 이전 글에서 Lucy-Hedgehog algorithm(이하 Lucy)에 대해 배웠다. Lucy는 다음의 점화식을 이용한다.

$\[S(v,p)=S(v,p-1)-(S(v/p,p-1)-S(p-1,p-1))\]

우리는 앞서 2개의 배열 A,B를 사용하여 목표 값인 $S(N, sqrt(N))$을 계산하였다. 여기서 확인해야할 것은 $S(v,p)$가 나타내는 값은 2부터 $v$까지의 값 중 $p$까지 에라토스테네스의 체를 돌렸을 때 걸러지지 않은 값들의 개수라는 것이다. 

최적화를 위해서 우리는 $\sqrt{N}\le y<N$을 만족하는 $y$를 일단 아무거나 선택한다. 우리의 배열 A, B들은 각각 $1, 2, 3, \cdots,\sqrt{N}$, $N/1,N/2,\cdots,N/\sqrt{N}$값을 지니고 있었다. 이들을 key들이라 부르자. 이 최적화는 $y$보다 큰 값을 가지는 key만 Lucy의 공식을 이용하고 나머지는 직접 계산한다는 것이다. 즉, 우리는 S(v, p)를 계산함과 동시에 에라토스테네스의 체를 돌리면서 key갑이 y보다 작은 값들에 대해서 직접 대입해준다. 우리가 $p$에 대한 에라토스테네스의 체를 돌리고 있다고 가정하자. 만약 $v\le y$가 성립하면 $S(v,p)$의 값을 직접 구해 바로 대입한다는 것이다. 이렇게 $S(v,p)$의 값을 얻어내기 위해서는 우리는 빠르게 $v$이하의 소수로 추측되는 값이 몇 개 있는지 확인해야 한다. 

이를 하기 위해서 Fenwick Tree를 사용한다. 자연수 각각에 1/0을 대응시키자. 1이 대응되었다는 것은 소수로 추측된다는 것이고 0이 대응되는 것은 합성수라는 것이다. 이들을 누적합으로 생각하면 $v$이하의 소수로 추측되는 값을 쉽게 계산할 수 있다. 따라서 에라토스테네스의 체에서 $j$가 걸러지는 경우 Fenwick Tree의 $j$번째 인덱스의 값이 0으로 변한다는 것이다. 본 글에서는 sieve라는 이름으로 Fenwick Tree를 구현하였다. 

특히, 계산의 편의를 위해서 sieveRaw라는 배열을 만들었으며 이는 $i$라는 자연수가 합성수이면 $i$번째 인덱스가 true이고 그렇지 않으면 false라고 마킹한다. 

sieve와 sieveRaw는 모두 에라토스테네스의 체가 반영되면서 점차 업데이트 되어간다. 

배열 A가 포함하는 인덱스는 모두 $y$보다 작거나 같으므로 배열 B만 업데이트 한다. 

 

따라서 알고리즘을 정리해보자면 다음과 같다. 

1. sieve라는 이름의 Fenwick Tree를 구현하는데 1을 제외한 $y$까지의 자연수에 모두 1을 대응시켜 놓는다.

2. sieveRaw라는 배열은 1을 True 나머지를 false라고 둔다. 

3. $p$를 2부터 $sqrt(N)$까지 돌리는 에라토스테네스의 체를 실행한다.

3a. $p$가 합성수이면 다음으로 넘어간다.

3b. $S(p-1, p-1)$의 값은 $p-1$까지의 소수의 개수이므로 Fenwick Tree로 계산해준다. 저번 Lucy의 포스팅에서도 언급했듯이 $v^2>p$인 경우만 계산하면 된다. 우리가 계산해야할 것은 $S(N/i, p)$이다. 따라서 $i\le N/p^{2}$를 만족한다. 또한, 우리는 $v> y$인 경우만을 고려하므로 $i\le N/y$가 성립한다. 따라서 $i$를 1부터 $\min(N/p^{2},N/y)$까지 iterate해주며 B[i]를 계산한다. 계산 과정에서 $N/i*p \le y$인 경우가 있다. 이 경우 Fenwick Tree를 이용해 계산해주고 그 외의 경우는 B[i * p]로 계산한다. 이 과정이 s0라는 함수에 나타난다. 

3c. $j=p*p$로 두고 에라토스테네스의 체를 돌려준다. 새롭게 합성수가 된 수를 Fenwick Tree에 반영해준다. 

 

이를 구현한 코드는 다음과 같다. 

코드 보기

#include<bits/stdc++.h>
#define MAX 10000001
using namespace std;
typedef long long ll;
 
int sieve(MAX);
bool sieveRaw[MAX];
ll B[MAX];
ll N, C, y;
ll num[MAX], pri[MAX];
 
void update(ll idx, ll value){
    while(idx <= y){
        sieve[idx] = sieve[idx] + value;
        idx = idx + (idx & -idx);
    }
}
 
ll sum(int idx){
    ll ret = 0;
    while(idx > 0){
        ret = ret + sieve[idx];
        idx = idx - (idx & -idx);
    }
    return ret;
}
 
void make_fenwick_tree(){
    for(ll i = 2;i <= y;i++){
        update(i, 1);
    }
}
 
ll s0(ll v){
    if(v <= y) return sum(v);
    return B[N / v]; 
}
 
void calc_primes(){
    ll i, j;
    for(ll p = 2;p <= C;p++){
        if(sieveRaw[p]) continue;
        ll sp = sum(p - 1);
        ll lim = min(N / y, N / (p * p));
        for(i = 1;i <= lim;i++){
            B[i] -= s0(N / (i * p)) - sp;
        }
        j = p * p;
        while(j <= y){
            if(sieveRaw[j]){
                j += p;
                continue;
            }
            sieveRaw[j] = true;
            update(j, -1);
            j += p;
        }
    }
}
 
void init(){
    ll i;
    for(i = 0;i < MAX;i++) sieve[i] = 0;
    make_fenwick_tree();
    sieveRaw[1] = true;
    for(i = 2;i <= y;i++) sieveRaw[i] = false;
    for(i = 1;i <= C;i++){
        B[i] = N / i - 1;
    }
}

 

여기에서 좀 더 최적화를 해보자. 우리는 B[i]의 모든 값이 필요한게 아니다. B[1]만 계산하면 N이하의 소수의 개수를 셀 수 있다. $p_{i}$를 $\sqrt{N}$보다 작거나 같은 최대의 소수라고 하자. 그럼 B[1]을 계산하기 위해서는 $S(N, p_{i-1}), S(N/p_{i}, p_{i-1})$을 계산해야한다. 이들을 계산하기 위해서는 $S(N/p_{i-1}, p_{i-2}), S(N/p_{i},p_{i-2}), S(N/p_{i}p_{i-1},p_{i-2})$들을 계산해야한다. 따라서 우리가 실제로 계산해야하는 값들은 $p$에서 에라토스테네스의 체를 하고 있다고 할 때 $p$보다 큰 소수 $q$들에 대해서 $S(N/q, p)$인 것이다. 따라서 calc_primes 부분은 다음과 같이 최적화할 수 있다. 

 

코드 보기

void calc_primes(){
    ll i, j;
    for(ll p = 2;p <= C;p++){
        if(sieveRaw[p]) continue;
        ll sp = sum(p - 1);
        ll lim = min(N / y, N / (p * p));
        B[1] -= s0(N / p) - sp;
        for(i = p + 1;i <= lim;i++){
            if(sieveRaw[i]) continue;
            B[i] -= s0(N / (i * p)) - sp;
        }
        j = p * p;
        while(j <= y){
            if(sieveRaw[j]){
                j += p;
                continue;
            }
            sieveRaw[j] = true;
            update(j, -1);
            j += p;
        }
    }
}