Lucy-Hedgehog Algorithm

이 알고리즘은 에라토스테네스의 체를 더 효율적으로 돌리기 위한 알고리즘으로 다음과 같은 점화식을 이용한다. 

$$dp[v][m]=dp[v][m-1]-(dp[\lfloor\frac{v}{m}\rfloor][m-1]-dp[m-1][m-1])$$

여기서 $dp[v][m]$은 $m$까지의 수로 체를 돌렸을 때 1부터 $v$까지 수 중 소수로 간주되는 것이 몇 개인지 말하고 있다. 우리가 $m-1$까지의 수를 보았다고 하자. 만약 $m$이 소수가 아니라면 체를 돌릴 필요가 없으니 $dp[v][m]=dp[v][m-1]$이다. 또 $m^{2}$이 $v$보다 크다면 역시 체를 돌릴 필요가 없어 $dp[v][m]=dp[v][m-1]$이다. 그럼 $m^{2}<v$이며 $m$이 소수인 경우를 보자. $dp[v][m-1]$에서 제외해야 할 수는 $m$의 배수들임이 자명하다. 더 나아가 $p$가 소수라고 간주되는 것들이고 $m\cdot p$ 꼴이어여 한다. 왜냐면 소수라 간주되지 않는 수들은 이미 앞에서 걸러졌기 때문이다. 또, $p>m-1$이어야 한다. 이것은 소수이지만 이미 앞에서 걸렀다. 따라서 $p>m-1$이고 $p$가 소수인 개수를 세야 한다. $dp[\lfloor\frac{v}{m}][m-1]$은 $p$ 중 소수인 것들을 가지고 있다. 여기서 $dp[m-1][m-1]$ 즉, $m-1$까지 소수의 개수를 뺀다면 $p>m-1$인 소수 $p$의 개수를 알 수 있다. 따라서 저러한 식이 유도되는 것이다. 

저 점화식에 dp를 적용하기 위해서는 총 $2\sqrt{N}$개의 값이 필요하다. 먼저 $m$의 가능한 값들을 생각하자. $v$는 최대 $N$일 것인데 $m^{2}<v$이니 $m$은 $\sqrt{N}$까지 만 생각하면 된다. 따라서 첫 번째 인자로 생각할 수 있는 것은 $1,\cdots,\sqrt{N}$과 $\lfloor\frac{N}{1}\rfloor,\cdots,\lfloor\frac{N}{\sqrt{N}}\rfloor$ 이다.

이제 2개의 배열을 만들고 하나는 $1,\cdots,\sqrt{N}$ 부분을 다른 하나는 $\lfloor\frac{N}{1}\rfloor,\cdots,\lfloor\frac{N}{\sqrt{N}}\rfloor$ 부분을 관리하게 하자. $m$을 1부터 $\sqrt{N}$까지 증가시키는 방향으로 구현한다. 각 $m$마다 새롭게 $dp[v][m]$을 갱신하는 경우는 $v$가 큰 것 먼저 값을 바꿔야 데이터 손실이 없다. $dp[m-1][m-1]$ 부분은 사실 $m-1$이하의 소수의 개수를 의미하기 때문에 미리 $\sqrt{N}$까지 에라토스테네스의 체를 돌려 전처리를 해 두자.

아래에는 구현한 코드가 있다.

 

보려면 클릭

#include<stdio.h>
#include<math.h>
#define MAX 333333
typedef long long ll;
ll num_pr[MAX], A[MAX], B[MAX]; 
// A는 dp[1,2,...\sqrt{n}]을 관리, B는 dp[N/1, N/2,...,\sqrt{n}]을 관리. 
// B는 항상 \sqrt{n}이상인 것을 포함한다.  
// A, B의 초기 값은 dp[v][1] 인데 1을 탐색했다고 보면 dp[v][1]=v-1이다. 
// dp[v][m] 에서 만약 m ^ 2 > v 거나 m이 소수가 아니면 더 이상 진전 될 것이 없다.  그래서 m을 \sqrt{N} 까지 돌린다.  v의 최댓값이 N이기 때문. 
// top-down으로 관리해야 데이터 손실이 없다. B먼저, 그 다음 A. 
// i = m 이고 j = v 라고 생각하면 편하다. 
ll ispr[MAX];
ll N, C;
 
ll get_val(ll x){
    if(x <= C) return A[x];
    return B[N / x];
}
 
void calc_primes(){
    ll i, j;
    for(i = 1;i <= C;i++) A[i] = i - 1;
    for(i = 1;i <= C;i++) B[i] = N / i - 1;
    for(i = 2;i <= C;i++){
        if(ispr[i]) continue;
        for(j = 1;j <= C;j++){
            if(N / j < i * i) break;
            B[j] = B[j] - (get_val(N / j / i) - num_pr[i - 1]); 
        }
        for(j = C;j >= 0;j--){
            if(j < i * i) break;
            A[j] = A[j] - (get_val(j / i) - num_pr[i - 1]);
        }
    }
}
int main(){
    ll i, j;
    N = 1e11; C = sqrt(N);
    ispr[1] = 1;
    for(i = 2;i <= C;i += 2){
        if(ispr[i]) continue;
        for(j = i * i;j <= C;j += i){
            ispr[j] = 1;
        }
        if(i == 2) i = 1;
    }
    for(i = 2;i <= C;i++){
        num_pr[i] = num_pr[i - 1];
        if(ispr[i]) continue;
        num_pr[i]++;
    }
    calc_primes();
    printf("%lld", B[1]);
}

 

이를 응용한 예시로 $4n+1$꼴 소수의 개수를 구하는 프로그램을 만들 수 있다.

알고리즘의 핵심 원리는 모든 $4n+3$ 꼴 수를 보존하고 나중에 한꺼번에 제외하는 것이다.

아래의 dp식을 보자.

$$dp4[v][m]=\begin{cases} dp4[v][m-1] - (dp4[\lfloor\frac{v}{m}\rfloor][m-1]-dp1[m-1][m-1]-(\lfloor\frac{v}{m}\rfloor-3)/4-1), &m\equiv 1\pmod{4}\\  dp4[v][m-1] - ((dp[\lfloor\frac{v}{m}\rfloor][m-1] - dp4[\lfloor\frac{v}{m}\rfloor][m-1] + (\lfloor\frac{v}{m}\rfloor-3)/4+1)-(dp[m-1][m-1]-dp1[m-1][m-1])), &m\equiv 3\pmod{4}\end{cases}$$

이때 $dp4$는 에라토스테네스의 체를 $m$까지 보았을 때 남아있는 $v$이하의 $4n+1$꼴 소수들과 모든 $4n+3$꼴 소수들의 개수이다. $dp1$은 $m-1$까지 $4n+1$꼴 소수의 개수이다. 

$m$이 $4n+1$ 꼴 소수인 경우를 생각하자. 같은 원리로 $m\cdot p$ 꼴 수를 걸러내야 하는데 $4n+3$꼴은 반드시 보존해야 하므로 $p$가 $4n+1$꼴 소수인 것만 생각한다. $dp4[\lfloor\frac{v}{m}\rfloor][m-1]$가 $\lfloor\frac{v}{m}\rfloor$ 까지의 $4n+3$꼴 수들을 포함하므로 정확히 $4n+1$꼴 수들만 모으기 위해 $-(\lfloor\frac{v}{m}\rfloor-3)/4-1$항을 추가한 것이다. 

$m$이 $4n+3$ 꼴 소수인 경우를 생각하자. 이번에도 $m\cdot p$ 꼴 수를 거르는데 $4n+1$ 꼴 합성수를 걸러야 겠다. 그럼 $p$가 $4n+3$ 꼴 소수이므로 $dp[\lfloor\frac{v}{m}\rfloor][m-1]$ 즉, $m-1$까지 보았을 때 소수라고 간주되는 것들 중에서 $dp4[\lfloor\frac{v}{m}\rfloor][m-1]$의 $4n+1$꼴 소수라고 간주되는 것들을 빼면 $4n+3$꼴 소수라고 간주되는 것들의 개수를 얻을 수 있다. $p$가 $m-1$ 초과이므로 그 뒤의 항이 추가된 것이다. 

이 알고리즘은 소수의 개수를 관리함과 동시에 $4n+1$꼴 소수의 개수도 관리하여 시간이 더 오래 걸린다. 더 효율적인 방법이 있으면 알려주길 바란다!^^ 아래에는 구현한 코드가 있다.

 

보려면 클릭

#include<stdio.h>
#include<math.h>
#define MAX 333333
typedef long long ll;
ll num_pr[MAX], A[MAX], B[MAX]; 
ll A4[MAX], B4[MAX];
// A는 dp[1,2,...\sqrt{n}]을 관리, B는 dp[N/1, N/2,...,\sqrt{n}]을 관리. 
// B는 항상 \sqrt{n}이상인 것을 포함한다.  
// A, B의 초기 값은 dp[v][1] 인데 1을 탐색했다고 보면 dp[v][1]=v-1이다. 
// dp[v][m] 에서 만약 m ^ 2 > v 거나 m이 소수가 아니면 더 이상 진전 될 것이 없다.  그래서 m을 \sqrt{N} 까지 돌린다.  v의 최댓값이 N이기 때문. 
// top-down으로 관리해야 데이터 손실이 없다. B먼저, 그 다음 A. 
// i = m 이고 j = v 라고 생각하면 편하다. 
ll ispr[MAX];
ll num_pr1[MAX];
ll N, C;
 
ll get_val_4(ll x){
    if(x <= C) return A4[x];
    return B4[N / x];
}
 
ll get_val(ll x){
    if(x <= C) return A[x];
    return B[N / x];
}
//4n+3을 보존하는 방식
//A4 또는 B4는 주어진 수 v이하의 모든 4n+3꼴 수와 4n+1 꼴 semi 소수(에라토스테네스의 체를 m까지 돌렸을 때 소수라 판정된 수)를 포함한다.  
//목표는 결과 가 나올 때까지 4n+3꼴 모든 수를 보존하는 것. 이후 한번에 제거하는 것이다.  
void calc_primes(){
    ll i, j;
    for(i = 1;i <= C;i++) A[i] = i - 1 - i / 2 + 1;
    for(i = 1;i <= C;i++) B[i] = N / i - 1 - N / i / 2 + 1; //2 배수 미리 처리함. 
    for(i = 1;i <= C;i++) A4[i] = i - 1 - i / 2;
    for(i = 1;i <= C;i++) B4[i] = N / i - 1 - N / i / 2;
    
    for(i = 3;i <= C;i++){
        if(ispr[i]) continue;
        for(j = 1;j <= C;j++){
            if(N / j < i * i) break;
            if(i % 4 == 1) B4[j] = B4[j] - (get_val_4(N / j / i) - num_pr1[i - 1] - (N / j / i - 3) / 4 - 1);
            //(4n+1)*(4n+1)꼴 합성수를 제외하는 과정이다.  
            //N/j/i까지에는 4n+3꼴의 모든 수와 4n+1꼴 소수가 들어 있다. num_pr1[i-1]은 i-1까지 4n+1꼴 소수가 들어 있다. 
            //따라서 N/j/i까지 4n+3꼴 수들을 빼면 i-1보다 큰 4n+1꼴 소수들을 얻을 수 있다. 
            else B4[j] = B4[j] - ((get_val(N / j / i) - get_val_4(N / j / i) + (N / j / i - 3) / 4 + 1) - (num_pr[i - 1] - num_pr1[i - 1]));
            // (4n+3)*(4n+3) 꼴 합성수를 제외하는 과정이다. 
            //N/j/i까지 소수들에서 N/j/i까지 4n+1꼴 소수들을 뺴면 4n+3꼴 소수들만 얻을 수 있다. get_val_4(N/j/i)에는 4n+3꼴 수들도 있으니 제외해주자. 
            //num_pr1[i-1]은 i-1까지 4n+1꼴 소수들만 모았다.
            //****주의 사항 : get_val_4(N / j / i)에서 아직 남아있는 4n+1꼴 합성수가 있다. 하지만 get_val(N / j / i)에도 남아 있는 합성수가 있기 때문에 상쇄된다.
            //                  이것이 A,B와 A4,B4를 같이 iterating 하는 이유이다.  
        }
        for(j = C;j >= 1;j--){
            if(j < i * i) break;
            if(i % 4 == 1) A4[j] = A4[j] - (get_val_4(j / i) - num_pr1[i - 1] - (j / i - 3) / 4 - 1);
            else A4[j] = A4[j] - ((get_val(j / i) - get_val_4(j / i) + (j / i - 3) / 4 + 1) - (num_pr[i - 1] - num_pr1[i - 1]));
        }
        
        for(j = 1;j <= C;j++){
            if(N / j < i * i) break;
            B[j] = B[j] - (get_val(N / j / i) - num_pr[i - 1]);
        }
        for(j = C;j >= 1;j--){
            if(j < i * i) break;
            A[j] = A[j] - (get_val(j / i) - num_pr[i - 1]);
        }
    }
}
int main(){
    ll i, j;
    scanf("%lld", &N); 
    C = sqrt(N);
    ispr[1] = 1;
    for(i = 2;i <= C;i += 2){
        if(ispr[i]) continue;
        for(j = i * i;j <= C;j += i){
            ispr[j] = 1;
        }
        if(i == 2) i = 1;
    }
    for(i = 2;i <= C;i++){
        num_pr[i] = num_pr[i - 1];
        num_pr1[i] = num_pr1[i - 1];
        if(ispr[i]) continue;
        num_pr[i]++;
        if(i % 4 == 1) num_pr1[i]++;
    }
    calc_primes();
    printf("%lld", B4[1] - (N - 3) / 4 - 1);
}