[정수론] 소수와 관련된 알고리즘 (소수 판정, 에라토스테네스의 체, 오일러 피 함수)

 

소수, 에라토스테네스의 체, 표준분해, 오일러 피 함수

Image from Wikipedia (English)

소수와 합성수 (primes and composites)

소수: $1$과 자기 자신만을 약수로 가지는 수이다. 영어로 prime number라 하기 때문에 보통 $p$로 표기한다.

합성수: $1$보다 큰 자연수 중에서 소수가 아닌 수. 약수의 개수가 3개 이상이다.

※ $1$은 소수도 아니고 합성수도 아닌 수이다. $1$을 소수로 분류하면 소인수분해가 유일하지 않기 때문인 것으로 보인다.

※ 소수(prime number)의 발음은 '소쑤'이고, 소수(decimal notation)의 발음은 '소:수'이다. 

 

$p\le \sqrt{n}$인 모든 소수 $p$에 대하여 $p$의 배수가 되는 합성수들을 모두 제거한 나머지는 모두 소수이다.

따라서 $n=46$의 소수는 $p\le \sqrt{46}=6.78\cdots$이므로 $46$의 소수는 $2,3,5$뿐이다.

 

소수 판정법 1, $O(n)$

어떤 수 $n$이 소수인지 아닌지 판별하는 기본적인 방법은 $2$부터 $n-1$까지 모든 숫자가 $n$을 나누는지 확인하는 방법이다. 따라서 시간복잡도는 $O(n)$ 이다.

def is_prime(n):
    if n <= 1:
        return False
    for i in range(2, n):
        if n % i == 0:
            return False
    return True

 

소수 판정법 2, $O(\sqrt{n})$

어떤수 $i$가 $n$을 나누면 $n//i$ 역시 $n$을 나눈다.

그러므로 $n-1$까지 확인하면 중복되는 판단을 한다.

그러므로 주어진 수의 제곱근인 $\sqrt{n}$까지만 확인하면 된다.

 이 경우의 시간복잡도는 $O(\sqrt{n})$ 이다.

참고로, int(n ** 0.5) + 1와 같이 쓸 수 있지만, ** 0.5 연산은 애초에 부동소수점 연산이다.

def is_prime(n):
    if n <= 1:
        return False
    if n == 2:
        return True
    if n % 2 == 0:
        return False

    i = 3
    while i * i <= n:
        if n % i == 0:
            return False
        i += 2
    return True

$n$이 $2$가 아니라면 짝수를 먼저 걸러낸다.

따라서 while 까지 내려온 수는 짝수가 아니므로 $i$를 $2$씩 증가시켜 찾을 수 있다.

(그렇다고해서 전체 시간복잡도가 줄어들지는 않는다.)

알고리즘대회가 아닌 기업 코딩테스트 수준에서 $O(\sqrt{n})$ 방법을 이용하면 대부분 소수 판정 시간초과를 받지 않는다.

※ 매우 큰 수의 경우 밀러-라빈(소수판정) 또는 폴라드 로(소인수분해) 등을 이용한다. 거의 대회에서나 사용되는 고급 알고리즘이다. 

 

에라토스테네스의 체, $O(n\log \log n)$

1은 소수가 아니므로 제외한다.

2부터 소수의 배수를 제외하는 방법이다.

2는 소수이므로 2보다 큰 2의 배수를 제외한다.

3은 소수이므로 3보다 큰 3의 배수를 제외한다.

4는 이미 지워졌다.

5는 소수이므로 5보다 큰 5의 배수를 제외한다.

이렇게 반복하면 상당히 빠른 시간내에 합성수를 제외할 수 있다. 

is_prime = [True] * (N + 1)
is_prime[0] = is_prime[1] = False

i = 2
while i * i <= N:
    if is_prime[i]:
        for j in range(i * i, N + 1, i):
            is_prime[j] = False
    i += 1

# if is_prime[x] == True, then x is prime (1 <= x <= N)

 

특정 구간의 수가 소수인지 판정하는 방법에 효율적이다.

주로 소수를 판정해야 하는 수가 쿼리로 계속 주어지는 경우, $O(q\sqrt{n})$의 시간이 소요되므로 미리 소수인지 아닌지를 배열에 저장했다가, 필요할때 해당 인덱스로 접근하여 해결한다.

쿼리 개수 $q$가 $\sqrt{n}$보다 많이 등장하면 에라토스테네스의 체를 이용하는 것이 좋다.

 

BOJ 1929, 소수 구하기, 실버 3, https://www.acmicpc.net/problem/1929

$M$ 이상 $N$ 이하의 소수를 모두 출력하는 문제이다. 

$M$과 $N$의 하한, 상한이 주어지므로 에라토스테네스의 체를 이용하는 것이 좋겠다.

(물론 이 문제는 위의 $O(\sqrt{n})$으로 해도 AC받을 수 있다.)

 

BOJ 4134, 다음 소수, 실버 4, https://www.acmicpc.net/problem/4134

주어진 $n(\le 4\cdot {10}^9$)과 크거나 같은 수 중 가장 작은 소수를 출력하는 문제이다.

$n$의 범위가 매우 크기 때문에 에라토스테네스의 체를 이용하면 메모리 이슈와 시간초과를 받을 있을 수 있다.

게다가 $n$ 이상의 범위에서 소수를 찾아야 하기 때문에 배열 크기의 상한을 쉽게 어림하기 어렵다.

따라서 처음에 소개한 $O(\sqrt{n})$의 알고리즘을 이용하여 첫 소수를 찾아야 한다.

참고로 소수의 개수는 무한함이 증명되어있기 때문에 무한루프에 걸리지 않고 반드시 멈춘다.

import sys
#sys.stdin = open('input.txt', 'r')
input = sys.stdin.readline

def is_prime(n):
    if n <= 1:
        return False
    if n == 2:
        return True
    if n % 2 == 0:
        return False

    i = 3
    while i * i <= n:
        if n % i == 0:
            return False
        i += 2
    return True


T = int(input())
for _ in range(T):
    n = int(input())
    i = n
    while not is_prime(i):
        i += 1
    print(i)

산술의 기본정리 (The Fundamental Theorem of Arithmetic; FTA)

(1) $1$보다 큰 소수는 소수이거나, 소수들의 곱으로 표현할 수 있다.
(2) 이 표현법은 (곱하는 순서를 정해주면) 유일하다. 

Corollary
$1$보다 큰 정수 $n$은 아래와 같이 표준형(canonicㅁa form)으로 나타낼 수 있고 이는 유일(unique)하다.
$$n=p_1^{k_1}{p}_2^{k_2}\cdots p_r^{k_r}$$
이때 $k_i$은 양의 정수이며 각 $p_i$의 대소관계는 $p_1<p_2<\cdots <p_r$ 이다.

 

증명은 생략한다.

유클리드 호제법과 유사하게 $n>n_1>n_2>\cdots \ge 1$에서 반드시 멈추게 된다.

 

$n$의 표준분해가 다음과 같다고 하자.
$$n=p_1^{k_1}{p}_2^{k_2}\cdots p_r^{k_r}$$

[양의 약수의 개수]
$n$의 모든 양의 약수들의 개수를 $\tau (n)$이라 하면
$$\tau (n)=(k_1+1)(k_2+1)\cdots (k_r+1)$$
이다.

[양의 약수의 합]
$n$의 모든 양의 약수들의 합은 $\sigma (n)$라 하면
$$\sigma (n)=\sum _{d\mid n}d=\prod _{i=1}^r\frac{p_i^{k_i+1}-1}{p_i-1}$$
이다.

 

소인수분해 알고리즘

관련 문제

BOJ 11653, 소인수분해, 브론즈 1, https://www.acmicpc.net/problem/11653

소수판정법과 유사하게 인수 $p$가 $N$을 나누면 $N//=p$를 한다.

$p$가 여러번 나눌 수 있으므로 while 조건에 N % p == 0을 통해 반복 출력한다.

가장 큰 소수는 $\sqrt{n}$ 이내에 있으므로 $\sqrt{n}$까지만 탐색하면 된다.

시간복잡도는 $O(\sqrt{n})$이다.

N = int(input())
i = 2
while i * i <= N:
    while N % i == 0:
        print(i)
        N //= i
    i += 1
if N > 1:
    print(N)

 

BOJ 2312, 수 복원하기, 실버 3, https://www.acmicpc.net/problem/2312

import sys
#sys.stdin = open('input.txt', 'r')
#sys.setrecursionlimit(10 ** 7)
input = sys.stdin.readline

T = int(input())
for _ in range(T):
    N = int(input())
    p = 2
    while p * p <= N:
        p_count = 0
        while N % p == 0:
            p_count += 1
            N //= p
        if p_count:
            print(p, p_count)
        p += 1
    
    # N이 1보다 큰 경우, 남은 N은 소수임
    if N > 1:
        print(N, 1)

 

BOJ 16563, 어려운 소인수분해, 골드 4, https://www.acmicpc.net/problem/16563

이번엔 소인수분해를 할 수가 쿼리로 주어진다.

쿼리 개수가 ${10}^6$이므로 $O(q\sqrt{n})$을 계산해보면 시간초과를 얻는다.

 

에라토스테네스의 체를 응용하여 1차원 배열 F[x]를 x의 가장 작은 인수라 하자.

예를 들어 F[10]=2, F[9]=3, F[49]=7, F[p]=p이다, (p는 소수)

참고로 F[1] = 0이다.

 

예를 들어, 30을 인수분해한다고 하자.

그러면 F[30]=2이므로 2를 출력하고 30을 2로 나눈 15를 찾는다.

F[15]=3이므로 3을 출력하고 15를 3으로 나눈 5를 찾아간다.

F[5]=5이므로 5를 출력하고 5를 5로 나눈 1을 찾아간다.

마지막으로 남은 수 1은 while 조건에 위배되므로 while loop가 끝나게 된다.

import sys
#sys.stdin = open('input.txt', 'r')
#sys.setrecursionlimit(10 ** 7)
input = sys.stdin.readline

MAX_N = 5000000
F = [i for i in range(MAX_N + 1)]
F[1] = 0
i = 2
while i * i <= MAX_N:
    if F[i] == i:
        for j in range(i * i, MAX_N + 1, i):
            if F[j] == j:
                F[j] = i
    i += 1

def factorize(F, x):
    while x > 1:
        print(F[x], end=' ')
        x //= F[x]
    print()

N = int(input())
k_list = list(map(int, input().split()))
for k in k_list:
    factorize(F, k)

오일러 피 함수 (Euler's phi function)

정수론에는 자주 등장하는 함수들이 있는데, 오일러 피 함수는 그중 하나이다.

Euler phi-function
$n\in \mathbb{N}$일 때, $\varphi (n)$은 $1\le a\le n$ 중에서 $gcd(a,n)=1$를 만족하는 $a$의 수로 정의한다.
즉, $n$이하의 양의 정수 중에서 $n$과 서로소인 정수의 개수를 $\varphi (n)$이라 한다.

Properties
(1) $\varphi (n)=n-1⟺n$ is prime.
(2)  $n\ge 3$에 대하여 $\varphi (n)$은 짝수이다.
(3) 소수 $p$에 대하여 $\varphi (p^k)=p^k-p^{k-1}=p^k(1-\frac{1}{p})$
(4) $\varphi (mn)=\varphi (m)\varphi (n)$ 이다. 즉 $\varphi$는 곱셈함수(multiplicative function)이다.
(5) $n=p_1^{k_1}{p}_2^{k_2}\cdots p_r^{k_r}$이라 하면 
$$\varphi (n)=\prod _{i=1}^r\varphi (p_i^{k_i})=n(1-\frac{1}{p_1})(1-\frac{1}{p_2})\cdots (1-\frac{1}{p_r})$$

 

(1) $\varphi (13)=12$

(4) $\varphi (4\cdot 9)=\varphi (4)\varphi (9)=2\cdot 6=12$

(5) $\varphi (20)=\varphi (2^2\cdot 5)=2\cdot 4=8$

 

BOJ 11689, GCD(n, k)=1, 골드 1, https://www.acmicpc.net/problem/11689

오일러 피 함수를 직접 물어보는 문제이다.

에라토스테네스의 체를 응용한 반복문이므로 시간복잡도는 $O(\sqrt{n})$이다.

$p$가 $n$을 나누면, 더이상 $p$가 인수가 되지 않도록 계속 나눈다. 

그리고 오일러 함수의 정의에 따라 $(1-\frac{1}{p})$를 곱한다.

그런데 분수를 이용한 곱셈은 부동소수점을 이용하므로 적절하지 않다.

$n(1-\frac{1}{p})=n-\frac{n}{p}$이고, 위에서 $p$는 $n$의 인수임이 확실하므로 몫연산과 뺄셈을 이용하여 정확한 정수연산으로 치환할 수 있다. 

if n % p == 0:
    while n % p == 0:
        n //= p
    result -= result // p

 

전체 코드는 다음과 같다.

def euler_phi(n):
    result = n
    p = 2
    while p * p <= n:
        if n % p == 0:
            while n % p == 0:
                n //= p
            result -= result // p
        p += 1
    if n > 1:
        result -= result // n
    return result

N = int(input())
print(euler_phi(N))

 

수학적으로는 $gcd(1,1)=1$이므로 $\varphi (1)=1$이다.

그러나 프로그래밍 문제 상황에 따라 $1$의 경우 다르게 처리해야할 수 있다.

 

BOJ 23832, 서로소 그래프, 골드 1, https://www.acmicpc.net/problem/23832

서로 다른 두 정점이 서로소인 경우에만 간선을 연결한다.

그러니까 $gcd(i,j)=1$인 순서쌍을 세면된다. 

Brute Force로 모든 정점 쌍을 찾으면 $\frac{N(N-1)}{2}$이므로 시간초과이다.

$i>1$인 경우에는 $\varphi (i)$를 계산하면 되지만, $\varphi (1)=gcd(1,1)=1$을 계산하므로 $\varphi (1)$은 계산하지 않도록 주의한다.

(예시 코드 생략)

 

BOJ 4355, 서로소, 골드 1, https://www.acmicpc.net/problem/4355

단순히 서로소를 판단하고 싶으면 $gcd(a,b)=1$를 계산하면 된다.

그러나 $n={10}^9$이고 테스트케이스마다 쿼리가 들어오므로 단순 계산은 힘들겠다.

$n$보다 작은 양의 정수 중에서 $n$과 서로소인 수의 개수이다.

$\varphi (n)$은 $n$ 이하의 양의 정수 중에서 $n$과 서로소인 수의 개수이다.

따라서 이 문제에서 $n=1$인 경우에만 $0$을 출력하도록 한다.

(예시 코드 생략)