Thuật toán tính gần đúng nghịch đảo của căn của giá trị x

Thuật toán này theo mình là một "điều kỳ diệu". Chỉ với đúng 10 dòng code, mà trong đó, thực chất chỉ có 2 dòng thực hiện tính toán, với các phép tính đơn giản gồm phép cộng, nhân và dịch bit, thuật toán đã đưa ra một phép tính gần đúng cực kỳ hiệu quả để tính \(\displaystyle \frac 1 {\sqrt{x}}\). Đây là thuật toán được xem là "chìa khóa" , dùng trong game Quake 3, để mở cánh cổng đến thế giới giả lập 3D, vào thời kỳ sơ khai khi sức mạnh tính toán còn chưa đủ mạnh.

English version here

Source Code

float Q_rsqrt(float number)
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;                       // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration  

  return y;
}

Chúng ta cùng tìm hiểu thuật đoạn code trên

Chuẩn IEE754

Theo như chuẩn này, một số float được dùng 32 bit để lưu trữ làm 3 phần: S - 1 bit dấu thể hiện số âm / số dương, E - 8 bit thể hiện số mũ, M - 23 bit lưu giá trị được chuẩn hóa về \(2^{23}\).

Như vậy, một biến \(x \) được tính bởi công thức: \(x=(1+\frac M {2^{23}})*2^{E-127} \quad (*)\)

Và, thể hiện trên giá trị bit của \(x\) sẽ là: \(M + 2^{23}*E \quad (**)\)

Tách thông tin nhị phân của \(x\) từ \(\log_2 x\)

Với \(x\) có giá trị nhỏ \(x \in [0,1)\), chúng ta có thể tính xấp xỉ của \(x\) thông qua công thức xấp xỉ: \(x \approx \log_2(1+x)\). Vậy:

$$\log_2(1+x) = x + \mu; \quad (\mu \text{ is empirically chosen at 0.043})$$

Với \(x\) bất kỳ, ta có theo (*):

$$\log_2 x = \log_2 \left((1+\frac M {2^{23}})*2^{E-127}\right) = \log_2(1+\frac M {2^{23}}) + E - 127$$

Vì \(\displaystyle \frac M {2^{23}} \in [0,1)\), ta áp dụng công thức xấp xỉ:

$$\log_2(1 + \frac M {2^{23}}) = \frac M {2^{23}} + \mu$$

Lúc đó:

$$\begin{align*} \log_2 x &= \frac M {2^{23}} + \mu + E - 127\\ &= \frac M {2^{23}} + \frac {2^{23} * E} {2^{23}} + \mu - 127 \\ &= \frac 1 {2^{23}}(\underbrace{M + 2^{23}E}_{\text{bit representation (**)}}) + \mu - 127 \end{align*}$$

Tính \(\frac 1 {\sqrt x}\) dùng \(\log_2\)

Ý nghĩa của dòng code

i  = 0x5f3759df - ( i >> 1 );

Ta có:

$$\log_2(\frac 1 {\sqrt{x}}) = \log_2(x^{-\frac 1 2}) = -\frac 1 2 \log_2(x)$$

Đặt \(\Gamma=\frac 1 {\sqrt{x}} \implies \log_2(\Gamma) = \log_2 (\frac 1 {\sqrt{x}}) = -\frac 1 2 \log_2(x)​\)

Ta thay vào 2 vế của biểu diễn bit:

$$\begin{align*} \frac 1 {2^{23}}(M_{\Gamma} + 2^{23}*E_{\Gamma}) + \mu - 127 = -\frac1 2 \left( \frac 1 {2^{23}}(M_x + 2^{23}*E_x) + \mu -127 \right) \\ \iff M_{\Gamma} + 2^{23}E_{\Gamma} = \underbrace{\frac 3 2 2^{23}(127-\mu)}_{\text{0x5f3759df}} - \underbrace{\frac 1 2(\underbrace{M_x + 2^{23}*E_x)}{x}}_{i>>1} \end{align*}$$

Đây là bí ẩn của dòng code tính toán đầu tiên.

Áp dụng Newton-Ralphson để tìm nghiệm gần đúng

$$\boxed{x_{new} = x_{old} - \frac {f(x)} {f' (x)}}$$

Đặt \(\displaystyle f(y)=\frac 1 {y^2} - x\), như vậy: \(\displaystyle f(y)=0 \iff y=\frac 1 {\sqrt{x}}\). Để giải nghiệm cho \(f(y)=0\), ta dùng phương pháp lặp Newton-Ralphson.

Đầu tiên ta đi tính đạo hàm của \(f(y)\)

$$f(y)=\frac 1 {y^2} - x \implies f'(y) = \frac {y-xy^3} {-2}$$

Vậy, áp dụng công thức, ta có:

$$y_{new} = y - \frac {f(y)}{f' (y)} = y - \frac {y-xy^3} {-2} = y(\frac 3 2 - \frac x 2 y^2)$$

Và đây cũng là ý nghĩa ẩn giấu của dùng code cuối cùng:

 y  = y * ( threehalfs - ( x2 * y * y ) );

Demo và JS code

JS Code

function fastInverseSqrt(x) {
            const threeHalfs = 1.5;

            let i = new Float32Array(1);
            let y = new Float32Array(1);

            y[0] = x;

            i = new Int32Array(y.buffer); // Reinterprets the bits as an integer
            i[0] = 0x5f3759df - (i[0] >> 1); // Magic number and bit shift

            y = new Float32Array(i.buffer); // Reinterpret the bits back to float
            y[0] = y[0] * (threeHalfs - (x * 0.5 * y[0] * y[0])); // First iteration of Newton's method
            y[0] = y[0] * (threeHalfs - (x * 0.5 * y[0] * y[0])); 

            return y[0];
        }