在上一篇隨筆“使用快速傅里葉變換計算大整數乘法”中，已經講述了使用快速傅里葉變換計算大整數乘法的原理。在這一篇隨筆中，我們就使用快速傅里葉變換來實現一個提供任意精度的算術運算的靜態類：BigArithmetic。

下面就是 BigArithmetic 類源程序代碼：?

??1?using?System;
??2?using?System.Diagnostics;
??3?
??4?namespace?Skyiv.Numeric
??5?{
??6???///?<summary>
??7???///?提供任意精度的算術運算的靜態類。使用快速傅里葉變換。
??8???///?本類對字節數組進行算術運算，字節數組以?100?為基。
??9???///?字節數組中第一個元素存儲的數字是最高有效位。
?10???///?</summary>
?11???static?class?BigArithmetic
?12???{
?13?????//=?C語言數值算法程序大全(第二版)，ISBN?7-5053-2931-6?/?TP?993
?14?????//=?Numerical?Recipes?in?C,?The?Art?of?Scientific?Computing,?Second?Edition
?15?????//=?Cambridge?University?Press?1988,?1992
?16?????//=?[美]?W.H.Press,?S.A.Teukolsky,?W.T.Vetterling,?B.P.Flannery?著
?17?????//=?傅祖蕓?趙梅娜?丁巖?等譯，傅祖蕓?校，電子工業出版社，1995年10月第一版
?18?
?19?????static?readonly?byte?Len?=?2;?//?字節數組的元素包含的十進制數字的個數
?20?????static?readonly?byte?Base?=?(byte)Math.Pow(10,?Len);?//?字節數組的基
?21?????static?readonly?byte?MaxValue?=?(byte)(Base?-?1);????//?字節數組的元素的最大值
?22?
?23?????//=?pp.431-432,?four1,?12.2?快速傅里葉變換(FFT)
?24?????///?<summary>
?25?????///?復函數的快速傅里葉變換
?26?????///?變量?nn?是復數據點的個數，實型數組?data[1..2*nn]?的實際界長是兩倍?nn，
?27?????///?而每個復數值占據了兩個相繼的存儲單元。nn?必須是?2?的整數冪
?28?????///?</summary>
?29?????///?<param?name="data">實型數組?data[1..2*nn]。注意，下標從?1?開始</param>
?30?????///?<param?name="isInverse">是否逆變換。注意:?逆變換未乘上歸一化因子?1/nn</param>
?31?????public?static?void?ComplexFFT(double[]?data,?bool?isInverse)
?32?????{
?33???????int?n?=?data.Length?-?1;?//?n?必須是?2?的正整數冪
?34???????int?nn?=?n?>>?1;?????????//?變量?nn?是復數據點的個數
?35???????for?(int?i?=?1,?j?=?1;?i?<?n;?i?+=?2)?//?這個循環實現位序顛倒
?36???????{
?37?????????if?(j?>?i)
?38?????????{
?39???????????Utility.Swap(ref?data[j],?ref?data[i]);
?40???????????Utility.Swap(ref?data[j?+?1],?ref?data[i?+?1]);
?41?????????}
?42?????????int?m?=?nn;
?43?????????for?(;?m?>=?2?&&?j?>?m;?m?>>=?1)?j?-=?m;
?44?????????j?+=?m;
?45???????}
?46???????for?(int?mmax?=?2,?istep?=?4;?n?>?mmax;?mmax?=?istep)?//?執行?log2(nn)?次外循環
?47???????{
?48?????????istep?=?mmax?<<?1;?//?下面是關于三角遞歸的初始賦值
?49?????????double?theta?=?(isInverse???-2?:?2)?*?Math.PI?/?mmax;
?50?????????double?wtemp?=?Math.Sin(0.5?*?theta);
?51?????????double?wpr?=?-2?*?wtemp?*?wtemp;
?52?????????double?wpi?=?Math.Sin(theta);
?53?????????double?wr?=?1;
?54?????????double?wi?=?0;
?55?????????for?(int?m?=?1;?m?<?mmax;?m?+=?2)
?56?????????{
?57???????????for?(int?i?=?m;?i?<=?n;?i?+=?istep)
?58???????????{
?59?????????????int?j?=?i?+?mmax;?//?下面是?Danielson-Lanczos?公式
?60?????????????double?tempr?=?wr?*?data[j]?-?wi?*?data[j?+?1];
?61?????????????double?tempi?=?wr?*?data[j?+?1]?+?wi?*?data[j];
?62?????????????data[j]?=?data[i]?-?tempr;
?63?????????????data[j?+?1]?=?data[i?+?1]?-?tempi;
?64?????????????data[i]?+=?tempr;
?65?????????????data[i?+?1]?+=?tempi;
?66???????????}
?67???????????wr?=?(wtemp?=?wr)?*?wpr?-?wi?*?wpi?+?wr;?//?三角遞歸
?68???????????wi?=?wi?*?wpr?+?wtemp?*?wpi?+?wi;
?69?????????}
?70???????}
?71?????}
?72?

該類的第一個方法是 ComplexFFT，計算復函數的快速傅里葉變換。注意，ComplexFFT 并沒有使用復數(不象 C++，C# 也沒有提供復數)，而是讓每個復數值占據實型數組的兩個相繼的存儲單元。還有，要求輸入的復數據點的個數必須是 2 的整數冪。該方法也能夠計算復函數的快速傅里葉逆變換。

該程序的算法是使用 1942 年 Danielson 和 Lanczos 證明的引理：一個界長為 N 的離散傅里葉變換可以重新寫成兩個界長各為 N/2 的離散傅里葉變換之和。在算法的第一部分，將數據整理成位序顛倒的次序。而在第二部分，有一個執行 log₂N 次的外循環。它依次計算界長為 2, 4, 8, ..., N 的變換。對于這一過程的每一步來說，為了履行 Danielson-Lanczos 引理，有兩個嵌套的內循環，其涉及到已計算的子變換和每個變換的元素。通過限制外部調用正弦和余弦到外層循環，可以使運算更有效，在外層循環中只要調用它們 log₂N 次。倍角的正弦和余弦的計算是通過內循環中簡單的遞歸關系進行的，如下所示：

cos(θ + δ) = cosθ - [ α cosθ? +? βsinθ ]
sin(θ + δ) = sinθ - [ α sinθ? -? βcosθ ]

其中 α, β 是預先計算的系數：α = 2 sin²(δ/2)，? β = sinδ 。

?

?73?????//=?pp.436,?realft,?12.3.2?單個實函數的?FFT
?74?????///?<summary>
?75?????///?單個實函數的快速傅里葉變換
?76?????///?計算一組?n?個實值數據點的傅里葉變換。用復傅里葉變換的正半頻率替換這些數據，
?77?????///?它存儲在數組?data[1..n]?中。復變換的第一個和最后一個分量的實數值分別返回
?78?????///?單元?data[1]?和?data[2]?中。n?必須是?2?的冪次。這個程序也能計算復數據數組
?79?????///?的逆變換，只要該數組是實值數據的變換(在這種情況下，其結果必須乘以?1/n)即可。
?80?????///?</summary>
?81?????///?<param?name="data">實型數組?data[1..n]。注意，下標從?1?開始</param>
?82?????///?<param?name="isInverse">是否逆變換。注意:?逆變換未乘上歸一化因子?1/n</param>
?83?????public?static?void?RealFFT(double[]?data,?bool?isInverse)
?84?????{
?85???????int?n?=?data.Length?-?1;?//?n?必須是?2?的整數冪
?86???????if?(!isInverse)?ComplexFFT(data,?isInverse);?//?此處是正向變換
?87???????double?theta?=?(isInverse???-2?:?2)?*?Math.PI?/?n;?//?遞歸的初始賦值
?88???????double?wtemp?=?Math.Sin(0.5?*?theta);
?89???????double?wpr?=?-2?*?wtemp?*?wtemp;
?90???????double?wpi?=?Math.Sin(theta);
?91???????double?wr?=?1?+?wpr;
?92???????double?wi?=?wpi;
?93???????double?c1?=?0.5;
?94???????double?c2?=?isInverse???0.5?:?-0.5;
?95???????int?n3?=?n?+?3;
?96???????int?n4?=?n?>>?2;
?97???????for?(int?i?=?2;?i?<=?n4;?i++)
?98???????{
?99?????????int?i1?=?i?+?i?-?1,?i2?=?i1?+?1,?i3?=?n3?-?i2,?i4?=?i3?+?1;
100?????????double?h1r?=?c1?*?(data[i1]?+?data[i3]);?//?兩個分離變換是
101?????????double?h1i?=?c1?*?(data[i2]?-?data[i4]);?//?從?data?中分離出來
102?????????double?h2r?=?-c2?*?(data[i2]?+?data[i4]);
103?????????double?h2i?=?c2?*?(data[i1]?-?data[i3]);
104?????????data[i1]?=?h1r?+?wr?*?h2r?-?wi?*?h2i;?//?此處重新組合以形成
105?????????data[i2]?=?h1i?+?wr?*?h2i?+?wi?*?h2r;?//?原始實型數據的真實變換
106?????????data[i3]?=?h1r?-?wr?*?h2r?+?wi?*?h2i;
107?????????data[i4]?=?-h1i?+?wr?*?h2i?+?wi?*?h2r;
108?????????wr?=?(wtemp?=?wr)?*?wpr?-?wi?*?wpi?+?wr;?//?遞歸式
109?????????wi?=?wi?*?wpr?+?wtemp?*?wpi?+?wi;
110???????}
111???????double?tmp?=?data[1];
112???????if?(!isInverse)
113???????{
114?????????data[1]?=?tmp?+?data[2];?//?同時擠壓第一個和最后一個數據
115?????????data[2]?=?tmp?-?data[2];?//?使它們都在原始數組中
116???????}
117???????else
118???????{
119?????????data[1]?=?c1?*?(tmp?+?data[2]);
120?????????data[2]?=?c1?*?(tmp?-?data[2]);
121?????????ComplexFFT(data,?isInverse);?//?此處是逆變換
122???????}
123?????}
124?

第二個方法是 RealFFT，計算單個實函數的快速傅里葉變換。因為在很多情況下期望求快速傅里葉變換的數據由實值樣本組成。如果將這些樣本放入一個復型數組中，并令其所有虛部為零的話，從執行時間和對存儲需求二者來看，其效率是很低的。所以，我們將原始數據放入一個界長只有一半的復型數組中，其偶數項放入該數組的實部，奇數項放入它的虛部。然后調用 ComplexFFT 來計算快速傅里葉變換。

125?????///?<summary>
126?????///?比較?x[0..n-1]?和?y[0..n-1]
127?????///?</summary>
128?????///?<param?name="x">第一操作數?x[0..n-1]</param>
129?????///?<param?name="y">第二操作數?y[0..n-1]</param>
130?????///?<param?name="n">兩個操作數?x?和?y?的精度</param>
131?????///?<returns>比較結果：-1:小于?1:大于?0:等于</returns>
132?????public?static?int?Compare(byte[]?x,?byte[]?y,?int?n)
133?????{
134???????Debug.Assert(x.Length?>=?n?&&?y.Length?>=?n);
135???????for?(int?i?=?0;?i?<?n;?i++)
136?????????if?(x[i]?!=?y[i])
137???????????return?(x[i]?<?y[i])???-1?:?1;
138???????return?0;
139?????}
140?
141?????//=?pp.775,?mpneg,?20.6?任意精度的運算
142?????///?<summary>
143?????///?求補碼。注意，被操作數被修改。
144?????///?</summary>
145?????///?<param?name="data">被操作數?data[0..n-1]</param>
146?????///?<param?name="n">被操作數?data?的精度</param>
147?????///?<returns>被操作數的補碼?data[0..n-1]</returns>
148?????public?static?byte[]?Negative(byte[]?data,?int?n)
149?????{
150???????Debug.Assert(data.Length?>=?n);
151???????for?(int?k?=?Base,?i?=?n?-?1;?i?>=?0;?i--)
152?????????data[i]?=?(byte)((k?=?MaxValue?+?k?/?Base?-?data[i])?%?Base);
153???????return?data;
154?????}
155?
156?????//=?pp.774,?mpsub,?20.6?任意精度的運算
157?????///?<summary>
158?????///?減法。從?minuend[0..n-1]?中減去?subtrahend[0..n-1]，得到?difference[0..n-1]
159?????///?</summary>
160?????///?<param?name="difference">差?difference[0..n-1]</param>
161?????///?<param?name="minuend">被減數?minuend[0..n-1]</param>
162?????///?<param?name="subtrahend">減數?subtrahend[0..n-1]</param>
163?????///?<param?name="n">被減數?minuend?和減數?subtrahend?的精度</param>
164?????///?<returns>差?difference[0..n-1]</returns>
165?????public?static?byte[]?Subtract(byte[]?difference,?byte[]?minuend,?byte[]?subtrahend,?int?n)
166?????{
167???????Debug.Assert(minuend.Length?>=?n?&&?subtrahend.Length?>=?n?&&?difference.Length?>=?n);
168???????for?(int?k?=?Base,?i?=?n?-?1;?i?>=?0;?i--)
169?????????difference[i]?=?(byte)((k?=?MaxValue?+?k?/?Base?+?minuend[i]?-?subtrahend[i])?%?Base);
170???????return?difference;
171?????}
172?
173?????//=?pp.774,?mpadd,?20.6?任意精度的運算
174?????///?<summary>
175?????///?加法。augend[0..n-1]?與?addend[0..n-1]?相加，得到?sum[0..n]
176?????///?</summary>
177?????///?<param?name="sum">和?sum[0..n]</param>
178?????///?<param?name="augend">被加數?augend[0..n-1]</param>
179?????///?<param?name="addend">加數?addend[0..n-1]</param>
180?????///?<param?name="n">被加數?augend?和加數?addend?的精度</param>
181?????///?<returns>和?sum[0..n]</returns>
182?????public?static?byte[]?Add(byte[]?sum,?byte[]?augend,?byte[]?addend,?int?n)
183?????{
184???????Debug.Assert(augend.Length?>=?n?&&?addend.Length?>=?n?&&?sum.Length?>=?n?+?1);
185???????int?k?=?0;
186???????for?(int?i?=?n?-?1;?i?>=?0;?i--)
187?????????sum[i?+?1]?=?(byte)((k?=?k?/?Base?+?augend[i]?+?addend[i])?%?Base);
188???????sum[0]?+=?(byte)(k?/?Base);
189???????return?sum;
190?????}
191?
192?????//=?pp.774,?mpadd,?20.6?任意精度的運算
193?????///?<summary>
194?????///?捷加法。augend[0..n-1]?與整數?addend?相加，得到?sum[0..n]
195?????///?</summary>
196?????///?<param?name="sum">和?sum[0..n]</param>
197?????///?<param?name="augend">被加數?augend[0..n-1]</param>
198?????///?<param?name="n">被加數?augend?的精度</param>
199?????///?<param?name="addend">加數?addend</param>
200?????///?<returns>和?sum[0..n]</returns>
201?????public?static?byte[]?Add(byte[]?sum,?byte[]?augend,?int?n,?byte?addend)
202?????{
203???????Debug.Assert(augend.Length?>=?n?&&?sum.Length?>=?n?+?1);
204???????int?k?=?Base?*?addend;
205???????for?(int?i?=?n?-?1;?i?>=?0;?i--)
206?????????sum[i?+?1]?=?(byte)((k?=?k?/?Base?+?augend[i])?%?Base);
207???????sum[0]?+=?(byte)(k?/?Base);
208???????return?sum;
209?????}
210?
211?????//=?pp.775,?mpsdv,?20.6?任意精度的運算
212?????///?<summary>
213?????///?捷除法。dividend[0..n-1]?除以整數?divisor，得到?quotient[0..n-1]
214?????///?</summary>
215?????///?<param?name="quotient">商?quotient[0..n-1]</param>
216?????///?<param?name="dividend">被除數?dividend[0..n-1]</param>
217?????///?<param?name="n">被除數?dividend?的精度</param>
218?????///?<param?name="divisor">除數?divisor</param>
219?????///?<returns>商?quotient[0..n-1]</returns>
220?????public?static?byte[]?Divide(byte[]?quotient,?byte[]?dividend,?int?n,?byte?divisor)
221?????{
222???????Debug.Assert(quotient.Length?>=?n?&&?dividend.Length?>=?n);
223???????for?(int?r?=?0,?k?=?0,?i?=?0;?i?<?n;?i++,?r?=?k?%?divisor)
224?????????quotient[i]?=?(byte)((k?=?Base?*?r?+?dividend[i])?/?divisor);
225???????return?quotient;
226?????}
227?

接下來是比較、求補碼、減法、加法、捷加法、捷除法，都相當簡單，程序中已經有詳細的注釋了。

228?????//=?pp.776-777,?mpmul,?20.6?任意精度的運算
229?????///?<summary>
230?????///?乘法。multiplicand[0..n-1]?與?multiplier[0..m-1]?相乘，得到?product[0..n+m-1]
231?????///?</summary>
232?????///?<param?name="product">積?product[0..n+m-1]</param>
233?????///?<param?name="multiplicand">被乘數?multiplicand[0..n-1]</param>
234?????///?<param?name="n">被乘數?multiplicand?的精度</param>
235?????///?<param?name="multiplier">乘數?multiplier[0..m-1]</param>
236?????///?<param?name="m">乘數?multiplier?的精度</param>
237?????///?<returns>積?product[0..n+m-1]</returns>
238?????public?static?byte[]?Multiply(byte[]?product,?byte[]?multiplicand,?int?n,?byte[]?multiplier,?int?m)
239?????{
240???????int?mn?=?m?+?n,?nn?=?1;
241???????Debug.Assert(product.Length?>=?mn?&&?multiplicand.Length?>=?n?&&?multiplier.Length?>=?m);
242???????while?(nn?<?mn)?nn?<<=?1;?//?為變換找出最小可用的?2?的冪次
243???????double[]?a?=?new?double[nn?+?1],?b?=?new?double[nn?+?1];
244???????for?(int?i?=?0;?i?<?n;?i++)?a[i?+?1]?=?multiplicand[i];
245???????for?(int?i?=?0;?i?<?m;?i++)?b[i?+?1]?=?multiplier[i];
246???????RealFFT(a,?false);?//?執行卷積，首先求出二個傅里葉變換
247???????RealFFT(b,?false);
248???????b[1]?*=?a[1];?//?復數相乘的結果(實部和虛部)
249???????b[2]?*=?a[2];
250???????for?(int?i?=?3;?i?<=?nn;?i?+=?2)
251???????{
252?????????double?t?=?b[i];
253?????????b[i]?=?t?*?a[i]?-?b[i?+?1]?*?a[i?+?1];
254?????????b[i?+?1]?=?t?*?a[i?+?1]?+?b[i?+?1]?*?a[i];
255???????}
256???????RealFFT(b,?true);?//?進行傅里葉逆變換
257???????byte[]?bs?=?new?byte[nn?+?1];
258???????long?cy?=?0;?//?執行最后完成所有進位的過程
259???????for?(int?i?=?nn,?n2?=?nn?/?2;?i?>=?1;?i--)
260???????{
261?????????long?t?=?(long)(b[i]?/?n2?+?cy?+?0.5);
262?????????bs[i]?=?(byte)(t?%?Base);?//?原書中這句使用循環，有嚴重的性能問題
263?????????cy?=?t?/?Base;
264???????}
265???????if?(cy?>=?Base)?throw?new?OverflowException("FFT?Multiply");
266???????bs[0]?=?(byte)cy;
267???????Array.Copy(bs,?product,?n?+?m);
268???????return?product;
269?????}
270?

接下來的方法是 Multiply，乘法。其算法是使用 RealFFT 求被乘數和乘數的快速傅里葉變換，將結果相乘，然后進行傅里葉逆變換得到卷積，最后執行適當的進位。其原理已經在上一篇隨筆“使用快速傅里葉變換計算大整數乘法”中講述得很清楚了。

271?????//=?pp.778,?mpdiv,?20.6?任意精度的運算
272?????///?<summary>
273?????///?除法。dividend[0..n-1]?除以?divisor[0..m-1]，m?≤?n，
274?????///?得到：商?quotient[0..n-m]，余數?remainder[0..m-1]。
275?????///?</summary>
276?????///?<param?name="quotient">商?quotient[0..n-m]</param>
277?????///?<param?name="remainder">余數?remainder[0..m-1]</param>
278?????///?<param?name="dividend">被除數?dividend[0..n-1]</param>
279?????///?<param?name="n">被除數?dividend?的精度</param>
280?????///?<param?name="divisor">除數?divisor[0..m-1]</param>
281?????///?<param?name="m">除數?divisor?的精度</param>
282?????///?<returns>商?quotient[0..n-m]</returns>
283?????public?static?byte[]?DivRem(byte[]?quotient,?byte[]?remainder,?byte[]?dividend,?int?n,?byte[]?divisor,?int?m)
264?????{
285???????Debug.Assert(m?<=?n?&&?dividend.Length?>=?n?&&?divisor.Length?>=?m?&&?quotient.Length?>=?n?-?m?+?1?&&?remainder.Length?>=?m);
286???????int?MACC?=?3;
287???????byte[]?s?=?new?byte[n?+?MACC],?t?=?new?byte[n?-?m?+?MACC?+?n];
288???????Inverse(s,?n?-?m?+?MACC,?divisor,?m);?//?s?=?1?/?divisor
289???????Array.Copy(Multiply(t,?s,?n?-?m?+?MACC,?dividend,?n),?1,?quotient,?0,?n?-?m?+?1);?//?quotient?=?dividend?/?divisor
290???????Array.Copy(Multiply(t,?quotient,?n?-?m?+?1,?divisor,?m),?1,?s,?0,?n);?//??s?=?quotient?*?divisor
291???????Subtract(s,?dividend,?s,?n);?//?s?=?dividend?-?quotient?*?divisor
292???????Array.Copy(s,?n?-?m,?remainder,?0,?m);
293???????if?(Compare(remainder,?divisor,?m)?>=?0)?//?調整商和余數
294???????{
295?????????Subtract(remainder,?remainder,?divisor,?m);
296?????????Add(s,?quotient,?n?-?m?+?1,?1);
297?????????Array.Copy(s,?1,?quotient,?0,?n?-?m?+?1);
298???????}
299???????return?quotient;
300?????}
301?
302?????//=?pp.777?-?778,?mpinv,?20.6?任意精度的運算
303?????///?<summary>
304?????///?求倒數。
305?????///?</summary>
306?????///?<param?name="inverse">倒數?inverse[0..n-1]，在?inverse[0]?后有基數的小數點</param>
307?????///?<param?name="n">倒數?inverse?的精度</param>
308?????///?<param?name="data">被操作數?data[0..m-1]，data[0]?>?0，在?data[0]?后有基數的小數點</param>
309?????///?<param?name="m">被操作數?data?的精度</param>
310?????///?<returns>倒數?inverse[0..n-1]，在?inverse[0]?后有基數的小數點</returns>
311?????public?static?byte[]?Inverse(byte[]?inverse,?int?n,?byte[]?data,?int?m)
312?????{
313???????Debug.Assert(inverse.Length?>=?n?&&?data.Length?>=?m);
314???????InitialValue(inverse,?n,?data,?m,?false);
315???????if?(n?==?1)?return?inverse;
316???????byte[]?s?=?new?byte[n],??t?=?new?byte[n?+?n];
317???????for?(;?;?)?//?牛頓迭代法:?inverse?=?inverse?*?(?2?-?data?*?inverse?)??=>??inverse?=?1?/?data
318?????? {
319???????? Array.Copy(Multiply(t,?inverse,?n,?data,?m),?1,?s,?0,?n);?//?s?=?data?*?inverse
320???????? Negative(s,?n);?????????????????????????????????????????//?s?=?-(data?*?inverse)
321???????? s[0]?-=?(byte)(Base?-?2);?????????????????????????????//?s?=?2?-?data?*?inverse
322???????? Array.Copy(Multiply(t,?s,?n,?inverse,?n),?1,?inverse,?0,?n);?//?inverse?=?inverse?*?s
323?????????int?i?=?1;
324?????????for?(;?i?<?n?-?1?&&?s[i]?==?0;?i++)?;?//?判斷?s?的小數部分是否為零
325?????????if?(i?==?n?-?1)?return?inverse;?//?若?s?收斂到?1?則返回?inverse?=?1?/?data
326?????? }
327???? }
328?

?

接下來的方法是 DivRem，除法，以及 Inverse，求倒數。

除法用除數的倒數乘以被除數來計算，倒數值用牛頓迭代法：

U_i+1?= U_i?(2 - VU_i)

來計算，這導致 U_∞二次收斂于 1/V。實際上，許多超級計算機和 RISC 機都是使用這種迭代法來實現除法的。

要注意在 DivRem 中，求得的余數可能大于等于除數，此時商比實際的值要小一。所以在程序的 293 到 298 行調整商和余數。

?

329?????//=?pp.778-779,?mpsqrt,?20.6?任意精度的運算
330?????///?<summary>
331?????///?求平方根?sqrt，以及平方根的倒數?invSqrt。invSqrt?也可設為?sqrt，此時，invSqrt?也是平方根。
332?????///?</summary>
333?????///?<param?name="sqrt">平方根?sqrt[0..n-1]，在?sqrt[0]?后有基數的小數點</param>
334?????///?<param?name="invSqrt">平方根的倒數?invSqrt[0..n-1]，在?invSqrt[0]?后有基數的小數點</param>
335?????///?<param?name="n">平方根的精度</param>
336?????///?<param?name="data">被操作數?data[0..m-1]，data[0]?>?0，在?data[0]?后有基數的小數點</param>
337?????///?<param?name="m">被操作數?data?的精度</param>
338?????///?<returns>平方根?sqrt[0..n-1]，在?sqrt[0]?后有基數的小數點</returns>
339?????public?static?byte[]?Sqrt(byte[]?sqrt,?byte[]?invSqrt,?int?n,?byte[]?data,?int?m)
340???? {
341?????? Debug.Assert(sqrt.Length?>=?n?&&?invSqrt.Length?>=?n?&&?data.Length?>=?m);
342???????if?(n?<=?1)?throw?new?ArgumentOutOfRangeException("n",?"must?greater?than?1");
343?????? InitialValue(invSqrt,?n,?data,?m,?true);
344???????byte[]?s?=?new?byte[n],?t?=?new?byte[n?+?Math.Max(m,?n)];
345???????for?(;?;?)?//?invSqrt?=?invSqrt?*?(3?-?data?*?invSqrt?*?invSqrt)?/?2?=>?invSqrt?=?1?/?sqrt(data)
346?????? {
347???????? Array.Copy(Multiply(t,?invSqrt,?n,?invSqrt,?n),?1,?s,?0,?n);?//?s?=?invSqrt?*?invSqrt
348???????? Array.Copy(Multiply(t,?s,?n,?data,?m),?1,?s,?0,?n);???//?s?=?data?*?invSqrt?*?invSqrt
349???????? Negative(s,?n);?????????????????????????????????????//?s?=?-(data?*?invSqrt?*?invSqrt)
350???????? s[0]?-=?(byte)(Base?-?3);?????????????????????????//?s?=?3?-?data?*?invSqrt?*?invSqrt
351???????? Divide(s,?s,?n,?2);??????????????????????????????//?s?=?(3?-?data?*?invSqrt?*?invSqrt)?/?2
352???????? Array.Copy(Multiply(t,?s,?n,?invSqrt,?n),?1,?invSqrt,?0,?n);???//?invSqrt?=?invSqrt?*?s
353?????????int?i?=?1;
354?????????for?(;?i?<?n?-?1?&&?s[i]?==?0;?i++)?;?//?判斷?s?的小數部分是否為零
355?????????if?(i?<?n?-?1)?continue;?//?若?s?沒有收斂到?1?則繼續迭代
356??????? Array.Copy(Multiply(t,?invSqrt,?n,?data,?m),?1,?sqrt,?0,?n);?//?sqrt?=?invSqrt?*?data?=?sqrt(data)
357?????????return?sqrt;
358?????? }
359???? }
360?
361?????///?<summary>
362?????///?采用浮點運算以得到一個初始近似值?u[0..n-1]:?u?=?1?/?data?或者?u?=?1?/?sqrt(data)
363?????///?</summary>
364?????///?<param?name="u">初始近似值?u[0..n-1]</param>
365?????///?<param?name="n">所需的精度</param>
366?????///?<param?name="data">被操作數?data[0..m-1]</param>
367?????///?<param?name="m">被操作數?data?的精度</param>
368?????///?<param?name="isSqrt">是否求平方根</param>
369?????///?<returns>初始近似值?u[0..n-1]</returns>
370?????static?byte[]?InitialValue(byte[]?u,?int?n,?byte[]?data,?int?m,?bool?isSqrt)
371???? {
372?????? Debug.Assert(u.Length?>=?n?&&?data.Length?>=?m);
373???????int?scale?=?16?/?Len;?//?double?可達到?16?位有效數字
374???????double?fu?=?0;
375???????for?(int?i?=?Math.Min(scale,?m)?-?1;?i?>=?0;?i--)?fu?=?fu?/?Base?+?data[i];
376?????? fu?=?1?/?(isSqrt???Math.Sqrt(fu)?:?fu);
377???????for?(int?i?=?0;?i?<?Math.Min(scale?+?1,?n);?i++)
378?????? {
379?????????int?k?=?(int)fu;
380???????? u[i]?=?(byte)k;
381???????? fu?=?Base?*?(fu?-?k);
382?????? }
383???????return?u;
384???? }
385?? }
386?}

最后的方法是 Sqrt，求平方根。用牛頓迭代法計算平方根和除法類似。若：

U_i+1?= U_i?(3 - VU_i²) / 2

則 U_∞二次收斂于 1/√V，最后乘以 V 就得到√V。

?

在上一篇隨筆“使用快速傅里葉變換計算大整數乘法”中提到：

由于快速傅里葉變換是采用了浮點運算，因此我們需要足夠的精度，以使在出現舍入誤差時，結果中每個組成部分的準確整數值仍是可辨認的。長度為 N 的 B 進制數可產生大到 B²N 階的卷積分量。我們知道，雙精度浮點數的尾數是 53 個二進位，所以：

2 x log₂B + log₂N + 幾個 x log₂log₂N < 53

上式中左邊最后一項是為了快速傅里葉變換的舍入誤差。

我們假設上式中的“幾個”為“三個”，那么，經過簡單的計算，得到以下結果：

基數 B長度 N十進制數字個數

256	約 107	約 2.4 x 107
100	約 108	約 2 x 108
10	約 109	約 109

注意，基數 B 選取得越小，計算速度就越慢。

轉自：http://www.cnblogs.com/skyivben/archive/2008/07/25/1250891.html

轉載于:https://www.cnblogs.com/freeopen/p/5482948.html

總結

以上是生活随笔為你收集整理的使用快速傅里叶变换计算大整数乘法-代码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。