Cấu trúc số thực, sai số hệ thống và việc làm tròn số

[ledoninh]

Tên bài viết: Cấu trúc số thực, sai số hệ thống và việc làm tròn số
Tác giả: ledoninh
Cấp độ bài viết: Chưa đánh giá
Tóm tắt: Thấy dân IT (kể cả sinh viên IT) im lặng trước nhiều câu hỏi về số thực, cách làm tròn số thực, mình không vui lắm, nên đặt ra topic này để cùng thảo luận (trả lời từng thắc mắc thì nhiều quá và dễ bị lặp lại câu hỏi). Hy vọng các bạn không “im ru bà rù” như ở các topic trước của mình (Cấu trúc BMP, giai thừa mở rộng).

---

Khái niệm “số thực” ở đây chỉ bao gồm kiểu Single (4 bytes) và kiểu Double (8 bytes) của VB, không bao gồm kiểu Real của Pascal (6 bytes, ai quan tâm thì trình bày sau).

1. Cấu trúc
Cả hai kiểu số thực này được xây dựng theo chuẩn IEEE:
x = (dấu) 2^n (1+a) với 1 > a >= 0, tức là gồm 1 bit dấu (cao nhất), các bit bậc nhị phân cao nhất (n biểu diễn qua 11 bit đối với Double hoặc 8 bit đối với Single) và các bit biểu diễn giá trị còn lại (a qua 52 bit đối với Double hoặc 23 bit đối với Single).

2. Sai số hệ thống và cách tạo thành các bytes
Như vậy, nếu số a không thể biểu diễn chính xác được qua một số hữu hạn (cụ thể là 52 và 23) các bit thì nó được làm tròn, gây ra sai số hệ thống, sai số này là cộng trừ 2^(-53) đối với Double và 2^(-24) đối với Single.
Sai số như thế là rất nhỏ, và mình xin nói trước rằng đó không phải là lý do khiến hàm Round cho kết quả sai. Cái này lại là một loại sai số hệ thống khác, khi hiển thị theo hệ thập phân. (1 chữ số thẩp phân tiêu tốn hơn 3 bít nhị phân, 15 chữ số hết gần 49 bit, còn dư hơn 3 bit nhưng không đủ cho 1 chữ số nữa). Minh sẽ đưa ra phương án Round sau cùng.

2.1. Tính n và a
Để biểu diễn cặp (n,a), số Double và số Single, ta có thể dùng kiểu biến sau :

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. Type thanhphan
2.     bac As Integer
3.     triso As Double
4. End Type
5. Type ieee64
6.     b(1 To 8) As Byte
7. End Type
8. Type ieee32
9.     b(1 To 4) As Byte
10. End Type

GeSHi © Codebox Plus

Xác định n và a theo số thực r bằng hàm sau :

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. Function expo(r As Double) As thanhphan 'Tìm lũy thừa cao nhất
2. Dim n%, t#
3. t = Abs(r)
4. If t = 0 Then Exit Function ' Khỏi xét trong hàm này
5. n = 0 'Với 2 < n <= 1 thì 2^0 là cao nhất, không cần While... sau đây
6. While t >= 2 ' Nếu t = 2^1 trở lên thì làm cho đến < 2
7. n = n + 1
8. t = t / 2
9. Wend
10. While t < 1 ' Nếu t < 2^0 thì làm cho đến >= 1
11. n = n - 1
12. t = t * 2
13. Wend
14. expo.bac = n
15. expo.triso = t - 1 'Trị số sau dấu phẩy
16. End Function

GeSHi © Codebox Plus

2.2. Cách tạo các bytes
Khi làm tròn số a thành 52 bit, ta có thể dùng kiểu “cắt” (hàm Int) hay kiểu “quy tròn” (hàm Round). Vì rằng dùng Int tổng quát hơn, có thể thay thế Round (Ví dụ Round(1.x, 0) có thể thay bằng Int(1.x + 0.5)), nên mình sẽ dùng hàm Int.

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. Function createDouble(r As Double) As ieee64
2. 'Gồm 1 bit dấu, 11 bit bậc nhị phân n, 52 bit trị số lẻ
3. Dim ii%, tt#
4. Dim bb As Byte, tp As thanhphan
5. With createDouble
6. For ii = 1 To 8
7. .b(ii) = 0
8. Next ii
9. End With
10. ‘ Với số 0 n = -(vô cực), nhưng được gán = 0
11. If r = 0# Then Exit Function
12. tp = expo(r)
13. ‘*** Chỗ này phải chèn thêm một đoạn nếu dùng kiểu “quy tròn”
14. 'Quy ước tịnh tiến n lên 1023 đơn vị sau đây để n > 0
15. ‘ và như vậy, nếu không muốn exit khi gặp số Double 0#,
16. ‘ ta có thể thay lệnh “If r = 0# Then Exit Function” ở trên
17. ‘ bằng lệnh “If r = 0# Then tp.bac = -1023”
18. tp.bac = tp.bac + 1023
19.  
20. 'Byte thứ 8 gồm 1 bit dấu và 7 bit cao của bậc n
21. bb = tp.bac \ 16 'Trong 11 bit, lấy 7, còn 4 bit thấp tính sau
22. If r < 0 Then bb = bb + 128 ' dấu "-" ở bit cao nhất
23. createDouble.b(8) = bb
24.  
25. 'Byte thứ 7 gồm 4 bit còn lại của bậc và 4 bit cao của trị số a
26. tt = tt + tp.triso
27. bb = Int(tt * 16)
28. tt = tt * 16 - bb 'Phần dư tiếp theo
29. bb = bb + (tp.bac Mod 16) * 16
30. createDouble.b(7) = bb
31.  
32. 'Sau đây là các bytes tiếp theo. Chú ý: bình thường là 6,
33. ‘nhưng có thể tăng nếu làm cho mình cách đọc Double riêng
34. For ii = 6 To 1 Step -1
35. bb = Int(tt * 256)
36. tt = tt * 256 - bb
37. createDouble.b(ii) = bb
38. Next ii
39. ‘*** Kiểu “cắt” hay kiểu “quy tròn” thì vấn đề nằm ở
40. ‘ bước cuối của vòng lặp (byte số 1). Như trên là “cắt”,
41. ‘ còn trong file gửi kèm là “quy tròn”. Sở dĩ phải thêm
42. ‘ 1 đoạn vào phần đầu (nếu muốn “quy tròn”) là do
43. ‘ không muốn điều chỉnh các bytes trước nếu gặp
44. ‘ trường hợp byte số 1 sau khi quy tròn bị “vượt khung”
45. End Function

GeSHi © Codebox Plus

Đoạn thêm 0.5 chữ số cuối để Int tương đương Round vào (chỗ cần chèn nêu trên) như sau:

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. tt = 1
2. For ii = 1 To 53
3. tt = tt / 2
4. Next ii
5. If tp.triso + tt >= 1 Then
6. tp.bac = tp.bac + 1024
7. For ii = 1 To 8
8. createDouble.b(ii) = 0
9. Next ii
10. Exit Function
11. End If

GeSHi © Codebox Plus

2.3. Khắc phục lỗi trong việc làm tròn số
Việc làm tròn số thực trong khuôn khổ VB hay Excel (đều của Microsoft) là việc khó hiểu, có lẽ là chỉ có Microsoft mới trả lời được.
Đối với kiểu Double, có thể thử trong 1 project bất kỳ bằng cách gõ

x# = 1.23456789012345

Đến “5” là vừa đủ 15 chữ số, nếu gõ thêm “5” thì bị cắt, nếu gõ thêm 6 thì cũng bị cắt, nhưng “5” cũ biến thành “6”.

Kiểu Single lại tệ hại theo hướng khác (có lẽ vì 7 chữ số thập phân cần 23.2535 bit, lớn hơn 23 mà nó có). Thử bằng CommandButton lệnh sau:

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. Private Sub Command1_Click()
2. x! = 1.23456749 ‘Dấu “!” = As Single
3. MsgBox x!
4. End Sub

GeSHi © Codebox Plus

Khi chạy lệnh, bạn sẽ thấy kết quả là 1.234568, mặc dù sau “7” là “4” không tới “5”. Có vẻ như VB làm tròn số Single từng bước từ phải qua trái ? Như vậy, mẹo làm tròn số Double của taykhongbatgiac (trung chuyển qua Single nếu có thể) là có cơ sở. Nhưng với số Double ngoài phạm vi Single và với chính số Single thì làm thế nào ?
Tại sao bit cuối cùng có giá trị rất nhỏ mà thường số Double có sai số lớn thế (thí dụ, 1.4465 được hiểu là 1.44649999999995) ? Đó là vì khi đổi ra hệ thập phân, ta phải “cắt đuôi” khá nhiều. 4 bit của byte 7 chỉ phải chia cho 16, byte 6 phải chia cho 16*256, byte 5 chia cho 16*256*256, vv… bắt đầu từ byte 5 trở xuống khi chuyển phải cắt bớt !
Tóm lại là phải dùng mẹo thôi, không nhờ gì được các bit và byte. Để ý 2 yếu tố sau:
a) Rắc rối chính là chữ số 5 cuối cùng mà ta lại muốn làm tròn đến chữ số trước nó.
b) Số 0.5 bao giờ cũng được hiển thị đúng giá trị vì nó là lũy thừa của 2: 0.5 = 2^(-1).
Từ hai nhận xét trên, ta có thể xây dựng hàm làm tròn MyRound như sau:

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. Function MyRound(x As Double, n As Integer) as Double
2. x = x * 10 ^ n ‘nếu x quá lớn, gần “đụng trần” số Double thì không được
3. x = x + 0.5
4. x = int(x)
5. MyRound = x / 10 ^ n
6. End Function

GeSHi © Codebox Plus

Hết.

[truongphu]

Tôi rất ngưỡng mộ anh ledoninh
về nhà (sau cuộc gặp) thấy bài của anh (hiếm!)
Đọc nhưng chưa nắm bắt gì

Từ hai nhận xét trên, ta có thể xây dựng hàm làm tròn MyRound như sau:

Tôi nhớ loáng thoáng.. khi tôi viết hàm làm tròn trước đây vài tháng
Cu T7 nói hàm tôi viết chả giá trị gì
vậy nên (tự giác viết bài khác thay thế.. )
Mời anh ledoninh xem code sau
(Riêng code của anh, ngày mai tôi sẽ đọc kỹ, tha lỗi )

Mã: Chọn hết

1. Function MyRound(x As Double, n As Integer) As Double
2. x = x * 10 ^ n 'n?u x quá l?n, g?n “d?ng tr?n” s? Double thì không du?c
3. x = x + 0.5
4. x = Int(x)
5. MyRound = x / 10 ^ n
6. End Function
7.  
8. Private Sub Command1_Click()
9. Print MyRound(99999999999.5555, 3) ' <-- sô' Ðung trâ`n theo viê't code, chua phai sô' Ðung trâ`n theo Double
10. End Sub
11.  
12. Private Sub Command2_Click()  ' code T7, hình như là thế, tôi cũng quên vì bài đã xóa
13. Print Format(99999999999.5555, "#.###")
14. End Sub

[ledoninh]

Xin chào, tất cả !
Tôi không ngờ gặp lại topic này ở đây, chắc các admin bốc qua từ diễn đàn cũ. Ngay bác truongphu còn chưa đọc ở đó thì có thể còn nhiều người chưa đọc. Tôi ngại chuyển vì rất bận rộn, nhưng đã có người giúp thì tôi sẽ tìm cách edit lại, vì bối cảnh khi đó (nhiều người thắc mắc về vấn đề này) khác bây giờ.
Riêng vấn đề mà bác truongphu đưa ra, tôi có biết, nhưng ở đây tôi chỉ giới hạn ở thực tế (sai số tương đối đủ nhỏ, cho dù sai số tuyệt đối có vẻ lớn). Có một tài liệu vật lý nói rằng con số 10^100 là giới hạn tối đa của vũ trụ, như kích thước vũ trụ so với hạt quark, khối lượng vũ trụ so với photon đều chưa tới 10^100 lần ...
VB giới hạn ở 16 chữ số, calculator mà kế toán hay dùng chỉ có 12 chữ số, đó là người ta cũng hài lòng ở tính "đủ xài". Tôi cũng vậy, không định dùng tới string để xử lý số lớn, vì tôi chỉ quan tâm đến cách khắc phục kết quả làm tròn ngoài ý muốn của ta.
Tôi nhấn mạnh là ý muốn của ta thôi, vì từ nhỏ chúng ta đã quen với cách làm tròn như thế - như thế. Các ngôn ngữ lập trình có lý do để dùng cách làm tròn của họ (xem bài "Đâu cũng nói chuyện tiền, điên đầu thiệt !" - http://www.caulacbovb.com/forum/viewtop ... =48&t=5362 ) . Quả thật, nếu ta là người phát tiền, mà 0.5 chỉ được làm tròn lên 1 thì lỗ to. Hoặc, tính trung bình nhiệu giá trị thực nghiệm (đã được làm tròn từng số theo kiểu của ta) thì có thể tự chuốc lấy sai lệch.
Tóm lại, kiểu của ta dùng theo thói quen của ta và áp dụng cho ít số đơn lẻ thôi.
Rất cảm ơn bác đã cùng suy nghĩ và mong bác giúp biên tập lại cho nó phù hợp với ý nghĩa khiêm tốn của chủ đề này. Còn dùng hàm Format thì tôi không thích, vì MS VB dùng nó để hiển thị number thành string thôi, nhưng không chặt chẽ nên ta mới lợi dụng được.

[thuanfun]

Theo tôi thì hàm expo của bác ledoninh chỉ đơn giản như sau là OK rồi.

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. Function expo(r As Double) As thanhphan 'Tìm lũy thừa cao nhất
2. Dim t As Double, n As Integer
3. t = Abs(r)
4. If t = 0 Then Exit Function
5. n = Int(Log(t) / Log(2))
6. expo.bac = n
7. expo.triso = t / Exp(n * Log(2)) - 1
8. End Function

GeSHi © Codebox Plus

Hơn nữa bác ledoninh nên trích dẫn đầy đủ nguồn gốc các phát biểu sau:

Có một tài liệu vật lý nói rằng con số 10^100 là giới hạn tối đa của vũ trụ, như kích thước vũ trụ so với hạt quark, khối lượng vũ trụ so với photon đều chưa tới 10^100 lần ...

Khối lượng vũ trụ: Cỡ 10^55 kg Xem http://hypertextbook.com/facts/2006/KristineMcPherson.shtml
Khối lượng nghỉ của photon: 0 kg Xem http://vi.wikipedia.org/wiki/Photon
Vũ trụ:http://vi.wikipedia.org/wiki/V%C5%A9_tr%E1%BB%A5
Quark: http://en.wikipedia.org/wiki/Quark

[thuanfun]

Ngoài ra nếu dựa vào kiểu String ta có thể làm được nhiều thứ tuyệt vời hơn:

Mã: [Chọn hết] [Mở rộng/Thu hẹp] [Tải về]

1. Private Function MyRound(ByVal x As String, n As Integer) As String
2. Dim  k As Integer
3. k = InStr(x, ".")
4. If (k = 0) Or (n <= 0) Or ((k > 0) And (n >= Len(x) - k)) Then
5.    MyRound = x
6. Else
7.    If CInt(Mid(x, k + n + 1, 1)) >= 5 Then
8.      MyRound = Left(x, k + n - 1) + CStr(CInt(Mid(x, k + n, 1)) + 1)
9.    Else
10.      MyRound = Left(x, k + n)
11.    End If
12. End If
13. End Function
14.  

GeSHi © Codebox Plus

[ledoninh]

Hoan nghênh bạn có ý kiến làm đơn giản hoá hàm làm tròn ! Nội dung chính của chủ đề là viết về dữ liệu số thực trong VB, hàm làm tròn chỉ mang tính định hướng (code không chi tiết).
Về con số 10^100, nó tròn tới mức khó tin, cũng chỉ có tính chất định hướng (thông thường thì đúng), có nghĩa là chuyện vật lý thường thức ấy mà (không nhớ tên tài liệu). Tuy nhiên, cần nói thêm rằng photon không tồn tại ở trạng thái nghỉ - cái gọi là "khối lượng nghỉ của photon" cũng lại là 1 ước lệ. Là hạt cơ bản chỉ tồn tại trong sóng điện từ với tốc độ ánh sáng, nó có khối lượng (công thức Einstein), và vì thế có động lượng. Nói một cách công bằng, ta có thể có sóng điện từ với tần số thật nhỏ, do đó năng lượng của photon gần bằng 0, suy ra khối lượng của nó gần bằng 0 ... Thật ra tôi chưa bao giờ thử tính khối lượng của photon ánh sáng đỏ hay tím bằng bao nhiêu !
Nếu bạn vẫn chưa hài lòng thì tôi đành chịu.