부동소수점 비교에서 무작정 epsilon을 쓰는 습관이 왜 종종 나쁜 선택인지, 그리고 더 나은 접근이 무엇인지 다양한 사례로 설명합니다.
2026 Apr 14 참고: 이 글의 제목은 의도적인 클릭베이트다. 나는 이 주장에 실제로 동의하긴 하지만, 더 솔직한 제목은 아마 이런 느낌일 것이다: epsilon을 사용해 부동소수점을 비교하는 것은 괜찮지 않다.
아마 “부동소수점 수는 절대로 정확한 같음으로 비교하면 안 되고, 반드시 어떤 식으로든 epsilon 비교를 사용해야 한다”는 주문을 들어본 적이 있을 것이다. 예를 들면 이런 식이다.
bool approxEqual(float x, float y)
{
return abs(x - y) < 1e-4f;
}
내가 15년 넘게 코드를 작성해 오면서 — 그리고 그 코드는 기하, 그래픽스, 물리, 시뮬레이션 등을 자주 다뤄서 부동소수점을 매일같이 써야 했다 — 이런 epsilon 비교가 실제로 좋은 해법인 경우는 한두 번밖에 만나지 못했다. 거의 언제나 더 나은 해법이 있었고, 그 해법은 코드를 어떤 방식으로든 다시 쓰는 것이거나, 아니면 그냥 x == y처럼 부동소수점을 그대로 비교하는 것이었다. 그리고 거의 언제나 epsilon 해법은 실제로는 가능한 선택지 중 최악에 가까운 것 중 하나였다.
어떤 종류의 epsilon을 추가하는 것이 첫 반사작용처럼 느껴질 수 있지만, 사실은 훨씬 더 나은 — 그리고 종종 훨씬 더 단순한 — 해법이 있는 예시를 여러 개 보여주겠다. 하지만 먼저, 부동소수점 수에 대해 이야기해 보자.
epsilon 비교라는 발상 전체는, 부동소수점 수를 어떤 무작위적인 블랙박스 기계처럼 보는 일반적인 인식에서 비롯된 듯하다. 이 기계는 컴퓨팅의 신들이 강요해서 가끔 부정확한 결과를 뱉는다고 여겨지곤 한다. 하지만 실제로는 이것은 꽤 결정론적인(컴파일러 옵션, CPU 플래그 등을 제외하면) 그리고 매우 표준화된 시스템이다.
부동소수점 수가 본질적으로 부정확한 이유는 가능한 모든 실수를 표현할 수 없기 때문이다. 사실 유한한 양의 메모리로는 어떤 방식으로도 그럴 수 없다. 수학이 원래 그렇기 때문이다 — 실수는 (심지어 유리수만 보더라도) 너무 많다. 그리고 우리는 아마도 각 수에 대해 고정된(단지 유한한 것이 아니라) 비트 수만 할당하고 싶을 것이므로, 표현 가능한 수는 유한한 집합뿐이라는 사실을 받아들여야 한다(구체적으로는 많아야 개다). 나머지 수들에 대해서는 근사치를 다뤄야 한다.
하지만 이 글을 부동소수점 수가 어떻게 표현되는지에 대한 강의로 만들고 싶지는 않다. 그건 wikipedia가 충분히 잘 설명한다고 생각한다. 더 깊이 있는 자료로는 David Goldberg의 고전적인 글이나, 같은 아이디어를 따르는 좀 더 최근의 글을 보라. 더 중요한 것은, 부동소수점 수의 이런 “부정확성”이 그 동작에 어떤 “불확실성”이나 “무작위성”이 있다는 뜻은 아니라는 점이다.
예를 들어, 두 부동소수점 수에 대한 단일 산술 연산(예: 덧셈, 곱셈 등)은, 입력이 정확하다고 가정했을 때 실제 정답에 가장 가까운 부동소수점 수를 결과로 내놓아야 한다(동률인 경우, 즉 표현 가능한 두 수가 실제 정답과 똑같이 가깝다면 반올림 규칙이 있다). 결과가 근사치이긴 해도, 여전히 가능한 한 진실에 가깝다는 보장이 있고, 매우 결정론적이라는 점에 주목하라.
하지만 이것은 우리가 익숙한 수학 공식이 부동소수점에서는 항상 성립하지 않는다는 뜻이기도 하다. 예를 들어 덧셈(그리고 곱셈)은 교환법칙()을 만족하는 것이 보장되지만, 반드시 _결합법칙_을 만족하는 것은 아니다. 즉 가 될 수 있다. 그런 예시는 찾기 꽤 쉽다. 여기에 32비트 부동소수점에서 동작하는 예가 있다.
// Outputs 0.89999998
std::cout << std::setprecision(8) << ((0.2f + 0.3f) + 0.4f) << '\n';
// Outputs 0.90000004
std::cout << std::setprecision(8) << (0.2f + (0.3f + 0.4f)) << '\n';
이 식들이 같지는 않지만, 매우 가깝다는 점에 주목하라(차이는 약 6e-8이다). 그리고 부동소수점 표준은 우리가 이 차이가 얼마나 클 수 있는지를 _예측_하는 데 쓸 수 있는 몇 가지 보장을 제공한다.
그렇다면 문제는 무엇일까? 어떤 계산의 결과가 근사치라는 것을 알고 있고, 어떤 기대 결과와도 근사적으로 비교한다. 얼핏 보면 꽤 그럴듯하다.
정말 그럴까?
나는 이런 epsilon 비교에 대해 크게 3가지 문제가 있다고 본다.
두 번째가 아마 가장 까다롭다. 프로그램의 한 부분은 2차원 점 두 개가 맨해튼 거리로 1e-4 이하이면 같다고 보고, 다른 부분은 점 두 개가 L-inf 거리(좌표 차이의 최댓값)로 1e-6 이하이면 같다고 보고, 입력 점들 자체는 또 다른 알고리즘으로 생성되고, 그러다 보면 입력-출력 불변식이 전부 꼬여 버린다. 그러면 선 렌더링이 망가지는데, 오직 이 특정 시나리오에서, 이 특정 데이터로, 밤에 보름달이 뜰 때만 재현된다. 디버깅 잘 해보길 바란다.
선 렌더링이 가끔 잘못되는 정도면 그나마 덜 심각하지만, 이것은 엄청나게 다양한 방식으로 나타날 수 있고, 심하면 프로그램이 크래시할 수도 있다. 특히 다른 기하 코드와 많이 결합되면 정말 끔찍해진다. 나는 이런 경우를 많이 겪었고, 서로 맞지 않는 epsilon 비교가 근본 원인인 경우가 아주 흔했다.
문제는 이런 epsilon들이 거의 언제나 그냥 찍어서 정해진다는 것이다. 수많은 epsilon 중에서 어떤 하나를 고르는 올바른 방법은 없다. 이런 비교가 여러 개 있다면, 어떤 epsilon 조합도 모든 것이 올바르게 동작하도록 보장해 주지 못한다.
epsilon의 또 다른 문제는, 이런 비교가 추이적이지 않다는 점이다. 이것은 기술적인 트집처럼 들릴지 모르지만, 실제로는 대부분의 알고리즘이 비교 같은 것이 실제로 추이적이라고 가정한다. 그리고 이런 알고리즘은 비추이적인 비교를 넣으면 그냥 깨질 수 있다(헛소리를 내거나 심지어 크래시할 수도 있다).
그럼 우리는 무엇을 해야 할까? 문제를 _생각_해야 한다! 충격적이겠지만 말이다. 우리는 왜 애초에 부동소수점을 비교하고 있는가? 알고리즘을 믿지 못해서인가? 데이터를 믿지 못해서인가? CPU를 믿지 못해서인가? 정답은 하나가 아니니, 예시를 몇 가지 보자.
턴제 게임을 만든다고 하자. 유닛은 격자 위를 이동하고, 이동 포인트가 있어서 한 턴에 여러 번 움직일 수 있다. 하지만 UX를 위해 이전 이동이 끝난 뒤에만 다음 이동을 실행할 수 있게 하고 싶다.
이제 아마 유닛의 위치를 어떤 방식으로든 보간하고 있을 것이고, 그냥 목표 격자 칸으로 순간이동시키지는 않을 것이다. 그러니 플레이어가 다음 목표 칸을 고를 수 있도록 허용하기 전에, 정확히 언제 이동이 끝났는지를 확인해야 한다.
그냥 일정 시간을 기다릴 수도 있다(그리고 많은 경우 그것도 충분히 좋은 해법이다!). 하지만 유닛마다 애니메이션이 다르고 따라서 시간도 다르며, 애니메이션을 줄이는 접근성 설정도 있다. 그래서 애니메이션 시간에 의존하는 것은 좋은 생각이 아니라고 판단한다.
그러다 보니, 이동이 끝나는 순간은 정확히 유닛의 위치가 목표 칸 중심과 일치할 때라는 사실을 떠올린다. 그래서 이런 코드를 쓴다.
void update() {
if (selectedUnit.position != targetCell.center)
return;
// Do the frame update
}
그런데 이동이 실행된 뒤에도 아무 일도 일어나지 않고, 게임은 사실상 멈춘다. 왜냐하면 selectedUnit.position != targetCell.center 조건이 절대로 참이 되지 않기 때문이다. 전형적인 easing이 들어간 선형 보간에서는 이런 코드가
vec2f getPosition(float time) {
return lerp(start, end, easing(time));
}
부동소수점 반올림을 충분히 많이 일으켜서, time == 1.f일 때도 결과가 end와 절대 같아지지 않을 수 있다. 심지어 어떤 보간 방식에서는 애초에 영원히 참이 되지 않는 것이 의도이기도 하다!
젠장, 멍청한 부동소수점! — 하고 투덜거리며, 부동소수점의 성배인 epsilon 비교를 집어넣는다.
void update() {
if (distance(selectedUnit.position, targetCell.center) > 1e-4)
return;
// Do the frame update
}
이렇게 하면 문제가 해결되니, 뭐가 그렇게 나쁘다는 걸까? 이유는 여러 가지다.
1e-4라는 구체적인 epsilon은 보간 방식을 바꾸면 깨질 수 있다.그럼 어떻게 해결할까? 한 가지 방법은 일종의 허용 반경을 쓰는 것이다. 유닛이 목표 칸 중심에서 예를 들어 0.25 이내로 들어오면 애니메이션을 멈추고 플레이어가 다음 명령을 내릴 수 있게 한다. 그리고 실제로 좋은 값을 많은 테스트로 찾아낸다.
이게 epsilon과 뭐가 다를까? 사실 크게 다르지 않다. 다만 이제는 최소한 랜덤한 값이 아니라 무언가 근거가 있다는 점이 다르다. 여기서 진짜 문제는 _데이터 모델_과 _표현_을 섞고 있다는 점이다. 게임 상태 머신의 내부 동작은 어떤 3D 모델이 어디에 있는지에 신경 쓰면 안 된다. 물론 이것은 생각보다 어려운 일인 경우가 많다(대략 UX라는 것이 존재하는 이유이기도 하다). 하지만 충분히 가능하다.
이 경우 가장 좋은 해법은(내 생각에는!) 애초에 유닛이 이동을 완전히 끝낼 때까지 기다리지 않고도 사용자가 명령을 내릴 수 있게 하는 것이다. 사용자가 어떤 격자 칸을 클릭하면, 렌더링 쪽에는 이동 애니메이션이 필요하다는 사실만 알려 주고, 게임 상태의 내부 모델은 유닛이 이미 다음 칸에 있다고 생각하게 하면 된다.
이것을 구현하는 방법은 많다. 예를 들어 요청된 애니메이션을 큐에 넣고 하나씩 재생할 수도 있고, 목표값이 갑자기 바뀌어도 튼튼한 애니메이션/보간 방식을 사용할 수도 있다. 중요한 것은, 이것이 충분히 구현 가능하고, 비교적 간단하며, epsilon이 필요 없다는 점이다.
구면 선형 보간은 구 위의 점(즉 단위 벡터)을 보간하는 방법으로, 보간된 벡터가 일정한 회전 속도를 유지하는 곡선을 따라가게 만든다. 단순한 normalize(lerp(a, b, t))로는 안 된다 — 그렇게 하면 벡터는 양 끝 근처에서는 느리고 가운데에서는 빠르게 움직인다. slerp의 멋진 시각화가 있다. 이 함수는 종종 쿼터니언에 대해서만 구현되지만, 사실 임의 차원의 모든 벡터에 유용하다(다만 쿼터니언의 경우에는 와 가 같은 회전을 나타내므로 조금 다르다).
아주 편리하게도, 입력 벡터가 정규화되어 있다면 식은 꽤 단순해진다.
vec3 slerp(vec3 a, vec3 b, float t) {
float angle = acos(dot(a, b));
return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}
그런데 이 코드는 가끔 NaN을 만들어 낸다. 이유는 몇 가지다.
acos는 NaN을 반환한다.angle이 0이 될 수 있고, 그러면 0을 0으로 나누게 되어 또다시 NaN이 나온다.첫 번째 문제는 다루기 쉽다. 그냥 acos 인자를 clamp로 감싸면 된다.
vec3 slerp(vec3 a, vec3 b, float t) {
float angle = acos(clamp(dot(a, b), -1.f, 1.f));
return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}
두 번째 문제는 좀 더 미묘하다. 다행히 angle이 충분히 작으면, 구면 선형 보간은 일반적인 선형 보간으로 바뀐다. 그래서 이 경우를 감지해서 lerp(a, b, t)로 전환할 수 있다. 하지만 얼마나 작아야 “충분히 작은” 걸까?
여기에도 그냥 epsilon을 던져 넣고 싶어진다.
vec3 slerp(vec3 a, vec3 b, float t) {
float angle = acos(dot(a, b));
if (angle < 1e-4) {
return lerp(a, b, t);
}
return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}
하지만 이렇게 하면, 결과 벡터가 기대한 것과 가깝더라도, 코드의 정밀도는 원래보다 낮아진다. 그리고 다시 말하지만, 랜덤한 epsilon은 늘 “왜 하필 이 값인가?”라는 질문을 낳는다. 더 잘할 수 있다.
glm과 Eigen은 둘 다 그럴듯한 검사를 사용한다.
vec3 slerp(vec3 a, vec3 b, float t) {
float d = dot(a, b);
if (d > 1.f - FLT_EPSILON) {
return lerp(a, b, t);
}
float angle = acos(dot(a, b));
return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}
여기서 FLT_EPSILON은 정확히 이다 — 단정밀도 부동소수점 값 중에서, 부동소수점 덧셈 기준으로 가 되는 가장 작은 값이다. 솔직히 말하면, 이것도 내게는 랜덤한 epsilon보다 그다지 나아 보이지 않는다. FLT_EPSILON을 던져 넣는다고 해서, 그것이 우리 경우에 정확히 맞는 임계값이라는 증명이 없다면 문제가 해결되지는 않는다.
코드를 분석해 보자. 주된 문제는 angle이 0이 되면 NaN이 생긴다는 것이다. 부차적인 문제는 angle이 너무 작으면 정밀도 오차가 생길 수 있다는 것이다.
먼저 정밀도를 생각해 보자. angle은 그냥 임의의 수가 아니라 acos 호출 결과이고, 우리는 acos의 인자가 1에 가까울 때를 걱정하고 있다. 사실 angle 자체는 다소 중요하지 않다. 우리는 주로 각도 그 자체보다 sin(angle)에 더 관심이 있다.
이제 는 그냥 이다. 우리의 경우 x = dot(a, b)다. 일 때,
항은 신경 쓸 만큼 크지 않고, 는 그냥 1에 가까운 상수다. 여기서 핵심은 이다. 이 FLT_EPSILON처럼 아주 작은 값이라 해도, sin(angle) 항은 그 제곱근 정도가 된다. 작은 수에서는 제곱근이 그 수를 꽤 크게 만든다. 예를 들어 FLT_EPSILON은 대략 1.2e-07이지만, 그 제곱근은 대략 3.5e-04로 약 2000배 더 크다.
여기서 내가 말하고 싶은 것은 이것이다. 우리가 신경 쓰는 경우에는 acos의 인자가 1에 가깝다. 그러므로 angle 자체가 0에 가깝다는 사실보다, acos 인자가 가질 수 있는 표현 가능한 값들 사이의 간격이 정밀도 문제의 훨씬 더 큰 원인이다.
실제로 1보다 작은 표현 가능한 최댓값은 1 - FLT_EPSILON이 아니라 1 - FLT_EPSILON / 2.f인데, 이것은 여기서 FLT_EPSILON 상수가 다소 임의적으로 선택되었다는 내 주장에 어느 정도 힘을 실어 준다.
그러니 우리가 특별한 무언가를 원하지 않는 한, 정밀도 측면에서 각도가 작을 수 있다는 점은 무시하고, NaN이 생기는 문제에 집중해도 된다. 이 코드에서 인자가 유한값이라고 가정하면, NaN이 생기는 유일한 방법은 0으로 나누는 것이다. 그리고 그건 angle이 0일 때만 발생하며, angle이 0이 되는 것은 dot(a, b) >= 1일 때뿐이다. 앞서 이야기했듯이, dot(a, b)의 바로 다음으로 작은 표현 가능한 값은 1 - FLT_EPSILON / 2.f이며, 이때의 angle은 대략 sqrt(FLT_EPSILON)이고 이는 약 3.5e-04이다. 전체 범위가 대략 [1e-38 .. 1e38]라는 것을 생각하면 꽤 큰 부동소수점 값이므로, 적어도 그럭저럭 괜찮은 acos 구현이라면 이 경우 0을 반환하지는 않을 것이라고 확신할 수 있다.
따라서 dot(a, b) < 1인 경우에는 실제로 아무 문제도 없다! 우리가 확인해야 할 것은 dot(a, b) == 1인지, 혹은 앞서 추가한 clamp와 합치면 dot(a, b) >= 1인지뿐이다.
vec3 slerp(vec3 a, vec3 b, float t) {
float d = dot(a, b);
if (d >= 1.f) {
return lerp(a, b, t);
}
float angle = acos(dot(a, b));
return (sin((1 - t) * angle) * a + sin(t * angle) * b) / sin(angle);
}
(내적이 에 가까운 경우는 의도적으로 빼놓았다. 그 경우에는 보간 문제의 해가 유일하지 않기 때문에 따로 다루는 편이 낫다.)
또는 sin(angle) 자체를 직접 검사할 수도 있다. 다만 이 경우 모서리 케이스에서도 비싼 acos 호출을 아낄 수는 없다.
vec3 slerp(vec3 a, vec3 b, float t) {
float angle = acos(clamp(dot(a, b), -1.f, 1.f));
float s = sin(angle);
if (s == 0.f) {
return lerp(a, b, t);
}
return (sin((1 - t) * angle) * a + sin(t * angle) * b) / s;
}
벡터의 (유클리드) 길이를 계산하는 것이 꽤 흔하고 유용한 연산이라는 말을 들어도 놀라지 않을 것이다. 구현도 꽤 단순하다.
float length(vec3 v) {
return sqrt(dot(v, v));
}
혹은 dot 호출을 펼치면,
float length(vec3 v) {
return sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}
아마 이런 코드를 수십 번은 봤을 것이다. 뭐가 문제일까? 사실 별문제는 없다. 다만 아주 작은 벡터(길이가 1e-19 정도인 벡터)에 대해서는 0을 반환하고, 조금 더 큰 벡터에 대해서도 정밀도가 크게 손실된다. 결과 자체는 완전히 유효한 부동소수점 값일 수 있지만, 그 제곱, 즉 제곱근 안의 식이 부동소수점으로 적절히 표현하기엔 너무 작을 수 있다.
그래서 뭐가 문제일까? 어차피 그런 벡터는 사실상 0과 구별할 일이 거의 없다. 하지만 “길이가 0인 벡터는 오직 영벡터뿐이다”라는 불변식이 있으면 매우 편리할 것이고, 특히 부동소수점을 조심스럽게 다루고 싶을 때 코드 작성이 훨씬 쉬워진다.
사실 이를 위한 매우 간단한 방법이 있다. 벡터 좌표 절댓값들의 최댓값 M을 구한 뒤, M * length(v / M)을 반환하면 된다. 부동소수점 보장 덕분에 좌표 중 하나는 적어도 1이 되므로, 제곱근 안의 식은 최소 1이고 최대 — 즉 차원 수(vec3라면 3)다.
그러면 코드는 이렇게 바뀐다.
float length(vec3 v) {
float M = max(abs(v.x), abs(v.y), abs(v.z));
vec3 u = v / M;
return M * sqrt(u.x * u.x + u.y * u.y + u.z * u.z);
}
이건 어셈블리 명령 수가 대략 두 배 정도 되지만, 이것이 핫 루프 안에서 일어나는 게 아니라면, 추가 정밀도와 안전성을 위해 충분히 가치 있다고 생각한다.
유일한 문제는 벡터가 영벡터일 때(그렇다, 그냥 작은 게 아니라 정말 0일 때) 이 함수가 NaN을 반환한다는 것이다! 하지만 벡터 좌표 중 하나라도 0이 아니면(서브노멀 값이라 해도), M은 엄밀히 0보다 크고 알고리즘은 올바르게 동작한다(그리고 순진한 버전보다 정밀도가 더 나빠지는 일도 없다). 그러므로 올바른 수정은 말 그대로 M을 0과 비교하는 것이다.
float length(vec3 v) {
float M = max(abs(v.x), abs(v.y), abs(v.z));
if (M == 0.f) return 0.f;
vec3 u = v / M;
return M * sqrt(u.x * u.x + u.y * u.y + u.z * u.z);
}
만약 여기서 if (M < 1e-4)라고 썼다면, 이 함수가 존재하는 이유 자체를 망가뜨리는 셈이 된다.
선형 연립방정식을 푸는 것이 꽤 흔하고 유용한 작업이라는 말을 들어도 놀라지 않을 것이다. 물리와 공학 문제의 절반은 결국 어떤 선형 시스템을 푸는 문제로 귀결된다(나머지 절반은 고유값 방정식이다).
일반적인 비희소 시스템을 푸는 표준적인 방법은 사실상 하나뿐이다. 바로 Gauss-Jordan 소거법이다. (사실 QR 분해도 쓸 수 있다 — 알기로는 더 느리지만 수치적으로는 더 안정적이다.) 전체 알고리즘은 꽤 길지만, 아주 복잡하지는 않다. 그저 for 루프 몇 개와 계수들에 대한 산술 연산의 묶음일 뿐이다.
중요한 것은 이 알고리즘이 상당히 불안정할 수 있다는 점이다. 어느 시점에서는 행렬의 남은 부분에서 좌상단 원소를 가져와 그 원소로 행 전체를 나눈다. 이 값이 작으면 소거 과정의 상쇄(cancellation)에서 오는 부동소수점 오차가 증폭되고, 0이면 NaN이 생긴다.
하지만 실제로는 _부분 피벗팅_이라는 것을 하면 이 알고리즘은 안정적으로 쓸 수 있다(불안정성은 극히 드물고 보통 이론적인 수준이다). 즉, 현재 열에서 절댓값이 가장 큰 계수를 가진 행을 찾아 그 행을 원래 맨 위 행 대신 사용한다. 알고리즘이 조금 복잡해지지만, 실전에서는 충분히 쓸 만해진다.
유일한 문제는 여전히 나눗셈을 하고 있고, 매우 작은 값 또는 심지어 0으로 나눌 수도 있다는 점이다! 이건 교묘한 트릭으로 피할 수 있는 문제가 아니다. 행렬은 특이(singular) 할 수 있고, 그런 경우 수학은 아예 해가 없다고 말한다.
그러므로 우선 이 루틴은 optional<vector>를 반환해야 하고, 둘째로 이 나눗셈을 검사해야 한다.
optional<vector> solve(matrix const& m, vector const& v) {
// Gaussian elimitation code...
for (int i = 0; i < m.columns(); ++i) {
// Find maximum m[j][i] among all remaining rows
float M = 0.f;
for (int j = i; j < m.rows(); ++j)
M = max(M, abs(m[j][i]));
if (M < 1e-4)
return std::nullopt;
// Proceed with the algorithm
}
}
물론 이 1e-4는 아무 데서도 나온 값이 아니다. 문제를 해결하는 대신 조용히 덮어버리려고 박아 넣은 전형적인 epsilon이다.
특이에 가까운 어떤 시스템에 대해서는, 우리 알고리즘이 거의 임의적인 임계값을 바탕으로 실패했다고 보고할 것이다. 더 나쁜 점은, 우리가 계산한 이 M 값이 행렬의 어떤 본질적인 성질도 아니라는 것이다. 그저 알고리즘의 중간 과정에서 얻은 값일 뿐이다.
이 임계값은 최소한 사용자 자신이 제공하도록 해야 한다(기본 매개변수를 줄 수는 있겠지만). 하지만 그것조차 꽤 의심스럽다. 행렬이 특이인지, 또는 특이에 가까운지를 제대로 판단하는 방법은 그 조건수를 계산하거나 특이값을 살펴보는 것이지, 어떤 임의의 중간값에 임의의 임계값을 두는 것이 아니다!
내 요점은 이렇다. 행렬이 얼마나 특이한지 판단하는 것은 우리의 일이 아니라 사용자의 일이다. 우리 일은, Gauss-Jordan 소거법의 구현자로서, 얻을 수 있는 한 최선의 답을 제공하거나, 그게 불가능하면 실패를 보고하는 것이다. 우리 코드가 진짜로 실패하는 유일한 경우는 0을 0으로 나누는 경우뿐이다. 그 외에는, 비록 매우 부정확할 수는 있어도(행렬에 따라), 어쨌든 어떤 그럴듯한 답은 나온다.
따라서 내 생각에 올바른 코드는 그냥 이것이다.
// ...
if (M == 0.f)
return std::nullopt;
// ...
솔직히 말해 이것도 큰 수를 작은 수로 나눠서 무한대를 얻는 상황까지 막아 주지는 않는다. 하지만 epsilon도 그건 막아 주지 못한다.
레이 트레이서를 만들든(복셀이든 아니든), 마우스로 오브젝트 선택을 구현하든, 또는 수많은 다른 일을 하든, 광선-박스 교차 루틴이 필요해진다.
알고리즘 자체는 꽤 단순하다. 광선이 박스에 들어오는 시간(즉 광선 를 따라가는 매개변수 )과 박스를 떠나는 시간을 계산한다. 앞의 시간이 뒤의 시간보다 작으면, 그것이 교차다. 들어오는 시간은 각 좌표축에서의 진입 시간들의 최댓값이고, 떠나는 시간은 각 좌표축에서의 이탈 시간들의 최솟값이다.
이 scratchapixel 글이 이를 꽤 잘 설명한다. 코드는 대략 이런 모양이다.
void sort(float & x, float & y) {
if (x > y) swap(x, y);
}
pair<float, float> intersect(ray r, box b) {
float txmin = (b.min.x - r.origin.x) / r.direction.x;
float txmax = (b.max.x - r.origin.x) / r.direction.x;
float tymin = (b.min.y - r.origin.y) / r.direction.y;
float tymax = (b.max.y - r.origin.y) / r.direction.y;
float tzmin = (b.min.z - r.origin.z) / r.direction.z;
float tzmax = (b.max.z - r.origin.z) / r.direction.z;
sort(txmin, txmax);
sort(tymin, tymax);
sort(tzmin, tzmax);
float tmax = min(txmax, min(tymax, tzmax));
float tmin = max(txmin, min(tymin, tzmin));
return {tmin, tmax};
}
참고: sort를 약간 다시 쓰면 minss/maxss SSE 명령을 쓰도록 강제할 수 있어서, intersect 함수를 분기 없는 형태로 만들 수 있다.
짧고 읽기 좋은 함수이고, 99.99%의 경우 완벽하게 잘 동작한다. 그런데 가끔 direction.x가 0인 광선을 만나게 되고, 그러면 나눗셈이 무한대나 NaN을 만들어 낸다. 이런!
그럼 epsilon을 넣자고? 아니다. 우선 abs(direction.x) < 1e-4라면 뭘 할 건가? 결국 direction.x로 나누지 않으면서도 교차를 제대로 찾는 코드를 추가해야 하고, 어쩌면 epsilon 없이도 그냥 잘 동작할 수 있다.
코드를 분석해 보자! direction.x == 0이면 무슨 일이 일어날까? 광선은 YZ 평면에 평행하므로 실제 교차는 거기서 일어난다. 하지만 X 좌표를 그냥 무시할 수는 없다. 광선 시작점이 [b.min.x, b.max.x] 구간 밖에 있으면 교차가 없고, 그 안에 있으면 Y와 Z 좌표를 살펴봐야 한다.
정확히 direction.x == 0일 때는 어떻게 될까? A / direction.x를 계산하면, A < 0이면 , A > 0이면 를 반환한다. 여기서 A는 b.min.x - r.origin.x 또는 b.max.x - r.origin.x다. 따라서 광선 시작점 r.origin.x가 [b.min.x, b.max.x] 구간 안에 있으면 txmin == -INF이고 txmax == INF가 된다. 반대로 원점이 구간 밖에 있으면 txmin과 txmax는 둘 다 -INF이거나 둘 다 INF가 된다.
교차가 존재하고 원점이 구간 안에 있다면, tmin 계산은 사실상 txmin == -INF를 무시한다. 왜냐하면 언제나 max(-INF, x) == x이기 때문이다. 마찬가지로 tmax는 txmax == INF를 사실상 무시한다. 그러면 이 경우 코드는 요구한 대로 YZ 평면에서의 올바른 교차를 계산한다.
반대로 원점이 구간 안에 없으면, txmin == txmax == INF가 되어 tmin == INF가 되지만 tmax는 어떤 유한값이어서 교차가 없게 된다(tmin > tmax이므로). 혹은 txmin == txmax == -INF가 되어 tmax == -INF가 되지만 tmin은 유한값이어서 역시 교차가 없게 된다.
이 말은, direction.x == 0인 경우에도 이 코드가 실제로는 정확하게 동작한다 는 뜻이다! epsilon은 필요 없다. 이미 무한대를 정확히 우리가 원하는 방식으로 처리하고 있다. 정말 그럴까?
우리는 IEEE754에는 0이 두 개 있다는 사실을 잊었다. 양의 +0과 음의 -0이다. 음의 0으로 나누면 역시 무한대가 나오지만, 부호도 뒤집힌다! 그래서 txmin == INF이고 txmax == -INF가 되어, 실제로는 교차가 있는데도 교차가 없다고 보고하는 경우가 생길 수 있다.
다행히도 sort(txmin, txmax) 호출이 이 문제를 해결한다. 이 둘을 바꿔서 txmin == -INF이고 txmax == INF인 경우로 만들어 주기 때문이다. 따라서 다시 한 번, 우리의 코드는 이미 이 경우를 올바르게 처리하고 있다.
위의 Scratch-a-pixel 링크는 이 문제에 대한 다른 해법도 다루는데, 그쪽은 분기가 아주 많이 들어가서 branchless는 아니지만 더 일찍 탈출할 수 있다.
사실 아직 고치지 못한 경우가 하나 더 있다. r.direction.x == 0이고 동시에 r.origin.x == b.min.x이면, 0을 0으로 나누게 되어 NaN이 나온다. 안타깝게도 이건 코드가 자동으로 해결해 주지 않는다. NaN은 어떤 것과 비교해도 항상 false를 반환하므로, sort, min, max의 구현 방식에 따라 tmin이나 tmax가 NaN이 될 수 있고, 사실상 교차가 없다는 뜻이 된다. 어떤 용도에서는 이것도 괜찮을지 모르지만, 이것을 진짜 교차로 보고하고 싶다면 간단한 해결책이 떠오르지 않는다. r.direction.x == 0 && r.origin.x == b.min.x를 검사해서 이 경우 tmin = -INF로 두는 방법은 있지만, 그러면 코드는 더 이상 branchless가 아니다. 그래도 실제로는 이런 경우가 극히 드물기 때문에 그냥 넘어갈 수도 있다 :)
대부분의 2차원 볼록 껍질 알고리즘은 결국 세 점 에 대해 마지막 점 가 광선 의 왼쪽에 있는지를 검사하는 문제로 귀결된다. 같은 말로, 에서 로 움직였다가 다시 에서 로 방향을 틀 때 왼쪽으로 도는지 오른쪽으로 도는지를 묻는 것이다. 그래서 이 판정은 흔히 left_turn이라고 불린다.
이것은 두 벡터 와 의 상대적 방향을 검사하는 문제로 바뀌며, 간단한 행렬식으로 계산할 수 있다.
float det(vec2 a, vec2 b) {
return a.x * b.y - a.y * b.x;
}
bool left_turn(vec2 a, vec2 b, vec2 c) {
return det(b - a, c - a) > 0.f;
}
멋진 점은, 부동소수점 관련 문제를 전부 이 판정 함수 안에 가둘 수 있다는 것이다. 볼록 껍질 알고리즘 자체는 좌표에 관심이 없고, 그저 left_turn 판정을 어떤 블랙박스로 제공받기만 하면 추상적인 알고리즘으로 작동할 수 있다.
하지만 계산기하 알고리즘의 대부분은 경계 상황에 극도로 민감하고, 그런 경계 상황이 이 알고리즘의 특정 불변식을 깨뜨리면 알고리즘 전체가 완전히 망가질 수 있다.
left_turn 판정에서 그런 불변식 중 하나는 점들을 순환시켜도 결과가 바뀌지 않는다는 것이다.
left_turn(A,B,C) == left_turn(B,C,A) == left_turn(C,A,B)
거의 일직선 위에 있는 세 점을 찾는 것은 아주 쉽고, 이런 경우 이 불변식이 깨진다. 이 현상은 너무 흔해서, 최소한 1000개 정도의 점을 다뤄야 하는 볼록 껍질 알고리즘이라면 이를 반드시 고려해야 한다.
우리 해법은 이미 알고 있지 않은가? — 그냥 여기에도 epsilon을 붙이면 된다.
bool left_turn(vec2 a, vec2 b, vec2 c) {
return det(b - a, c - a) > 1e-4f;
}
물론 전혀 해결되지 않는다. 부동소수점 오차는 이 판정이 위의 순환 불변식을 어기게 만들기에 충분하다.
이 문제를 푸는 방법은 많다.
이 가운데 첫 번째 방법(고정 격자로 반올림)이 아마 가장 단순하고 가장 실용적이다. 하지만 그렇게 할 수 없는 경우도 있고, 입력 데이터를 교란하지 않은 채 그대로 다뤄야 하는 경우도 있으며, 그때는 더 복잡한 방법이 필요하다. 어쨌든 epsilon은 당신을 구해 주지 못한다. 결국 알고리즘은 헛소리를 내거나, 그냥 크래시할 것이다.
어쨌든 이런 판정은 두 값이 아니라 세 값, 즉 left, right, collinear(혹은 1990년대 감성으로 -1, 0, 1)를 반환한다고 생각하는 편이 더 낫다. 이것은 삼방 비교와 매우 비슷하다.
마지막 두 사례는 epsilon이 실제로 좋은 해법이라고 내가 생각하는 경우들이다!
기하 시각화 라이브러리를 작성하고 있다고 하자. 혹은 일종의 지도 엔진을 만들고 있을 수도 있다. 사용자는 우리가 통제할 수 없는 폴리라인을 제공한다. 우리의 일은 그 위에서 계산을 수행하고, 어쩌면 그것을 렌더링하는 것이다.
그런데 폴리라인 렌더링은 꽤 복잡한 주제다. 특히 우리는 보통 다양한 종류의 join과 cap을 사용해서 폴리라인을 삼각형 메쉬로 분해한다. 그리고 이를 계산하려면 특정 선분의 법선이나, 연속된 두 선분 사이의 각도 같은 것이 필요하다. 그런데 물론 연속된 두 점이 같거나 너무 가까우면 모든 것이 깨진다. 법선과 각도 계산이 심하게 틀어지고, 아름다운 선 위에 보기 흉한 아티팩트가 생긴다.
물론 연속된 같은 점은 허용하지 않는다고 선언할 수도 있고, 어떤 라이브러리나 알고리즘에는 그것이 합리적인 요구사항일 수 있다. 하지만 사용자 대상 엔진이나 프레임워크에서는 그다지 적절하지 않다. 어쨌든 같은 점을 걸러내는 것은 쉽다 — 말 그대로 std::unique 한 번 호출하면 된다. 하지만 같은 점은 아니되 너무 가까운 점들은 어떻게 할 것인가?
수학적으로는, 점들이 너무 가까운 것은 문제가 아니다. 출력에는 매우 가느다란 삼각형이 생기겠지만, 결과 메시의 위상은 깨지지 않는다. 하지만 실제로는 그런 가느다란 삼각형을 정확히 계산할 만큼의 부동소수점 정밀도가 부족하다(그리고 그것이 바로 보기 흉한 아티팩트를 만드는 원인이다). 반면 이번에는 가능한 한 정확하게 무언가를 계산하는 것이 목표가 아니라, 데이터를 시각화하는 것이 목표다. 실제로 아무도 이런 가느다란 삼각형을 보지 못한다. 그냥 너무 얇기 때문이다.
즉, 그냥 렌더링하지 않아도 괜찮다. 다시 말해, 애초에 생성하지 않아도 괜찮다. 다시 말해, 너무 가까운 입력 점을 그냥 걸러내도 괜찮다. 그런데 얼마나 가까워야 “너무 가까운” 걸까?
다시 한 번 말하지만, 임의의 epsilon을 박아 넣는 것은 대개 나쁜 해법이다. 다만 이 경우에는 아주 나쁘다기보다 그냥 평범하게 별로인 정도다. 이 용도에 맞는 좋은 epsilon을 찾는 방법은 있다.
마지막 방법을 택하지 않는 이상, 올바른 값은 1e-4나 1e-6가 아닐 것이라고 장담할 수 있다.
아마 가장 자명한 경우일 것이다. 내가 수학 라이브러리를 작성한다면, 폭넓은 테스트 커버리지를 원한다. 테스트는 어떤 함수가 기대한 값을 반환하는지 확인해야 하므로, 결과를 기준값과 비교해야 한다. 하지만 부동소수점 오차 때문에, 결과는 거의 언제나 약간 다를 것이다.
이제 여기에는 두 가지 방법이 있다. 첫 번째는 정확한 IEEE754 준수를 가정하고, 테스트를 동등 비교로 작성하는 것이다. 반올림이 있다 해도, 올바르게 설정된 CPU 플래그 아래에서는 모든 표준 준수 하드웨어에서 같아야 하므로, == 비교를 안전하게 쓸 수 있다. 많은 경우 입력 데이터를 잘 고르면 아예 반올림이 전혀 일어나지 않게 만들 수도 있다!
많은 함수(예: 벡터 덧셈, 뺄셈, 곱셈)에서는 그냥 정수 좌표를 쓸 수 있고(이들 중 많은 값은 float로 정확히 표현 가능하다), 혹은 2의 거듭제곱으로 나눈 정수를 쓸 수도 있다. 어떤 경우에는 이것만으로는 충분하지 않다. 예를 들어 벡터 길이를 계산할 때는 좌표가 정수여도 결과는 무리수일 수 있다.
그런 경우에도 약간 요령을 부릴 수 있다. 피타고라스 수를 이용하면, 길이가 부동소수점으로 정확히 표현된다는 것을 아는 테스트 케이스를 만들 수 있다. 예를 들어
length(vec2(3.f, 4.f)) == 5.f
length(vec2(5.f, 12.f)) == 13.f
length(vec2(0.4375f, 1.5f)) == 1.5625f
하지만 이것은 꽤 번거롭다. 우리 코드 전체를 위해 점점 더 복잡한 특수 테스트를 직접 만들어야 하기 때문이다. 그럼 두 번째 방법은 뭘까?
epsilon을 써라! 꽤 작은 epsilon을 쓰는 것이 좋다. 예를 들어 FLT_EPSILON 같은 것 말이다. 하지만 이 경우에는 솔직히 좋은 해법이다. 대체로 수학 라이브러리를 epsilon 비교가 잡아내지 못할 방식으로 망가뜨리기는 매우 어렵다.
다만 어떤 함수가 추가 보장(예: length는 정확히 0인 벡터에 대해서만 0을 반환한다든가, left_turn은 인자를 순환시켜도 불변이라든가)을 제공한다면, 그런 보장에 대해서는 별도의 특화된 테스트를 작성하는 것이 가장 좋다.
그렇다면 epsilon은 좋은가 나쁜가? 대체로는 나쁘지만, 가끔은 괜찮다. 이 글을 읽었다고 해서 당신의 epsilon 범벅 코드가 갑자기 전부 망가질까? 아마 그렇지는 않을 것이다.
현실의 많은 경우에는 이 모든 것이 별로 중요하지 않다. 범용 고품질 수학 라이브러리를 작성한다면 이런 것을 신경 써야 할 의무가 있다. 말랑한 플랫폼 게임용 장난감 물리 엔진을 대충 만들고 있다면, epsilon으로도 아마 충분히 잘 굴러갈 것이다. 그리고 깨지면 그냥 1e-4를 2e-4로 바꾸고 계속 가면 된다.
어떤 공학 문제에 대한 진짜 답은, 어떤 집단이나 인터넷의 아무 사람 말만 따르는 것이 아니라, 자기 머리로 생각하는 것이다. 충격적이지만 말이다. 그래도 이 글에서 뭔가 하나쯤은 배웠기를 바란다!
혹시 내 글이 마음에 든다면, 내가 하는 다른 작업도 후원해 주는 것을 고려해 달라!
예를 들어 내 YouTube devlogs를 봐 달라. 이런 영상도 있다: