◆愛知県 Y.M.Ojisan さんからの解答。
【はじめに】
次のＨＰによれば、「ゴールドバッハ予想」同様本問題も未解決問題のようです。
http://www80.sakura.ne.jp/~aozora/suuron/node16.html
あるＨＰによれば「ゴールドバッハ予想」は４０兆まで計算機により確認されているそうです（年代不詳）。
そこでここではパソコンで何処まで確認できるか（目標は４０兆）に挑戦しました。
【解答】
０〜４０００億 および ４０兆〜４０兆１０００億の範囲で反例なし。
【予備調査】
まずパソコンで凝らずに計算できる１０億までを確認しました。
その結果、素数の存在確率が素数定理π（ｘ）＝ｘ／log(ｘ)に従うとして統計処理した結果と大差ないことを確認しました。
下記は減算側素数の最大値、平均値、標準偏差値の統計理論値との比較結果です。
２素数をＰ１，Ｐ２とし 凵≠o２−Ｐ１とします。
凾ｪ大きいとき それぞれの値は次のように近似できます。
（P1素数列に対しポアソン分布です）。
なお最大値はあと一個程度存在する可能性のある位置としました。
平均値 P1m=E[P1] は
P1m／log(P1m)=log()2/2(log()-1)
最大値 P1max は
P1max／log(P1max)= log(/2)log()2/2(log()-1)
標準偏差値 P1σ は
(P1σ+P1m)／log(P1σ+P1m)= log(/2)log()2/(log()-1)−２
から得られる値です。
以上の結果から４０兆まで調べる場合でも、Ｐ１は高々４１２４程度まで調べれば、例外が１個出る程度と予測されます。
【数値実験法】
１億(偶数5000万個)につき１０秒程度で検査できるソフトを作成することができました。
（PC=１．７Ｇ Ｎote）
方法は、
（１）４byte整数の配列上に素数判定ビットマップMap[i]を考えます。
Map[i]のjビットは奇数3+2*(j+32i)に対応しており、値が１なら素数であるとします。
このビットマップを４０兆までに対応して√（４０兆）までエラトステネスのふるい法で作成します。
（これはCPU２秒ほど）
(2)ここからＰ１に対応するビット列並びが逆順のP1mapを作成します。
前の予測から８０９５まで分(１２８バイト)だけで十分です。
(下図例は13まで)
(3)検査する偶数をＧ１〜Ｇ２とするとき、（Ｇ１＋３）〜（Ｇ２＋８０９５）の間の素数のビットマップP2mapをMapを用いてやはりエラトステネス法で作成します。
(4)偶数のマップGmapをクリアします。
Gmap[i]のjビットは偶数2*(j+32i)に対応します。
(5)P2mapをＰ２＝(Ｇ１+3)から順に調べて、１（素数）であればP1mapその位置を先頭にしてGmapにOrします。
これはGmapの(Ｐ２−Ｐ１)の位置をマークしていると言うことです。
(6)そしてＧ１〜Ｇ２に対応するbit位置が全部１なら例外なしになります。
(7)実際にはP1mapは0〜31bitずらしたものを予め全部用意しておき、計算時間を短縮しました。
この問題では検査すべきＰ１mapの巾が非常に狭くて済むので、この方法におけるＧ１を任意の位置から始めても、計算時間に無駄はほとんどありません。
そこで目標の４０兆あたりも容易に検査することができ、計算時間もさほど延びないことが確認できました。
従って、パソコンでも５０日ほど回せば確認することが可能です。
【数値実験結果】
今回は一夜分で、目標の１％の４０００億まで例外がなく、また目標の４０兆付近の０．２５％でも例外がないこと、およびＰ１の最大が４０００も超えなかったことを報告して終了します。
結果のうち、小区間におけるＰ１の最大値の分布を追加拡大したものを下図に示します。
予測を否定できない結果です。
なお、分かりにくいですが一応ソースを添付しておきます。
Ｃ＋＋ソース
【感想】
最初ＰＣでは無理かなと思いましたが、意外と高速に検査することができました。
このプログラムの場合、検査領域を分割して並列に計算することができるので、最近はやりのグリッドコンピューティング（ネット上につながるＰＣを使って並列に高速計算する方法）を、数学の部屋でボランティアを募って行えば、１日でもできそうな程度でした。