費馬數乘積形式與平方和形式轉換研究

周利榮 胡天磊

摘要：該文分析了費馬數素因子乘積形式轉換為平方和形式的算法，費馬數平方和形式轉換為素因子乘積形式的算法。利用轉換算法2在極短的時間內（0.02 s）將62位素數934616397153579777691635581996068965840512375416381885802 80321分解成平方和形式。綜合利用兩種算法部分分解費馬數F12。

關鍵詞：費馬數;費馬平方和定理;呂卡定理

中圖分類號：TP312? ? ? 文獻標識碼：A

文章編號：1009-3044（2021）28-0114-03

開放科學（資源服務）標識碼（OSID）：

1 背景

表達式Fn=[22n]+1稱為費馬數，當n取0、l、2、3、4時的值分別為3、5、17、257、65537，并均為素數。由此，費馬提出一個猜想：形如Fn=[22n]+l的數一定為素數。數學家歐拉發現，當n取5時，F5=4294967297=641*6700417不為素數。費馬素數除了被費馬本人所證實的那五個外竟然沒有再發現一個。費馬數問題促進素性判定算法和整數的分解算法的研究。

2 費馬數性質

有關定理及性質：

1）費馬平方和定理：任何一個4p+1型素數可表成兩個整數的平方和，且表示方法是唯一的。

2）呂卡定理：費馬數Fn=[22n]+1的素因數必具p=k.[2n+2]+l的形式，其中k為正整數[1]。

3）性質1：當n≥2時，費馬數Fn=[22n]+1的末位數為7[2]。

4）性質2：（[u2+v2]）*（[a2+b2]）=[（ua+vb）2] +（[va-ub）2]。

5）性質3：（[u2+v2]）*（[a2+b2]）=[（ua-vb）2] +（[va+ub）2]。

3 費馬數素因子乘積形式轉換成平方和形式

3.1 轉換的依據

由呂卡定理可知費馬數的素因數必具k.[2n+2]+l的形式，即為4p+1型的素數。由費馬平方和定理可知任何一個4p+1型素數可表成兩個整數的平方和，且表示方法是唯一的。由性質2：（[u2+v2]）*（[a2+b2]）= [（ua+vb）2] +（[va-ub）2]，可知整數的平方和的乘積仍為平方和。因此，費馬數一定可以轉換成整數的平方和形式。

3.2 轉換算法1（窮舉法）

輸入：費馬數的素因數pi（1≤i≤k）。

輸出：費馬數的平方和形式。

1）i=1;

2）tmp=pi ，t=[pi]-1;

3）tm=tmp-[t2];

4）如果tm 是完全平方數，轉（5），否則 t--，轉（3）;

5）保存pi的平方和形式;

6）i++，如果 i≤k，轉（2），否則轉（7）;

7）由性質2輸出費馬數的平方和形式，程序結束。

設（x，y）是滿足pi=x2 + y2的格點，則（y，x）也是滿足pi=x2 + y2的格點，即格點在第一象限的分布滿足對稱性，因此，t的搜索區間為[[pi/2]，[pi]-1]。

3.3 轉換算法2（利用素因子平方和之間的比例關系）

性質2、性質3的成立是顯而易見的。由性質2、性質3可知，兩因子費馬數Fn=p1*p2=（[u2+v2]）*（[a2+b2]）=[（ua+vb）2] +（[va-ub）2]=[（vb-ua）2] +（[va+ub）2]。即兩因子費馬數有兩種平方和形式。由于Fn=[22n]+l=[（22n-1）2+12]，因此，[（22n-1）2+12]必定是其中一種平方和形式。設u

推論 1：vb-ua =1。

推論 2：[uv]≈[ba]。

推論 3：[uu2+v2]≈[ba2+b2]。

推論 4：[vu2+v2]≈[aa2+b2]。

對于兩因子費馬數，有了上述推論，如果已將較小因子p1分解成平方和形式p1=（[u2+v2]）。利用推論4可得a的近似值，從a的近似值開始遞增將p2分解成平方和形式[a2+b2]，這樣極大地減少了搜索范圍。

對于兩因子費馬數Fn=p1*p2， p1

1）用窮舉法將較小因子p1分解成平方和形式p1=（[u2+v2]）;

2）計算c=v/[p1];

3）由p2=（[a2+b2]）及a/[p2]≈c將較大因子p2分解成平方和形式p2=（[a2+b2]）。

4）由性質2輸出費馬數的平方和形式。

例：F5=4294967297=641*6700417。p1=641，p2=6700417。窮舉法分解p1=641=（[42]+[252]），設u=4，v=25，計算c=25 /[641 ]=0.987，由推論4知：a/[p2]≈c=0.987得a≈2555，循環兩次即可分解6700417 =（[25562]+[4092]），設 a=2556，b=409。ua+vb= 4* 2556+25 * 409 =20449，va-ub=25*2556-4*409=62264， F5=（[204492]+[622642]）。

例：F8= 115792089237316195423570985008687907853269984665640564039457584007913129639937=1238926361552897 *93461639715357977769163558199606896584051237541638188580280321。p1=1238926361552897，p2=93461639715357977769163558199606896584051237541638188580280321。窮舉法分解p1=1238926361552897=（[242465592+255153042]）， u=24246559， v=25515304，計算c=25515304 /[1238926361552897 ]=0.7248998337342965585004679297962，由a/[p2]≈c得a≈7008009763344264886811252498144，循環兩次即可分解p2 =（[70080097633442648868112524981452]+[66595374368066614409021877834642]）。設a=7008009763344264886811252498145，b=6659537436806661440902187783464。ua+vb= 2424655 *7008009763344264886811252498145+25515304 * 6659537436806661440902187783464 =33984024439900551177939471112034026611，va-ub=25515304*7008009763344264886811252498145-24246559*6659537436806661440902187783464=17340632172455487023654788790090010704， F8=（[173406321724554870236547887900900107042]+[339840244399005511779394711120340266112]）。

3.4 算法的效率比較

算法2由于利用推論4，算法的時間主要用于第（1）步：用窮舉法將較小因子p1分解成平方和形式，（2）（3）（4）步所用時間極少（分解F5的第二個素因子p2 為平方和形式僅用0.0001s;分解F6的第二個素因子p2 為平方和形式僅用0.001s;分解F7的第二個素因子p2 為平方和形式僅用0.01s;分解F8的第二個素因子p2 （62位素數）為平方和形式僅用0.02s，是目前被分解的最大素數）。因此，算法2的效率遠優于算法1。上表中，算法1計算F7、F8較大因子p2平方和形式所需要時間=循環一次時間*循環次數估算得到。可見當n≥7時，算法1的效率極低，幾乎不可能。

4 費馬數平方和形式轉換成乘積形式（適用于兩因子費馬數）

4.1 轉換算法3（解方程組法）

主要利用性質2：（[u2+v2]）*（[a2+b2]）= [（ua+vb）2] +（[va-ub）2]。如果已知費馬數的分解形式Fn=x2 + y2，由性質2可知，存在u、v、a、b，使得[（ua+vb）2] +（[va-ub）2]= x2 + y2= （[u2+v2]）*（[a2+b2]），只要能求出u、v、a、b，則可將n分解成素因子乘積形式。

要求出u、v、a、b的值，必須有四個方程組成的方程組。而條件只有兩個方程：x = ua+vb，y = va-ub。但x，y具有x= x1+x2，y = y2-y1的 形式，因此，可設x= x1+x2，y = y2-y1，則x2 + y2= （x1+x2）2 + （y2 -y1）2。如果x1，x2，y1，y2滿足x1=[ ua，]x2= vb，y1=[ va]，y2=[ ub]，則計算u=gcd（x1，y2），v= gcd（x2，y1），則得n的一個因子（[u2+v2]），分解成功。

如果x1，x2，y1，y2滿足x1=[ ua，]x2= vb，y1=[ va]，y2=[ ub]，則x1x2=uavb，y1y2=vaub，因此y1y2=x1x2，則必有方程組：

y1y2=x1x2……①

y1-y2=y ……②

由于y已知，只要知道x1，x2，則可求y1，y2，進而求出u、v、a、b，則可將n分解成素因子乘積形式。由②得：y2=y+y1，代入①得：[y12]+y.y1-x1.x2=0，是一個一元二次方程，解方程得：y1=（-y+[y2-4x1x2]）/2，y2= y+（-y+[y2-4x1x2]）/2。

推論 5：當n≥2時，費馬數Fn= x2 + y2中x，y必有一個為奇數。

由性質1不難推出此結論。設x為奇數，如果x1=[x/2];x2=x-x1，則有x2-x1=1。由3.3節的分析可知：（ua-vb）=1，（[va+ub]）=[22n-1]，又有：ua+vb=x ，va-ub=y 。因此，x1，x2，y1，y2必滿足x1=[ ua，]x2= vb，y1=[ va]，y2=[ ub]，[y2-4x1x2]=[（va-ub）2]+4uavb=[（va+ub）2]是完全平方數。y1=（-y+[y2-4x1x2]）/2，y2= y+（-y+[y2-4x1x2]）/2必為整數。因此，算法只要循環一次即可。

算法3描述如下：

1）求出x1=[x/2];

2）計算x2=x-x1;

3）解方程組y1y2=x1x2，y1-y2=y;

4）y1、y2必為整數，由u=gcd（x1，y2），v= gcd（x2，y1）求出u，v，將n分解，程序結束。

例：F6= 18446744073709551617= x2 + y2=[14387937592+40468032562]，x=1438793759，y= 4046803256，由于x為奇數，令x1=719396879，x2=719396880， 解方程組y1y2=x1x2，y2-y1=y得y1=（-y+[y2-4x1x2]）/2=124082020， y2= y+y1=4170885276為整數。u= gcd（y2，x1）=89，v=gcd（y1，x2）=516。（[u2+v2]）=274177，18446744073709551617=274177*67280421310721。

4.2 轉換算法4（求最大公約數法）

由4.1節的分析可得到如下四個方程：

vb-ua =1? ? ?……①

va+ub=[22n-1]? ?……②

ua+vb=x? ? ?……③

va-ub=y? ? ?……④

①+③得：2vb= x+1，vb= （x+1）/2;②+④得：2va=[22n-1]+y，va=（[22n-1]+y）/2，v=gcd（（x+1）/2，（[22n-1]+y）/2）。

③-①得：2ua= x-1，ua= （x-1）/2;②-④得：2ub=[22n-1-]y，ub=（[22n-1-]y）/2，u=gcd（（x-1）/2，（[22n-1]-y）/2）。

算法4描述如下：

1）求出v=gcd（（x+1）/2，（[22n-1]+y）/2）;

2）求出u=gcd（（x-1）/2，（[22n-1]-y）/2）;

3）得Fn的一個因子（[u2+v2]）;

4）將Fn分解n=（[u2+v2]）*[Fn/（[u2+v2]）]，程序結束。

例：F7=340282366920938463463374607431768211457=

x2 + y2=[163823502215354644792+84794438579364025042]，

設x=16382350221535464479，y=[8479443857936402504]，

v=gcd（（x+1）/2，（[22n-1]+y）/2）=208648999，u= gcd（（x-1）/2，（[22n-1]-y）/2）=126945596。（[u2+v2]）=59649589127497217，340 2823669209384634633746074317682114577=59649589127497 217 * 5704689200685129054721。

4.3 算法的效率比較

算法3由于第3）步解方程組需要時間比較多，效率低于算法4，但效率仍然是極高的。從上表可知，對于兩因子費馬數，在已知費馬數平方和形式的情況下，將其轉換成素因子乘積形式的效率是極高的。

5 兩種轉換算法在分解多因子費馬數F12中的應用

5.1 素因子乘積形式轉換成平方和形式

第一步：利用呂卡定理得到費馬數F12的兩個小素因子。

由呂卡定理可知費馬數的素因數必具k.[2n+2]+l的形式，k從1開始搜索，當k=7時， p1=7*[214]+1=114689且p1| F12，得到F12的第1個素因子p1。將F12除以p1，k又從1開始搜索，當k=1588時，p2= 1588*[214]+1=26017793且p2|（F12/ p1），得到F12的第2個素因子p2。總耗時約為0.03s。

第二步：利用算法1將費馬數F12的兩個小素因子分解成平方和形式，p1=114689=[2602]+[2172]，p2=26017793=[20472]+[46722]。耗時約為1.9s。F12是多因子費馬數，推論1-4并不成立，算法2失效。

第三步：利用性質2、性質3將費馬數F12的兩個小素因子乘積C13=p1* p2=2983954661377分解成兩種平方和形式，2983954661377=[15460442]+[7705212]，x=770521，y=1546044;2983954661377=[16589192]+[4816042]，n=[481604]，m=1658919。耗時約為0.02s。

5.2 平方和形式轉換成素因子乘積形式

由性質2、性質3可得到如下四個方程：

ua-vb =n? ?……①

va+ub=m? ?……②

ua+vb=y? ? ……③

va-ub=x? ? ……④

③-①得：2vb= y-n，vb=（y-n）/2;②+④得：2va=m+x，va=（m+x）/2，v=gcd（（y-n）/2，（m+x）/2）。

①+③得：2ua= y+n，ua=（y+n）/2;②-④得：2ub=m[-]x，ub=（m[-]x）/2，u=gcd（（y+n）/2，（m-x）/2）。

將C13轉換成素因子乘積形式，可用如下算法。

算法5描述如下：

1）求出v=gcd（（y-n）/2，（m+x）/2）;

2）求出u=gcd（（y+n）/2，（m-x）/2）;

3）得Fn的一個因子（[u2+v2]）;

4）將Fn分解n=（[u2+v2]）*[Fn/（[u2+v2]）]，程序結束。

在已知2983954661377=[15460442]+[7705212]，x=770521，y=1546044;2983954661377=[16589192]+[4816042]，n=[481604]，m=1658919的情況下，v=gcd（（y-n）/2，（m +x）/2）=260，u=gcd（（y+n）/2，（m -x）/2）=217，p1=[u2+v2] =114689。

∴2983954661377=114689*26017793。

即將C13=2983954661377分解成素因子乘積形式，耗時約為0.02s。

5.3 研究兩種形式轉換的意義

費馬數的分解算法有：橢圓曲線分解算法（分解F10、F11）、連分數算法（分解F7）、數域篩算法（分解F9）、Pollard's rho分解算法（分解F8）[3-4]。由文獻[5]F12已經分解出六個素因子，還有一個1133位的末尾合數C1133未分解。即F12 = 114689 * 26017793 * 63766529 * 190274191361 *1256132134125569 * 568630647535356955169033410940867804839360742060818433 * C1133。

上述費馬數的分解算法：橢圓曲線分解算法、連分數算法、數域篩算法、Pollard's rho分解算法對C1133顯然是無效的。如果知道C1133的兩種平方和形式，由算法5可求出C1133的一個素因子p7，且效率極高，利用Miller–Rabin算法[6]、ECPP、APRCL、AKS[7]等算法判定C1133/ p7的素性，如果為C1133/ p7素數，則F12完全分解，否則求出C1133/ p7的兩種平方和形式，由算法5可求出C1133/ p7的一個素因子p8，如此循環直至C1133完全分解。

將C1133這樣的大整數分解為平方和形式尚無有效算法，如能像算法2一樣找出多因子費馬數的素因子平方和變量（u、v、a、b）之間的比例關系，由p1、p2、p3、p4、p5、p6求出C1133的平方和形式，是一種可能途徑。這為分解費馬數提供了新的思路和可能性。

C++中用unsigned _int64表示整數范圍為[0，[264]]，即最大整數為18446744073709551616，費馬F6= 18446744073709551617大于此最大整數。NTL庫是一個用于大整數運算的C++庫，可以用于任意長度的整數計算及整數的多項式算法，以上實驗結果均在win10系統、VS 2017軟件+NTL、intel 酷睿處理器i3? @3.7G hz /8GB內存條件下取得。

參考文獻：

[1] 賈耿華.關于費馬數的研究[D].成都：成都理工大學，2006.

[2] 劉培杰.費馬數[J].自然雜志，1991，13（8）：608-612.

[3] 周利榮，胡天磊.基于萊梅素數判定定理的安全素數構造算法[J].計算機工程與應用，2016，52（13）：152-156，182.

[4] 周利榮，胡天磊.Demytko素數構造算法優化及應用研究[J].電腦編程技巧與維護，2018（6）：54-59.

[5] Wilfrid Keller.Fermat numbers website data[EB/OL].[2020-12-20].http：//www.prothsearch.com/fermat.html.

[6] Ishmukhametov S T，Rubtsova R，Savelyev N.The error probability of the miller-Rabin primality test[J].Lobachevskii Journal of Mathematics，2018，39（7）：1010-1015.

[7] 魏成行.素性檢測算法研究及其在現代密碼學中的應用[D].濟南：山東大學，2009.

【通聯編輯：謝媛媛】