费马数乘积形式与平方和形式转换研究

周利荣 胡天磊

摘要：该文分析了费马数素因子乘积形式转换为平方和形式的算法，费马数平方和形式转换为素因子乘积形式的算法。利用转换算法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.

【通联编辑：谢媛媛】