費爾馬“二平方”素數
問題的提出除２這個非凡的素數外，所有的素數都可以分成兩類：第一類是被4除余1的素數，如5,13,17,2937,41;第二類是被4除余3的素數，如3,7,11,19,23,31。第一類素數都能表示成兩個整數的平方和（第二類則不能）。
例如：5=1-1+2*213=2*2+3*317=1*1+4*4 29=2*2+5*5這就是聞名的費爾馬“二平方”定理。有趣的是：上述等式右側的數有的又恰恰是兩個素數的平方，如13,29的表達式，我們起名叫作費爾馬“二平方”素數，即假如一個素數能夠表示成兩個素數的平方和的形式：Ｆ＝Ｘ＊Ｘ＋Ｙ ＊Ｙ （１）其中Ｆ、Ｘ、Ｙ 都是素數，它就是費爾馬“二平方”素數。
編程思路本文擬用ｃ 語言編程，求４２億之內的費爾馬“二平方”素數。假如按定義從左向右，先求一個素數Ｆ，然後再去找相應的素數Ｘ、Ｙ ，工作量重復太大。我們可以對上述公式進行分析：
１、左側Ｆ 是素數，它肯定是奇數，那麼右側兩式的和也應該是奇數，這樣Ｘ 和Ｙ 為一奇一偶，因為奇數的平方還是奇數，偶數的平方還是偶數。Ｘ、Ｙ 又要求是素數，而既是偶數又是素數的數只有一個，就是２。我們假定Ｘ＝２。所以（１）式可以簡化為：Ｆ＝２＊２＋Ｙ ＊Ｙ（２）也就是說，費爾馬“二平方”素數的表示形式是惟一的。
２、按式（２）由右向左，由小到大找素數Ｙ ，再計算出相應的Ｆ，判定其是否素數。
３、求出素數Ｙ 後將其保存起來，在判定其它數是否素數時可直接用已求出的素數去除，如此反復。
源程序
#include<math.h>
void main()
{
unsigned long i,j,a[10000],m,m1=3,m2=7,b=1,n=0,d=1,x=4000000000;
a[1]=2;
10:for(i=m1;i<=m2;i++,i++)
{
if(i%a[1]==0) goto 13;
for(j=2;j<=d-1;j++)
if(i%a[j]==0) goto 13;
a[b++]=i; m=i*i+4;
if(m>x) goto 14;
for(j=2;j<=b-2;j++)
if(m%a[j]==0) goto 13;
printf("%20lu=2*2+%5lu*%5lu",m,i,i);
if(++n%2==0) printf("
");
13:m1=m2+4; m2=a[++d]*a[d]-2;
goto 10;
14:printf("
total=%lu
",n);
}
結論
運行程序會發現，除“２９＝２＊２＋５＊５”以外，所有的費爾馬“二平方”素數個位數字都是３，相應Ｙ 的個位數字都是３或７。費爾馬“二平方”素數分布（修改程序中變量ｘ 的值得到）也很耐人尋味，請看下表（表中１０萬以內包含１萬以內，下同）：
范圍個數最大的一個的表達式
１萬１０９４１３＝２＊２＋９７＊９７
１０萬２０９７９７３＝２＊２＋３１３＊３１３
１００萬４２９９４０１３＝２＊２＋９９７＊９９７
１０００萬７６９２２３３７３＝２＊２＋３０３７＊３０３７
１億１８３９７７５２７７３＝２＊２＋９８８７＊９８８７
１０億４２７９９９００２４５３＝２＊２＋３１６０７＊３１６０７
２０億５５１１９８３１８８０９３＝２＊２＋４４５３３＊４４５３３
３０億６４１２９９３５１２３７３＝２＊２＋５４７１３＊５４７１３
４０億７１８３９７７４４６４９３＝２＊２＋６３０６７＊６３０６７
費爾馬“二平方”素數太少了，４０億內才７１８個，千萬分之二還不到呢。
隨著數的范圍的增大，似乎越來越稀少，但再往後永遠是這樣嗎？會不會在某個范圍內反而又稠密起來呢？
費爾馬“二平方”素數是無窮多個呢，還是有限多個呢？假如是有限個，又是多少個呢？最大的一個又是什麼數呢？
這些問題的證實可能很簡單，也許很復雜，真說不定會成為像“哥德巴赫猜想”那樣的謎呢！
----------------------------------------------------------------------
下面是作者原程序,因為是中文全角,上面的就改了一下原文放在下面對照:
＃ｉｎｃｌｕｄｅ″ｍａｔｈ ．ｈ″
ｍａｉｎ（）
｛ｕｎｓｉｇｎｅｄ ｌｏｎｇｉ ，ｊ ，ａ［１００００］，ｍ，ｍ１＝３，ｍ２＝７，ｂ ＝１，ｎ ＝０，ｄ ＝１，ｘ＝４０００００００００；
ａ［１］＝２；
ｌ０：ｆｏｒ （ｉ ＝ｍ１；ｉ ＜＝ｍ２；ｉ ＋＋，ｉ ＋＋）
｛ｉｆ （ｉ ％ａ［１］＝＝０）ｇｏｔｏ ｌ３；
ｆｏｒ （ｊ ＝２；ｊ ＜＝ｄ －１；ｊ ＋＋）
ｉｆ （ｉ ％ａ［ｊ］＝＝０）ｇｏｔｏ ｌ３；
ａ［ｂ ＋＋］＝ｉ ；ｍ＝ｉ ＊ｉ ＋４；
ｉｆ （ｍ＞ｘ）ｇｏｔｏ ｌ４；
ｆｏｒ （ｊ ＝２；ｊ ＜＝ｂ －２；ｊ ＋＋）
ｉｆ （ｍ％ａ［ｊ］＝＝０）ｇｏｔｏ ｌ３；
ｐｒｉｎｔｆ（″％２０ｌｕ ＝２＊２＋％５ｌｕ ＊％５ｌｕ″，ｍ，ｉ ，ｉ）；
ｉｆ （＋＋ｎ ％２＝＝０）ｐｒｉｎｔｆ（″＼ｎ″）；
ｌ３：；｝ｍ１＝ｍ２＋４；ｍ２＝ａ［＋＋ｄ］＊ａ［ｄ］－２；
ｇｏｔｏ ｌ０；
ｌ４：ｐｒｉｎｔｆ（″＼ｎｔｏｔａｌ ＝％ｌｕ＼ｎ″，ｎ）；
｝
