天天看點

淺談浮點數運算的誤差

測試程式

我們知道,浮點數運算存在舍入誤差。在某些特殊的情況下,舍入誤差還可以累計到非常大的地步。讓我們來看一下測試程式吧:

1 using System;
 2  
 3 static class DecimalSumTester
 4 {
 5   static void Main(string[] args)
 6   {
 7     try
 8     {
 9       var n = (args.Length > 0) ? int.Parse(args[0]) : 10;
10       for (var i = 0; i < 3; i++) Console.WriteLine(F(n, i));
11     }
12     catch (Exception ex) { Console.WriteLine(ex.Message); }
13   }
14     
15   static decimal F(int n, int k)
16   {
17     var z = 0.1m + k * 1000000000000000000000000000m;
18     var w = decimal.Round(z) / 2;
19     while (n-- > 0) z += z / 2 - w;
20     return z - w * 2;
21   }
22 }      

在這個程式中:

  • 第 19 行通過 while 循環不斷進行累加: z += z / 2 - w; 。w 是不變的,而 z 是通過不斷累加而增大的。
  • 第 9 行讀取指令行參數作為 n 的值來決定 while 循環次數。
  • 第 15 至 21 行的 F 方法第一個參數 n 就是循環次數。第二個參數 k 決定累加初值的整數部分的大小。
  • 在這個程式中參數 k 僅取 0、1 和 2 三種情況。如果算術運算沒有誤差的話,這三種情況下 F 方法的傳回值應該一樣。

在 Linux 中編譯和運作

在 Arch Linux 64-bit 作業系統的 Mono 3.0.4 環境下編譯和運作:

work$ dmcs DecimalSumTester.cs
work$ mono DecimalSumTester.exe 24
1683.4112196028232574462890625
2730.9
2730.9
      

上述結果中第一行是計算出來的準确值。第二行和第三行的值理論上應該等于第一行。但是由于浮點數運算的舍入誤差累計的結果,最終答案的誤差相當大。在我的上一篇随筆“淺談 System.Decimal 結構”中提到,decimal 的算術運算在 Linux 環境下的舍入規則是四舍五入,這可能導緻誤差累計得比較大。而在 Windows 環境下的舍入規則是四舍六入五向偶。那麼我們就接着往下看吧。

在 Windows 中編譯和運作

在 Windows 7 SP1 32-bit 作業系統的 Microsoft .NET Framework 4.5 環境下編譯和運作:

D:\work> csc DecimalSumTester.cs
Microsoft(R) Visual C# 編譯器版本 4.0.30319.17929
用于 Microsoft(R) .NET Framework 4.5
版權所有 (C) Microsoft Corporation。保留所有權利。
D:\work> DecimalSumTester 24
1683.4112196028232574462890625
2097.9
0.1
      

這個運作結果和在 Linux 環境下大不相同。但是誤差也是相當的大。

調試程式

讓我們來看看在運算過程中發生了什麼吧。在上述測試程式中插入一些調試語句:

1 using System;
 2  
 3 static class DecimalSumDebug
 4 {
 5   static void Main(string[] args)
 6   {
 7     var n = (args.Length > 0) ? int.Parse(args[0]) : 10;
 8     for (var i = 0; i < 3; i++)
 9       Console.WriteLine(F(n, i).ToString().PadRight(77, '-'));
10   }
11     
12   static decimal F(int n, int k)
13   {
14     var z = 0.1m + k * 1000000000000000000000000000m;
15     var w = decimal.Round(z) / 2;
16     for (decimal x, y; n-- > 0; z += x)
17     {
18       x = (y = z / 2) - w;
19       Console.WriteLine("{0,-30}: {1,-30}: {2}", z, y, x);
20     }
21     return z - w * 2;
22   }
23 }      

這個程式和前面的測試程式的功能是相同,僅僅是運算過程中增加輸出中間變量的值的語句。

在 Linux 中調試

在 Arch Linux 的 Mono 環境下編譯和運作,輸出一大堆調試資訊:

work$ dmcs DecimalSumDebug.cs
work$ mono DecimalSumDebug.exe 10
0.1                           : 0.05                          : 0.05
0.15                          : 0.075                         : 0.075
0.225                         : 0.1125                        : 0.1125
0.3375                        : 0.16875                       : 0.16875
0.50625                       : 0.253125                      : 0.253125
0.759375                      : 0.3796875                     : 0.3796875
1.1390625                     : 0.56953125                    : 0.56953125
1.70859375                    : 0.854296875                   : 0.854296875
2.562890625                   : 1.2814453125                  : 1.2814453125
3.8443359375                  : 1.92216796875                 : 1.92216796875
5.76650390625----------------------------------------------------------------
1000000000000000000000000000.1: 500000000000000000000000000.05: 0.05
1000000000000000000000000000.2: 500000000000000000000000000.1 : 0.1
1000000000000000000000000000.3: 500000000000000000000000000.15: 0.15
1000000000000000000000000000.5: 500000000000000000000000000.25: 0.25
1000000000000000000000000000.8: 500000000000000000000000000.4 : 0.4
1000000000000000000000000001.2: 500000000000000000000000000.6 : 0.6
1000000000000000000000000001.8: 500000000000000000000000000.9 : 0.9
1000000000000000000000000002.7: 500000000000000000000000001.35: 1.35
1000000000000000000000000004.1: 500000000000000000000000002.05: 2.05
1000000000000000000000000006.2: 500000000000000000000000003.1 : 3.1
9.3--------------------------------------------------------------------------
2000000000000000000000000000.1: 1000000000000000000000000000.1: 0.1
2000000000000000000000000000.2: 1000000000000000000000000000.1: 0.1
2000000000000000000000000000.3: 1000000000000000000000000000.2: 0.2
2000000000000000000000000000.5: 1000000000000000000000000000.3: 0.3
2000000000000000000000000000.8: 1000000000000000000000000000.4: 0.4
2000000000000000000000000001.2: 1000000000000000000000000000.6: 0.6
2000000000000000000000000001.8: 1000000000000000000000000000.9: 0.9
2000000000000000000000000002.7: 1000000000000000000000000001.4: 1.4
2000000000000000000000000004.1: 1000000000000000000000000002.1: 2.1
2000000000000000000000000006.2: 1000000000000000000000000003.1: 3.1
9.3--------------------------------------------------------------------------
      

上述結果中第一組是準确值,沒有發生舍入誤差。後兩組的舍入情況不同,但最終結果居然一樣。至于為什麼進行這樣的舍入,請參閱我的上一篇随筆。

在 Windows 中調試

在 Windows 的 .NET Framework 中編譯和運作,也輸出一大堆調試資訊:

D:\work> DecimalSumDebug 10
0.1                           : 0.05                          : 0.05
0.15                          : 0.075                         : 0.075
0.225                         : 0.1125                        : 0.1125
0.3375                        : 0.16875                       : 0.16875
0.50625                       : 0.253125                      : 0.253125
0.759375                      : 0.3796875                     : 0.3796875
1.1390625                     : 0.56953125                    : 0.56953125
1.70859375                    : 0.854296875                   : 0.854296875
2.562890625                   : 1.2814453125                  : 1.2814453125
3.8443359375                  : 1.92216796875                 : 1.92216796875
5.76650390625----------------------------------------------------------------
1000000000000000000000000000.1: 500000000000000000000000000.05: 0.05
1000000000000000000000000000.2: 500000000000000000000000000.1 : 0.1
1000000000000000000000000000.3: 500000000000000000000000000.15: 0.15
1000000000000000000000000000.4: 500000000000000000000000000.2 : 0.2
1000000000000000000000000000.6: 500000000000000000000000000.3 : 0.3
1000000000000000000000000000.9: 500000000000000000000000000.45: 0.45
1000000000000000000000000001.4: 500000000000000000000000000.7 : 0.7
1000000000000000000000000002.1: 500000000000000000000000001.05: 1.05
1000000000000000000000000003.2: 500000000000000000000000001.6 : 1.6
1000000000000000000000000004.8: 500000000000000000000000002.4 : 2.4
7.2--------------------------------------------------------------------------
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
2000000000000000000000000000.1: 1000000000000000000000000000  : 0
0.1--------------------------------------------------------------------------
      

同樣,第一組的結果是準确值,後兩組有着不同的舍入誤差。第二組的舍入誤差比 Linux 中的要好。第三組就很糟糕了,被累加的值直接被舍入到 0 了。

System.Double 資料類型的情況

把前面的測試程式稍加修改,就可用于測試 double 資料類型的舍入誤差:

1 using System;
 2  
 3 static class DoubleSumTester
 4 {
 5   static void Main(string[] args)
 6   {
 7     try
 8     {
 9       var e = (args.Length > 0) ? int.Parse(args[0]) : 15;
10       for (var i = 0; i < 3; i++) Console.WriteLine(F(10, i, e));
11     }
12     catch (Exception ex) { Console.WriteLine(ex.Message); }
13   }
14     
15   static double F(int n, int k, int e)
16   {
17     var z = 0.1 + k * Math.Pow(10, e);
18     var w = (long)z / 2.0;
19     while (n-- > 0) z += z / 2 - w;
20     return z - w * 2;
21   }
22 }      

這時把循環次數固定為 10,指令行參數指定影響舍入誤差的整數部分的 10 的幂指數。運作結果如下所示:

work$ dmcs DoubleSumTester.cs
work$ mono DoubleSumTester.exe 0
5.76650390625
5.76650390625
5.76650390625002
work$ mono DoubleSumTester.exe 1
5.76650390625
5.76650390624993
5.76650390625
work$ mono DoubleSumTester.exe 4
5.76650390625
5.76650390631767
5.76650390624854
work$ mono DoubleSumTester.exe 10
5.76650390625
5.76657485961914
5.76650238037109
work$ mono DoubleSumTester.exe 12
5.76650390625
5.761962890625
5.7666015625
work$ mono DoubleSumTester.exe 14
5.76650390625
5.6875
5.0625
work$ mono DoubleSumTester.exe 15
5.76650390625
9
0
work$ mono DoubleSumTester.exe 16
5.76650390625
0
0
      

這次,Linux 和 Windows 中的運作結果是相同的。同樣,每次運作的第一行是準确值,其餘兩行是不同的舍入誤差形成的結果。

參考資料

  1. 部落格園:淺談 System.Decimal 結構

繼續閱讀