윤년을 체크하는 방법

컴공 학생이라면 1~2학년에 한번쯤은 터미널에 달력 표시하는 프로그램을 만들어 봤을 것이다. 나도 1학년 때 순수 C로 구현했던 기억이 있는데, 그 때 윤년의 정의를 처음으로 제대로 알았다.

4년마다 한번씩 윤년
(근데) 100년마다는 윤년 아님
(근데) 400년마다는 윤년

아마 조건문 가르칠 때 다중 조건 가르치기 용이한 예시라 단골로 나오지 않았나 싶다. 아무튼 이걸 단순히 구현하면 아래가 된다.

bool is_leap_year(uint32_t y) {
	if ((y % 4) != 0) return false;
	if ((y % 100) != 0) return true;
	if ((y % 400) == 0) return true;
	return false;
}

근데 아까 이런 코드를 봤다. libstdc++에도 기존의 코드 대신 포함된 것 같다.

bool is_leap_year_fast(uint32_t y) {
    return ((y * 1073750999) & 3221352463) <= 126976;
}

이게 어떻게 가능한지, 찾는 건 어떻게 찾은 건지 궁금해서 설명하는 글을 읽어봤다.

가장 기초적인 최적화

우선 가장 기초적인 최적화로 (y % 100) != 0을 (y % 25) != 0으로 바꿀 수 있다.
- (y % 4) != 0을 확인하는 단계에서 이미 y는 4의 배수이기 때문에(= $4|y$ )
- $4|y\ \wedge\ 25|y \iff 100|y$ 을 $25|y \iff 100|y$ 로 축약할 수 있다.
이어서, (y % 400) == 0를 (y % 16) == 0으로 바꿀 수 있다.
- (y % 100) != 0이 거짓인 단계에서 이미 y는 25의 배수이기 때문에(= $25|y$ )
- $25|y\ \wedge\ 16|y \iff 400|y$ 을 $16|y \iff 400|y$ 로 축약할 수 있다.
위의 최적화는 4와 16이라는 점 때문에 (2의 배수이기 때문에) 나머지 연산을 비트마스킹으로 바꿀 수 있어서 유용하다.
- if ((y % 4) != 0) -> if ((y & 3) != 0)
  - 2진수 형태에서 가장 아래 두 자리가 0이면: 4의 배수
- if ((y % 16) == 0) -> if ((y & 15) == 0)
  - 2진수 형태에서 가장 아래 네 자리가 0이면: 16의 배수
이 부분을 처음엔 당연하다고 생각하고 넘어갔는데, 다시 생각해보니까 제대로 이해가 안되어서 찾아보니 Chinese remainder theorem; CRT까지 흘러갔다. 조합식 개념을 몰라서 그랬던 거 같고, 지금 내용과는 크게 관련 없어서 위 부분까지만 이해하고 다시 넘어갔다.

임의 홀수의 배수인지 확인하기

4와 16은 비트마스킹으로 처리하고, 남은 25는 곱셈 연산 및 비교연산으로 바꾼다.
이 최적화는 gcc에 기본으로 적용되어 있다.
- gcc로 컴파일하고 다시 C로 역번역하면, 아래처럼 바뀐다.
- if ((y % 25) != 0) -> if (y * 3264175145u > 171798691u)
- 왜 이런 숫자가 나올까?
컴퓨터에게 나눗셈은 곱셈보다 훨씬 느리다.
- 나눗셈은 곱셈처럼 병렬처리가 불가능하고, 앞단계 연산을 기다려야 하기 때문
- 그래서, 보통 컴파일러는 최적화를 위해 역수를 곱하는 방식으로 계산한다.
if (y * 3264175145u > 171798691u)의 의미
- $y \div 25$ $y \div 25$ 를 $y \times \dfrac{1}{25}$ $y \times \frac{1}{25}$ 로 계산하고 싶은데, $\dfrac{1}{25}$ $\frac{1}{25}$ 를 2진수로 나타내기 위해 32비트 고정소수점으로 근사하려고 한다.
  - 다만, 지금 알고 싶은 것은 $y\div 25$ 의 값이 아니라, 실제로 $y$ 가 25의 배수인지 아닌지만 궁금하다. 그래서 $y$ 를 25의 배수로 가정하고 다른 방식으로 계산한다.
- $y$ 가 25의 배수라면, $y \times \dfrac{k}{25}$ 는 k가 어떤 정수든 항상 정수가 된다(소수부가 정확히 0이 된다). 반대로 25의 배수가 아니라면, 소수부가 0이 아닌 값이다.
- 계산하기 편한 $k$ $k$ ( $2^{32}\times \dfrac{k}{25}$ $2^{32} \times \frac{k}{25}$ 가 정수에 최대한 가깝게 나오는 k)를 구해보면 $k=19$ $k = 19$ 일 때 $2^{32}\times \dfrac{19}{25}=3264175144.96$ $2^{32} \times \frac{19}{25} = 3264175144.96$ 이 된다.
  - $k$ 를 구하는 과정은 정수론 혹은 수치해석적 과정이 필요하고 좀 복잡한 거 같아서 여기서는 넘어갔다.
  - 만약 이 수를 현실처럼 곱할 수 있어서 곱한다면, y가 25의 배수라면 소수부(여기선 32비트 고정소수점이므로 하위 32비트)가 정확히 0일 것이다. (== 0)
- 실제로는 소수를 곱할 수 없기 때문에, 저 값을 올린 3264175145u를 곱하면, 소수부에 오차가 생길 것이다. 그 오차는 올린 오차가 0.04이기 때문에 $0.04 \times (2^{32}-1) = 171798691.8$ 이하일 것이다. (<= 171798691.8)
- 즉, 정수로 고치면 171798691u를 초과하면 소수부가 0이 아니고 곧 25의 배수가 아닌 것이다. 이를 식으로 정리하면 if (y * 3264175145u > 171798691u)
- y가 uint32_t이므로, 코드 상에서는 3264175145u를 곱했을 때 32비트 레지스터에 남아있는 수는 전부 32비트 소수부이다.
위 내용들을 정리하면 코드가 아래처럼 정리된다.
- (사실 최근의 컴파일러는 가장 위 단순 구현으로 적어도 알아서 아래처럼 최적화 해준다.)

bool is_leap_year2(uint32_t y) {
	if ((y & 3) != 0) return false;
	if (y * 3264175145u > 171798691u) return true;
	if ((y & 15) == 0) return true;
	return false;
}

비트 조작 접근법

위 코드도 충분히 잘 최적화된 좋은 코드지만, 분기가 있기 때문에 아예 분기 없는 형태로 (연산만으로) 만들 수 있는 형태를 블로그 저자가 시도해봤다.
- 블로그 저자가 후보로 잡은 ((y*f) & m) <= t가 z3를 시도해보기 좋은 형태였다.
z3은 Microsoft Research에서 만든 ...(작성 중)

출처:

https://github.com/gcc-mirror/gcc/commit/28b1b55c60464da23c65c6a8a9515c29f8d67ad8
https://hueffner.de/falk/blog/a-leap-year-check-in-three-instructions.html
https://arxiv.org/abs/1902.01961