Fast Float Sqrt

 이미지 처리 연산 과정에서 제곱근 값을 계산해야 경우가 많이 생긴다. 이 경우 math-라이브러리의 sqrt() 함수를 사용하는 것보다 원하는 정밀도까지만 계산을 하는 근사적인 제곱근 함수를 만들어서 사용하는 것이 계산 비용 측면에서 유리할 수 있다. 

주어진 양수 $x$의 제곱근 $y = \sqrt{x}$를 구하는 것은 $F(y) = y^2 - x$의 근을 구하는 것과 같다. 프로그램적으로 근은 반복적인 방법에 의해서 구해지는데 많이 사용하는 알고리즘 중 하나가 Newton 방법이다. Newton 방법은 한 단계에서 근의 후보가 정해지면 그 값에서 F에 접하는 직선을 구한 후 그 직선이 가로축과 만나는 점을 다시 새로운 근으로 잡는 반복 과정을 일정한 에러 범위 내에 들어올 때까지 수행한다.

 

함수 $F(y)$의 근은 아래의 점화식에 의해서 구해지게 된다:

$$\tt y(n + 1) = y(n) - F(y(n)) / F'(y(n)), n = 0,1, 2,...;$$

제곱근의 경우에는 ($F(y) = y^2 - x$) 점화식은 

$$\tt y(n + 1) = (y(n) + x / y(n)) / 2;$$

의 형태로 주어진다. 그러나, 이 점화식에는 계산 부하가 상대적으로 큰 나누기 연산이 들어 있다. 그런데 $\sqrt{x} = x \times (1 / \sqrt{x})$로 표현할 수 있으므로 제곱근의 역수를 쉽게 계산하는 알고리즘을 찾으면 된다. $x$의 제곱근의 역수를 구하는 점화식은 (이 경우 $F(y) = 1 / y^2 - x$)

$$\tt y(n + 1) = y(n) * (3 - y(n) * y(n) * x) / 2;$$

이므로 전체적인 제곱근을 구하는 pseudo 코드는

float fast_sqrt(float x) {
    xhalf = 0.5F * x;
    y = initial_guess_of (1/sqrt(x));
    while (error is not small)
        // find approx inverse sqrt root of (x);
        y = y * 1.5 - y * y * y * xhalf;
    end while

    // finally return sqrt(x) = x * (1 / sqrt(x));
    return (x * y);
}

보통 iteration 횟수는 고정이 되어 있다. 남은 문제는 초기값 추정뿐이다. Newton 방법에서 잘 추정된 초기값은 몇 번의 반복으로도 원하는 정밀도를 얻을 수 있게 하지만 잘못된 추정은 수렴 속도를 더디게 한다. 또한 초기값의 추정 과정도 매우 간단해야 하는데, 좋은 추정을 하기 위해서는 IEEE에서 규정하는 double/float 포맷에 대한 정보와 수치해석 지식이 필요하다.

 

IEEE의 32비트 float 상수의 비트 구성은 $\tt 0 \sim 22$번째 비트까지는 소수점 아래의 유효숫자(Mantissa: $\tt M$)를 표기하고, $\tt 23$번째부터 $\tt 30$번째까지는 지수 부분(Exponent: $\tt E$)을 표시하는데 $\tt 127$만큼 더해진 값이다(+-지수를 표현하기 위한 것으로 실제 지수는 $\tt e = E - 127$). 마지막 $\tt 31$비트는 부호를 결정한다.

임의의 양의 float 상수 $\tt x >0$가 

$$\tt x = (-1)^0\times 2^{(b_{30}b_{29}....b_{23})_2 -127} \times (1.b_{22} b_{21}...b_{0})_2\\= 2^{E-127} \left( 1+ \sum_{i=1}^{23} b_{23-i} 2^{-i}\right)\\= (1 + M) 2^{E-127} = (1 + M) 2^{e};$$

로 쓰여지면, 정수로 변환된 비트 데이터는 $\tt ix = *(int~*)\&x = (E + M) 2^{23}$의 값을 갖는다. $\tt x$ 제곱근의 역수는

$$\tt \frac{1}{\sqrt{x}} = \frac{1}{\sqrt{1 + M}} 2^{- e / 2} ;$$

로 표현되므로 비트 데이터에서 지수부는 $\tt (- e / 2) + 127 = 190 - E / 2$로 표현된다. 

 

이제 필요한 것은 잘 선택된 $\tt 1/\sqrt{x}$의 초기값을 찾는 것으로, 첫 번째로 시도할 근사는 유효숫자 부분 $(\tt 1/\sqrt{1 + M})$을 제거하고, 단순히 지수 부분($\tt 2^{- e / 2}$)만 살리면 된다. 정수 표현으로 바꾸면 $\tt 190$이 $\tt 23$에서 $\tt 30$번째 비트를 채우도록 $\tt 23$을 시프트 하고, $\tt - E / 2$ 부분은 원래 float 상수의 정수 표현을 나누면 유효숫자 부분에 해당하는 정도만 에러가 발생한다:

   정수 표현 = (190 - E / 2) << 23;

                  = (190 << 23) - (x의 정수형 변환 & (0xFF << 23)) >> 1;

                  ~ (190 << 23) - (x의 정수형 변환) >> 1;    // 유효 숫자 / 2 만큼 에러 발생;

                  = 0x5F000000 - (*(int *)&x) >> 1;

유효숫자 부분은 좀 더 높은 초기 정밀도를 갖도록 근사를 할 수 있는데 결과만 쓰면 $\tt 0x5F000000$ 대신에 $\tt 0x5F3759DF$를 이용하면 된다.

//빠른 float 제곱근의 역수;
float fast_invsqrt(float x) {
    float xhalf = 0.5F * x;
    // initial guess of 1 / sqrt(x);
    int i = 0x5F3759DF - (*(int *)&x) >> 1;
    x = *(float *)&i;    
    //iterate
    x = x * (1.5F - x * x * xhalf);
    x = x * (1.5F - x * x * xhalf);
    // end of calculation of 1 / sqrt(x);
    return (x);
};

// 빠른 float 제곱근;
float fast_sqrt(float x) {
    float xhalf = 0.5F * x;
    //initial guess of 1 / sqrt(x);
    int i = 0x5F3759DF - (*(int *)&x) >> 1;
    float y = *(float *)&i;    
    //iterate
    y = y * (1.5F - y * y * xhalf);
    y = y * (1.5F - y * y * xhalf);
    //end of calculation of 1/sqrt(x);
    //return sqrt(x)=x * 1/sqrt(x);
    return (x * y)   
}; 

$ \log_2 (x) = e + \log_2 (1 + M)$이므로 x의 근삿값은 $\log_2 (1+M)$의 근삿값을 잘 선택하면 된다. 함수 $f(M)=\log_2 (1+M) - M$는 $M=1/\ln(2)-1$에서 최댓값 $\tt 0.0860713$을 가지므로

$$0\le \log_2 (1+M)-M \le 0.0860713$$

이다. 따라서 오차를 최댓값의 절반으로 잡아서 근사하면

$$\log_2 (1+M) \approx M+ \delta,\quad \delta = \frac{0.0860713}{2}=0.0430385$$ 

그리고 $\delta$는 $x$에 무관한 값이 된다. $x$의 비트 데이터를 정수로 표현하면

$$ i_x= 2^{23} (E+M) = 2^{23}(e + \delta +M) + 2^{23}(127-\delta) \approx 2^{23}\log_2(x) + 2^{23}(127-\delta).$$ 로 근사할 수 있다. 두 번째 항은 $x$에 무관한 값이다. 따라서, 제곱근의 역수에 해당하는 비트 데이터의 근사적인 정수 표현은

$$i_{\frac{1}{\sqrt{x}}} \approx 2^{23}\log_2 (1/\sqrt{x}) + 2^{23}(127-\delta)\approx \frac{3}{2} 2^{23} (127-\delta ) -\frac{1}{2} i_x$$

로 쓸 수 있고, 위에서 구한 $\delta$을 대입하면

$$\tt  i_{\frac{1}{\sqrt{x}}} \approx  0x5F37BCB6 - i_x>>1$$ 

을 얻을 수 있다. $\delta = 0.0450466$을 대입하면 위에 적힌 코드에서 사용한 $\tt 0x5F3759DF$가 됨을 알 수 있다. 그러나 어떤 근거로 이 값을 선택했는지는...