Linear-time construction of the sieve of Eratosthenes

에라토스테네스의 체는 소수를 구할 때 흔히 쓰는 방식이다. 구현도 쉽기 때문에 정수론과 관련된 PS 문제에서 자주 만날 수 있다.

 

일반적으로는 아래와 같이 구현한다.

vector<int> sieve_of_eratosthenes(int n) {
    vector<bool> sieve(n + 1);
    vector<int> prime;
    for (int k = 2; k <= n; k++) {
        if (sieve[k]) {
            continue;
        }
        prime.push_back(k);
        for (int i = 2; i * k <= n; i++) {
            sieve[i * k] = true;
        }
    }
    return prime;
}

이 알고리즘은 얼마나 효율적일까? 간단하게 각 정수를 평균적으로 몇 번 방문하는지 보면 된다. 이 방식은 $n$까지의 모든 양의 정수의 방문 횟수가 각 정수의 소인수 개수이다. 어떤 정수의 소인수 개수를 알려주는 함수로 prime omega function이라는 것이 있는데, $n$의 소인수 개수는 평균적으로 $O(\log\log n)$개이다.

 

여기서 방법을 살짝 바꿔서 각 정수를 정확히 한 번만 방문하게 할 수 있다. 각 정수의 소인수 하나당 방문을 하는 대신 최소 소인수에 대해서만 방문을 하는 것이다. 아이디어는 안쪽 루프에서 현재 보고 있는 메인 루프의 정수 $k$의 배수를 방문하되 각 소수 $p$에 대해 $pk$꼴의 수만 방문하고, 작은 $p$부터 방문하면서 $k$가 $p$의 배수가 되면 멈추는 것이다. 이렇게 하면 $p$가 $pk$의 최소 소인수일 때만 방문을 하게 되므로 우리의 목적을 달성할 수 있다.

 

주의할 점은 일반적인 에라토스테네스의 체와 방문하는 방식이 다르기 때문에, $k$가 합성수일 때도 방문을 해주어야 한다.

vector<int> sieve_of_eratosthenes(int n) {
    vector<bool> sieve(n + 1);
    vector<int> prime;
    for (int k = 2; k <= n; k++) {
        if (!sieve[k]) {
            prime.push_back(k);
        }
        for (int i = 0; prime[i] * k <= n; i++) {
            int p = prime[i];
            sieve[p * k] = true;
            if (k % p == 0) {
                break;
            }
        }
    }
    return prime;
}

냉정히 말해서 이렇게 한다 해도 동작 시간을 크게 줄여주지는 못한다. 기존 체에서 평균 방문 횟수인 $\log\log n$이 너무나도 작은 값이기 때문이다. 하지만 이 방식이 기존 방식과 차별화되는 가장 큰 특징은 각 정수 $n$을 정확히 한 번 방문하고, 방문할 때 $n$의 최소 소인수를 즉시 알 수 있으며 $n$을 방문하기 전 $n$의 약수를 모두 방문한다는 것이다. (마지막 사실은 증명하지는 않았으나 귀납법으로 쉽게 확인 가능하므로 생략한다.) 이 특징 때문에 특정 정수 이하의 모든 양의 정수에 대해 곱셈 함수의 값을 구하는 데 유용하다.

 

어떤 산술 함수 $f:\mathbb{Z}^+\to\mathbb{C}$가 $f(1)=1$이고 서로소인 두 양의 정수 $a$와 $b$에 대해 $f(ab)=f(a)f(b)$를 만족하면 곱셈 함수라고 한다. 곱셈 함수의 특징은 모든 소수 $p$와 양의 정수 $\ell$에 대해 $f(p^\ell)$의 값이 정해지면 함수가 유일하게 결정된다는 것이다.

 

이를 이용하면 dynamic programming으로 곱셈 함수를 계산할 수 있다. 어떤 양의 정수 $n$의 최소 소인수 $p$와 $p$로 나누어 떨어지지 않는 $m$에 대해 $n=p^\ell m$으로 나타날 때, $f(n)=f(p^\ell)f(m)$이기 때문이다. 우리는 $p$의 값만 알고 $p^\ell$의 값은 모르지 않느냐고 생각할 수 있지만, 이 값 역시 (곱셈 함수는 아니지만) dynamic programming으로 계산할 수 있다.

 

대표적인 곱셈 함수인 Euler's totient function의 값을 구하는 코드를 마지막으로 글을 마친다. 위키피디아 링크에서 더 다양한 곱셈 함수에 대해 알아볼 수 있다.

vector<int> euler_totient(int n) {
    vector<bool> sieve(n + 1);
    vector<int> prime;
    vector<int> phi(n + 1);
    phi[1] = 1;
    for (int k = 2; k <= n; k++) {
        if (!sieve[k]) {
            prime.push_back(k);
            phi[k] = k - 1;
        }
        for (int i = 0; prime[i] * k <= n; i++) {
            int p = prime[i];
            sieve[p * k] = true;
            if (k % p == 0) {
                phi[p * k] = phi[k] * p;
                break;
            } else {
                phi[p * k] = phi[k] * (p - 1);
            }
        }
    }
    return phi;
}