Một phương pháp tính nhanh căn bậc 2

Như các bạn đã biết để tính căn bậc 2 của một số, chúng ta thường dùng các phương pháp lặp để tính giá trị. Các phương pháp này thường có nhược điểm là khá chậm và thật là tệ nếu chúng ta cần tính căn bậc 2 của một số có giá trị lớn. Trong bài viết này chúng ta sẽ cùng tìm hiểu một phương pháp được dùng để tính căn bậc 2 trong đồ hoạ máy tính, nơi mà các giá trị cần tính cần độ chính xác tương đối không quá cao. Kỹ thuật này được áp dụng rất tốt trong trường hợp chúng ta chỉ cần một vài chữ số thập phân chính xác.

Trong kỹ thuật này, chúng ta sẽ sử dụng các bit quan trọng nhất của mantissa để làm chỉ mục cho một bảng tìm kiếm căn bậc 2. Bằng cách sử dụng bảng tra cứu này và giá trị một nửa ở phần số mũ (exponent), chúng ta có thể tìm được giá trị căn bậc 2. Với một số thập phân bất kì chúng ta có thể biểu diễn dưới dạng nhị phân: $$2^{\pm e_1e_2e_3*} \pm m_1m_2…m_n $$

Trong đó m và e đều được biểu diễn dưới dạng bit. Tất cả các số thập phân đều có thể biểu diễn dưới dạng chuẩn hoá có dạng: $$eeeeeeee \pm 1.mmmmmmm$$

Do tất cả các số thập phân đều có thể chuẩn hoá nên hầu hết cả hệ thống chỉ lưu phần phân số. Để tiện cho việc minh hoạ chúng ta sẽ biểu diễn số thập phân với 8 bit phần số mũ (exponent) và 7 bit mantissa, số sẽ được lưu dưới dạng như sau:

$$[2\text{^}]e_7e_6e_5e_4e_3e_2e_1e_0[*] \pm [1.] m_6m_5m_4m_3m_2m_1m_0$$

$$\sqrt{2^{m}n} = 2^{\frac m 2}\sqrt{n} \text{ (1)}$$

Từ công thức 1 ta có thể thấy được việc tìm căn bậc 2 của một số đã được đơn giản hoá bằng cách chia đôi phần exponent và tìm căn bậc 2 của phần mantissa. Khi phần exponent là một số lẻ chúng ta không thể chia được cho 2, chúng ta cần tách phần exponent thành một số chẵn thành 2 phần: $$e_7e_6e_5e_4e_3e_2e_10$$ và một phần được gọi là “quaternary mantissa” $$e_0m_6m_5m_4m_3m_2m_1m_0$$ trong nửa khoảng $$[1..4)$$

Đáng lưu ý ở đây, chúng ta có thể thấy được do bit cuối phần exponent là 0 nên phần exponent chắc chắn là số chẵn. Còn ở phần mantissa bit dấu được bỏ qua. Tất cả các giá trị được lưu trong nửa khoảng $$[1..4)$$ mà không đánh mất dữ liệu không tin nào ở số thập phân.

Bảng truy xuất được tạo ra khi chương trình được khởi tạo và lưu các giá trị căn bậc 2 của số

$$2^{e_0}*1.m_6m_5m_4m_3m_2m_1m_0$$

Tiếp đến chúng ta hãy xem thử mã giả của 2 chương trình tạo bảng tra cứu và tính căn bậc 2 của một số:

1. Cách tính giá trị U căn bậc 2 của một số V

e ← exponent(V): Lấy phần mũ của số V.
i ← mantissa(V): Lấy phần mantissa của số V.
if (e bit-and 1): Kiểm tra xem bit cuối cùng của phần mũ có bằng 1 hay không (nghĩa là phần mũ có lẻ hay chẵn).
set the high bit of i: Nếu phần mũ lẻ, đặt bit cao nhất của phần mantissa bằng 1.
e ← e/2: Chia phần mũ cho 2 (lưu ý là phép chia này phải giữ nguyên dấu).
j ← T[i]: Tìm giá trị tương ứng của i trong bảng tra T để lấy phần mantissa của căn bậc hai.
U ← 2^e * j: Lặp lại kết quả bằng cách nhân phần mantissa vừa tìm được với 2 mũ phần mũ mới.

Ý tưởng chính của thuật toán:

Chia nhỏ vấn đề: Thuật toán chia việc tính căn bậc hai thành hai phần: tính căn bậc hai của phần mantissa và điều chỉnh phần mũ.
Bảng tra cứu: Sử dụng một bảng tra (T) để lưu trữ trước các giá trị căn bậc hai của các phần mantissa thường gặp. Điều này giúp tăng tốc độ tính toán.
Điều chỉnh phần mũ: Chia đôi phần mũ để giảm một nửa giá trị căn bậc hai.

2. Cách tạo bảng tra cứu

build_table(precision: integer, table: array [0..2^precision] of integer): Đây là hàm để xây dựng bảng tra.
i, j: integer; f, sf: real: Khai báo các biến.
for i: integer 0, 2^precision – 1: Lặp qua tất cả các giá trị có thể của phần mantissa.
f ← 1.i; table[i] ← mantissa(sqrt(f)): Tính căn bậc hai của số 1.i và lưu phần mantissa vào bảng tra.
f ← 2^1 * 1.i; table[i + 2^precision] ← mantissa(sqrt(f)): Tính căn bậc hai của số 2^1 * 1.i và lưu phần mantissa vào bảng tra cứu.

Ý tưởng của thuật toán:

Xây dựng bảng tra cứu: Hàm này tạo ra bảng tra cứu T chứa các giá trị căn bậc hai của các phần mantissa có thể xảy ra. Bảng tra cứu này được sử dụng trong quá trình tính toán chính.
Độ chính xác: Độ chính xác của kết quả phụ thuộc vào kích thước của bảng tra (precision). Càng nhiều phần tử trong bảng tra cứu, kết quả càng chính xác.

Sau đây là đoạn code đã được tối ưu được trình bày thông qua ngôn ngữ Javascript. Thuật toán và kết quả của phương pháp được tham khảo từ sách “Graphics Gems” của các tác giả Paul Lalonde and Robert Dawson.

Thuật toán gốc được tác giả trình bày tại đây: https://github.com/erich666/GraphicsGems/blob/master/gems/SquareRoot.c

				
					// Initialize the square root table
const t = new Uint16Array(256); // 0x100 = 256

function buildTable() {
    let f = new Float32Array(1); // Use a typed array to access bits directly
    let fi = new Uint32Array(f.buffer); // View the same data as unsigned integers

    for (let i = 0; i <= 0x7f; i++) {
        fi[0] = 0;

        // Build a float with the bit pattern i as mantissa and exponent 0 (stored as 127)
        fi[0] = (i << 16) | (127 << 23); 
        f[0] = Math.sqrt(f[0]); 

        // Take the square root, then strip the first 7 bits of the mantissa into the table
        t[i] = (fi[0] & 0x7fffff) >> 16; 

        // Repeat the process, this time with an exponent of 1 (stored as 128)
        fi[0] = 0;
        fi[0] = (i << 16) | (128 << 23);
        f[0] = Math.sqrt(f[0]);
        t[i + 0x80] = (fi[0] & 0x7fffff) >> 16;
    }
}

// Fast square root by table lookup
function fsqrt(n) {
    if (n === 0) return 0; 

    let f = new Float32Array(1);
    f[0] = n;
    let num = new Uint32Array(f.buffer);

    let e = (num[0] >> 23) - 127; 
    num[0] &= 0x7fffff; 

    if (e & 0x01) num[0] |= 0x800000; 
    e >>= 1; 

    // Table lookup and reconstruction
    num[0] = ((t[num[0] >> 16]) << 16) | ((e + 127) << 23);
    return f[0];
}

// Call buildTable() once to initialize the lookup table
buildTable(); 

console.log(fsqrt(10), Math.sqrt(10))
// Output:

// 3.15625 3.1622776601683795

Một phương pháp tính nhanh căn bậc 2

Subscribe to SkyGLab