-----------------
title : 컴파일러는 정수 나눗셈을 어떻게 최적화 할까?
cat_title : 컴파일러는 정수 나눗셈을 어떻게 최적화 할까?
path : /프로그래밍
publish_date : 2020-07-06
-----------------

얼마전에 어셈블리를 살펴보다가, 정수 나눗셈과 나머지 연산에 대해서 최적화 옵션을 키게 되면 컴파일러가 매우 흥미로운 코드를 생성한다는 점을 발견하였습니다. 이 점에 대해서 간략하게 글을 적어보고자 합니다.

### 일반적인 경우

일반적으로 CPU 에는 정수 나눗셈을 하는 명령어인 `div` 가 있습니다. 예를 들어서 아래와 같이 두 개의 인자를 나누는 함수를 살펴봅시다.

```cpp-formatted
int func(unsigned int a, unsigned int b) { return a / b; }
```

위 코드를 컴파일 하게 되면

```asm
func(unsigned int, unsigned int):
        mov     eax, edi  # eax <- a
        xor     edx, edx  # edx <- 0
        div     esi # edx:eax 를 esi (= b) 로 나눈다.
        ret
```

와 같은 어셈블리를 생성합니다. `div` 어셈블리는 조금 특이한데, 위와 같이 인자로 32 비트 레지스터를 전달하게 되면 `edx:eax` 를 하나의 64 비트 짜리 정수로 생각한 뒤, 해당 값을 `esi` 로 나누게 됩니다. 그리고 그 몫은 `eax` 에 들어가게 되죠.

매우 간단한 어셈블리 입니다. 위와 같이 피제수(= a) 와 제수(= b) 모두 컴파일 타임에 알 수 없는 경우라면 위와 같이 제너릭한 어셈블리를 생성하는 것이 최선입니다.

### 제수를 컴파일 타임에 알 수 있다면?

반면에 아래와 같은 코드를 살펴봅시다.

```cpp-formatted
int func(unsigned int d) { return d / 100; }
```

이 경우 컴파일 시에 사뭇 다른 어셈블리를 생성합니다.

```asm
func(unsigned int):
        mov     eax, edi
        imul    rax, rax, 1374389535
        shr     rax, 37
        ret
```

띠용? `1374389535` 이라는 숫자는 도대체 어디서 온 것일까요? 그리고 왜 이번에는 `div` 명령어를 사용하지 않을까요?

#### div 는 느리다

일단 무엇보다도 CPU 에서 나눗셈 연산은 **무지하게** 느립니다. 기본적으로 덧셈과 뺄셈 연산이 보통 1 사이클 안에 되고 곱셈도 이 보다 살짝 느린 3 ~ 4 사이클 안에 수행됩니다. 그런데 (32 비트) 정수 나눗셈의 경우 20 ~ 26 사이클이 걸리죠. 이 시간 동안 곱셈을 무려 5 번이나 더 할 수 있습니다.

물론 예외는 있습니다. 2 의 지수승으로 나누는 연산은 CPU 가 굉장히 빨리 수행할 수 있습니다. 예를 들어서 $$2^5$$ 으로 나눈 다는 것은 사실상 **오른쪽으로 쉬프트 5 번** 하는 것과 마찬가지 이기 때문입니다. 그리고 CPU 상 쉬프트 연산과 같은 간단한 비트 연산들은 모두 1 사이클 안에 수행됩니다. 엄청나게 빠르죠!

따라서 컴파일러는 되도록이면 나눗셈을 곱셈으로 바꿀 수 있으면 바꾸고, 불가피 하게 나눗셈을 하게 된다면 되도록이면 2 의 지수승으로 나눗셈을 해버리려고 하게 됩니다.

위 경우도 마찬가지 입니다. $$ [\frac{d}{100}] $$ 은 사실 $$ [d \times \frac{1}{100}]$$ 과 똑같죠. \sidenote{참고로 [ ] 는 가우스 연산자로 해당 실수를 정수로 버림하라는 의미 입니다. 예를 들어서 [2.3] 은 2 고 [1.9999] 는 1 이 되죠. 컴퓨터에서 정수 나눗셈을 수행하는 이치와 같습니다.}

문제는 정수로 $$\frac{1}{100}$$ 을 표현할 길이 없다는 점입니다. 하지만 만약에 분자 분모에 적당히 큰 $$2^n$$ 을 곱해준다면 어떨까요? 예를 들어서 $$2^{10}$$ 을 곱해준다고 생각해봅시다.

$$[\frac{d}{100}] = [\frac{d \cdot 2^{10}}{100 \cdot 2^{10}}] = [d \cdot \frac{2^{10}}{100} \cdot \frac{1}{2^{10}}]$$

마지막 $$\cdot \frac{1}{2^{10}}$$ 는 2 의 지수승 나눗셈 이므로 쉬프트 연산을 통해 아주 빠르게 수행할 수 있습니다. 그런다면 $$\frac{2^{10}}{100}$$ 은 어떨까요? $$2^{10}$$ 이 `1024` 이니까, 대충 $$\frac{2^{10}}{100} \approx 11$$ 로 근사해봅시다.

그렇다면 $$[\frac{d}{100}]$$ 을 대충 $$[d \cdot 11 \frac{1}{2^{10}}]$$ 으로 근사할 수 있을까요?  예를 들어서 $$d = 100$$ 이면 전자의 경우 1 이 되지만 후자의 경우도 마찬가지로 $$[\frac{1100}{1024}]$$ 가 되서 1 이 됩니다.

하지만 조금만 생각해보면 작은 $$d$$ 에 대해서도 오차가 발생한다는 사실을 알 수 있습니다. 에를 들어서 $$d$$ 가 5000 이라면 전자는 50 이 되지만 후자는 $$[\frac{55000}{1024}]$$ 인 53 이 됩니다.

이렇게 오차가 발생하는 이유는 곱해주는 값인 $$2^{10}$$ 이 너무 작기 때문입니다. 그렇다면 얼마나 큰 2 의 지수승을 구해줘야지 $$d$$ 가 `unsigned int` 의 범위인 0 부터 4,294,967,295 까지의 모든 정수들에 대해서 올바른 결과값을 돌려줄 수 있을까요?

편의상 $$2^n$$ 을 곱해주었다고 생각해봅시다.

이 경우 우리가 계산하려는 식을 엄밀하게 쓰면

$$[[\frac{2^n + M}{100} \cdot d] \cdot \frac{1}{2^n}]$$

이 되는데, $$M$$ 을 $$2^n + M$$ 이 100 으로 나누어 떨어지게 하는 최소한의 값으로 잡으면 ($$M \equiv -2^n \pmod{100}$$) 첫 번째 $$[]$$ 안에 항상 정수가 되서 그냥

$$[\frac{2^n + M}{100} \cdot d \cdot \frac{1}{2^n}]$$

와 같이 간단히 쓸 수 있습니다. 위 식을 다시 써보면

$$[\frac{1}{100}(d + \frac{M \cdot d}{2^n})]$$

라 쓸 수 있습니다.

자, 그러면 원래 계산하려던

$$[\frac{d}{100}]$$

과 비교해봅시다. 전자와 후자가 `int` 범위 안의 $$d$$ 에 대해서 모두 같으려면 $$n$$ 이 얼마가 되야 할까요? 만약에, $$\frac{M\cdot d}{2^n}$$ 이 1 보다 항상 작게 되면 위 두 값은 항상 일치하게 됩니다.

왜냐하면 $$[99/100]$$ 이나 $$[99.9999/100]$$ 이나 0 이 되는 것은 매한가지 일테니까요. 반면에 그 차이가 1 이상이 된다면 더이상 같다고 보장할 수 없게 됩니다. 왜냐하면 $$[99/100]$$ 은 0 이지만 $$[100/100]$$ 은 1 이니까요!

그런데 $$d$$ 의 최대값은 $$2^{32} - 1$$ 이므로, 우리가 쓰고자 하는 조건은

$$\frac{M \cdot (2^{32} - 1)}{2^n} < 1$$

이 되겠습니다. 자, 그럼 $$n$$ 이 32 이면 충분할까요? $$2^{32}$$ 는 4294967296 이므로 $$M$$ 은 4 가 될 것입니다. 하지만 $$\frac{4 \cdot (2^{32} - 1)}{2^32} \approx 4$$ 가 되어서 조건을 만족시키지 않죠.

조건을 만족하는 $$n$$ 을 하나씩 차근차근 찾아보면

* $$n = 33$$ : $$2^{33}$$ 이 `8589934592` 니까 $$M=8$$ 이 되서 $$\frac{8 \cdot (2^{32} - 1)}{2^{33}} \approx 4 > 1$$
* $$n = 34$$ : $$2^{34}$$ 이 `17179869184` 니까 $$M=16$$ 이 되서 $$\frac{16 \cdot (2^{32} - 1)}{2^{34}} \approx 4 > 1$$
* $$n = 35$$ : $$2^{35}$$ 이 `34359738368` 니까 $$M=32$$ 이 되서 $$\frac{32 \cdot (2^{32} - 1)}{2^{35}} \approx 4 > 1$$
* $$n = 36$$ : $$2^{36}$$ 이 `68719476736` 니까 $$M=64$$ 이 되서 $$\frac{64 \cdot (2^{32} - 1)}{2^{36}} \approx 4 > 1$$
* $$n = 37$$ : $$2^{37}$$ 이 `137438953472` 니까 $$M=28$$ 이 되서 $$\frac{28 \cdot (2^{32} - 1)}{2^{37}} \approx 0.88 < 1$$

따라서 위를 만족하는 최소의 $$n$$ 은 37 임을 알 수 있습니다. 그리고 곱할 수는 실제로 $$\frac{137438953472 + 28}{100}$$ 인 1374389535 임을 확인할 수 있죠!

### 위 방법이 언제나 가능한 것은 아니다.

그렇다면 아래 코드는 어떻게 바꿀 수 있을까요?

```cpp-formatted
int func(unsigned int d) { return d / 7; }
```

일단 조건은 $$M \equiv -2^{n} \pmod{7}$$ 이어야 하고, $$\frac{M \cdot (2^{32} - 1)}{2^n} < 1$$ 여야 하죠. 이를 만족하는 최소의 $$n$$ 은 계산해보면 $$n = 35$$ 임을 알 수 있습니다. 따라서 곱해줘야 하는 값은 $$[\frac{2^{35}}{7}] + 1$$ 인 4908534053 가 됩니다.

한 가지 문제가 있습니다. 4908534053 는 `unsigned int` 의 최대 범위를 벗어나는 크기 입니다. 즉 $$2^{32}$$ 보다 큰 수이죠. 그런데 이 수를 예를 들어서 $$2^{31}$$ 과 곱하게 되면 그 결과는 $$2^{64}$$ 를 벗어나게 됩니다.

```asm
        imul    rax, rax, 4908534053
```

따라서 이 부분에서 `rax` 에 제대로 된 결과값이 들어갈 수 없습니다.

그렇다면 이 경우 컴파일러는 어떤 코드를 생성할까요?

```asm
func(unsigned int):
        mov     eax, edi
        imul    rax, rax, 613566757
        shr     rax, 32
        sub     edi, eax
        shr     edi
        add     eax, edi
        shr     eax, 2
        ret
```

흠, 조금 읽기 쉬운 C 코드로 바꿔보면 아래와 같습니다.

```cpp-formatted
int func(unsgiend int d) {
  return ((d - ((d * 613566757) >> 32)) / 2 + ((d * 613566757) >> 32)) / 4
}
```

수식으로 표현해보자면

$$[\frac{[\frac{d - [\frac{613566757 d}{2^{32}}]}{2}] + [\frac{613566757 d}{2^{32}}]}{4}]$$

가 됩니다.

처음에 저는 이 식을 보았을 때 도대체 이게 어떠한 방식으로 유도되었는지 도무지 감을 잡을 수 없었습니다. 다른 수들을 가지고 컴파일 해 보니 공통적으로 위와 같은 꼴의 식으로 컴파일 된다는 사실을 알 수 있었습니다.

아쉽게도 위 형태의 식이 **어떠한 방식으로** 유추되었는지는 모르겠습니다. 하지만 며칠간의 삽질 끝에 위와 같은 형태의 식이 $$n < 2^{32}$$ 인 모든 양의 정수들의 정수 나눗셈을 정확히 계산하는 데에는 충분하다는 사실을 증명할 수 있었습니다. 아래 그 사실을 간단히 소개하고보자 합니다.

#### 정리

적당한 $$C < 2^{32}$$ 와 $$X < 32$$ 에 대해서

$$[\frac{[\frac{n - [\frac{C n}{2^{32}}]}{2}] + [\frac{C n}{2^{32}}]}{2^X}] = [\frac{n}{d}]$$

가 $$2^{32}$$ 보다 작은 양의 정수 $$n$$ 에 대해서 항상 성립합니다.

위를 증명하기 위해 먼저 적당한 $$C$$ 와 $$X$$ 를 어떻게 찾을 수 있는지 살펴봅시다.

#### 보조 정리

모든 $$D < 2^{32}$$ 에 대해서

$$0 \le \frac{2^{32} + C}{2^{33 + X}} - \frac{1}{D} \le \frac{1}{D} \frac{1}{2^{32}}$$

를 만족하는 $$C < 2^{32}$$ 와 $$X < 32$$ 가 항상 존재한다. (즉 적당한 $$C$$ 와 $$X$$ 가 있어서 $$\frac{1}{D}$$ 를 충분히 근사할 수 있다는 의미 입니다.)

증명은 다음과 같습니다. 먼저 $$D < 2^{32}$$ 이므로 $$2^{X} \le D < 2^{X + 1}$$ 이 성립하도록 $$X$$ 를 잡을 수 있습니다. 편의상 $$D = 2^X + t$$ 라고 합시다. 대입을 하자면 $$C$$ 는 아래 식을 만족해야 하겠죠.

$$0 \le \frac{2^{32} + C}{2^{33 + X}} - \frac{1}{2^X + t} \le \frac{1}{2^X + t} \frac{1}{2^{32}}$$

양변에 $$2^X$$ 를 곱해주면

$$0 \le \frac{2^{32} + C}{2^{33}} - \frac{1}{1+ \frac{t}{2^X}} \le \frac{1}{1 + \frac{t}{2^X}} \frac{1}{2^{32}}$$

위 식을 정리하면 아래와 같은 식을 얻을 수 있습니다.

$$\frac{2^{33}}{1+ \frac{t}{2^X}} \le 2^{32} + C \le \frac{2^{33} + 2}{1+ \frac{t}{2^X}} = \frac{2^{33}}{1+ \frac{t}{2^X}} + \frac{2}{1+ \frac{t}{2^X}}$$

자 근데, 우리는 $$t < 2^X$$ 라는 점을 알고 있습니다. 따라서 $$\frac{2}{1+ \frac{t}{2^X}} \in [1, 2]$$ 이겠죠. 따라서 일단 $$\frac{2^{33}}{1+ \frac{t}{2^X}}$$ 은 최소한 $$2^{32}$$ 보다 크므로 $$C > 0$$ 을 만족할 수 있습니다. 문제는 정수 $$C$$ 가 존재하기 위해서는 위 구간의 크기가 1 이상이어야 하는데, $$\frac{2}{1+ \frac{t}{2^X}} \ge 1$$ 이므로, 위 범위 내에 들어가는 정수 $$C$$ 는 반드시 존재하게 됩니다.

보조 정리를 조금 정리해보자면;

$$\frac{n}{D}2^{1+X} - n \le \frac{Cn}{2^{32}} \le \frac{n}{D}2^{1+X} (1 + \frac{1}{2^{32}}) - n$$

와 같은 관계식을 얻을 수 있습니다. 따라서

$$[\frac{n}{D}2^{1+X}] - n \le [\frac{Cn}{2^{32}}] \le [\frac{n}{D}2^{1+X} (1 + \frac{1}{2^{32}})] - n$$

로 쓸 수 있죠.

자 이제 그렇다면 원래의 정리로 다시 돌아가보겠습니다. 먼저 $$[x] \le x$$ 가 항상 성립하므로,

$$[\frac{[\frac{n - [\frac{C n}{2^{32}}]}{2}] + [\frac{C n}{2^{32}}]}{2^X}]  \le [\frac{\frac{n - [\frac{C n}{2^{32}}]}{2} + [\frac{C n}{2^{32}}] }{2^X}] = [\frac{n + [\frac{C n}{2^{32}}]}{2^{X+1}}]$$

가 되는데, 보조정리를 이용하면

$$ [\frac{n + [\frac{C n}{2^{32}}]}{2^{X+1}}] \le [\frac{[\frac{n}{D}2^{1+X} (1 + \frac{1}{2^{32}})]}{2^{X+1}}] \le [\frac{n}{D} (1 + \frac{1}{2^{32}})] \le \frac{n}{D}$$

가 됩니다. (이 때 마지막 조건은 $$n < 2^{32}$$ 인 가정 하에 성립)

이번에는 최소값을 계산해봅시다. 마찬가지로 보조정리를 대입해보면;

$$ [\frac{[\frac{n - [\frac{C n}{2^{32}}]}{2}] + [\frac{C n}{2^{32}}]}{2^X}] \ge [\frac{1}{2^X}(-[\frac{1}{2}[\frac{n}{D}2^{1+X}(1+\frac{1}{2^{32}})]] + [\frac{n}{D}2^{1+X}])]$$

의 관계식을 얻을 수 있습니다. 여기서

$$[\frac{1}{2^X}(-[\frac{1}{2}[\frac{n}{D}2^{1+X}(1+\frac{1}{2^{32}})]] + [\frac{n}{D}2^{1+X}])] \ge [\frac{[\frac{n}{D}2^{1+X}]}{2^{1+X}}] \ge [\frac{n}{D}]$$

를 얻을 수 있습니다.

따라서

$$ [\frac{n}{D}] \le [\frac{[\frac{n - [\frac{C n}{2^{32}}]}{2}] + [\frac{C n}{2^{32}}]}{2^X}] \le \frac{n}{D}$$

인데 $$[\frac{[\frac{n - [\frac{C n}{2^{32}}]}{2}] + [\frac{C n}{2^{32}}]}{2^X}]$$ 는 항상 정수니까, 결국

$$[\frac{[\frac{n - [\frac{C n}{2^{32}}]}{2}] + [\frac{C n}{2^{32}}]}{2^X}] = [\frac{n}{d}]$$

이 $$n < 2^{32}$$ 인 모든 양의 정수 $$n$$ 에 대해 언제나 성립하게 됩니다.

정말 신기합니다!

실제로 

```py
def func(D):
    max_C = 0
    X = 0
    for i in range(32):
        C = int((2 ** (33 + i)) / D) - (2 ** 32) + 1
        if C < (2 ** 32):
            max_C = C
            X = i
        else:
            return (max_C, X)
```

위와 같은 파이썬 코드를 통해 컴파일러가 사용하는 $$C$$ 와 $$X$$ 의 값을 찾을 수 있습니다. 예를 들어서 `func(117)` 을 하면 403800345 와 6 을 리턴하는데, 실제로

```cpp
int func(unsigned int d) {
    return d / 117;
}
```

을 컴파일 한다면

```asm
func(unsigned int):
        mov     eax, edi
        imul    rax, rax, 403800345
        shr     rax, 32
        sub     edi, eax
        shr     edi
        add     eax, edi
        shr     eax, 6
        ret
```

로 컴파일 됨을 알 수 있습니다.

### 더 궁금 하다면..

만일 정수 나눗셈 연산 최적화에 더 관심이 있으신 분들은 [libdivide](https://libdivide.com/) 라이브러리를 살펴보는 것을 추천드립니다. 이 라이브러리는 일반적인 명령어 말고도, `SSDE2` 와 같은 벡터 연산 명령어들을 적극 활용해서 정수 나눗셈 최적화를 더욱 수월하게 수행할 수 있습니다!

덧붙여 [Hacker's Delight](https://www.amazon.com/Hackers-Delight-2nd-Henry-Warren/dp/0321842685) 의 10 장을 보면 정수 나눗셈을 수행하는 여러가지 방법을 소개하고 있으니 이 부분이 더 궁금하신 분들은 참고하시면 도움이 될 것 같습니다.


