선형 체를 통한 곱셈적 함수의 계산

곱셈적 함수

정수를 정의역으로 갖는 함수 $f$가 서로소인 두 정수 $a, b$에 대해 $f(a)f(b) = f(ab)$를 만족할 때, 함수 $f$를 곱셈적이라고 합니다. 대표적인 곱셈적 함수로는 다음과 같은 함수들이 있습니다.

• $\mu(n)$: 뫼비우스 함수
• $\tau(n)$: 약수 개수 함수
• $\phi(n)$: 오일러 피 함수
• $\sigma(n)$: 약수 합 함수

$1$은 자기 자신과 서로소기에, $f(1) = f(1)f(1)$이 성립해야 합니다. $f(1) = 0$이라면 모든 자연수 $a$에 대해 $f(a) = f(1)f(a) = 0$이 되기에, 보통 고려하지 않습니다. 따라서 우리는 앞으로 항상 곱셈적 함수 $f$에 대해 $f(1) = 1$이라고 생각합니다. 곱셈적 함수의 곱셈적 성질에 의해, 소수 p에 대해 $f$의 함숫값 $f(p^k)$만을 알면 다른 모든 정수에 대해 $f$의 함숫값을 구할 수 있습니다.

그 전에 $f(p^k)$는 $O(1)$에 준해 계산할 수 있다고 가정하겠습니다. 예를 들어, $\mu(p^{k}) = -[k = 1]$, $\tau(p^k) = k+1$, $\phi(p^k) = p^k - p^{k-1}$, $\sigma(p) = \frac{p^k-1}{p-1}$ 등 $p$와 $k$에 대한 간단한 식으로 나타내어집니다.

위의 내용을 바탕으로 우리는 특정한 함숫값 $n$에 대해 $f(n)$을 구하기 위해서 먼저 소인수분해를 진행한 후, 각각의 소인수에 대해 함숫값을 구해 곱하는 방식을 사용할 수 있습니다. 하지만, 만약 1부터 N까지 모든 $n$에 대해 $f(n)$의 값이 궁금하다면 어떻게 할까요? 위의 방법을 모든 수에 대해 사용하면, 폴라드 로 알고리즘을 통해 ${\cal O}(N^{1/4})$에 소인수분해를 할 수 있다 하더라도 ${\cal O}(N^{5/4})$만큼의 연산이 필요합니다. 이때 선형 체를 통해 이를 ${\cal O}(N)$에 계산할 수 있습니다.

선형 체의 동작 과정

먼저 각 함수의 값을 저장할 테이블이 필요하겠죠. 기존 코드에서의 소수 여부를 저장하는 테이블 sieve와 크기가 같은 테이블 f를 하나 더 선언해줍시다. 소수에 대해서는 이미 ${\cal O}(1)$에 구할 수 있다고 가정했기 때문에, 합성수 $n$에 대해서만 고려하겠습니다.

선형 체에서는 각 $n$을 딱 한 번씩만 체에서 지웁니다. $n$은 $q$를 2부터 증가시키는 과정에서, $n = p \cdot q$인 $q$를 순회할 때 지워집니다. 여기서 $p = {\rm lpf}(n)$이고요. 이전 글에서 말했듯 항상 $p \leq {\rm lpf}(q)$는 성립합니다. 우리는 각 $n$을 지우는 과정에서 $f[n]$을 채울 것입니다.

먼저 $p < {\rm lpf}(q)$라면 항상 $\gcd(p, q) = 1$이 됩니다. ${\rm lpf}(q)$는 $q$를 나누는 가장 작은 소수이고, $p$는 그 소수보다 더 작은 소수이기 때문이죠. $q$가 소수라면 $f[q]$는 ${\cal O}(1)$에 구할 수 있고, 합성수라면 $q/{\rm lpf}(q)$를 순회할 때 $f[q]$의 값을 이미 구해놓았기에 $f[p]$와 $f[q]$의 값은 이미 알고 있다고 할 수 있습니다. 그러므로, 이 경우에는 $f[n] = f(p)f[q]$로 간단히 $f[n]$을 채울 수 있습니다.

조금 더 복잡한 경우는 $p = {\rm lpf}(q)$입니다. 이를 구하기 위해 $p^k \mid q$이고 $p^{k+1} \nmid q$라고 가정합시다. 즉, $k$는 $p$가 $q$를 나누는 최대 횟수입니다. 그러면 $\gcd{p, q/p^k} = 1$이 성립하기에 $f(n) = f(pq) = f(p^{k+1} \dfrac{q}{p^k}) = f(p^{k+1}) f(\dfrac{q}{p^k})$라는 식을 세울 수 있습니다. 또한 $f(q) = f(p^k \dfrac{q}{p^k})$이므로 $f(q) = f(p^k) f(\dfrac{q}{p^k})$인 것도 알 수 있습니다. 위 두 식을 종합하면 $f(n) = f(pq) = f(p^{k+1}) f(\dfrac{q}{p^k}) = \dfrac{f(p^{k+1})}{f(p^k)} f(q)$가 됩니다.

일부 함수의 경우 $\frac{f(p^{k+1})}{f(p^k)}$의 값이 매우 간단히 나타납니다. 예를 들어, 오일러 피 함수의 경우는 $\frac{\phi(p^{k+1})}{\phi(p^k)} = \frac{p^{k+1}-p^k}{p^k-p^{k-1}} = p$입니다. 예외적으로 뫼비우스 함수의 경우는 $k \geq 2$일 때 $\frac{\mu(p^{k+1})}{\mu(p^k)}$가 $0/0$ 꼴로 나타내어지기 때문에, 앞선 식을 정리하여 보면 $k \geq 1$일 때 $\mu(p^{k+1}) = 0$이라 $f(n) = f(p^{k+1}) f(\frac{q}{p^k}) = 0$이 됩니다. 이러한 경우에는 매우 간단히 식을 정리할 수 있습니다. 위에서 예시를 든 오일러 피 함수의 경우는 $f[n] = p \cdot f[q]$이고, 뫼비우스 함수의 경우는 $f[n] = 0$으로 정리할 수 있습니다.

하지만 일반적으로는 $\frac{f(p^{k+1})}{f(p^k)}$가 $k$에도 의존하는 식으로 나타내어집니다. 이를 구하기 위해 우리는 $k$의 값을 저장하는 테이블을 하나 더 만들 수 있습니다. 즉, ${\rm lpf}(n)$이 $n$을 나누는 횟수를 저장하는 $k[n]$ 배열을 도입할 수 있습니다. 여전히 소수 $p$에 대해 $k[p] = 1$이므로 합성수에 대해서만 생각하겠습니다.

다시 위의 $p < {\rm lpf}(q)$의 경우로 돌아가서, $k[n]$의 값을 어떻게 채울지 생각해봅시다. 그런데, $n = pq$이고 ${\rm lpf}[n] = p \neq {\rm lpf}(q)$이므로 $k[n]$은 항상 1이 됩니다. 또 $p = {\rm lpf}(q)$라면 $n$은 단순히 $q$에 ${\rm lpf}(q)$가 한 번 더 곱해진 것이기 때문에 $k[n] = k[q]+1$이 성립합니다. 간단히 $k$의 값도 구할 수 있게 되었으므로, 우리는 $\frac{f(p^{k+1})}{f(p^k)}$도 구할 수 있고, 결과적으로 $f[n] = \frac{f(p^{k+1})}{f(p^k)} \cdot f[q]$로 테이블을 채울 수 있습니다.

정리하자면, 우리는 $f$와 $k$, 두 개의 테이블을 채워야 합니다. 각 테이블은 $n$이 체에서 지워질 때 값이 채워집니다. $q$를 2부터 증가시키며 순회하는 과정에서, ${\rm lpf}(n) = p$이고 $n = pq$일 때,

• 만약 $p < {\rm lpf}(q)$라면, $k[n] = 1$, $f[n] = f(p)f[q]$로 채웁니다.
• 만약 $p = {\rm lpf}(q)$라면, $k[n] = k[q]+1$, $f[n] = \frac{f(p^{k+1})}{f(p^k)} \cdot f[q]$로 채웁니다.

추가적으로, $f(p)$를 매번 계산하는 대신, $p$를 순회할 때 $f[p] = f(p)$로 채워주어 불필요한 연산도 줄이고, 겸사겸사 $f$를 소수와 합성수를 포함해서 모든 수에 대해 $f$의 함숫값을 가지는 테이블로 사용할 수 있습니다. 이를 구현하면 다음과 같습니다.

구현

# 실제 구현에서는 함수로 빼는 대신 직접 식을 대입해 사용하는 게 더 간결합니다.
def get_f(p, k): return k+1 

N = 1000

sieve = [0] * N
f = [0] * N
k = [1] * N
primes = []

f[1] = 1 # f(1) = 1

for i in range(2, N):
    if sieve[i] == 0:
        primes.append(i)
        f[i] = get_f(i, 1) # 소수에 대해 f 채우기기
    
    for p in primes:
        if i*p >= N:
            break
        
        sieve[i*p] = 1
        if (i % p == 0): # lpf(i) = p
            f[i*p] = f[i] * get_f(p, k[i]+1) // get_f(p, k[i])
            k[i*p] = k[i] + 1
            break
        
        # lpf(i) < p
        f[i*p] = f[i]*f[p]
        k[i*p] = 1

print(f)

예시로 약수 개수 함수를 사용했습니다. 오일러 피 함수와 같이 $f(p^{k+1})/f(p^k)$의 값이 k에 의존하지 않는 경우에는 k 배열을 빼고 구현해도 무방합니다.

마치며

에라토스테네스의 체를 활용해 ${\cal O}(n \log \log n)$에 테이블을 채우는 방법도 존재합니다. 나중에 다른 코드를 읽거나 할 때 참고하시면 좋겠습니다.

또한, 저 코드에서 간단한 변형을 통해 f 대신 f의 누적 합을 구하도록 만들 수도 있습니다. 다만 f를 구한 후 다시 누적 합 배열을 만들어도 시간 복잡도는 같기에 굳이 자세한 내용을 적지는 않겠습니다. 읽어주셔서 감사합니다!