再谈 BigInteger - 使用快速傅里叶变换
我在“淺談 BigInteger”的隨筆中實現了一個 Skyiv.Numeric.BigInteger 類，那時乘法是使用常規的 O(N2) 的算法，所以比 .NET Framework 3.5 Base Class Library 中的 System.Numeric.BigInteger 類稍慢，后者的乘法是使用 Karatsuba 算法，其時間復雜度約為 O(N1.585)。
現在，我已經實現了一個提供任意精度的算術運算的靜態類: BigArithmetic，該類使用快速傅里葉變換計算大整數乘法，其時間復雜度降低到 O(N logN loglogN)。所以，Skyiv.Numeric.BigInteger 類也重新改寫調用 BigArithmetic 類的靜態方法來進行算術運算。
下面就是改寫后的 Skyiv.Numeric.BigInteger 類的源程序代碼：
??1?using?System;??2?using?System.Text;
??3?
??4?namespace?Skyiv.Numeric
??5?{
??6???sealed?class?BigInteger?:?IEquatable<BigInteger>,?IComparable<BigInteger>
??7???{
??8?????static?readonly?byte?Len?=?2;
??9?????static?readonly?byte?Base?=?(byte)Math.Pow(10,?Len);
?10?
?11?????sbyte?sign;??//?符號，取值：-1,?0,?1。
?12?????byte[]?data;?//?字節數組以?100?為基，字節數組中第一個元素存儲的數字是最高有效位。
?13?
因為 BigArithmetic 類是對字節數組進行算術運算，字節數組以 100 為基，為了方便調用 BigArithmetic 類的靜態方法進行算術運算，BigInteger 也改為使用以 100 為基的字節數組來存儲。原來是使用 109 為基的 List<int> 來存儲的。
我在改寫過程中曾經使用過以 108 為基的 List<int> 來存儲，在需要調用 BigArithmetic 類的靜態方法進行算術運算時將參數轉換為以 100 為基的字節數組，運算完成后再將運算結果轉換回來。這樣做的好處是加法和減法運算比較快(還使用原來的方法，不調用 BigArithmetic 類的靜態方法)，并且除了 operator *、DivRem 和 Sqrt 方法以外，其他所有的方法都不用改寫。不過經過綜合考慮，還是采用現在的方案。
?
?14?????BigInteger()?15?????{
?16?????}
?17?
?18?????BigInteger(long?x)
?19?????{
?20???????sign?=?(sbyte)((x?==?0)???0?:?((x?>?0)???1?:?-1));
?21???????data?=?new?byte[10];?//?long.MinValue?=?-9,223,372,036,854,775,808
?22???????ulong?z?=?(x?<?0)???(ulong)-x?:?(ulong)x;
?23???????for?(int?i?=?data.Length?-?1;?z?!=?0;?i--,?z?/=?Base)?data[i]?=?(byte)(z?%?Base);
?24???????Shrink();
?25?????}
?26?
?27?????BigInteger(BigInteger?x)
?28?????{
?29???????sign?=?x.sign;?//?x?!=?null
?30???????data?=?new?byte[x.data.Length];
?31???????Array.Copy(x.data,?data,?data.Length);
?32?????}
?33?
?34?????public?static?implicit?operator?BigInteger(long?x)
?35?????{
?36???????return?new?BigInteger(x);
?37?????}
?38?
私有的構造函數，和原來的大同小異。
?
?39?????public?static?BigInteger?Parse(string?s)?40?????{
?41???????if?(s?==?null)?return?null;
?42???????s?=?s.Trim().Replace(",",?null);
?43???????if?(s.Length?==?0)?return?0;
?44???????BigInteger?z?=?new?BigInteger();
?45???????z.sign?=?(sbyte)((s[0]?==?'-')???-1?:?1);
?46???????if?(s[0]?==?'-'?||?s[0]?==?'+')?s?=?s.Substring(1);
?47???????int?r?=?s.Length?%?Len;
?48???????z.data?=?new?byte[s.Length?/?Len?+?((r?!=?0)???1?:?0)];
?49???????int?i?=?0;
?50???????if?(r?!=?0)?z.data[i++]?=?byte.Parse(s.Substring(0,?r));
?51???????for?(;?i?<?z.data.Length;?i++,?r?+=?Len)?z.data[i]?=?byte.Parse(s.Substring(r,?Len));
?52???????z.Shrink();
?53???????return?z;
?54?????}
?55?
靜態的 Parse 方法，也和原來的大同小異。
程序的第 56 到 95 行以及 111 到 116 行是 Abs、Pow 方法以及單目 +、-、++ 和 -- 運算符，以及減法(-)運算符，和原來的相同。
?
?96?????public?static?BigInteger?operator?+(BigInteger?x,?BigInteger?y)?97 ??? {
?98?????? if?(x?==?null?||?y?==?null)?return?null;
?99?????? if?(x.AbsCompareTo(y)?<?0)?Utility.Swap(ref?x,?ref?y);
100?????? BigInteger?z?=?new?BigInteger();
101?????? z.sign?=?x.sign;
102?????? byte[]?bs?=?Utility.Expand(y.data,?x.data.Length);
103?????? bool?isAdd?=?x.sign?*?y.sign?==?1;
104 ????? z.data?=?new?byte[x.data.Length?+?(isAdd???1?:?0)];
105 ????? if?(isAdd)?BigArithmetic.Add(z.data,?x.data,?bs,?bs.Length);
106 ????? else?BigArithmetic.Subtract(z.data,?x.data,?bs,?bs.Length);
107 ????? z.Shrink();
108 ????? return?z;
109 ??? }
110
加法(+)運算符，調用 BigArithmetic 類的 Add 和 Subtract 方法。
118?????{
119???????if?(x?==?null?||?y?==?null)?return?null;
120???????if?(x.sign?*?y.sign?==?0)?return?0;
121???????BigInteger?z?=?new?BigInteger();
122???????z.sign?=?(sbyte)(x.sign?*?y.sign);
123???????z.data?=?new?byte[x.data.Length?+?y.data.Length];
124???????BigArithmetic.Multiply(z.data,?x.data,?x.data.Length,?y.data,?y.data.Length);
125???????z.Shrink();
126???????return?z;
127?????}
128?
乘法(*)運算，直接調用 BigArithmetic 類的 Multiply 方法。
程序的第 129 到 141 行是 / 和 % 運算符，和原來的相同。
?
142?????public?static?BigInteger?DivRem(BigInteger?dividend,?BigInteger?divisor,?out?BigInteger?remainder)143?????{
144???????remainder?=?null;
145???????if?(dividend?==?null?||?divisor?==?null)?return?null;
146???????if?(divisor.sign?==?0)?throw?new?DivideByZeroException();
147???????if?(dividend.AbsCompareTo(divisor)?<?0)
148???????{
149?????????remainder?=?new?BigInteger(dividend);
150?????????return?0;
151???????}
152???????BigInteger?quotient?=?new?BigInteger();
153???????remainder?=?new?BigInteger();
154???????quotient.data?=?new?byte[dividend.data.Length?-?divisor.data.Length?+?1];
155???????remainder.data?=?new?byte[divisor.data.Length];
156???????BigArithmetic.DivRem(quotient.data,?remainder.data,?dividend.data,
157?????????dividend.data.Length,?divisor.data,?divisor.data.Length);
158???????quotient.sign?=?(sbyte)(dividend.sign?*?divisor.sign);
159???????remainder.sign?=?dividend.sign;
160???????quotient.Shrink();
161???????remainder.Shrink();
162???????return?quotient;
163?????}
164?
靜態的 DevRem 方法，也是調用 BigArithmetic 類的 DivRem 方法。
?
165?????public?static?BigInteger?Sqrt(BigInteger?x)166?????{
167???????if?(x?==?null?||?x.sign?<?0)?return?null;
168???????if?(x.sign?==?0)?return?0;
169???????if?(x.data.Length?==?1)?return?new?BigInteger((long)Math.Sqrt(x.data[0]));
170???????BigInteger?z?=?new?BigInteger();
171???????z.sign?=?1;
172???????z.data?=?new?byte[x.data.Length?/?2?+?3];
173???????z.data?=?Adjust(BigArithmetic.Sqrt(z.data,?z.data,?z.data.Length,?x.data,?x.data.Length),?x.data.Length);
174???????z.Shrink();
175???????BigInteger?z1?=?z?+?1;?//?平方根有可能比實際小?1。
176???????return?(z1?*?z1?<=?x)???z1?:?z;
177?????}
178?
179?????static?byte[]?Adjust(byte[]?bs,?int?digits)
180?????{
181???????if?(bs[0]?>=?10)?throw?new?OverflowException("sqrt?adjust");
182???????byte[]?nbs?=?new?byte[(digits?+?1)?/?2];
183???????if?(digits?%?2?==?0)
184?????????for?(int?k?=?bs[0],?i?=?0;?i?<?nbs.Length;?i++,?k?=?bs[i]?%?10)
185???????????nbs[i]?=?(byte)(k?*?10?+?bs[i?+?1]?/?10);
186???????else?Array.Copy(bs,?nbs,?nbs.Length);
187???????return?nbs;
188?????}
189?
求平方根(Sqrt)，調用 BigArithmetic 類的 Sqrt 方法。但是要注意，BigArithmetic 類的 Sqrt 方法求得的平方根是小數，需要使用 Adjust 方法調整為整數。并且平方根有可能比實際的小 1，也需要判斷和調整。
?
190?????void?Shrink()191?????{
192???????int?i;
193???????for?(i?=?0;?i?<?data.Length;?i++)?if?(data[i]?!=?0)?break;
194???????if?(i?!=?0)
195???????{
196?????????byte[]?bs?=?new?byte[data.Length?-?i];
197?????????Array.Copy(data,?i,?bs,?0,?bs.Length);
198?????????data?=?bs;
199???????}
200???????if?(data.Length?==?0)?sign?=?0;
201?????}
202?
私有的 Shrink 方法用于在需要時收縮存儲數據的字節數組。
程序的第 203 到 238 行是各種邏輯運算符，和原來的完全一樣。
239?????public?override?string?ToString()240?????{
241???????StringBuilder?sb?=?new?StringBuilder();
242???????if?(sign?<?0)?sb.Append('-');
243???????sb.Append((data.Length?==?0)???0?:?(int)data[0]);
244???????for?(int?i?=?1;?i?<?data.Length;?i++)?sb.Append(data[i].ToString("D"?+?Len));
245???????return?sb.ToString();
246?????}
247?
從基類繼承的 ToString 方法，和原來的大同小異。
程序的第 248 到 274 行是從基類繼承的 GetHashCode、Equals 方法，以及用于實現接口的 Equals 和 CompareTo 方法，和原來的相同。
?
275?????int?AbsCompareTo(BigInteger?other)276?????{
277???????if?(data.Length?<?other.data.Length)?return?-1;
278???????if?(data.Length?>?other.data.Length)?return?1;
279???????return?BigArithmetic.Compare(data,?other.data,?data.Length);
280?????}
281???}
282?}
最后是私有的 AbsCompareTo 方法，調用 BigArithmetic 類的 Compare 方法。
下面是程序中用到的輔助類：
?1?using?System;?2?
?3?namespace?Skyiv
?4?{
?5???static?class?Utility
?6???{
?7?????public?static?T[]?Expand<T>(T[]?x,?int?n)
?8?????{
?9???????T[]?z?=?new?T[n];?//?assume?n?>=?x.Length
10???????Array.Copy(x,?0,?z,?n?-?x.Length,?x.Length);
11???????return?z;
12?????}
13?
14?????public?static?void?Swap<T>(ref?T?x,?ref?T?y)
15?????{
16???????T?z?=?x;
17???????x?=?y;
18???????y?=?z;
19?????}
20???}
21?}
?
經過運行測試程序，發現計算 256,142 位的乘積的運行時間從原來的 34.2497761 秒降低到 0.1749114 秒，速度提高了將近二百倍。如果計算的數字更大的話，速度將提高得更多。因為原來乘法的時間復雜度是 O(N2)。而現在乘法使用 FFT，時間復雜度是 O(N logN loglogN)。
我想，這個 BigInteger 類的運算速度已經是非常快的了。不知道有沒有其他 C# 的 BigInteger 類的運算速度更快。
?
本文中的所有源代碼可以到 https://bitbucket.org/ben.skyiv/biginteger?頁面下載。
也可以使用 hg clone http://bitbucket.org/ben.skyiv/biginteger/ 命令下載。
關于 hg ，請參閱 Mercurial 備忘錄。
2010-04-22 更新: 對 Skyiv.Numeric.BigInteger 的 GetHashCode、Parse 和 ToString 方法進行了優化。請參閱：再談 BigInteger - 優化
《新程序員》：云原生和全面數字化實踐50位技術專家共同創作，文字、視頻、音頻交互閱讀總結
以上是生活随笔為你收集整理的再谈 BigInteger - 使用快速傅里叶变换的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 三种运动让身高增长4-10cm
- 下一篇： 浮出水面