有限域上稀疏多元多項式插值算法*

唐 敏，鄧國強

桂林電子科技大學 數學與計算科學學院 廣西密碼學與信息安全重點實驗室，廣西 桂林 541004

1 引言

插值是函數逼近的重要方法，利用它可通過函數在有限個點處的取值情況，估算出函數在其他點處的近似值。在信號處理[1]、不確定性量化[2-4]、壓縮感知[5]、圖像處理[6-7]等領域需要有效地利用較少的函數值恢復一個函數，而這個函數往往具有稀疏的表示形式，稀疏插值能有效地處理此類問題。稀疏多元多項式插值是一種降低計算機代數算法時間復雜度的有效方法，是稀疏有理函數插值[8]及稀疏隱函數插值[9]的基本步驟，應用在結式計算[9]、多元多項式GCD計算[10-12]、方程求根、數值微分、數值積分、曲面的外形設計和有限元法及上述提及的諸多應用領域。因此提高稀疏多元多項式插值的效率十分有意義。

對于多元多項式，稠密插值所需的插值點個數為(d+1)n，其中d是變元的次數，n是變元的個數。研究工作表明，多元多項式插值算法的計算復雜度與變元個數可以不呈指數關系，而是與目標多項式的稀疏性相關，即是關于n、d、t和lbp的多項式函數，其中t是目標多項式的項數，素數p是有限域的特征。

1979年Zippel提出了第一個具有多項式時間復雜度的稀疏多元多項式插值算法[12]。該算法基于這樣的假設：如果在一個隨機的賦值點處多項式的取值為零，那么它就是一個零多項式。由于這個結論的成立是高概率的，因此Zippel算法是一個概率性的方法。在執行上，Zippel算法順序地對每個變元逐一插值，不能并行化，插值點個數為ndt，時間復雜度為O(ndt3)。1990年，Zippel改進了他的方法[13]，使用了形如的插值點，原算法中需要計算的一般線性系統變成了轉置的Vandermonde系統，求解時間從O(t3)降到O(t2)。由于Zippel算法是一種較好的具有多項式時間復雜度的方法，大部分計算機代數系統都內置了基于Zippel算法的多元GCD計算的實現，比如Macsyma、Magma和Mathematica。

1988年，Ben-Or和Tiwari提出了一個確定性的算法[14]，可以對以整數、有理數、實數或復數為系數的多元多項式進行稀疏插值。給定項數為t的多項式f的項數界T≥t，算法需要2t個插值點，即前n個素數 的 冪 次，0≤i≤2t-1。該方法與Zippel算法不同，不是逐個變元的插值，而是一次性地對f的各個單項式進行插值，而且可以并行化。主要的不足是插值點的存儲空間為tlbn位，當多項式規模較大時，由于中間表達式膨脹而導致計算速度非常慢。

2000年，Kaltofen等設計了一個Zippel和Tiwari方法相結合的混合算法，稱之為“racing算法”[15-16]。為了減少插值點個數，在Zippel算法中插值下一個變元的時候，同時執行Newton插值和單變元的Ben-Or/Tiwari插值，取其快者。

2010年，Javadi和Monagan提出了一個有限域上的稀疏多元多項式插值算法[17]，它也是一個概率算法。該方法是Ben-Or/Tiwari算法的一個改進，使用O(t)個插值點獨立完成每個變元的插值，共需O(nt)個插值點，并且可以并行化。為了準確計算每個變元在各個單項式中的準確次數，使用了二部圖完美匹配進行測試和驗證。該算法和Zippel算法及Kaltofen的racing算法進行了比較，在運行時間和插值點個數等方面具有一定的優勢。

2011年，Cuyt和Lee給出了稀疏有理函數插值算法[8]，其中涉及到單變元有理函數插值和基于Zippel算法及Tiwari算法的稀疏多元多項式插值。2016年，基于Cuyt的稀疏有理函數插值算法，Tang給出了稀疏隱函數插值算法[9]。在稀疏有理函數插值和隱函數插值算法中，稀疏多元多項式插值作為基本的子函數，對整個算法的效率起到關鍵作用。

本文提出了一個有限域上的稀疏多元多項式插值算法，它是Javadi/Monagan算法的改進，改進方案主要有兩點：一是消除了必須給定項數界T的限制，通過計算預先設計好的特定矩陣的行列式，得到f的準確項數；二是消除了必須給定次數界D的限制，通過構造輔助函數，利用概率法結合提前終止技術的Cauchy插值法，計算出f的準確次數，解決了Javadi和Monagan論文中提到的次數界D過高而導致高計算復雜度的問題。更進一步，改進算法無需給定T和D，對于實際問題在利用插值恢復或近似時更具實用性。

本文后續組織結構如下：第2章介紹了Javadi/Monagan算法的基本思想；第3章給出了改進的Javadi/Monagan算法；第4章分析了改進算法的時間復雜度；第5章中設計了三組數值實驗，對改進算法和Javadi/Monagan算法進行了比較，給出了實驗結果；第6章給出了一個應用實例。

2 Javadi/Monagan算法

本章給出稀疏多元多項式插值問題的形式化定義及Javadi/Monagan算法的思想。

2.1 稀疏多元多項式插值問題的形式化描述

令p是一個素數，f∈Zp[x1，x2，…，xn]是一個用黑盒表示的多元多項式：

稀疏多元多項式插值的目標是用盡可能少的插值點及較低的多項式時間復雜度算法還原多項式f。

2.2 Javadi/Monagan算法的思想

Javadi/Monagan算法可分成兩個步驟，首先確定各個單項式，然后確定系數。在第一個步驟中，逐一計算每個變元在各個單項式中的次數。第二個步驟需要求解一個Vandermonde線性代數系統。同時要求給定f的項數界T≥t及次數界D≥degf。

（1）確定單項式Mi，i=1，2，…，t。

③使用Berlekamp/Massey算法[18]生成n+1個關于z的單變元多項式，使得，其中i=0，1，…，2t-1,j=0，1，…，n。

④確定變元xk，k=1，2，…，n在f中的單項式Mi，i=1，2，…，t中的次數eij。詳見2.3節。

⑤由④可確定eij，i=1，2，…，t，j=1，2，…，n，則單項式

（2）計算系數。

通過求解線性方程組計算系數ai，i=1，2，…，t。令r1，r2，…，rt是Λ1(z)的根，且ri=Mi(α1，α2，…，αn)，回顧是輸入為時黑盒的輸出，因此有：

此為Vandermonde系統，有唯一解。

2.3 確定變元xk的次數

令R1={r1，r2，…，rt}和分別是Λ1和Λk+1的全部根構成的集合。令Dj包含單項式Mj中xk的所有可能的次數，即：

二部圖Gk定義如下：U和V是二部圖Gk中頂點個數為t的互補頂點子集，U和V中的結點表示R1和Rk中的元素，即ui∈U，vj∈V，分別用ri和標記，ui和vj之間有一條權值為dij的邊當且僅當

引理1[17]變元xk在所有單項式中的次數可唯一確定當且僅當二部圖Gk有唯一的完美匹配。

以下給出為計算xk的次數，構造的二部圖Gk沒有唯一的完美匹配時的解決方案。

隨機選擇不同于α1，α2，…，αn+1的元素令，使用 Berlekamp/Massey算法，以v′i=f(β′i)，0≤i≤ 2t-1為輸入，生成多項式Λ′k+1(z)。令是Λ′k+1(z)的根構成的集合，G′k是通過Λ1(z)和Λ′k+1(z)構造的二部圖。

定義1[17]二部圖G′k與Gk的交集定義如下：與G′k具有相同的結點，ri和之間有權值為dij的邊當且僅當在二部圖Gk中，ri連接且在二部圖G′k中，ri連接權值均為dij。

引理2[17]令eij=degxj(Mi)，在二部圖中，結點ri和之間有邊相連，且權值為eij。

定理1[17]令的交集有唯一的完美匹配。

引理1、引理2和定理1的證明參考文獻[17]。由定理1，如果Gk沒有唯一的完美匹配，可通過增加2t個插值點，構造二部圖來確定變元xk在單項式Mj(1≤j≤t)中的次數。

3 改進的Javadi/Monagan算法

本章給出對Javadi/Monagan算法進行改進的兩個主要策略。3.1節和3.2節分別描述了如何只利用黑盒的輸入輸出（即有限個插值點），計算多項式f的準確項數t和準確次數d，從而消除了原算法必須預先給定這兩個指標的上界的限制，同時提高了實用性及降低了計算時間復雜度(分析過程見第4章)。

3.1 f的項數t的計算

除了Javadi/Monagan算法，許多插值算法都要求給定目標多項式f的項數界T，本文給出一個準確計算多項式f的項數t的方法。該方法由定理2保證。

定理2[14]令V是元素為(V)ij=vi+j-2的矩陣。Vl是由V的前l行前l列組成的方陣。如果t是多項式f的準確項數，即f中非零單項式的個數，那么：

根據定理2，通過下述過程可計算f的準確項數：

例如：黑盒多元多項式f定義為：令p=1 009，隨機生成中的3個整數，比如α1=66，α2=12，α3=3。

因此多項式f的準確項數為5。

3.2 f的次數d的計算

在2.3節構造二部圖的過程中，計算單項式Mj中xk的所有可能的次數Dj時，需要檢測的次數為D+1。因此Javadi/Monagan算法在給定的次數界過高時，會導致高計算復雜度。本節利用單變元有理函數插值計算f的準確次數d，使得在確定每個變元的次數時，時間復雜度達到最低。方法描述如下：

（1）構造輔助有理函數H。

通過在f中引入齊次變元z，構造f為分子，g為分母的輔助有理函數H。其中f(x1，x2，…，xn)為黑盒多項式，g為隨機生成的關于變元z的一次多項式。若視F為關于變元z的單變元多項式，則系數為關于x1，x2，…，xn的多元多項式，注意到系數多項式與z是齊次的。

易證f的全次數與F中z的最高次相同。

（2）對H進行單變元有理函數插值。

任取n元組 (α1，α2，…，αn)，對H進行關于變元z的單變元有理函數插值，可得到：

（3）H(α1，α2，…，αn，z)的分子f(α1z，α2z，…，αnz)關于z的最高次即為黑盒多項式f的全次數。

假設f(α1z，α2z，…，αnz)和g(z)有非平凡的公因子，即g(z)|f(α1z，α2z，…，αnz)，也就是說，g(z)的根也是f(α1z，α2z，…，αnz)的根。因為g(z)是隨機生成的，而f(α1z，α2z，…，αnz)的根是有窮的，所以f(α1z，α2z，…，αnz)和g(z)高概率互素，從而f(α1z，α2z，…，αnz)關于z的次數即為黑盒多項式f的全次數。

令 (α1，α2)=(1，1)，g=5z+3 ，則：

利用單變元有理函數插值可得到H(1，1，z)=(z4+2z2)/(5z+3)，則f的準確次數為H的分子中z的最高次，即為4。

接下來給出單變元有理函數插值的實現過程。對于單變元有理函數f/g∈Κ(Z)，應用文獻[19]提出的概率方法，并結合提前終止技術的Cauchy插值來恢復f/g。來源于文獻[19-20]的引理給出了單變元有理函數插值的理論基礎。

引理3[20]令Κ是任意域，F(Z)，G(Z)，H(Z)∈Κ(Z)且 gcd(F，G)=1。令是非負整數，且1。令ik，f/g∈Κ(Z)不必是Κ中互不相同的元素，滿足：

對于l≥ 2 ，令hl(Z)，ql(Z)∈Κ(Z)分別是Euclidean多項式余項序列中的第l次余項和商。

對于l≥ 2 ，令wl(Z)，gl(Z)∈Κ(Z)是擴展Euclidean方法中的乘子wlh0+glh1=hl，即：

那么存在下標j≥ 1滿足，并且對下標j有：

根據引理3，如果給定單變元有理函數的次數上界，使用擴展Euclidean算法可計算有理函數的分子和分母。在文獻[19]中，提出了基于引理3的單變元有理函數插值算法。為確定準確的f/g的次數，可對k從1迭代到d+e+1，其中d、e分別是f和g的次數。

基于概率法和提前終止技術的Cauchy插值法恢復單變元有理函數，算法描述如下：

算法1單變元有理函數插值算法

其 中 Interpolate(i1，i2，…，ik;H(i1)，H(i2)，…，H(ik)) 表示給定k個點 (i1，H(i1))，(i2，H(i2))，…，(ik，H(ik))的自變量和函數值，利用Newton插值或Lagrange插值得到次數為k-1的單變元多項式；Rem表示余式，Quo表示商。

黑盒多項式f的次數即為deg(h1)。

3.3 改進算法

本節給出改進的Javadi/Monagan算法的偽代碼。

算法2改進的Javadi/Monagan算法(有限域上稀疏多元多項式插值算法)

以下給出改進的稀疏多元多項式插值算法的一個實例。黑盒多元多項式f定義為：

令p=101，變元個數n=3。重構f的過程如下：

（2）計算出多項式f的準確項數t=5。

（4）計算Λk(z)，k=1，2，…，4:

令vi=B(βi)，i=0，1，…，2t-1，當k＞1時，用替換使用Berlekamp/Massey算法生成由序列v0，v1，…，v2t-1得到的Λk∈Zp(z)：

（5）計算出多項式f的次數d=5。

（6）確定degxk(Mi),1≤i≤t,k=1，2，…，n。

以第1個變元x1為例，計算Λ1的根：

計算Λ2的根：

測試rj×(αn+1/α1)i，0≤i≤d，j=1，2，…，t,得到：

構造二部圖G1，如圖1所示。

Fig.1 Bipartite graphG1圖1 二部圖G1

因為二部圖G1沒有唯一的完美匹配，所以生成α5=54。

計算Λ5(z)=z5+100z4+73z3+85z2+51z+94及Λ5的根。構造二部圖G1′，如圖2所示。

Fig.2 Bipartite graphG1′圖2 二部圖G1′

構造Gk和Gk′的交集，如圖3所示。因此，x1在Mi,1≤i≤t中的次數分別為0,0,2,2,0。同理，計算x2和x3在Mi中的次數，分別為0,0,1,2,1和0,5,1,1,2。則

Fig.3 Bipartite graph圖3 二部圖

4 時間復雜度分析

本章討論有限域上改進的稀疏多元多項式插值算法（算法2）的時間復雜度。需要考慮的主要因素有：

（1）計算多項式f的項數t，時間復雜度為O(t2)。

（2）計算多項式f的次數d，時間復雜度為O(d)。

（3）調用Berlekamp/Massey算法計算Λk至多n+2次，每次調用的復雜度為O(t2)。

（4）使用文獻[13]中的技術求解Vandermonde系統，時間復雜度為O(t2)。

（5）求解Λk的根，時間復雜度為O(t2lbp)。

（6）計算變元xk的次數，時間復雜度為O(tlbt+dtlbt)。

（7）構造二部圖Gk，時間復雜度為O(dt2)，求Gk和G′k的交集，時間復雜度為O(dtlbd)。

綜上所述，改進算法的時間復雜度為：

從時間復雜度分析可以看出，如果給定的項數界T和次數界D過高，會導致高計算復雜度，而改進算法花費較少的代價計算出多項式f的準確項數t及次數d，達到了此類型算法在后續操作上的最低時間復雜度。

5 數值實驗

本章對改進算法和Javadi/Monagan算法進行性能比較。兩種算法的編程環境均為Maple15，程序運行的硬件環境為Intel?CoreTMi7 2.20 GHz處理器和4.00 GB內存，操作系統為Windows 7。注意兩個程序都是順序執行，在確定變元x1，x2，…，xn的次數時未并行化。

本章給出了3組測試集的性能比較結果，使用的多項式都是隨機生成的，比較的對象為CPU時間。黑盒中的多元多項式系數取自Zp，其中p=100 003。

測試集1本組測試集為6個包含3個變元的多元多項式。第i個多項式（1≤i≤6）使用如下的Maple命令隨機生成：

其中，第i個多項式包含2i個非零項，d=20是多項式的準確次數。該組測試分別執行了改進算法和Javadi/Monagan算法，記錄每個多項式在兩種算法下的運行時間（單位：s），如表1所示。表頭第1列Ex.表示例子的編號；第2列t表示項數；第3列括號中的百分比表示改進算法計算項數和次數的總時間與整個算法的運行時間的比值；改進算法無需給定項數界和次數界，Javadi/Monagan算法執行時，分別給定的次數界為D=20，D=30和D=30，項數界為多項式的準確項數T=2i(1≤i≤6)。表2和表3的設置同表1。

Table 1 Performance comparison of improved algorithm and Javadi/Monagan algorithm(n=3)表1 改進算法與Javadi/Monagan算法的比較(n=3)

從表1可以看出，隨著i的增加，兩種算法的執行時間也隨之增加。在給定準確項數界T=2i和準確次數界D=20的情況下，Javadi/Monagan算法的執行時間優于改進算法，原因是改進算法利用3.1節和3.2節的方法分別計算出準確的項數和次數，不需預先給定這兩個數據。顯然，這兩者的計算需要耗費一定的時間，但從表1中的第3列括號中的百分比可以看出，這兩者的計算時間占整體運行時間的比值隨著次數的增加越來越小，因而優勢越來越明顯。另一方面，隨著次數界D的增加，Javadi/Monagan算法需要的執行時間也隨之增加，當D較大時，算法的時間復雜度較高，這就是Javadi和Monagan在論文中提到的壞次數界問題（bad degree bound）。改進算法則有效地避免了這一問題。

測試集2本組測試集為6個包含6個變元的多元多項式。第i個多項式(1≤i≤6)使用如下的Maple命令隨機生成：

表2給出了測試集2下改進算法和Javadi/Monagan算法的執行時間，表格的各項含義與測試集1相同。

測試集3本組測試集為6個包含12個變元的多元多項式。第i個多項式(1≤i≤6)使用如下的Maple命令隨機生成：

Table 2 Performance comparison of improved algorithm and Javadi/Monagan algorithm(n=6)表2 改進算法與Javadi/Monagan算法的比較(n=6)

表3給出了測試集3下改進算法和Javadi/Monagan算法的執行時間，表格的各項含義與測試集1相同。

Table 3 Performance comparison of improved algorithm and Javadi/Monagan algorithm(n=12)表3 改進算法與Javadi/Monagan算法的比較(n=12)

從3組測試集可以看出，改進算法在確定插值多項式f的項數和次數上，花費的時間占整個算法運行時間的比例非常小，多項式規模越大，比例越小。如果Javadi/Monagan算法執行時，給定的次數界D越高，運行時間越長，而改進算法不受此影響。

6 應用實例

本章給出改進的Javadi/Monagan稀疏插值算法在幾何問題上的一個應用。著名的Morley三等分定理描述如下：在任意三角形ABC中，對3個內角∠CAB、∠ABC和∠ACB分別進行三等分，相鄰邊交于三點P、Q、R，則三角形PQR必為一個正三角形，如圖4所示。令a=BC,b=CA,c=AB,x=PQ=QR=RP,求x和a、b、c之間的關系。

Fig.4 Morley triangle圖4 Morley三角形

根據題設中的幾何關系和相關定理可知（推導過程略）：

其中，R是三角形ABC的外接圓半徑，有如下方程組：

令x和a、b、c之間的關系為f(x，a，b，c)=0 ，如果利用結式等符號消元法，消去方程組（1）中的變元p、q、r，那么可得到x和a、b、c之間的關系f(x，a，b，c)。實際上，如果把a、b、c視為符號參數，消元過程因中間表達式膨脹而無法完成，若將a、b、c用實數α1、α2、α3進行替換，消元過程可順利進行，結果是實例化的關系式f(x，α1，α2，α3)，例如（假定模p=105）：

如果視f(x，α1，α2，α3)為關于變元x的單變元多項式，則x的各次冪的系數是關于變元a、b、c的多元多項式，使用稀疏多元多項式插值恢復系數多項式，即可得到x和a、b、c之間的關系。在此例中，x的各個冪次的系數多項式的項數和次數都無法估計，如果采用原始的Javadi/Monagan算法，可能因為項數界和次數界設置不正確而導致無法得出結果或計算復雜度過高。若采用改進的Javadi/Monagan算法，可準確計算出各個系數多項式的項數和次數，最終得到的x和a、b、c之間的關系式：

此例中，f(x，a，b，c)共有804項，其中x27的系數多項式的項數為45，次數為32。x25的系數多項式的項數為28，次數為34。其余系數多項式的指標不再贅述。

在結式消元、GCD計算、幾何組合優化、信號處理等問題上，很多情況下無法預知目標多項式的次數和項數，故改進算法針對這些問題更具實用性。

7 結束語

稀疏多元多項式插值是很多計算機代數算法的子函數，也廣泛應用于信號處理、壓縮感知、圖像處理等領域。在給定項數界T和次數界D較高的情況下，改進算法可有效降低Javadi/Monagan算法的計算復雜度。進一步，改進算法消除了傳統插值算法中必須預先給定插值多項式項數T和次數D的必要條件，更具實用性。

3.3節的算法2給出了有限域上改進的Javadi/Monagan稀疏多元多項式插值算法，其中兩個循環體（算法2中步驟4到步驟8及步驟10到步驟19）的計算任務可并行化處理，一是生成關于變元z的單變元多項式Λ1(z)，Λ2(z)，…，Λn+1(z)，可分為相互獨立的n+1個子任務，二是確定變元x1，x2，…，xn在各個單項式中的次數，可分為相互獨立的n個子任務。假設有q個計算內核的集群系統，每個內核可獨立承擔上述計算Λj(z)，j=1，2，…，n+1 的 (n+1)/q個子任務，以及計算變元xj，j=1，2，…，n在各個單項式Mi，i=1，2，…，t中的次數eij的n q個子任務，計算結果可保存在一個文件中，最終可確定黑盒多項式f中的各個單項式。在后續工作中，擬將改進算法進行并行化處理，以進一步提高稀疏多元多項式插值算法的效率。