Vấn Đề Ẩn Chứa Sự Không Nhất Quán Trong Kết Quả Phép Toán Dấu Phẩy Động

Hôm qua, A Nam phát hiện một lỗi bug trong dự án do sự không nhất quán của phép toán dấu phẩy động. Điều kỳ lạ là cùng một đoạn mã C hoàn toàn giống nhau, tham số đầu vào cũng hoàn toàn tương đồng, nhưng lại cho ra hai kết quả khác biệt. Hiện tượng này khiến tôi vô cùng hứng thú và quyết định nghiên cứu kỹ lưỡng nguyên nhân cốt lõi.

Đoạn mã gốc tương đối phức tạp. Sau khi nắm rõ bản chất vấn đề, tôi đã đơn giản hóa đoạn mã gây lỗi để tái hiện hiện tượng này như sau:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

static void
foo(float x) {
  float xx = x * 0.01f;
  printf("%d\n", (int)(x * 0.01f));
  printf("%d\n", (int)xx);
}
int
main() {
  foo(2000.0f);
  return 0;
}

Khi biên dịch và chạy chương trình bằng gcc 4.9.2 với tùy chọn bắt buộc sử dụng toán học dấu phẩy động x87, bạn sẽ thấy kết quả đầy bất ngờ:

1
2
3

gcc a.c -mfpmath=387
19
20

Lần đầu xuất ra 19, lần sau lại là 20. Vì sao lại như vậy? Hãy cùng xem đoạn mã hợp ngữ do gcc tạo ra (đoạn mã rút gọn liên quan):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

  flds  16(%rbp)
  flds  .LC0(%rip)
  fmulp  %st, %st(1)
  fstps  -4(%rbp)     ; 1. Lưu kết quả x * 0.01f vào biến float trong bộ nhớ
  flds  16(%rbp)
  flds  .LC0(%rip)
  fmulp  %st, %st(1)
  fisttpl -20(%rbp)    ; 2. Chuyển trực tiếp kết quả x * 0.01f sang kiểu nguyên
  movl  -20(%rbp), %eax
  movl  %eax, %edx
  leaq  .LC1(%rip), %rcx
  call  printf
  flds  -4(%rbp)         ; 3. Đọc lại kết quả từ bước 1 đã lưu trong bộ nhớ
  fisttpl -20(%rbp)
  movl  -20(%rbp), %eax
  movl  %eax, %edx
  leaq  .LC1(%rip), %rcx
  call  printf

Chúng ta hãy phân tích ba điểm chú ý trong đoạn mã này:

Đầu tiên, giá trị 0.01 không thể biểu diễn chính xác trong hệ nhị phân, do đó phép nhân với 0.01 luôn chứa sai số nhất định. Mặc dù cả hai phép tính đều là x * 0.01f và theo quy tắc chuyển đổi của C, khi các toán hạng đều là float thì phép tính sẽ được thực hiện với độ chính xác float. Tuy nhiên, gcc không thiết lập chính xác chế độ kiểm soát độ chính xác của FPU.

Tại bước 2, kết quả phép nhân được chuyển trực tiếp từ thanh ghi dấu phẩy động sang kiểu nguyên. Trong khi đó, tại bước 1, kết quả bị ép xuống độ chính xác thấp hơn khi lưu vào bộ nhớ thông qua lệnh fstps. Đến bước 3, khi đọc lại giá trị này từ bộ nhớ bằng lệnh flds, giá trị trong thanh ghi st đã bị thay đổi so với bản gốc.

Sự khác biệt tinh tế này giữa hai lần lưu trữ và đọc giá trị từ thanh ghi dấu phẩy động chính là nguyên nhân khiến lệnh fisttpl chuyển đổi sang kiểu nguyên cho ra hai kết quả khác nhau.

Bổ sung ngày 14/6/2018: Trong một bài phỏng vấn đăng trên DDJ năm 1997, vấn đề tương tự đã được đề cập. Trích dẫn như sau:

“Tôi xin đưa ra một ví dụ xảy ra ngay hôm qua. Một sinh viên đang thực hiện phép tính mô phỏng plasma, trong đó xuất hiện các “bức tường phản xạ”. Khi một hạt vượt qua bức tường này, nó không nên đi tiếp mà phải bị phản xạ lại. Có những đoạn mã đại loại:

1 2 3 4 5 6 7

float x, y, z; int j; ... x = y + z; if (x >= j) replace (x); y = x; ...

Khi bật chế độ tối ưu hóa, giá trị x được tính toán trong thanh ghi với độ rộng bổ sung. Điều này xảy ra trên các máy Intel do thanh ghi có độ rộng lớn hơn, và cũng xuất hiện trên một số kiến trúc khác. Giá trị x thực tế không còn là kiểu float thuần túy mà được lưu trong thanh ghi rộng hơn. Khi máy lưu trữ x, nó ép giá trị này về độ chính xác float. Tuy nhiên, giá trị dùng trong phép so sánh lại là giá trị gốc trong thanh ghi rộng. Điều kiện “x trong thanh ghi >= j” không đảm bảo rằng x khi lưu vào y sẽ nhỏ hơn j. Đôi khi y=j xảy ra trong khi lẽ ra không bao giờ như vậy. Sinh viên của tôi đã kỳ vọng vào tính “minh bạch tham chiếu” mà các lập trình viên compiler thường đề cập, nhưng thất vọng vì các kỹ sư compiler đã tối ưu hóa hiệu suất bằng cách không tải lại giá trị từ bộ nhớ vào thanh ghi.”