Haskell - 소수 구하기

최근에 스칼라의 유명한 프레임워크인 Play Framework를 관리하던 Lightbend의 직원들의 퇴사 소식이 들려옴에 따라, 함수형 언어를 좋아하는 저로써는 장차 스칼라의 미래가 걱정되기 시작했습니다. 그래서 얼마전부터 하스켈 공부를 시작하게 되었는데요, 그 과정에서 최근에 배운 것 중, 인상 깊었던 것을 적어보고자 합니다.

하스켈로 소수 구하기

하스켈 책(Programming in Haskell)을 뒤적뒤적 하다보면, 리스트 챕터에서 소수를 구하는 코드를 작성하는 내용이 나옵니다. 다른 하스켈 책들도 리스트와 high-order function을 논할 때 항상 이 소수가 나옵니다.

이 포스팅에서는 하스켈을 이용하여 소수를 구하는 코드를 작성하고, 개선하는 작업을 해보겠습니다.

타 언어 살펴보기

우선, 소수를 구하는 코드를 다른 언어들로부터 살펴보겠습니다. 단순하게, n보다 같거나 작은 수의 소수를 구하는 코드라 하면, 1 부터 n까지 반복하면서 각 수의 약수들을 구하고, 약수가 1과 자기자신밖에 없다면 소수로 인식하는 코드입니다.

vector<int> primes(int until) {
    vector<int> ret;
    for (int i = 2; i <= until; i++) {
        bool isPrime = true;
        double sq = sqrt(i);
        for (int j = 2; j <= sq; j++) { // 1과 자기 자신은 생략합니다.
            if (i % j == 0) {
                isPrime = false;
                break;
            }
        }
        if (isPrime) ret.push_back(i);
    }
    return ret;
}

한편, 함수형 언어라 불리우는 스칼라는 아래와 같이 작성할 수 있습니다.

def primes(n: Int) = 2 +: (3 to n).filter (x => !(2 until Math.sqrt(x).toInt).exists(x % _ == 0))

스칼라도 함수형 언어이므로 표현을 간결하게 둘 수 있습니다.

이제 찐 함수형 언어는 어떤지 살펴봅시다.

하스켈

하스켈은 Pure functional language라고 부르죠.. 상태가 전혀 없는 함수형 언어입니다. 특히, 수학의 집합을 표시하는 기호와 많이 닮은 점이 정말로 신기했습니다.

간단하게 예를 들어, 짝수를 구하는 식으로 여기에서 쓰일 문법들을 살펴보겠습니다.

evens :: Int -> [Int]
evens n = [x | x <- [1..n], x `mod` 2 == 0]

• [1..n]: 하스켈에서 Ranged 배열을 기술하는 방식으로, [1..5]이라 하면 [1,2,3,4,5]가 나오고, [1,3..10]이라 하면, [1,3,5,7,9]가 나옵니다. 한편, 위 evens는 [2,4..n]이라 표현하여도 무방합니다.

• x <- [1..n]: 이 표현은 1부터 n까지의 수가 있을 때, 하나씩 꺼내보는 Iterator입니다. 이제 저 x에 1부터 n까지의 수가 binding될 것입니다.

• x <- [1..n], x `mod` 2 == 0: 이 표현은 [1..n] 리스트에서 각 x마다, 2로 나눈 나머지가 참인 값들을 전달합니다.

• [x | ...] : 이 표현은 가드 ‘|’ 의 오른쪽 조건을 충족하는 모든 x라고 이해할 수 있습니다. 한편, 여기에서도 변형을 줄 수 있습니다.

  [x * 2 | x <- [1..n]]

  그리고 이 식은 Functor의 map 함수와 성질이 유사합니다.

  map (*2) [1..n]

이상으로 한가지 문제에 대한 여러 표현법을 살펴봤습니다. 물론 하스켈의 단점도 있습니다. 순수 함수형이라는 이유로 선언 후 변경이 가능한 변수(Mutable variable)의 사용이 일절 없다보니 어렵지만, 일단 만들어지면 매우 Robust한 특징이 있습니다. 여기에서 가변성에 대해 다루지는 않겠습니다.

이제 본격적으로 하스켈로 소수를 구해봅시다.

Prime Numbers with Haskell

소수를 구하는 방법은 여러가지가 있는데, 우선 가장 직관적인 것에서 시작해보겠습니다. 우선 아래의 명제로 시작해봅시다.

소수는 1과 자기 자신으로만 나누어 떨어지는 수이다.

이 명제를 기준으로, 이 포스팅에서는 n이 주어졌을 때, n이하의 모든 소수를 구하는 코드를 구현해보겠습니다.

Initial

우선, 문장을 분석해봅시다.

함수명	설명
primes n	n이 주어졌을 때, n이하의 isPrime을 충족하는 수들의 집합
isPrime n	divisors n이 1과 자기 자신으로만 구성되어 있을 때 이를 소수라 부른다.
divisors n	n이하의 자연수 중, 이를 나눴을 때 나머지가 0인 수를 약수라 부른다.

하스켈은 먼저, primes n을 코드로 표현해봅시다. 자연수는 1 ~ n 의 정수입니다. 이는 하스켈로 [1..n]로 표현합니다.

한편 하스켈에서, 리스트를 정의할 때 수학하고 비슷한 표기법으로 정의한다 했는데, 그 원리는 아래와 같습니다.

[(실제 원소로 남을 값) | (원소 들의 집합), (원소 값에 대한 규칙)]

primes n = [x | x <- [1..n], isPrime x]

그리고 소수는 앞서 말한대로, 1과 자기자신만 존재하는지 체크하면 됩니다.

isPrime n = [1, n] == divisors n

참고로 순서가 중요합니다. 만약 divisors n이 [n, 1]을 돌려준다면 소수임에도 불구하고 False가 리턴될 것입니다. 그렇지 않도록 출력을 제한하면 됩니다.

한편, 어떤 수가 주어졌을 때, 이 수보다 작은 수들 중, 약수를 구하려면 어떤 수들의 집합에서, 나누어 떨어지는지 체크하는 조건문이 포함되어 있어야 하는데, 즉, 나머지를 체크하는 식이 포함되어야 합니다. 일반적인 프로그래밍 언어에서는 a % b를 많이 쓰지만, 하스켈에서는 a `mod` b 혹은, mod a b 이렇게 사용합니다.

divisors n = [x | x <- [1..n], n `mod` x == 0]

그렇지만, 1과 n의 비교는 여기에서는 불필요합니다.

divisors n = [x | x <- [2..n - 1], n `mod` x == 0]

이 코드들을 조합하면, 네.. 3줄만에 소수를 구할 수 있습니다.

이제 위 코드를 정리하면 아래와 같습니다.

divisors x = [x' | x' <- [2..x - 1], x `mod` x' = 0]
isPrime x = null $ divisors x
primes x = [x | x <- [2..x], prime x]

조금 반칙 같으니, 함수의 타입들도 붙여줘봅시다.

divisors :: Int -> [Int]
divisors x = [x' | x' <- [2..x-1], x `mod` x' = 0]

isPrime :: Int -> Bool
isPrime x = null $ divisors x

primes :: Int -> [Int]
primes x = [x | x <- [2..x], prime x]

사실 하스켈을 처음 배울 때 정말로 이렇게 짧아질 수 있는지 싶어 놀라웠습니다.

이제 실제 세계에서 쓸 수 있으려면, 제약조건이 하나 필요합니다. 시간입니다. 하스켈이 우선 native한 환경에서 동작하기 때문에, 당연히 왠만한 스크립트 코드나 JVM계열의 언어보다는 빠르겠습니다만, primes x의 시간 복잡도는 n까지의 수를 계산하라고 할 때, 우선 O(n^2)입니다.(divisors가 O(n), isPrime은 divisors로 인해 O(n), primes는 O(n * O(isPrimes)) = O(n^2)) 한마디로 최적화가 필요하죠..

Improvement

두가지 방법을 살펴보겠습니다. 약수를 개선하는 방법과, 에라토스테네스의 채를 이용하는 방법입니다.

약수의 개선

위의 코드에서 간단하게 개선해볼 수 있는 방법으로는, 반복자를 도는 횟수를 줄인다는 것에 착안합니다. 결국, 1 ~ x를 모두 비교하지 않고, 어떤 중간의 수에서 멈추는 방법입니다. 간단하게, 8이라는 수를 생각해봅시다. 정수가 아닌 실수들을 모두 고려하여 연산을 해봅시다.

xs = [(x, 8 / x) | x <- [1..8]] 

1 -> 8.00
2 -> 4.00  -- 
3 -> 2.67  --
4 -> 2.00
5 -> 1.60
6 -> 1.33
7 -> 1.14
8 -> 1.00

어느 지점에서 나누는 수와 몫의 크기가 반전되는지를 찾아야 하는데, 일단 반절은 아닌 것 같습니다. 9로 키워봅시다.

1 -> 9.00
2 -> 4.50
3 -> 3.00 --
4 -> 2.25
5 -> 1.80
6 -> 1.50
7 -> 1.29
8 -> 1.12
9 -> 1.00

숫자를 점점 키워가다보면, 제곱근 수와 일치하는 지점에서 수의 크기가 반전됨을 알 수 있습니다. 따라서, sqrt x만큼만 돌아도 나머지를 구할 수 있다는 결론이 됩니다. (저렇게 반전이 되는 지점을 찾는 방법을 이용해서 제곱근을 구할 수 있습니다. 뉴튼의 근사법이 그런 방법을 이용합니다.)

-- 제곱근 까지만 약수를 찾는 코드
divisorsHalf x = [x' | x' <- [2..n], x `mod` x' == 0]
    where n = round $ sqrt $ fromIntegral x / 1
-- divisorsHalf 9 = [1, 3]

-- 주어진 집합을 토대로 나머지 약수를 찾는 코드
divisorsRest x xs = [x `div` x' | x' <- xs, x'^2 /= x]

-- divisorsRest 9 $ divisorsHalf 9 = [9]

두 결과물을 합쳐서 [1, 3, 9]를 만들어야 합니다.

이제 이것을 가지고 divisors’ 함수를 다시 구성해봅시다.

divisors' x = l ++ r
    where l = divisorsHalf x
          r = divisorsRest x l

그리고 아까 구해둔 prime x와 primes x를 이용하면 처음 버전의 primes x보다는 빠른 성능으로 구할 수 있습니다.

에라토스테네스의 채

한편, 위의 방법들과 다르게 Memoization 기법을 이용하여 구할 수도 있습니다. 동적 프로그래밍… 만일, 이미 구해둔 소수로 나누어지는 값들은 그냥 살피지 않고 지나칠 수 있다면? 이 방법이 사실 에라토스테네스의 체입니다. 아시다시피, 미리 숫자들을 나열해놓고, 1 지우고, 2가 아닌 2의 배수를 지우고, 3이 아닌 3의 배수를 지워나가면서 남는 수들만으로 소수를 구하는 방법입니다. 이는 사실 상태변수를 가질 수 있는 일반적인 방법으로도 충분히 구해낼 수 있습니다.

아래는 C++로 간단하게 짜본 코드입니다.

vector<bool> xs(n + 1, true);
vector<int> ret;

for (int i = 2; i < n; i++) {
    for (int j = i * 2; j < n; j += i) {
        xs[j] = false;
    }
}

for (int i = 2; i < xs.size(); i++) {
    if (xs[i]) ret.push_back(i);
}

이처럼 Mutable한 자료구조가 있다면 손쉽게 구할 수 있습니다만, 아시다시피.. 함수형 언어는 전혀 그렇지가 않습니다.

나름대로 이 방법 저 방법을 생각해본 결과, 반복적으로 filter를 적용하는 것이 좋을 것 같았습니다. 우선 전체 수가 나열된 집합을 구하는 함수와, 자기 자신을 제외한 배수들을 걸러내는 함수를 만들어봅시다.

totalNumbers n = [2..n]
-- totalNumbers 10 = [2,3,4,5,6,7,8,9,10]

sieve :: Int -> [Int] -> [Int]
sieve xs n = filter (\x -> x == n || x `mod` n == 0) xs
-- sieve 2 (totalNumbers 10) = [3,5,7,9]

다음 반복자를 정의해봅시다. 2부터 시작하여 n 까지, 배수의 배열들을 이용해서 계속 제외해나가는 방식입니다.

primes' :: Int -> [Int]
primes' n = foldl sieve (totalNumbers n) [2..n]

동작시켜봤는데… 정말로 느립니다. 그냥 primes가 훨씬 빠르네요… 왜그런지 문제를 고민해봤습니다… 아무래도 이미 찾은 소수에 대해서도 계속해서 비교 및 연산을 하다보니 이런 문제가 생긴 것 같습니다.

에라토스테네스 채 기법 개선하기

도대체 뭐가 문제일까 고민해본 결과, 연산에 불필요한 수행이 보였습니다. [2,3,5,7,...10007]이 있다고 가정해봅시다. 10006일 때, 2 ~ 10006까지는 더이상 살펴볼 필요가 없습니다. 그리고 10006이 지워졌다면, 이미 그 아래의 약수로 인해 지워졌을 것입니다.

그러므로, 그 불필요한 연산을 줄이는 방법을 생각해봅시다.

• 우선, 어떤 수 n을 검사할 때, 자기 자신을 포함한 그 이하의 수는 검사할 필요가 없습니다.
• 그리고, 어떤 수 n이 없다면, 이미 그보다 작은 약수에 의해 지워진 것입니다. 그리고 이 n의 배수 또한 없습니다. 예를들어, 6이 없다면, 그 다음 수인 12는 이미 2에 의해 지워져있는 셈이지요.

여기서는 재미있는 특징이 있었는데, 주어진 리스트의 첫번째 원소는 검사하는 기준이 됩니다. 그렇다면, head와 tail로 구분하는 패턴매칭을 적용하면 됩니다. 그 다음, 남은 꼬리 부분에만 filter를 적용합시다. 한편, 꼬리 부분에는 head 원소가 없으므로, x == x' 라는 조건문이 없어도 됩니다.

primesLoop (h:xs) = h : primesLoop (filter p xs)
    where p x = x `mod` h /= 0

만약 더이상 검색할 수가 남아있지 않다면 비어있는 리스트를 되돌려줍니다.

primesLoop [] _ = []

최초에 이 함수를 호출할 때는 아래와 같이 호출합니다.

primesLoop [2..n]

이상을 종합해보면, 아래와 같은 코드가 나옵니다.

primes' :: Int -> [Int]
primes' n = primeLoop [2..n]
 where
  primeLoop [] = []
  primeLoop (h:xs) = h : primeLoop (filter p xs)
    where p x = x `mod` h /= 0

문제는, 이렇게 해도 primes보다 성능이 안나옵니다. 그렇지만, 여기에서 n을 지운다면 무한 소수 스트림을 얻을 수 있습니다.

primes' :: [Int]
primes' = primeLoop [2..]
 where
  primeLoop [] = []
  primeLoop (h:xs) = h : primeLoop (filter p xs)
    where p x = x `mod` h /= 0

이상을 갖고, take나 drop을 활용하여 수를 얻어올 수 있습니다.

• 소수 1000개 추출

  take 1000 primes' -- 소수 1000개 추출
  

• 소수 1000번째부터 1000개 추출

  take 1000 $ drop 1000 $ primes'
  

• 1000 이하의 소수들 추출

  takeWhile (<1000) primes'
  

참고로, primes’ 에서, 이미 구해진 값들은 캐싱이 이루어져서 첫 호출 후, 두번째 호출부터는 속도가 개선됩니다.

> :set +s -- 타이머 켜기
> takeWhile (<10000) primes'
...,9907,9923,9929,9931,9941,9949,9967,9973]
(0.59 secs, 268,324,608 bytes)
> takeWhile (<10000) primes'
...,9907,9923,9929,9931,9941,9949,9967,9973]
(0.04 secs, 4,712,216 bytes)

허허.. 저도 이거 하다가 놀랐네요.ㅋㅋㅋ 이렇게 DP가 쉽게 구현되었습니다.

컴파일러가 DP를 해주네…??

한편, take n primes'를 호출한 후, take (n + 1) primes'를 호출하면, 실제 어떤 일이 일어날까요? n + 1번째 원소는 take n이 호출 될 시점엔 아무런 filter를 거치지 않고 있습니다. 하지만, take (n + 1)이 불리운다면, n번째 소수와 n + 1번째 소수 사이의 모든 수들이 앞에서 구해둔 소수들과의 나눗셈 작업을 시작하게 됩니다.

동작 과정을 지켜보기 위해 trace라는 함수를 사용해봅시다. trace는 String -> a -> a 의 함수 타입으로, a를 그대로 되돌려주되, String을 화면에 출력하는 함수입니다. 예를들어, trace "안녕!" 3은 3을 그대로 되돌려주지만, 화면에 “안녕”을 출력합니다.

import Debug.Trace

primes' :: [Int]
primes' = primeLoop [2..]
 where
  primeLoop [] = []
  primeLoop (h:xs) = h : primeLoop (filter p xs)
    where p x = trace ("{h = " ++ (show h) ++ ", x = " ++ (show x) ++ "}") (x `mod` h /= 0)

이렇게 바꾸면 take 5 primes' 라는 연산에 아래와 같은 결과물을 출력합니다.

> take 5 primes'

[2{h = 2, x = 3}
,3{h = 2, x = 4} -- 4는 2로 나누어져서 제거됩니다.
{h = 2, x = 5}
{h = 3, x = 5}   -- 5는 소수
,5{h = 2, x = 6} -- 6은 2에 제거
{h = 2, x = 7}
{h = 3, x = 7}
{h = 5, x = 7}   -- 7은 소수
,7{h = 2, x = 8} -- 8은 2에 의해 제거..
{h = 2, x = 9}   
{h = 3, x = 9}   -- 9는 3에 의해 제거..
{h = 2, x = 10}  -- 10은 2에 의해 제거..
{h = 2, x = 11}
{h = 3, x = 11}
{h = 5, x = 11}
{h = 7, x = 11}  -- 11은 소수
,11]

그리고 다시 호출할 때는, trace 메시지는 보이지 않고, 그대로 결과를 돌려줍니다.

> take 5 primes'
[2,3,5,7,11]

그렇습니다! 이 값은 이미 연산이 되어있기 때문에 두번 하지 않는 것입니다! JVM도 이런식으로 최적화를 한다는 것을 본적은 있습니다. 함수가 항상 같은 값을 주면 JVM단에서 캐싱을 해둡니다. 예를들어, String의 해쉬코드가 실제 구현은 일일히 문자열을 돌아보는거지만, 막상 디버거를 찍어보면 두번째 방문시 값을 빠르게 돌려줍니다. 마찬가지로, 하스켈에서도 이정도 최적화는 지원해주나 봅니다…

이 시점에서 take 6 primes' 를 호출해봅시다.

> take 6 primes'
[2,3,5,7,11{h = 2, x = 12} -- 12는 2에 의해 제거
{h = 2, x = 13}
{h = 3, x = 13}
{h = 5, x = 13}
{h = 7, x = 13}
{h = 11, x = 13}           -- 13은 소수
,13]

앞에서 했던 계산은 두번 하지 않고, 6번째 수를 위해서만 연산을 수행합니다. 그리고, 신기한 점은, 13이라는 원소에 대해서 앞서 적용하지 않았던 모든 연산을 적용합니다. 그리고 더 재미있는 점은, h값이 소수로만 구성되어 있는 것을 알 수 있네요. 13과 다음 소수의 사이는 조금 더 길기 때문에, 한번만 더 수행해봅시다.

> take 7 primes'
[2,3,5,7,11,13{h = 2, x = 14}
{h = 2, x = 15}
{h = 3, x = 15}
{h = 2, x = 16}
{h = 2, x = 17}
{h = 3, x = 17}
{h = 5, x = 17}
{h = 7, x = 17}
{h = 11, x = 17}
{h = 13, x = 17}
,17]

14, 15, 16 각 숫자들은 앞에서 구해둔 소수들에 의해 제거됩니다. 즉, 새 소수가 추가되었을 때 다시 처음부터 연산을 수행하지 않고, 앞에서 연산하다 말았던 작업을 더 진행하게 됩니다. 즉, 미뤄두었던 필터가 실행되는 경우입니다. h의 값들이 모두 소수입니다. 즉, 4, 6과 같은 소수가 아닌 수들은 이미 앞에서 제거되고 없어 더이상 비교 작업을 수행하지 않고 있습니다.

앞에서 구해둔 소수가 10000개를 넘어간다면 10001번째 수는 10000번의 연산 작업을 통해 비교될 것입니다.

추가로, 소수가 짝수인 경우는 2 말고는 없습니다. 따라서, 홀수들만 갖고 연산하도록 변경한다면 시간이 조금 더 개선될 것입니다.

primes' :: [Int]
primes' = 2 : primeLoop [3, 5 ..]
 where
  primeLoop :: [Int] -> [Int]
  primeLoop []       = []
  primeLoop (h : xs) = h : primeLoop (filter (\x -> x `mod` h /= 0) xs)

결론

이번에 소수 구하기를 다양한 방법으로 구현해보면서 얻은 결론은, 일단 함수가 Pure하다는 것은 결국 Side effect가 없이 똑같이 재현이 가능하다는 것이고, 그 말은 연산 값을 다시 재활용해도 실제 동작에 영향이 없다는 것을 의미합니다. Side effect가 있을 때는, 그것이 Side effect가 있다는 것을 알게 하는 기법을 써야한다고 합니다. 일명 하스켈이 정말로 좋아하는 IO 모나드…

이제 이번화는 여기에서 끝내겠습니다. 다음화에서는 다른 주제로 예시를 만들어보겠습니다.