[c언어] 수열의 부분 합(Prefix Sum) 구하기 - 어떤 방법이 더 빠르고 정확할까

수열의 부분합 구하기! 중고등학교 수학시간에 많이 했던 일이죠. 어떤 수열 X가 a, b, c, d ... 와 같은 식으로 있다면 부분합(Prefix Sum, 혹은 Scan) S는 다음과 같이 계산됩니다.

S₁: a
S₂: a + b
S₃: a + b + c
S₄: a + b + c + d

수학시간에서는 이 부분합 수열의 일반항을 구하는 일을 주로 했지만, 컴퓨터 과학에서는 이 수열의 각 항을 빠르고 효율적으로 (또 정확히) 계산하는 방법에 대해 논하게 됩니다. 이 부분합을 구해서 어디에 쓰나 싶지만, 의외로 여러 분야에서 널리 쓰이고 있습니다. 

대표적인 사례로 누적분포(cumulative distribution)을 구하는 작업이 있겠습니다. 이는 특정 임의 분포에서 표본을 추출하는데 자주 사용됩니다. 예를 들어 A, B, C, D가 각각 뽑힐 확률이 4.0 : 1.0 : 0.4 : 3.7인 분포가 있다고 합시다. 여기서 알파벳 몇개를 뽑아내는 작업을 컴퓨터로 수행하려면 어떻게 해야할까요? 답부터 말하자면, 4.0, 1.0, 0.4, 3.7의 누적분포를 구한 뒤, 0~최댓값 사이의 난수를 생성하여 이진 탐색을 실시하는 것입니다. 즉, 누적 분포를 구하면 4.0, 5.0, 5.4, 9.1이 됩니다. 이제 여기서 0~9.1 사이의 난수를 뽑았는데, 4.158가 나왔다고 칩시다. 그러면 4.158은 4.0과 5.0 사이에 있으므로, 알파벳 B를 생성하는 식입니다. 누적 분포의 각 값들은 항상 증가하는 순서대로 위치하므로 이진 탐색을 이용하면 O(log N) 시간에 표본을 추출할 수 있습니다. 

이외에도 카운팅 정렬이나 List Ranking 등의 작업에도 쓰이고 있구요. 단순한 작업이지만 여러 곳에 쓰이고 있으니 이를 빠르고 정확하게 계산하는게 중요한 일이겠지요? 이번 포스팅에서는 c언어를 이용해 다양한 부분합 계산법을 구현해보고 속도와 정확도를 비교해보는 일을 해보고자 합니다. 여기에 사용된 코드 및 실험 결과는 https://github.com/bab2min/prefix-sum 에 정리해두었으니 궁금하신 분들은 직접 예제를 받아 돌려보셔도 좋습니다.

제일 간단한 방법

차례대로 더하면 되는 것이기 때문에 그냥 그대로 구현해봅시다.

void prefix_sum_simple(float* arr, int n) { float sum = 0; for(int i = 0; i < n; ++i) { sum += arr[i]; arr[i] = sum; } }

설명이 필요없을 정도로 간단합니다. 반복문을 돌면서 arr[0]부터 arr[n - 1]까지 차례로 더하고, 그 값을 다시 arr[0]부터 arr[n - 1]까지 써넣는 단순한 코드입니다. 이 코드가 우리의 베이스라인이 될겁니다. 너무 단순해 보여서 개선할 부분이 안 보이는데 과연 코드를 더 빠르게 할 수 있을까요? 

루프 풀기(Loop Unrolling)

사실 위의 코드는 float 덧셈 n번뿐만 아니라 int 덧셈도 n번 수행해야합니다. 반복 횟수를 카운트하기 위해서 i를 1씩 증가시키기 때문이죠. 다음과 같이 반복문을 풀어버리면 어떨까요?

void prefix_sum_unroll4(float* arr, int n) { int n4 = n & ~3; float sum = 0; for (int i = 0; i < n4; i += 4) { arr[i] += sum; arr[i + 1] += arr[i]; arr[i + 2] += arr[i + 1]; arr[i + 3] += arr[i + 2]; sum = arr[i + 3]; } for (int i = n4; i < n; ++i) { sum += arr[i]; arr[i] = sum; } }

먼저 float 덧셈 4번이 루프 한 번에 들어가도록 묶어보았습니다. 4의 배수인 지점까지는 첫번째 루프로 반복하고, 4로 나눠떨어지지 않는 나머지는 두번째 루프에서 처리합니다. 그럼 대략적으로 사용되는 연산량은 float 덧셈 n번에 int덧셈 n/4번이 되겠군요. int 덧셈을 꽤 절약했습니다. 이 방법은 실험결과 베이스라인보다 약 9%정도 빠른 것으로 나타났습니다.

4 말고 8이나 16의 배수로 푸는것도 가능하겠죠? 각각 실험해보면 8개씩 풀어낸 경우는 약 7%, 16개씩 풀어낸 경우는 약 8% 정도 베이스라인보다 빨랐습니다. 오히려 너무 많이 풀면 속도 향상에는 도움이 안된다는 걸 알 수 있습니다.

순서를 바꾸고 덧셈을 늘려보기

여기서부터는 직관에서 벗어나는 결과가 나오니 집중하세요! 다음은 연산 순서를 바꾸는 대신 덧셈을 한 번 늘린 코드입니다.

void prefix_sum_unroll4_reorder(float* arr, int n) { int n4 = n & ~3; float sum= 0; for (int i = 0; i < n4; i += 4) { arr[i] += sum; arr[i + 1] += arr[i]; arr[i + 3] += arr[i + 2]; arr[i + 2] += arr[i + 1]; arr[i + 3] += arr[i + 1]; sum = arr[i + 3]; } for (int i = n4; i < n; ++i) { sum += arr[i]; arr[i] = sum; } }

a[0] = 1, a[1] = 2, a[2] = 3, a[3] = 4라고 해봅시다. a[1] += a[0]을 하면 a[1] = 1+2가 됩니다. 그리고 a[3] += a[2]를 하면 a[3] = 3+4가 되겠죠. 그 다음 a[2] += a[1], a[3] += a[1]을 수행해 a[2] = 1+2+3, a[3] = 1+2+3+4가 되도록 하는 식입니다. 이렇게 하면 float 덧셈은 5번으로 늘어납니다. 그런데 속도는 베이스라인보다 36% 빨라집니다. 이상한 일이죠? 왜 이런 일이 발생하는지 이해하려면 현대 CPU 구조에 대해서 살펴보아야 합니다.

현대 CPU는 다양하고 복잡한 연산들을 빠르게 수행하기 위해서 명령어 파이프라인이라는 기법을 사용합니다. 이는 커다란 작업 하나를 단순한 작업으로 쪼개 각 작업들을 동시에 수행하는 것입니다. CPU 구조를 직접 설명하려면 어렵고 지루할테니, CPU가 하는 일과 매우 비슷한 '김밥싸기'에 비유해 파이프라인에 대해 설명해보도록 하겠습니다.

명령어 파이프라인

김밥집에서는 주문이 들어오는 대로 김밥을 싸야합니다. 재료가 무한히 준비되어 있다고 가정하고, 김밥을 싸는 과정은 다음과 같이 총 4단계로 구성됩니다. 

1) 김을 펼쳐놓고
2) 밥을 고르게 펴바른뒤
3) 채소와 햄을 넣고
4) 이쁘게 말아준다. 

불행히도 네 가지 작업을 다 잘하는 인재는 많지 않기 때문에 이런 사람을 고용하려면 돈이 많이 듭니다. 그래서 김밥집 사장은 네 가지 작업을 모두 수행하는 능력자 대신, 한 가지 작업만 잘할 수 있는 사람을 각각 4명 고용합니다. 1)~4)번 작업을 수행하는 각 일꾼들이 차례로 늘어선 다음, 자기 일이 끝나면 옆사람에게 일을 넘겨주는 식으로 분업하면 1명이 한 가지 작업만 하면서도 빠르게 김밥을 쌀 수 있겠지요. 이렇게 분업화하는 것을 파이프라이닝이라 합니다. CPU가 수행하는 연산들도 생각보다 복잡하기 때문에 한 장치가 모든 연산을 다 수행하도록 하지 않고, 단순한 장치 여러 개를 만들어 차례로 일을 하도록 설계해놓았습니다.

(CPU 파이프라인의 개념: 4개의 stage가 김밥 싸는 4단계라고 보면됩니다. 윗쪽의 0~8 cycle은 각 초라고 보면 되구요)

4명의 일꾼이 1초에 하나씩 작업을 처리한다고 가정해봅시다. 그러면 주문이 들어가서 모든 파이프라인을 통과해 김밥이 나오기까지는 4초가 걸릴겁니다. 주문이 매초에 하나씩 들어온다면 첫 주문 결과는 4초 뒤에 나오겠지만, 그 다음부터는 쉬지 않고 1초마다 주문 결과가 나오겠죠. 왜냐면 각 일꾼들이 1초에 하나씩 김을 펼치고, 밥을 펴바르고, 채소와 햄을 넣고, 말아줄테니깐요. 이와 같이 주문이 들어간 뒤 결과가 나오기 까지의 시간을 latency라 하며, 이론적으로 결과 하나씩이 나오는데 걸리는 평균 시간을 throughput이라 합니다. 따라서 이 김밥집의 latency = 4이고, throughput = 1이 됩니다.

주문이 매 초 하나씩 계속 들어가면 이 김밥집은 최고 효율로 일할 수 있습니다. 그런데 그렇지 못한 경우도 있겠죠. 입맛이 까다로운 손님이 왔습니다. 이 손님은 먼저 김밥 1개를 주문한뒤, 이 김밥의 맛을 보고 두번째로 주문할 김밥의 종류를 결정한다고 합니다. 그리고 또 두번째 김밥의 맛을 본뒤 세번째 김밥 종류를 결정하구요. 이 경우 김밥집의 효율은 어떻게 될까요? 먼저 첫 주문이 들어간 뒤 4초 뒤에 김밥이 1개 나올거고, 그걸 맛본뒤 주문이 들어갈테니 두번째 주문은 4초뒤에 들어갈 겁니다. 따라서 두번째 김밥은 8초뒤에 나오겠네요. 마찬가지로 세번째 김밥이 나오려면 12초가 걸리겠구요. 이 김밥집 일꾼들은 1초에 한개씩 김밥을 쌀 수 있는 능력이 있음에도 4초에 1개 밖에 못 싸는 상황이 발생됩니다. 왜냐면 까다로운 손님이 생산된 김밥 결과를 보고 주문을 넣기 때문이지요.

이제 CPU 이야기로 돌아와봅시다. 매 시간 들어가는 주문은 CPU로 입력되는 명령어, 그리고 실제로 만들어지는 김밥은 그 명령어 연산의 결과라고 보시면 됩니다. 모든 명령어는 latency가 있습니다. 간단한 작업이라면 latency가 짧겠지만, 복잡한 작업이라면 매우 길어질 수도 있죠. latency가 길더라도 throughput은 보통 0.5~1 정도되기 때문에 명령어들이 차례대로 들어온다면 CPU의 각 일꾼들은 쉼없이 일하여 최고의 효율을 낼 수 있습니다. 그런데 김밥집의 까다로운 손님처럼, 앞 명령어의 결과에 의존하여 다음 명령어를 결정해야 하는 코드가 있다면 그 효율이 떨어질 수 밖에 없습니다. 이런 코드로 대표적인 것이 if, for 등의 분기문, 그리고 바로 이전의 연산 결과를 이용해야하는 연산입니다.

일반적으로 현대 CPU의 float 덧셈 연산의 latency는 3정도입니다(throughput은 1입니다). 이 말인즉슨, float 덧셈 명령어가 CPU에 입력되고 3사이클이 지나야만 그 결과가 나온다는 것이죠. 덧셈 연산이 차례로 들어간다면 1사이클에 1개꼴로 결과를 낼 수 있지만, 그렇지 않은 경우라면 3사이클에 1개밖에 결과를 낼 수 밖에 없습니다.

float sum = 0; sum += arr[0]; arr[0] = sum; sum += arr[1]; arr[1] = sum; sum += arr[2]; arr[2] = sum; sum += arr[3]; arr[3] = sum; // ...

처음에 나왔던 단순한 반복문을 풀어보면 위와 같습니다. sum += arr[0]을 수행한 결과가 나와야지만 arr[0] = sum을 수행하고, sum += arr[1]을 수행할 수 있습니다. 마찬가지로 sum += arr[1]을 수행한 결과가 나와야지만 arr[2] = sum을 수행하고 sum += arr[2]를 수행할 수 있습니다. 이렇게 앞 연산 결과로 변경된 값을 뒤에서 사용해야하는 경우 데이터 의존성(data dependency)이 있다고 합니다. 따라서 위 코드를 실행할때는 CPU는 매 float 덧셈마다 앞 연산의 결과가 나올때까지 3사이클을 기다려야합니다.

arr[i] += sum; // W1 arr[i + 1] += arr[i]; // R1, W2 arr[i + 3] += arr[i + 2]; //W3 arr[i + 2] += arr[i + 1]; // R2 arr[i + 3] += arr[i + 1]; //R2, R3 sum = arr[i + 3];

반면 위의 코드를 봅시다. 서로 데이터 의존성이 있는 부분을 R과 W로 표시해봤습니다. (Rx는 Wx가 수행된 다음에야 수행할 수 있다는 뜻입니다.). 이 경우 덧셈 연산이 한 번 늘었지만, 의존성이 없어서 바로 수행될 수 있는 연산들이 꽤 있습니다. arr[i+3] += arr[i+2]는 arr[i] += sum, arr[i + 1] += arr[i]과 동시에 수행될 수 있습니다. 마찬가지로 arr[i + 2] += arr[i + 1]와 arr[i + 3] += arr[i + 1] 도 동시에 수행될 수 있습니다. 최신 CPU는 비순차적 실행(Out-of-order Execution)을 지원하기 때문에 순서를 바꿔서 미리 수행하거나 동시에 수행해도 괜찮은 연산들은 미리(혹은 동시에) 수행합니다. 결과적으로 덧셈 횟수가 1회 늘었음에도 낭비되는 사이클은 줄어들어 연산속도가 빨라지는 것이죠.

(여기서 동시에 수행한다는 것은 멀티스레딩을 의미하는 게 아닙니다. 하나의 코어 내에도 동일한 연산을 수행하는 장치가 여러 개 있는 경우도 있습니다. 이런 경우 하나의 코어가 여러 개의 명령어를 비순차적으로 처리하여 1보다 빠른 throughput을 내는 것, 즉 한 사이클에 둘 이상의 결과를 내는 것을 가리킵니다.)

더 공격적으로 의존성을 없애보기

덧셈을 늘려도 속도를 빠르게 할 수 있다는 것을 알았으니, 좀 더 공격적으로 데이터 의존성을 해소해보는 방향으로 가봅시다.

void prefix_sum_unroll4_shift(float* arr, int n) { int n4 = n & ~3; float sum = 0; for (int i = 0; i < n4; i += 4) { // 첫번째 덧셈 arr[i + 3] += arr[i + 2]; arr[i + 2] += arr[i + 1]; arr[i + 1] += arr[i]; // 두번째 덧셈 arr[i + 3] += arr[i + 1]; arr[i + 2] += arr[i]; // 마지막 덧셈 arr[i] += sum; arr[i + 1] += sum; arr[i + 2] += sum; arr[i + 3] += sum; sum = arr[i + 3]; } for (int i = n4; i < n; ++i) { arr[i] += sum; sum = arr[i]; } }

이번엔 float 덧셈이 총 9개입니다. 베이스라인보다 덧셈 횟수가 2배 이상 더 많아졌네요. 그런데 속도는 오히려 약 57% 정도 빨라집니다. 첫번째 덧셈 그룹에 있는 3개의 덧셈은 서로 데이터 의존성이 없습니다. 마찬가지로 두번째 그룹에 있는 2개의 덧셈도 서로 의존성이 없구요. 마지막에 있는 4개의 덧셈도 의존성이 없습니다. 따라서 각 그룹의 연산들은 모두 동시에 수행될 수 있습니다.

SIMD를 활용하기

최신 CPU에서는 여러 개의 float를 한번에 묶어서 연산할 수 있는 SIMD 명령어가 포함되어 있습니다. 이를 활용하면 더욱더 속도를 높일 수 있습니다.

inline __m128 scan_sse(__m128 x) { // x: [0, 1, 2, 3] // 첫번째 덧셈 x = _mm_add_ps(x, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(x), 4))); // x: [0, 0+1, 1+2, 2+3] // 두번째 덧셈 x = _mm_add_ps(x, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(x), 8))); // x: [0, 0+1, 0+1+2, 0+1+2+3] return x; } inline void prefix_sum_sse(float* arr, int n) { int n4 = n & ~3; __m128 sum = _mm_setzero_ps(); for (int i = 0; i < n4; i += 4) { __m128 x = _mm_load_ps(&arr[i]); __m128 out = scan_sse(x); // 마지막 덧셈 out = _mm_add_ps(out, sum); _mm_store_ps(&arr[i], out); sum = _mm_shuffle_ps(out, out, _MM_SHUFFLE(3, 3, 3, 3)); } if (!n4) n4 = 1; for (size_t i = n4; i < n; ++i) { arr[i] += arr[i - 1]; } }

2001년에 발표된 SSE2는 float 4개(총 128bit)를 동시에 처리할 수 있습니다. 위의 코드는 prefix_sum_unroll4_shift 함수의 연산을 SSE2 명령어를 이용해 대체한 것입니다. prefix_sum_unroll4_shift 의 첫번째, 두번째, 마지막 덧셈과 동일한 연산을 수행하는 코드에 똑같이 주석을 달아놓았습니다. SSE2가 지원되는 CPU(2020년 현재 현역인 모든 Intel, AMD CPU에서 지원됨)에서 이 코드는 베이스라인보다 170% 빠르게 실행됩니다.

와 부분합 구하는게 간단한 코드인줄만 알았는데, 생각보다 최적화할 거리가 아주 많습니다. 요즘 컴파일러가 강력하다고 다들 말하는데, 왜 컴파일러는 자동으로 이렇게 최적화를 해주지 않을까요? 컴파일러가 멍청하기 때문일까요? 컴파일러가 prefix sum 코드를 자동으로 최적화하지 못한 이유는 포스팅의 끝부분에서 설명하도록 하겠습니다.

정확도는 어떨까

지금까지 속도를 끌어올리는 방법에 대해 살펴봤는데요, 여기서부터는 연산의 정확도에 대해 고민해봅시다.  사실 수학적으로 따지면 덧셈을 어떻게 하든지 그 결과는 같습니다만, 컴퓨터의 부동소수점 덧셈은 그렇지 않습니다. 부동소수점은 유효자리수를 정하여 근사치를 저장하는 방식이기 때문입니다. float의 경우 24비트, double의 경우 53비트만을 정확하게 저장할 수 있죠. 근사값을 저장하기 때문에 연산 순서에 따라 그 결과가 달라질 수 있습니다. 

이해하기 쉽게 십진법으로 예를 들자면, 유효자리가 3자리인 상황에서 101 + 1.23 + 6.78를 계산한다고 해봅시다. 먼저 (101 + 1.43) + 6.48 순으로 계산해보면, 101 + 1.43 = 102.43이지만 우리의 암기력은 유효자리 3개 밖에 안되므로, 소수점 이하 .43은 반올림해버릴수 밖에 없습니다. 따라서 101 + 1.43 = 102가 됩니다. 이 상황에서 6.48을 더하면 102 + 6.48 = 108.48이지만 또 소수점 아래를 반올림해버리면 108을 얻게 됩니다. 반면 101 + (1.43 + 6.48)순으로 계산하면, 1.43 + 6.48 = 7.91로 계산이 되고(유효숫자 3자리 안에 들어오기 때문에 손실되는 부분이 없습니다), 101 + 7.91 = 108.91에서 소수점 아래를 반올림하면 109를 얻게 됩니다. 분명 같은 덧셈이지만 연산 순서에 따라 결과가 달라졌습니다. (그리고 후자의 경우가 참값에 더 가까운 결과를 얻었다는 것도 알 수 있습니다.) 

정리하자면, 부동소수점 방식은 한정된 비트로 유효숫자를 표시하기 때문에 연산 과정에서 오차가 발생하고 이 때문에 연산 순서를 바꾸면 결과가 달라질 수 있다는 것입니다. 특히 연산의 정밀도를 많이 잃어버리는 경우는 큰 값과 작은 값을 더하는 경우입니다. 이 경우 큰 값의 자리수에 맞춰서 계산을 해야하기 때문에 작은 값의 아래자리들은 대부분 버려지게 됩니다. (101 + 1.43과 1.43 + 6.48의 경우를 비교해보시면 이해하기 쉽습니다.)

이론은 충분히 습득한 거 같으니, 실제로 각 방법이 어느 정도 오차를 낼지 상상해봅시다. 먼저 베이스라인인 순서대로 합쳐나가는 방법을 따져봅시다.

void prefix_sum_simple(float* arr, int n) { float sum = 0; for(int i = 0; i < n; ++i) { sum += arr[i]; arr[i] = sum; } }

arr값이 전부 0 이상이라고 가정하면 sum은 뒤로 갈수록 점점 커질겁니다. 따라서 뒷쪽의 sum += arr[i]는 큰 값에 작은 값을 더하는 모양이 되어서 정확도 손실이 클 겁니다. 반면 arr이 양수와 음수가 고르게 섞여 있다면 뒤로 가도 sum이 대단히 크게 증가하거나 감소하지는 않을 것이므로 상대적으로 정확도 손실이 적을 거구요. 다음 함수의 오차가 어쩔지도 추정해봅시다.

void prefix_sum_unroll4_shift(float* arr, int n) { int n4 = n & ~3; float sum = 0; for (int i = 0; i < n4; i += 4) { // 첫번째 덧셈 arr[i + 3] += arr[i + 2]; arr[i + 2] += arr[i + 1]; arr[i + 1] += arr[i]; // 두번째 덧셈 arr[i + 3] += arr[i + 1]; arr[i + 2] += arr[i]; // 마지막 덧셈 arr[i] += sum; arr[i + 1] += sum; arr[i + 2] += sum; arr[i + 3] += sum; sum = arr[i + 3]; } for (int i = n4; i < n; ++i) { arr[i] += sum; sum = arr[i]; } }

arr이 전부 0 이상인 경우만 따져보자면, 마찬가지로 sum은 뒤로 갈수록 점점 커집니다. 다만 sum을 arr에 합하는 순서가 조금 다릅니다. arr[i+3]을 생각해보면 arr[i] + arr[i + 1] + arr[i + 2] + arr[i + 3]을 다 합한 뒤에 마지막으로 sum을 더하게 됩니다. sum이 매우 큰 경우라면 sum에 arr[i]~arr[i+3]를 각각 더하는 경우 아래 자리를 대부분 잃어버릴겁니다. 반면 arr[i]~arr[i+3]을 먼저 더하면 아래 자리를 보존한 상태에서 덧셈이 가능합니다. 그 다음 sum을 더해서 아래 자리수를 잃어버린다 하더라도 오차가 상대적으로 적겠죠.

정말로 그런지 코드를 돌려서 결과를 확인해봅시다. 백만개 float 배열을 만들어서 0~1사이의 수를 임의로 채워넣읍시다. 그리고 double 타입으로 합계를 구해 참값을 계산해두고, 위에서 소개한 방법들로 prefix sum을 구해 참값과의 오차를 계산해보도록 하지요.

	속도 상승	평균 절대 오차(MAE)
simple	1.00	2.037
unroll4	1.09	2.037
unroll8	1.07	2.037
unroll16	1.08	2.037
unroll4_reorder	1.36	0.768
unroll8_reorder	1.44	1.160
unroll16_reorder	1.52	1.198
unroll4_shift	1.57	0.683
unroll8_shift	1.23	0.344
sse	2.72	0.683
avx	2.08	0.344
kahan	0.37	0.001

표로 정리하는 김에 속도 상승분도 같이 넣었습니다. 참고로 속도는 총 5종류의 컴파일러(gcc 5.4, gcc 7.5, gcc 9.3, clang 11.0.3, msvc 19.16)로 컴파일하여 백만 float 배열의 prefix sum을 구하는 작업을 50번 반복한 평균을 가지고 계산하였습니다.

예상했던것처럼 베이스라인이 제일 큰 오차를 냈습니다(루프 풀기 역시 연산 순서가 바뀌진 않았으니 동일한 결과가 나왔구요). 반면 unroll4_shift는 오차가 0.683으로 베이스라인의 1/3 정도 밖에 안됩니다. 동일한 방법을 8개 단위로 확장한 unroll8_shift의 경우는 오차가 절반으로 더 줄어드는걸 볼 수 있습니다. 그런데 정확한 합을 구하는 제일 좋은 방법은 따로 있습니다. Kahan의 합 알고리즘이라는 방법인데요, 이는 다음과 같이 구현할 수 있습니다.

void kahan(float* arr, int n) { float sum = 0, c = 0; for (int i = 0; i < n; ++i) { float y = arr[i] - c; float t = sum + y; c = (t - sum) - y; // 수학적으로 따지면 c는 원래 0이 되어야합니다만 // 부동소수점 연산의 오차 때문에 0이 아닌 값을 가지게 됩니다. // 이 값이 sum을 구할때 발생한 오차의 근사값입니다. arr[i] = sum = t; } }

부동소수점 연산의 오차를 역이용하여, sum 계산시 발생하는 오차를 추정하여 c에 저장하는 식으로 전체 합계의 오차를 추정해나가는 알고리즘입니다. 평균 오차가 사실상 0에 가깝지만 속도는 많이 느리군요.

컴파일러는 왜 prefix_sum_simple을 자동으로 최적화하지 못했나

이쯤에서 컴파일러의 변명을 들어보면 좋을 것 같습니다. 왜 prefix_sum_simple 함수를 자동으로 최적화해서 prefix_sum_unroll4_shift 나 prefix_sum_sse처럼 바꾸지 못했는지. 그 이유는 부동소수점 float의 연산 순서를 바꾸면 그 결과가 바뀌기 때문입니다. 즉, 수학적으로 따지면 동일한 작업을 하는 코드지만, 실제로는 서로 다른 결과를 냅니다(MAE가 2가 넘을 수도 있고 0에 가까울수도 있죠). 컴파일러 최적화의 기본 원칙은 프로그램의 결과를 동일하게 유지하는 상태에서 속도를 빠르게 하거나 용량을 줄이는 것이기 때문에, 결과를 바꾸면서 최적화할 수는 없었던 겁니다. 

그러니 마냥 컴파일러를 비난할 수는 없는 상황이지요. 대신 CPU구조와 부동소수점 연산에 대해서 가벼이 생각했던 제 자신을 돌아보기로 하였습니다. 반성하는 차원에서 깃헙 레포지토리를 파서 결과를 공유했으니(https://github.com/bab2min/prefix-sum) 관심 있으신 분들은 참고하시면 되겠습니다.

참고로 이번에 알게된 prefix sum 최적화 기법은 다음 tomotopy 업데이트에 반영될 예정입니다. 배열 길이가 백 만처럼 큰 수가 아니다 보니 오차가 많이 발생하지는 않으므로 평균 오차도 안정적이고 속도도 빠른 sse 버전(혹은 sse 미지원 환경에서는 unroll4_shift 버전)를 사용할 예정입니다. 전체 연산의 대부분이 random distribution sampling인지라 prefix sum을 구할 일이 잦은데, 속도 향상이 클것으로 예상됩니다.