高精度運算)

高精度運算(轉)

目錄

整數

整數表示

進位轉換

四則運算

快速乘法

一般原理

Karatsuba乘法

Toom-Cook乘法

FFT乘法

我們所熟知的科學計算一般就是指數值計算,數值計算是計算數學的一個主要部分,它研究用電腦求解各種數學問題的數值計算方法及其理論與軟體實現[1]

,
關於數值計算的研究在電腦被發明之前就已經有了相當的基礎,它涉及到的內容包括函數的數值逼近,數值微分與數值積分,非線性方程數值解,數值線性代數,
常微分方程與偏微分方程數值解等.數值計算中處理的對象並不僅僅是數值,還包括由數值構成的簡單資料結構,例如一般的多項式,無窮級數,矩陣等,數值計算
處理問題的一般方法是通過數學推導將問題化歸到這些數學對象的運算上.

作為應用數學,數值計算的主要目標是解決來自於生產實踐的工程學問題.與此同時,數學工作者做數學研究本身也是一種生產實踐,數學研究過程中同樣會
產生許多問題,與工程學問題不同,這些問題往往只能用抽象的符號來表達,僅用數值計算的方法是不易解決的,對於這類問題解決方案的研究逐漸形成了應用數學
的一個新的分支,為了與數值計算相區別,常常稱之為符號計算.類似的,我們可以給符號計算下一個簡單的定義:符號計算是一門研究用電腦求解各種數學問題
的符號計算方法及其理論與軟體實現的科學,它是數學(家)的計算數學.符號計算中處理的資料和結果都是符號,這種符號可以是字母,公式,也可以是數.與數
值計算不同的是,數是作為一種符號出現在符號計算中的,這就要求關於數的運算應該是絕對精確的,我們接下來就要討論關於數的高精度運算.

整數

在基於硬體的整數指令中,電腦能夠處理的整數是有界的,在目前典型的電腦中整數的溢出界都不超過
,而符號計算中常常需要處理更大的整數,例如階乘,斐波拉契數列這樣簡單的數論Function Compute,另一個不平凡的例子是所謂的中間表示膨脹[2]

,例如採用Euclid演算法計算兩個整係數多項式的最大公約數時,即使輸入的兩個多項式和輸出的最大公約數多項式都具有較小的係數,計算過程中的中間結果仍然很可能出現非常大的係數.設
F&=7x^7+2x^6-3x^5-3x^3+x+5,//
G&=9x^5-3x^4-4x^2+7x+7,
/end{align*}" width="580" height="50">
在計算過程中將有理數化為整數,我們將得到如下的多項式序列
&1890x^4-4572x^3-6930x^2-846x+4527,//
&294168996x^3+257191200x^2-20614662x-142937946,//
&-103685278369841305200x^2-32576054233115610000x//
&+122463167842311670000,//
&2956790833503849546789342057565207098291763520000x//
&+555325261806247996966034784074025291687620160000,//
&1092074685733031219201041602791259862659169966184593803518//
&602418777140682884334769647063543607737698426880000000000,
/end{align*}" width="580" height="235">
最後的那個整數達到了118位.除此之外,高精度浮點數的表示和運算也是直接依賴於高精度整數的.

整數表示

為了簡單起見,在下面的討論中,如果沒有明確指出,我們所說的整數都是正整數.關於高精度整數的表示方法,一個簡單的想法是選定進位
,用整數的
進位表示(向量)來表示整數,為了更好的理解這一點,首先對整數的進位表示作一個嚴格的定義.

定義1
設
,
,稱記號
為
進位整數,如果
具有的整數值由如下規則確定:

1. 
   ;
2. 
   .

即

事實上,給定進位
後,整數的進位表示可以看成由有限集
的直和之並形成的無窮維向量空間(類似於一維多項式空間)到整數集的同構.也就是說,每個整數都可以唯一地表示成一個
進位整數,這就是進位轉換的基礎,具體的計算公式見下面的定理.

定理1
給定整數
,存在非負整數
使得
,並且存在唯一的
位
進位整數滿足
,其中

從電腦科學的角度來看,整數的進位表示包含著資料結構的思想,它將高精度整數所包含的資訊用一列較小的單元來共同儲存和表示.在現有的符號計算軟體中,高精度整數典型的儲存格式[3]

如下:
/begin{tabular}{|c|p{8em}|}
/hline
$/pm e$ & /hfil $d$ /hfil//
/hline
/end{tabular}
/end{center}" width="580" height="27">
其中
就是整數的
進位表示(向量),
是
進位表示的位元,即這個整數在
進位下是多少位元.進位
在軟體設計之初就要確定下來,高精度整數四則運算的時間隨著整數的位元增加而增加,所以在不超過機器精度整數溢出界的範圍內,應該選取盡量大的基數
,以便充分利用基於硬體的整數指令.在現有的軟體中,常常取
為10或2的方冪,前者不必在輸入輸出時進行進位轉換,適用於對輸入輸出速度要求特別高的場合;後者用於計算時有著天然的優勢,在後面將會看到,當進位
為2的方冪時,很容易通過移位完成一部分乘除法運算.

進位轉換

高精度整數輸入輸出時,常常需要做進位轉換,即給定正整數
的
進位表示
,需要得到
在
進位下的表示.計算方法有如下兩種:

• 在
  進位下除以
  .設
  的
  進位表示為
  ,根據進位表示的定義,在
  進位下依次計算出
• 在
  進位下乘以
  .設
  進位整數
  ,利用Horner法則,在
  進位下依次計算
  ,
  從
  到0.最後得到
  即為
  在
  進位下的表示.

當轉換極長的數時,最方便的是由轉換數字塊開始,這些數字塊可以用單精確度技術來處理,而後用簡單的多精度技術把這些塊組合在一起[3]

.以第一種方法為例,反覆地用
除它,於是得到
的
進位表示,然後對於
進位表示的每一位通過單精確度操作給出
位
進位數字.

四則運算

本章的開頭時曾經強調過,數是作為一種符號出現在符號計算中的,這就要求關於數的運算應該是絕對精確的.在機器精度的
範圍內,現有的電腦可以輕鬆的完成整數的四則運算,這主要得益於電腦設計師們的工作.而接下來要討論的高精度整數四則運算,其基本原理其實和基於硬體
的整數指令是一致的,事實上,這就是阿拉伯計數法的發明者們很早就發展出的一整套藉助紙筆進行的四則運算理論,我們從小就開始學習熟練地利用它們來操作整
數.需要特別指出的是,這裡將要提到的加減乘除與作為群運算的加減乘除是不同的,因為加減乘除作為群運算是一個抽象的二元映射,而我們現在要做的其實是基
於這些群運算和上一節中進位表示所做的"表面"工作,說得具體一些,就是輸入兩個數的
進位表示,輸出群運算結果的
進位表示,關鍵點集中在求出
進位表示上.

整數加法

首先來看高精度整數的加法,為了簡單起見,假定相加的兩個數的位元相同,如果位元不同首先通過添零補齊.注意,這樣的假定不會影響演算法的效率.

演算法1
(整數加法)

輸入:整數
,
.

輸出:
的和
.

定義輔助序列
按遞推公式依次求出

容易看出,這裡定義的輔助序列
其實就是加到第
位的進位,由於不用輸出各位上的進位,演算法實際實現時可以只用一個變數來儲存
.儘管我們對這種筆式加法早已爛熟於胸,但是演算法的有效性還是嚴格地證明一下比較好.

證明
根據進位表示和
的定義,我們有
a_k+b_k+d_k&=([{a/over B^k}]-B[{a/over B^{k+1}}])+([{b/over B^k}]-B[{b/over B^{k+1}}])//
&+([{a+b/over B^k}]-[{a/over B^k}]-[{b/over B^k}])//
&=[{a+b/over B^k}]-B[{a/over B^{k+1}}]-B[{b/over B^{k+1}}],
/end{align*}" width="580" height="116">
由此推出
(a_k+b_k+d_k)/bmod B&=[{a+b/over B^k}]/mod B=c_k,//
[{a_k+b_k+d_k/over B}]&=[{a+b/over B^{k+1}}]-[{a/over B^{k+1}}]-[{b/over B^{k+1}}]=d_{k+1}.
/end{align*}" width="580" height="75">
證畢.
□

現在我們將這個演算法改寫成適合於編程的更加緊湊的形式.

演算法2
(整數加法)

輸入:兩個
位整數
,
.

輸出:
的和
,
或
.

1. 令
   .
2. 令
   ,
   .
3. 
   .若
   ,轉到
   ;否則令
   ,返回
   .

整數減法

高精度整數減法的方法與加法完全類似,我們直接給出演算法.

演算法3
(整數減法)

輸入:整數
,
.

輸出:
的差
.

定義輔助序列
按遞推公式依次求出

輔助序列
對應於到第
位的借位.

同樣地,也可以將這個演算法改寫成適合於編程的更加緊湊的形式.

演算法4
(整數減法)

輸入:兩個
位整數
,
,
.

輸出:
的差
.

1. 令
   .
2. 令
   ,
   .
3. 
   .若
   ,轉到2;否則返回
   .

如果高精度整數有求補操作,那麼仿照機器精度整數補碼錶示的思想,還可以通過補碼加法來做減法:用
來表示
位
進位數
的補碼,則

高精度整數的加減法演算法實際上就是對我們熟悉的十進位筆算加減法的精確化描述並推廣到
進位記數系統後得到的,用
表示用它們計算
位
進位整數和或差時所需要的
進位個位元加減法的次數,可以證明
,
這和必要的輸入輸出操作相比是都是線性層級的,因此有理由認為這就是理想的演算法.儘管如此,在後面有關多項式模算術的討論中我們將看到,作為多項式模算術
的一個直接應用,利用有限域上的整數模算術可以把高精度整數加減法分解成在多個處理器上並行的獨立任務,這樣一來它就會比經典演算法要快.

整數乘法

接下來討論一位元乘多位元的乘法,這裡用到的方法和加減法仍然是相似的,同樣可以直接給出演算法.

演算法5
(一位元乘多位元)

輸入:整數
,
.

輸出:
的積
.

定義輔助序列
按遞推公式依次求出

對於多位元乘多位元的情形,如
,則
每一個
都是一位元乘多位元,而乘上
只需在右邊添上
個0,或者說左移
位,於是可以通過
次一位元乘多位元,
次加法和添零(移位)求出
來.除此之外,在後面的"快速乘法"一節中將會看到,對於多位元乘多位元,我們還有更快的演算法.

同樣地,也可以將這個演算法改寫成適合於編程的更加緊湊的形式.

演算法6
(整數乘法)

輸入:
位整數
,
位整數
.

輸出:
的積
.

1. 令
   ,
   .
2. 若
   ,則
   ,轉到6.
3. 令
   ,
   .
4. 令
   ,
   ,
   .
5. 
   .若
   ,轉到4;否則令
   .
6. 
   .若
   ,轉到2;否則返回
   .

注1
如果
不清零,將有
.

商為一位元的除法

在四則運算的筆算方法中,除法可能是最複雜的,因為除法需要試商,試商包含著猜測的成分,而不能形成有效演算法.機器
精度整數除法也要試商,但得益於二進位表示,對它們而言商只有0,1兩種情況,只需要比較大小就可以求出.為了有效地計算高精度整數的除法,我們需要將試
商的過程演算法化[3]

,首先來看商為一位元的除法,即假定商
滿足
,筆算除法的經驗告訴我們,被除數和除數的最高几位元對試商是很重要的,我們常常憑藉目測(口算)最高几位元就能基本上把商給確定下來.

定理2
設整數
,則商
滿足不等式
證明
首先
故
,不等式的前半部分得證.

而由
可以得到
/left[/frac{a_{s+1}B+a_s}{b_s}/right]+1&/ge/frac{a_{s+1}B+a_s+1}{b_s}//
&=/frac{a_{s+1}B^{s+1}+(a_s+1)B^s}{b_sB^s}//
&>/frac{a}{b},
/end{align*}" width="580" height="122">
故
,證畢.
□

記
則
是商
一個很好的上界,下面將看到它已經很接近商的真值:
比商的真值至多大2.

定理3
若
的最高位不小於
,即
,則證明
用反證法.假設
,則
因為
所以
,推出
因此
,與題設條件矛盾,證畢.
□

上面只考慮了被除數的頭兩位與除數的首位,如果允許做更精細的考察,譬如說考慮被除數的頭三位與除數的頭兩位,我們還可以進一步接近商的真值.

定理4
設
,若
,則證明
將
代入不等式中得
化簡得到
a&/ge a_{s+1}B^2+a_sB+a_{s-1}//
&/ge q(b_sB+b_{s-1})//
&>q(b-B^{s-1})
/end{align*}" width="580" height="81">
因為
,所以
,證畢.
□

如果定理中的條件不滿足,將
減一,這時我們總可以說
比商的真值至多大1.綜合以上這些想法,現在可以寫出商為一位元的除法的詳細過程了.

演算法7
(商為一位元的除法)

輸入:整數
,
.

輸出:
的商
.

1. 若
   ,
   同乘以
   .
2. 計算
3. 計算
   ,如果
   ,
   減一.
4. 計算餘數
   ,如果
   ,返回
   ,否則返回
   .

整數除法

注意到
在演算法的最後一步也同時被算出,以此為基礎可以寫出一般的商為多位元的除法演算法.

演算法8
(整數除法)

輸入:整數
,
.

輸出:
的商
.

定義輔助序列
按遞推公式依次求出

為了利用快速乘法,還可以利用牛頓迭代來計算整數除法,具體來說就是先利用數值計算中浮點數的牛頓迭代求出除數倒數具有一定精度的近似值,然後利用整數乘法將被除數乘上去,當除數的位元較多時,這種方法是很有效.

快速乘法

一般原理

高精度整數是電腦代數系統的內建基本類型,四則運算的快慢對系統的效能好壞有著決定性的影響,目前最著名的高精度整數運算庫是GNU的GMP[4]

,
許多著名的電腦代數系統都是基於GMP構建的.加法和減法的複雜度關於整數位元是線性,考慮到輸入輸出的複雜度關於整數位元也是線性,因此從演算法上
來看加減法已經達到了複雜度的下界.高精度除法總可以通過Newton迭代歸結為高精度乘法,所以四則運算的主要問題集中在乘法上,設
為乘數的位元,就目前已知的情況而言,不同乘法演算法的複雜度可以從平凡的
(普通乘法),
(Karatsuba乘法),
(Toom-3乘法),
(複數域上的FFT)直到
(有
限域上的FFT),最後一個已經相當接近線性複雜度了.複雜度較低的演算法往往有較大的常數因子,例如當兩個相乘數都很小時,普通乘法反而是最快的,實用中
常常將這些不同的乘法演算法結合起來使用,每次做乘法時都根據相乘兩數的大小動態地選擇具體採用哪一種演算法,選擇時用到的參數一般通過實驗來確定.

多項式點值表示

在前面的利用一位元乘多位元來做整數乘法時,已經提到過整數乘法還有更快的演算法.在介紹這些進階演算法之前,先來分析一下,筆算乘法的演算法是如此地直截了當,以至於如果單從整數的角度考慮,應該不可能有更快的演算法了.我們現在換一個角度,把目光投向多項式.

對於一個次數界為
的多項式
其係數標記法就是一個由係數組成的向量
,這是我們習以為常的一種表示形式,並且是唯一的;而點值標記法[5]

是多項式在
個不同點處的點值對構成的集合
其中
.一個多項式可以有很多種不同的點值表示,這是因為每一組
都決定著一種點值表示.係數表示與點值表示之間是可以相互轉換的,係數表示到點值表示稱為多項式多點求值,點值表示到係數表示稱為多項式插值,關於這兩個主題的進階討論將留到與多項式相關的章節中,下面關於整數快速乘法的討論只需要這幾個概念就足夠了.

整數乘法與多項式

已知兩個次數界為
的多項式的係數表示分別為
設乘積的係數表示為
,那麼
這樣計算乘積時總共需要
次係數乘法.現在我們首先取定
個不同點
,然後利用多項式多點求值計算出
與
在這組點上的點值表示,如果
,那麼顯然有
,只需要
次係數乘法就可以計算出乘積