Gói float chuẩn VHDL tính sai sqrt(2) và nhiều phép toán khác

Gói thư viện số thực dấu phẩy động ieee.float_pkg — một phần của chuẩn VHDL IEEE 1076-2008 và 1076-2019 — trả về kết quả không được làm tròn chính xác ở nhiều phép toán cơ bản. Nghiêm trọng nhất: sqrt(2.0) ở định dạng float64 cho ra 3FF6A09E667F3BCC, trong khi giá trị đúng phải là 3FF6A09E667F3BCD — thấp hơn một đơn vị ở vị trí cuối cùng (ULP, đơn vị nhỏ nhất có thể biểu diễn).

Vấn đề không chỉ nằm ở hàm căn bậc hai. Một cuộc kiểm tra có hệ thống đã phát hiện hàng trăm nghìn lời gọi sai trên nhiều phép toán. Tác giả đã công bố toàn bộ bằng chứng, nguyên nhân theo số dòng và bản vá trong một kho mã nguồn mở.

Phạm vi ảnh hưởng

Mã nguồn được kiểm tra là tệp chính thức của IEEE: ieee/float_generic_pkg-body.vhdl trên nhánh phát hành của kho opensource.ieee.org/vasg/Packages, dài 5.712 dòng, ghi $Revision: 1220. Nhánh phát triển của nhóm làm việc IEEE P1076 trên GitLab có nội dung y hệt.

Đáng chú ý là GHDL và NVC — hai trình mô phỏng VHDL mã nguồn mở phổ biến nhất — đóng gói nguyên văn tệp này mà không sửa gì. Các phép thử được chạy trên GHDL 7.0.0-dev và NVC 1.23.1; khi cả hai chạy cùng bộ vector, chúng ghi ra tệp kết quả giống nhau.

Cách kiểm chứng

Tác giả đặt gói thư viện cạnh một bộ tham chiếu số học chính xác (exact.py): mỗi toán hạng là một dấu và một phân số của hai số nguyên, phép toán được thực hiện trên phân số rồi chỉ làm tròn một lần vào định dạng đích theo kiểu làm tròn mà lời gọi chỉ định.

Với cấu hình mặc định của gói (làm tròn về gần nhất, 3 bit bảo vệ, bật số không chuẩn hóa), nhiều kết quả sai rõ rệt:

• mac(-6.0, -1.0, -5.0) trả về -15.0 thay vì 1.0
• mac(1.0, 1.0, -infinity) trả về +infinity thay vì -infinity
• reciprocal(2.0**127) ở float32 trả về +infinity thay vì 2.0**-127
• resize(2.0**-512 từ float64 sang float32) trả về 1.0 thay vì 0.0
• scalb(1.0, 300) trả về 0.0 thay vì +infinity
• to_signed(0.50000006, 32) trả về 0 thay vì 1

Các toán hạng này không phải trường hợp hiếm. Trong dãy số nguyên từ 1 đến 4.096, sqrt trả về float32 lệch một đơn vị ở 256 trường hợp và float64 lệch ở 264 trường hợp. Trong 4.913 bộ ba toán hạng của mac, có 290 bộ ba cho kết quả float32 sai. Việc thu hẹp từ float64 xuống float32 trả về +infinity cho mọi lũy thừa của 2 từ 2^-1074 đến 2^-1024, và làm dừng trình mô phỏng ở khoảng 2^-662 đến 2^-639.

Thống kê sai số

Mỗi tổ hợp phép toán, định dạng và kiểu làm tròn được chạy 10.000 lời gọi, với toán hạng nhắm vào những điểm dễ hỏng: số không chuẩn hóa, hai đầu dải số mũ, phân số làm tròn tràn sang binade kế tiếp, điểm hòa chính xác, số không có dấu, vô cực.

Ở float32, tổng cộng 81.068 kết quả sai và 272 lần dừng trên 859.771 lời gọi. Ở float64 là 86.081 sai và 445 lần dừng trên 857.234 lời gọi. float16 có 113.980 sai trên 866.782 lời gọi; float128 có 90.744 sai và 585 lần dừng trên 855.777 lời gọi.

Với định dạng nhỏ, có thể thử toàn bộ toán hạng: định dạng 7 bit (3 bit số mũ, 3 bit phần lẻ) chạy 8.692.326 lời gọi và có 1.271.619 kết quả sai. Định dạng 10 bit chạy 17.340.682 lời gọi với 883.716 kết quả sai. Riêng mac ở chế độ làm tròn về gần nhất sai ở 243.520 trong 2.097.152 bộ ba của định dạng 7 bit.

Nguyên nhân theo từng dòng mã

Báo cáo chỉ ra lỗi cụ thể ở từng hàm (số dòng theo tệp của IEEE):

• normalize, dòng 422: số mũ bị cắt xuống chỉ còn hai bit nhiều hơn định dạng đích, dẫn đến phép lấy modulo. Từ float64 sang float32, số mũ còn 10 bit, nên bị tính theo modulo 1.024: 2^-512 thành 1.0, 2^600 thành 2^88. Khi số mũ bị cắt rơi vào vùng số không chuẩn hóa, chỉ số dịch vượt ra ngoài phần lẻ và làm dừng trình mô phỏng. Hàm này cũng trả về vô cực hoặc số 0 bất kể kiểu làm tròn.
• divide (dòng 1243) và reciprocal (dòng 1099): bit dính (sticky bit) đưa vào bước làm tròn bị gán cứng là '1', nên thương chính xác bị làm tròn như thể còn dư.
• reciprocal, dòng 1053–1093: số mũ kết quả -e - 3 được giữ trong số bit bằng với e, nên binade cao nhất bị tràn. Ở float32, mọi đối số từ 2^127 trở lên trả về vô cực.
• sqrt, dòng 1787: phương pháp lặp Newton từ (1 + x)/2 với số bước cố định, không có bit dính. Với số bình thường, kết quả lệch tối đa một đơn vị; với số không chuẩn hóa, điểm khởi đầu quá xa nên 128 trong 359 trường hợp float32 sai, có trường hợp sai hơn hai đơn vị.
• mac, dòng 1438–1536: tích bị cắt về độ rộng một toán hạng trước khi cộng với c, và có nhánh trừ số lớn cho số nhỏ khi dấu khác nhau.
• add, dòng 839: số 0 chính xác luôn mang dấu dương, trái với quy định của IEEE 854 6.3 khi làm tròn về âm vô cực.
• to_signed, to_unsigned, to_integer: ba bit dưới phần nguyên được giữ, phần còn lại bị bỏ mà không có bit dính.
• dividebyp2 và scalb: đưa kết quả vào normalize không có bit bảo vệ, nên kết quả không chuẩn hóa bị cắt thay vì được làm tròn.

Bản vá và kết quả

Bản vá float_generic_pkg-body.patch chỉ sửa phần thân gói, 33 đoạn, không đổi phần khai báo. Các thay đổi chính:

• normalize xét số mũ ở độ rộng đầy đủ để phân loại số 0, số không chuẩn hóa, số bình thường và vô cực, đồng thời trả về số lớn nhất hoặc nhỏ nhất khi kiểu làm tròn yêu cầu.
• divide và reciprocal lấy bit dính từ phần dư của phép chia; reciprocal giữ số mũ rộng thêm hai bit và gán dấu cho số 0, vô cực theo đối số.
• sqrt lấy căn bậc hai nguyên của phần định trị một bit mỗi bước, dùng phần dư làm bit dính.
• mac đặt toàn bộ tích và c vào một cửa sổ 3 * fraction_width + 9 bit, cộng trừ tại đó và làm tròn một lần — đúng chuẩn fused multiply-add.
• add, resize và các hàm chuyển đổi giữ dấu của số 0 và làm tròn đúng theo kiểu làm tròn được chỉ định.

Trên cùng bộ vector, gói đã vá cho 0 kết quả sai và 0 lần dừng ở tất cả các nhóm thử, kể cả 8.692.326 lời gọi của định dạng 7 bit và 17.340.682 lời gọi của định dạng 10 bit.

Gói số dấu phẩy tĩnh cũng có lỗi

ieee.fixed_pkg — cũng lấy từ nguồn IEEE và được NVC đóng gói nguyên văn — được chạy toàn bộ toán hạng của sáu cặp dải chỉ số nhỏ, tổng 751.418 lời gọi. Phép cộng, trừ, nhân, rem, mod, so sánh, add_carry, resize, to_integer, scalb, đảo dấu và abs đều đúng. Riêng phép chia sai 41.026 lời gọi.

Ví dụ: 0.125 / 1.875 trong ufixed (3 downto -3) trả về 0.0625, trong khi thương thật là 0.0667 và số gần nhất định dạng chứa được là 0.0703. Nguyên nhân là thương được chia nguyên đến các bit bảo vệ, và bước làm tròn chỉ thấy các bit đó chứ không thấy phần dư.

Cách tự kiểm chứng

Chỉ cần GHDL, có thể tự thấy lỗi bằng hai lệnh:

ghdl -a --std=08 plain.vhd
ghdl --elab-run --std=08 plain

Tệp plain.vhd thực hiện 25 lời gọi, trong đó 24 lời gọi in ra WRONG và lời gọi thứ 25 làm dừng trình mô phỏng với lỗi chỉ số ngoài biên tại dòng 358 của phần thân gói.

Lỗi đã được báo cáo cho nhóm làm việc IEEE P1076 vào ngày 13 tháng 11 năm 2026 với mã VHDL-Issues#388, đồng thời gửi báo cáo cho các dự án ghdl/ghdl#3486, nickg/nvc#1758 và FPHDL/fphdl#4.

Góc nhìn

Với các kỹ sư thiết kế phần cứng và kỹ sư FPGA dùng VHDL, đây là lời nhắc rằng "có trong chuẩn" không đồng nghĩa với "đúng". Gói float_pkg thường được dùng để mô phỏng thuật toán dấu phẩy động trước khi tổng hợp, nên sai số một ULP ở bước mô phỏng có thể khiến việc đối chiếu với mô hình tham chiếu bằng C hoặc Python trở nên khó hiểu. Các phép mac, sqrt và chuyển đổi định dạng là những chỗ cần cẩn trọng nhất, đặc biệt khi kết quả tiếp tục bị chuyển đổi hoặc tích lũy qua nhiều bước.

Về phía công cụ, người dùng GHDL và NVC hiện phải tự chịu trách nhiệm nếu cần kết quả chính xác tuyệt đối — hoặc áp bản vá, hoặc thay bằng một thư viện khác. Việc cả hai trình mô phỏng cùng ghi ra một kết quả sai cũng cho thấy nguy cơ của việc coi sự đồng thuận giữa hai công cụ là bằng chứng về tính đúng đắn.