数学の最近のブログ記事

数値積分

user-pic
0

発散級数の和(その2)のところで、


の積分の数値計算を、台形公式により粗く、0.596364028767647ととりあえず計算してみたのだが、その後、黒川信重「ラマヌジャン ζの衝撃」を眺めていて、その88ページで以下の記載をみつけて唖然とした。私が計算した値と、誤差を考慮してもなにか本質的に違う。
オイラーは、"1-1!+2!-3!+4!-5!+..." = 0.59637164と求めています(全集の編者は0.59637255と修正)

それで、Wikiで、ボレルの総和法の記事を調べて、上記の値が0.59634736232319407...とあって、何が正しい値に近いのかいよいよわからなくなった。

像てい

そこで、手許にある、数値解析入門[増訂版]、山本哲朗著、サイエンス社を参考にして、台形公式より精度がよいとされるシンプソンの3/8公式を用いて上記積分の数値計算を再度行うことにした。

シンプソンの3/8公式では、以下のようにして、関数f(x)の数値計算を行う。
S = 3h(f(x0) + 3f(x1) + 3f(x2) + f(x3)/8 + 3h(f(x3) + 3f(x4) + 3f(x5) + f(x6)/8 +
...+ 3h(f(xN-3) + 3f(xN-2) + 3f(xN-1) + f(xN)/8

これのソースコードは、以下のとおりであり、また計算結果は、0.59634721063974127642で、0.596347までしか正しくない。ただ、「ラマヌジャン ζの衝撃」にあった上記数字は、おそらく正しくないとはいえる。それにしても、シンプソンの公式といいつつ、台形公式よりも必ずしも精度が上がらず、仕方なく刻み目hを小さくすると異様に計算時間がかかってしまうのだった。また、蛇足であるがいくつかの注意点も記す:

  • 積分範囲Wを仮にセットしたものの、実際は0から∞までの積分であり、計算ループの脱出条件を、一回前のループの計算値との差が、10-23以内とした。これを小さくするほど精度は上がる。
  • 精度により効いてくるのは、積分範囲の刻み目hの幅であるが、これをあまりに小さくすると計算時間が膨大になってしまう。そこでとりあえず、h = 10-10とした。計算時間は、だいたい30分くらいだった。
  • 上記のような小さい数のセットのため、TH_set(int n)という関数を用意した。これは、10-nの実数値を返す。
// シンプソンの3/8公式による積分計算

#include <stdio.h>
#include <math.h>

double f(double x);
double TH_set(int n);

int main()        {

int i = 0;
double W = 100.0;  // 積分範囲
double S = 0.0; // 積分値
double b = 0.0;  // 変数の値
double Sx; // 1つ前の値
double Threshold  = TH_set(23);  // 差分の閾値 10のマイナス23乗にセット
double h =  TH_set(10);  // 刻み幅10のマイナス10乗にセット


 S += 3.0*h*(f(b) + 3.0*f(b+h) + 3.0*f(b + 2.0*h) + f(b + 3.0*h))/8.0;
 b += 3.0*h;

do                      {
 Sx = S;
 S += 3.0*h*(f(b) + 3.0*f(b+h) + 3.0*f(b + 2.0*h) + f(b + 3.0*h))/8.0;
 b += 3.0*h;
 
  i++;
//  if ( i%100000 == 0 ) { printf("%6.20f\t%5.20f\n",S, S - Sx); i = 0; }
 
  if ( S - Sx < Threshold ) break;

} while ( b + 2*h < W );


 printf("%6.20f\n",S);
 return 1;
}

double f(double x)  {
 return exp(-x)/(1.0 + x);
}

double TH_set(int n) {
 int i;
 double TH = 1.0;

 for (i = 0; i < n; i++)  {
  TH /= 10.0;
 }
 
 return TH;

}

黒川信重「ラマヌジャン ζの衝撃」の発散級数の和の続きで、次のような式があった。

今更驚かないが、これもトリッキーで、左辺は普通に考えたら左辺は-∞になるはずであるが、右辺は意味不明な小数になっている。そこで、「ラマヌジャン ζの衝撃」に従い、説明を述べる。

まず、次のような関数Fnを考える。

そして、F0は、下記の計算で1になる。

また、F1も、下記の計算で1になる。

以上から帰納的に、下記の式が成立する。

よって、最初の式は、以下のように、Fnを使って書きあらわされる。このような和の式は、ボレルの総和法と呼ばれている。ラマヌジャンが大変好んだ式でもある。

そこで問題は、最後の積分の式であるが、これは台形公式により、数値積分してみる。すなわち、積分範囲は、[0,∞)であるが、exp(-x)/(1.0 + x)の分子が極めて急速に減衰するので、[0,100)で考え、これを10000刻みする。
これは直ちに結果が得られて、0.596364028767647となるので、一応所望の値である。

#include <stdio.h>
#include <math.h>

double f(double x);

int main()        {

double W = 100.0;  // 積分範囲
double h = W/10000.0;  // 刻み幅
double S = 0.0; // 積分値
double b = 0.0;  // 変数の値

// 普通に台形公式を使う

 S = h*f(b)/2.0;
 b += h;

while ( b < W - h )   {
  S += h*f(b);
  b += h;
}

 S += h*f(b)/2.0;
 
 printf("%6.15f\n",S);
 
 return 1;
}

double f(double x)  {
 return exp(-x)/(1.0 + x);
}

初等関数

user-pic
0

数学セミナーの2018年6月号は、「複素関数の質問箱」という特集で、最近リーマン予想のごく初歩のところから勉強し始めた私には願ったりの特集である。なぜなら、リーマン予想は、リーマンのゼータ関数の解析接続を前提とした理論なので、複素関数論を理解していることが必須であるところ、私は複素関数論なんてほとんどまともに勉強したことがないからである。

その「複素関数の質問箱」は、以下の七編からなる。

  1. 複素数ってそもそも何なの?
  2. 複素関数論中の実関数
  3. 複素微分とはなんですか、ふつうの微分とはどう違うのですか
  4. 複素積分とはなんですか、ふつうの積分とはどう違うのですか
  5. 正則関数の"綺麗"な性質はなぜ
  6. 複素対数とは何か/パラドックスをこえて
  7. 解析接続の意味と意義について

その「複素関数論中の実関数」と題する、中村弥生氏の記事を読んでいるうちに、以下のような初等関数の定義が書かれているのをみて、新鮮な驚きを感じた。

さて、代数関数、指数関数、三角関数、対数関数、およぴそれらの線形結合、積、商、合成関数、逆関数をとる操作によって得られる関数を初等関数といいます。初等関数のうち、特に代数関数でないものを初等超越関数といいます。

私は初等関数というのは、高校までで習う関数くらいに曖昧に理解していたので、この明快な定義に感動さえした。すなわち、代数関数、指数関数、三角関数、対数関数に積、商、合成関数、逆関数をとる操作を繰り返していくと、見かけ上相当に複雑な関数が得られる。それでも、それらは初等関数と呼ばれるのだな、と思った。

で、そこからさらに気づいたことがある。微分操作可能性である。それは、解析的な意味での微分可能性とは違う。要するに、決まった公式により式変形できるのである。なので、別の言い方をすると、微分操作的式変形を許すなら、初等関数なのである。

ところが初等関数を含む、常微分方程式、偏微分方程式、積分となると急に様子が変わってくる。それらの解は最早、一般的には初等関数とはならない。例えばアーベルは、代数関数の積分を考察して、楕円関数論を切り開いた。

厳密な定義は知らないが、初等関数でない関数の一部、すなわち理工学に応用される関数は特殊関数と呼ばれるのだと思う。例えば手許にある、井上正雄著「応用関数論」共立全書、昭和29年刊には、楕円関数、Γ関数、ルジャンドル関数、ベッセル関数が論じられており、これらが特殊関数の典型なのだろう。

この「応用関数論」で、最初の方に複素関数論の基礎の章がある。すなわち、関数というものは複素関数論で基礎づけられなくてはならず、なので「複素関数論中の実関数」に初等関数の定義も出てくるのだな、と思った。

黒川信重「ラマヌジャン ζの衝撃」現代数学社の発散級数の和のところを眺めていて、これは決して初めてみたわけではないが、何度みてもトリッキーである。典型的には、以下の式である。

普通に考えたら、左辺は-∞になるので、ありえないのだが、ともかくこの級数は、アーベルの総和法という技法を使って、以下のような多項式の級数に変換される。それにしても、アーベル大先生が、こんなまやかし?に加担していたとは(笑)

プロの数学者には、この式変形は自明なのだろうが、私のような弱小アマチュアには理解できないので、自分にもわかるように段階を追って行間を埋めてみた。まず、以下のように、1+xの逆数の級数展開を得る。

すると、(1+x)2の逆数の級数展開は、以下のとおりとなり、さらに変形される。

最後の式変形は、xkの係数に着目すると、(左側の括弧からのn=0の項)×(右側の括弧からのn=kの項)+(左側の括弧からのn=1の項)×(右側の括弧からのn=k-1の項)+・・・+(左側の括弧からのn=kの項)×(右側の括弧からのn=0の項)となって、そのような項の各々はどれも(-1)kであって、この項の数が(k+1)個あることに基づく。

よって、以下の式変形が成り立つ。この式の左辺は確かに、x → 1-0の極限で、1/4になる。

と、一応式は追えたものの、トリックの種明かしは分からないままである。

ところが、上記の式の応用で、以下の、よりトリッキーな式が導き出される。

これの導出は、まずc= 1 + 2 + 3 + 4 + ・・・とおく。 そして、c - 4cを以下のようにつくる。4cを、ずらして一つ飛ばしで引くのがミソである。

c = 1 + 2 + 3 + 4 + ・・・
4c =    4     + 8 + ・・・

すると、c - 4c = -3c = 1 - 2 + 3 - 4 + ・・・
なので、c = 1 + 2 + 3 + 4 + ・・・は、1 - 2 + 3 - 4 + ・・・の結果を-3で割ったものになるという次第である。

これを単なる数学的トリックと嗤ってはいけない。なぜなら、例えばこれらの式は、超弦理論の空間次元にも関与してくるのである。もしかしてゼータ関数の解析接続などに関係があるのかもしれないが、私は不勉強でそのあたりの意味が全くわからない。他のより複雑な発散級数の和と併せて、勉強してみる。

数学の本

user-pic
0

数学の本というのはなんともベタなタイトルだが、それは数学セミナー2018年5月号が「数学の本」特集だからである。それが代数や解析などの専門書の特集なら、なにを今更なのだが、今回は数学の一般書、啓蒙書、雑誌の特集なのだった。それで私は、例によって本ノートに、気に入った本をメモしていくことにした。

すると、いきなりガーンとショックを受けたことがあった。それは、フェルマーの定理に関するベストセラー、サイモン・シン、青木薫訳「フェルマーの最終定理」新潮文庫をみて、私がまだこの本を読んでいないことを確認したからである。

私はそもそも数学者でなく、せいせい数楽者なのだが、いかなる意味でも数学の一般書を軽視する立場にはない。また私は、ここにもリストされていて、リーマン予想に関する一般書の名著とされるマーカス・デュ・ソートイ、富永星訳「素数の音楽」新潮文庫をかなり前に読んで感動したことがあるので、サイモン・シン「フェルマーの最終定理」を読まない理由を見つけることは困難である。

にも拘わらず「フェルマーの最終定理」を読まないできたのは、そのことについて書かれた本を何冊も読んできたので、その背景についてはもうだいたい知っているだろうという思い込みがあったようである。しかし、冷静に顧みて私は、ワイルズの私生活をどれくらい知っているだろうか。ほとんど知りはしないのである。それで、なるべく早めにサイモン・シン「フェルマーの最終定理」を読むことにした。

あとは、リストされている本でそのうち読みたいと思ったものをリストしてみる。

  • G.H.ハーディ、 C.P.スノー「ある数学者の生涯と弁明 」(シュプリンガー数学クラブ) 柳生 孝昭 (訳)
    ハーディは、ラマヌジャンにも関与していて、興味深い。
  • マーカス・デュ・ソートイ、富永星訳「シンメトリーの地図帳」新潮文庫
    19世紀のパリ、決闘の果てに夭折した天才ガロアは新しい数学の言語を生み出していた。古代ギリシャから続く対称性探求の旅に挑む数学界の探検家たちは、その言葉を駆使して "シンメトリーの素数"を網羅した「地図帳」を完成させる。最後に発見されたのは19万6883次元の空間に潜む巨大結晶「モンスター」だった――。『素数の音楽』の著者と旅する、美しくも奇妙な数学の世界。
    という内容紹介をみても私には想像がつかない。面白そうである。
  • エドワード・フレンケル、青木薫訳「数学の大統一に挑む」文藝春秋
    これは確か、ラングランズ・プログラムの話だったか。そのうち読もうと思っていて、忘れてしまっていた。
  • 森口 繁一、応用数学夜話 (ちくま学芸文庫)
    私は、コンピュータの計算が好きで、なので森口 繁一のファンである。なので、これも私には必読なのだった。

最大の素数

user-pic
0

素数が無数にあるということは太古の昔から既知であるので、最大の素数なんて、ない。それで、ここでいう最大の素数というのは、現時点で素数であることが知られている最大の数という意味である。

それは、GIMPSによれば、277232917-1だということである。私は、Lucas-Lehmer判定法を自分のパソコン上で実装してみるというささやかな努力を通じて、それがいかに大変かを認識し、自分でも新たな大きい素数を探究してみようという努力を諦めた。

そこでせめて、277232917-1を計算してみようという、相当に妥協した考えに至った。しかし、ともう一つの自分がつぶやく。そんなの、単なる力業で、なんの工夫もいらないのではないか、と。

で、とりあえず去年作成した多倍長ライブラリを参照してまず考えたのは、int MUL(int n, int *a_num, int a)をナイーブに使用して、2の77232917乗を計算するものである。これの前準備として用意する多倍長計算のための整数配列の長さを決める必要がある。2の77232917乗は桁数としては、23249425桁であり、Radix = 100000とすると、必要な配列数は、これを5で割って、4649885なので、余裕をみて4650000個の配列とする。そして、2を1000回掛けるというナイーブな計算をさせてみると、処理時間は約48秒だった。すると、2を77232917回掛けると、見積の処理時間は、1029時間、すなわちなんと約43日もかかってしまう。

これはいくらなんでも、と思いしばし思考して、そうだ、2でなく2の14乗すなわち、16384ずつ掛ければよいのだと気付いた。これでもなんとか100000進法の枠組みに収まっている。しかも、2でなく2の14乗を1000回掛けても、処理時間はほとんど同じで、これにより処理時間を43日の1/14すなわち、3日程度に短縮できる。なお、欲を出して1000000進法にし、2の15乗を掛けようとすると、4バイト整数の32ビットの枠組みを壊してしまうので、無理である。

一方、77232917 = 5516636×14 + 13なので、2の14乗を5516636回掛け、次に2の13乗を掛け、最後に1を引くことで求めるメルセンヌ素数が得られる。それをプログラミングしたのが以下のとおりである。

といっても、まだ実行していないので、計算結果は得られていない。

#include <stdio.h>
#include <time.h>

int Radix = 100000;
int n = 4650000;
int *a_num;
int CLEAR(int n, int *a_num, int a);
int DISP(FILE *fo, int n, int *a_num);
int MUL(int n, int *a_num, int a);
int SUB(int n, int *a_num, int a);

int main()       {
int i;
clock_t t_begin, t_end;
FILE *fo;

if ((fo = fopen("Max_Prime.txt","wt")) == NULL )
 return -1;

t_begin = clock();  // 計測開始

a_num = new int[n];

// 277232917 - 1の計算
// 77232917 = 5516636×14 + 13

CLEAR(n,a_num,1);
for ( i = 1; i <= 5516636; i++ )
 MUL(n, a_num, 16384); // 2の14乗を5516636回掛ける
 
MUL(n, a_num, 8192); // 2の13乗を掛ける
 
SUB(n,a_num,1);   // 1を引く
 
DISP(fo, n, a_num);

t_end = clock();  // 計測終了

printf("Elapsed time : %7.3f秒\n",(double) (t_end - t_begin)/1000);

delete[] a_num;
return 1;
}

int CLEAR(int n, int *a_num, int a)   {
int i;
for ( i = 0; i < n; i++ )  {
    a_num[i] = 0;	
}
a_num[0] = a;
return 1;
}

int DISP(FILE *fo, int n, int *a_num)   {

int i;
int disp = 0;
int first = 1;

for ( i = n - 1; i >= 0; i-- )     {
 if ( a_num[i] != 0 )  disp = 1;
 if ( disp == 1 )     {
	if ( first == 1 )  {
		fprintf(fo,"%d",a_num[i]);
		first = 0;
    } else             {
		fprintf(fo,"%05d",a_num[i]);
	}
  }
}

if ( disp == 0 ) fprintf(fo,"0");

fprintf(fo,"\n");	

fclose(fo);

return 1;
}

int MUL(int n, int *a_num, int a)     {
int i,j;
	
  for ( i = 0; i < n; i++ )
   a_num[i] = (a * a_num[i]);  // a倍する
  // 桁揃え
  for ( j = 0; j < n - 1; j++ )  {
	a_num[j+1] += (a_num[j] / Radix); // 桁上がり
	a_num[j] = (a_num[j] % Radix);  // 10万進法	
  }	

return 1;

}

int SUB(int n, int *a_num, int a)      {
int i,j;

a_num[0] -= a;

 for (i = 0; i < n-1; i++ )  {
   if ( a_num[i] < 0 )        {
	a_num[i+1]--;
	a_num[i] += Radix;
   }
 }

  // 桁揃え
  for ( j = 0; j < n - 1; j++ )  {
	a_num[j+1] += (a_num[j] / Radix); // 桁上がり
	a_num[j] = (a_num[j] % Radix);  // 10万進法	
  }

return 1;
}

素数間距離

user-pic
0

エラトステネスの篩などのアルゴリズムで生成された素数を、小さい順からみていると、1つの素数と次の素数までの距離が気になってきた。特に、その距離が2である素数対は、双子素数と呼ばれて、古来より有名である。私も、駄文ながら、双子素数などの記事を書いてみた。

さて、双子素数を一般化して、素数間距離がnの素数対をかりに、n-SP(n-spaced prime number pair)と呼ぶことにする。nが1の素数対は、(2,3)しかないので、ここでは3以上の素数だけを考慮することにする。すると当然ながら、nは偶数となる。

そこで、素数定理で数えた、1000000000までの50847534個の素数について、素数間距離の分布を計算してみた。なお、素数定理では1000000000までの素数の個数が50847533個と示されているが、これは私が数え上げのために書いたPerlプログラムに単純ミスがあったためで、すなわち最後に$n++をしておくのを忘れたためであった。

その分布のグラフをExcelで描いてみると、以下のとおりである。

これだと全体の分布の傾向は分かるが細かい数字がみえないので、念のため具体的な数字を示す表も示す。

DistanceFrequencyDistanceFrequencyDistanceFrequencyDistanceFrequencyDistanceFrequencyDistanceFrequency
23424506434246796608979182695109103484767124468957
142464565161846097183351032201824043221569679242367474
2611195852812192433021779913268389634719531361171524
385487464064878042954456443896344633472048577247
503280665224580454410754562114625818194860371839
62115558641189516621678768883967012556472126663
7462526765511378105313805352282379828478077
862779388288789058057921928294176699631078
9816175100169001022239310410310106871910815459
1109065112713911410892116471011845021209621
122297512431361265863128204313028131323510
13414871361297138258914015161429531441555
1466681487241501486152506154583156798
158306160360162534164247166182168443
170198172158174255176128178106180200
1828018489186101188331906319274
194321964119873200282022320446
206132081821048212122141121615
21862201222211224622662283
2303232123413236523812405
242424422464248425042521
254025602580260126202640
266026802700272027402761
278028002821------

この表から見て取れるのは、意外にも、素数間距離が6の倍数のときの頻度が特異的に高いことである。もし双子素数が無限にあるなら、少なくとも素数間距離が6の素数対(6-SP)も無限にあるかもしれない。

また、当然に一般にnが大きくなるほど頻度は下がるのであるが、nが6の倍数のとき頻度が急に増える。どうも6という数はやはり素数間距離で特殊である。

そもそも、6-SPは実質的に2-SPよりも多いのかもしれない。もしどちらも加算無限なら大小なんて論じても仕方ないようであるが、任意の有限の正の数Nを選んだとき、常にN以下で6-SPが2-SPより大きいなら、6-SPは実質的に2-SPよりも多いと呼んでもいいのだと思われる。

あと、素数間距離の概念は数論で意外に本質的らしくて、こちらのJames Maynardの論文によれば、以下が成立するらしい。すなわち、n→∞とするとき、

  • pn+1 - pn ≦ 12
  • pn+2 - pn ≦ 600

これは、私が誤解しているのでなれば、素数間距離が12以下になる素数対は無限に存在するという意味である。これは素数間距離が2になる素数対は無限に存在することまでは主張するものではないので、双子素数の無限性は証明されていないが、かなり肉薄しているといえる。

最後に蛇足ながら、上記1000000000までの素数間距離分布を計算したプログラムを示しておく。

open FI, "<sieve.txt" or die "sieve.txt not found.\n";
open FO, ">SP.txt";

$p = 3;  # セットする素数の初期値
$max = 0;

while (<FI>)      {
 chomp;
 next if ( $_ < 5 ); # 5になるまでスキップ
 
 $s = $_ - $p;  # 素数間ギャップ
 $s /= 2;       # 配列に入れるので2で割る
 $k = sprintf "%04d",$s;
 $SP{$k} = $k;
 $SP{$k}->{N}++;
 $p = $_;
 if ( $max < $k )  {
   $max = $k;
 }
}

for ( $i = 1; $i <= $max; $i++ )   {
  $k = sprintf "%04d",$i;
  printf FO "%d\t%d\n",$i*2,$SP{$k}->{N};
}

close FI;
close FO;

メビウス関数の続きであるが、下記の級数を実際に1000000項くらいまで計算することにした。



下記のリストで、sは一番目の級数の値で、tは二番目の級数の値。
100000:s = -0.0004872276170375616322は、100000項までの和が-0.0004872276170375616322になることを示す。
このリストをみると、sは確かにほぼ0で、tはほぼ-1なのだが、値は微妙に脈動しており、本当に0と-1にそれぞれ収束しそうな気がしない。それにしてもメビウス関数の計算自体、素因数分解が関与するせいか、意外に重いのだった。

100000:s = -0.0004872276170375616322
100000:t = -1.0055846798890000215000
200000:s = 0.0000250654400056779753
200000:t = -0.9997135290446533817000
300000:s = 0.0007458916171313174480
300000:t = -0.9905970977815699019000
400000:s = 0.0000131148976087111140
400000:t = -0.9998111928133626147000
500000:s = 0.0000071327527987717798
500000:t = -0.9999215377875856348000
600000:s = -0.0003862460112389568424
600000:t = -1.0051335596563819585000
700000:s = 0.0003270415264529654515
700000:t = -0.9955985480298048529000
800000:s = -0.0000076252662109503673
800000:t = -1.0001187859693974058000
900000:s = -0.0002518487853556570176
900000:t = -1.0034498162961291712000
999900:s = 0.0002146054804472317858
999900:t = -0.9970212678878801471000

なお、念のため計算に使用したプログラムのソースコードを示す。

#include <stdio.h>
#include <math.h>

int IsPrime(int n);
int Moebius(int n);
int *Prime;

int main(void)       {
 int i,m,n, LIM, P_max;
 double s, t;


// 素数配列の用意
 Prime = new int[78500]; // 1000000までの素数は、78498個なので、これで間に合う。
 
 Prime[0] = 1;
 i = 1;
 
 for ( n = 2; n < 1000000; n++ )  {
	if ( IsPrime(n) == 1 )
	  Prime[i++] = n;
 }

 P_max = i;
 
 s = 0.0;
 t = 0.0;
 LIM = 100000;
 for ( i = 1; i < 999900; i++ )    {  // 999900まで足す
	m = Moebius(i);
	s += (double) m/i;
	t += (double) m*log(i)/i;
	if ( i >= LIM )    {
      printf("%d:s = %5.22f\n",i,s);
      printf("%d:t = %5.22f\n",i,t);
      LIM += 100000;
	}
 }

 printf("%d:s = %5.22f\n",i,s);
 printf("%d:t = %5.22f\n",i,t);
 
 delete[] Prime;
  return 1;

}

int IsPrime(int n)  {
 int i = 2;
 
 if ( n == 2 ) return 1;
 
 while ( i*i <= n )  {
  if ( (n%i) == 0 )
   return 0;
  i++;
 }
 return 1;
}

int Moebius(int n)     {
 int k,m,s;

 if ( n == 1 ) return 1;
 m = n;
 k = 0;
 s = 0;
 
 while ( Prime[k] <= n )  {
  
  if ( k > 0 && (m%Prime[k]) == 0 )  {
	m /= Prime[k];
	if ((m%Prime[k]) == 0 ) return 0;
	s++;
  }
   k++;	
 }
 if ((s%2) == 0 ) return 1;
 return -1;
}

ほぼ9年前に、オイラーの無限和という記事を書いたが、最近エラトステネスの篩のプログラムで、10億までの50847535個の素数について、オイラーの無限和、すなわち逆数の和を計算してみた。このため、エラトステネスの篩のプログラムの出力を、sieve.txtとして予め保存した。

その計算結果は、以下のとおりである。

         11: 1.26709956709956680
        101: 1.81271819114788050
       1009: 2.19907120745259020
      10007: 2.48315987728253030
     100003: 2.70528217874728980
    1000003: 2.88732909956469360
   10000019: 3.04144948127937240
  100000007: 3.17497523992056240
  999999937: 3.29275571880031490

このリストは例えば、素数10000019までの逆数和が3.04144948127937240である、と読む。最後の桁が0になる訳がないので、丸められているのだろう。オイラーの無限和に書いたように、「ある大きい素数である、1801241230056600467の逆数まで足したとしても、素数の逆数和は4を超えない。」のであるから、4を超えるために気の遠くなるどころではない長い道のりが必要である。3を超えるためには、ほんの一千万までの素数の逆数和でよかったのであるが。

なお、参考までに、単なる自然数の逆数和の計算も付しておく。こちらが発散することの証明を理解するのは容易だが、それでも、オイラーの無限和よりもそれほど韋駄天とは思えない。

         10: 2.92896825396825380
        100: 5.18737751763962060
       1000: 7.48547086055034330
      10000: 9.78760603604434820
     100000:12.09014612986333500
    1000000:14.39272672286498900
   10000000:16.69531136585727200
  100000000:18.99789641385255500
 1000000000:21.30048150234855000

なお、小山信也「素数とゼータ関数」共立講座 数学の輝き6の9ページによれば、素数の逆数の和が無限大になるというオイラーの定理は、数学史上で最も偉大な定理の一つとされるという。ううん、そうなのか。私のようなアマチュアにはその価値は到底分からない。しかし、その偉大さの一つの例証として、かのエルデシュが別証明を与えている。

あと、蛇足ながら上の計算に使ったPerlのプログラムを示しておく。

open FI, "<sieve.txt" or "sieve.txt not found.";

$t = 0;
$lim = 10;

while (<FI>)    {
 chomp;
 if ( $_ >= 2 )   {
  $t += 1/$_;
 }
 if ( $_ > $lim )       {
  printf "%11d:%20.17f\n",$_,$t;
  $lim *= 10;
 }
 $e = $_;
}

  printf "%11d:%20.17f\n",$e,$t;

$t = 0;
$lim = 10;

for ( $i = 1; $i <= 1000000000; $i++ ) {
   $t += 1/$i;
 if ( $i == $lim )       {
  printf "%11d:%20.17f\n",$i,$t;
  $lim *= 10;
 }
}

素数の積

user-pic
0

黒川信重「ラマヌジャン ζの衝撃」現代数学社の関連で、小山信也「素数とゼータ関数」共立出版を読んでいて、最初の方で、素数の個数は無限にあるという命題が述べられていて、ユークリッドによる証明、すなわち、素数が有限個しかないと仮定して、それをp1、p2、・・・pnとすると、p1p2・・・pn+1という数を作ってみると、これはp1、p2、・・・pnのどの数でも割り切れないので、新たな素数となる。すると素数がp1、p2、・・・pnしかないという仮定が矛盾するので、素数が有限個しかないという仮定が正しくない、ということになるのであるが、小山氏はこれは背理法ではないとする。なぜならこれは新しい素数の構成方法を示したに過ぎないからだという。

その記述を見ていて私はふと、当たり前のことを改めて考えた。すなわち、素数は結局無数にあるのだから、p1p2・・・pn+1という数を作ってみても、ここにリストしていない素数がそれを割り切る可能性は小さくないのでは、ということである。

それで私は次のようなことを考えた。それは二つの素数p,qがあったとして、積pqが偶数なら1を足し、積pqが奇数なら2を足す。このとき、p = qであってもよいとする。そうして得られた数が素数になるなら、素数p,qのペアを積素数(product prime)と呼ぶことにする。すると、積素数は無数にあるか、という問いかけである。

そこで、以下のような計算をしてみた。

素数の上限素数の組み合わせの個数積素数の個数割合
101590.6
100351810.23
10001436517730.123
10000757065629000.083

この表から分かることは、積素数は、素数が大きくなるにつれて次第に割合が下がってくることである。もしかして、ある上限の値を超えるとパタリと積素数はなくなるのだろうか。そうとも思えない。

ちなみに、素数n以下の積素数の個数を(log n)2/nと見積もってみた。この式は、n = 1000まではまあまあ合うのだが、10000になると、実際よりはるかに小さい数となってしまう。

あと、積素数且つ双子素数というのはほとんどないのだよね、と思って調べてみると、10000以下で以下の24個しかなかった。62900個の積素数中なので、相当に少ない。 (5,3) (7,5) (151,149) (181,179) (241,239) (271,269) (421,419)
(571,569) (1291,1289) (1321,1319) (2311,2309) (2551,2549) (2731,2729)
(3361,3359) (3391,3389) (4261,4259) (4651,4649) (5851,5849) (5881,5879)
(6361,6359) (6781,6779) (9001,8999) (9241,9239) (9631,9629)

ということで、積素数はもしかして無限にあるかもしれないが、積素数且つ双子素数は有限なのだろうという曖昧なことしか言えないのだった。

以下、恐ろしく成分未調整だが、計算に使ったプログラムのソースコードを掲げる。

#include <stdio.h>

int IsPrime(int n);
int *Prime;

int main(void)       {
 int i,j,n,N, P_max1, P_max2;
 int s;
 int t, T;


// 素数配列の用意
 Prime = new int[9600];
 
 i = 0;
 
 for ( n = 2; n < 100000; n++ )  {
	if ( IsPrime(n) == 1 )
	  Prime[i++] = n;
 }

 P_max1 = i;
 
  N = 10000;
  
  i = 0;
  
  while ( Prime[i] <= N )  {
	if ( i >= P_max1 ) return -1;
	i++;
  }
  
  P_max2 = i;
  
  T = t = 0;
  
  for ( i = 0; i <= P_max2; i++ )     {
	for ( j = 0; j <= P_max2; j++ )
	 if ( j <= i )        {
		t++;
		s = Prime[i]*Prime[j];
		if ((s%2) == 0 ) 
		  s++;
		else
		  s += 2;
        if ( IsPrime(s) == 1 )   {
			T++;
			if ( Prime[i] - Prime[j] == 2 )
		    printf("(%d,%d) ",Prime[i],Prime[j]);
		}
	}

  }
  
  printf("\n-----------\n%d/%d for %d\n",T,t,N);
  return 1;

}

int IsPrime(int n)  {
 int i = 2;
 
 if ( n == 2 ) return 1;
 
 while ( i*i <= n )  {
  if ( (n%i) == 0 )
   return 0;
  i++;
 }
 return 1;
}

月別 アーカイブ

このアーカイブについて

このページには、過去に書かれたブログ記事のうち数学カテゴリに属しているものが含まれています。

前のカテゴリは教育です。

次のカテゴリは日常です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。