コンピュータの最近のブログ記事

読売新聞をみていて、ええっ、と目を瞠る記事があった。それは、正義のハッカー女子と題する記事である。

正義のハッカー女子とは誰の事かというと、武蔵野市にあるNTTセキュアプラットフォーム研究所の中島明日香さん。コンピュータ・システムの脆弱性を見出す、ホワイト・ハッカーという役割の仕事をしているという。

中島さんは高校時代はIT音痴だったという。でも、無類の読書好きで、ネットで読んだ、女性プログラマー七瀬晶さんのProject Sevenという本に影響を受けたという。

父親の仕事の関係でニューヨークにいた時期があり、帰国して慶応大の環境情報学部で真剣にコンピュータを学んだ。まだ二十代後半なのに、講談社ブルーバックスで、「サイバー攻撃」という著書もあるという。これは知らなかったので、早速注文した。

NTTは、日本でコンピュータセキュリティの研究のメッカのような場所で、それは一つには、暗号理論の大家、岡本龍明博士がいらっしゃるからだろう。NTTセキュアプラットフォーム研究所も、いわゆるセキュリティのグループと、暗号のグループの二つに分かれるはずである。

日本にはもちろんNTT以外にもコンピュータ・セキュリティを研究している会社は沢山あるが、中島明日香さんの著書以外にも主だったコンピュータ・セキュリティの専門書の著者はたいていNTT所属の人が著者だったりする。例えば、私の愛読書の一つの「コンピュータネットワークセキュリティ」八木・秋山・村山著、コロナ社もそうである。

私は実は仕事の関係でときどき、NTTの武蔵野の研究所を訪れるので、たまたま中島さんにお会いできれば幸いなのだが、そもそもホワイト・ハッカーという人々と気軽に会えると思ったら大間違いで、マーク・グリーニーの暗殺者の飛躍にでてくる范講も天才ホワイト・ハッカーで、そのスキルゆえに西欧の国々から熱望され、一方でそれでは困るので中国当局が暗殺しようと狙い、ジェントリーが必死で護衛したのだった。

最大の素数の続きであるが、もしかして一気に2nで、nが数千万のときの計算を高速で実行する方法として有効なのではとふと、思いついたことがある。それは次のような方法である。

  • STEP1: x = 2と代入し、y ← x*x
  • STEP2: x = yと代入し、y ← x*x
  • STEP3: yが所定の値に達したらyを出力して終わる。そうでなければ、STEP2に戻る。
という極めて素朴なものであるが、2 → 22 → 24 → 28 → 216 → 232 → ・・・と極めて高速に数が増加して、ループがn回まわると、結果は2の2n乗になっている。

この方式で277232917を計算するには、
77232917 = 226 + 222 + 219 + 218 + 217 + 213 + 211 + 210 + 29 + 28 + 27 + 35と分解して計算する。

この分解は、例えば次のようなPerlプログラムを使用する。

IND(72232917);

sub IND     {
 my $n = shift;
 my $m = 1;
 my $i = 0;

 if ( $n < 100 )     {
   printf "reminder = %d\n",$n;
   return;
 }
 
 while ( $n >= $m )     {
   $m *= 2;
   $i++;
 }

 printf "%d\n",$i - 1;
 IND($n - $m/2);
 
 return;

}

これは速いのではと思ったら、大きなボトルネックがあった。それは、多倍長×多倍長という極めてコストの大きい計算を要することで、結果的に却って時間がかかるのだった。

具体的には、以下のとおりである。
  1. 2の215乗までの計算 : 115.875秒
  2. 2の217乗までの計算 : 461.625秒
  3. 2の218乗までの計算 : 924.688秒

このように、2の2n+1乗は、2の2n乗の計算のほぼ2倍の時間がかかると見積もれる。すると、2の227乗までの計算は、924.688×29秒 ≒ 131.5時間 ≒ 5.48日で、結局ストレートに214乗をどんどん掛けていった方が、3日ちょっとで終わる。多倍長×多倍長を使わないので、却って速い。残念。ともかく、ソースコードを以下に示す。


#include <time.h>
#include "MLA4.h"

int n = 5000000;

int main(void)           {
 int i;
 int *a_num, *b_num;
 FILE *fo;
 clock_t t_begin, t_end;
 
 if ((fo = fopen("MMUL.txt","wt")) == NULL ) return -1;
 
 t_begin = clock();  // 計測開始

 
 a_num = new int[n];
 b_num = new int[n];
 
 CLEAR(n,a_num,2);
 CLEAR(n,b_num,2);
 
 for ( i = 0; i < 18; i++ )   {  // 最大の素数の場合、この18が、26程度になる
	MMUL(n, a_num, a_num, b_num);
	COPY(n,a_num, b_num);
 }
 
 DISPF(fo,n,b_num);

 delete[] a_num;
 delete[] b_num;

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

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

 fclose(fo);
 return 1;
}

最近のC++のコンパイラが64ビット整数を実装しているとは知っていたのだが、どうも使いこなせないできた。

すなわち、知る限りで、64ビット整数はlong long intで宣言するということになっていて、例えば
long long int *a_num; とポインタは宣言できるのだが、
a_num = new long long int[1000];のように実際に配列をアロケートする式が文法エラーになってしまう。

また、long long int a;と宣言して、a = 12399871717171;と大きい数を代入しようとするとコンパイル時に警告が出て、aには実際は想定外の値が代入されてしまう。

さらに、printf("%d\n",a);も想定外の値を表示する。

そこでいろいろウェブで調べていて、long long intと等価な宣言として、__int64があることを知った。これだと、__int64 *a_num;も、a_num = new __int64[1000];も問題ない。

a = 12399871717171;は、a = 12399871717171LL;と、64ビット整数であることを示すLLを末尾に付ければいいことが分かった。

printf("%d\n",a);は、printf("%I64d\n",a);と書けばいいことが分かった。これはprintf("%||d\n",a);と書いても等価なようなのだが、少なくともEmbacadero C++ Buidler XE4だと、正しく動いたり、そうでなかったりすることがあるので、専らprintf("%I64d\n",a);を使うことにした。

以下に、正しく動作したサンプルプログラムのソースコードを掲げる。


#include <stdio.h>

int main()     {

int i;
__int64 a,b,c,d;
__int64 *c_num;

c_num = new __int64[10];

a = 10010277776676688LL;
b = 5555552291LL;
c = a%b;

printf("%I64d\n",a);
printf("%I64d\n",c);

for ( i = 0; i < 10; i++ ) {
 c_num[i] = a++;
}
 
for ( i = 0; i < 10; i++ )  {
 printf("%I64d\n",c_num[i]);
}

delete[] c_num;
return 1;
}

最大の素数

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;
}

BIOS update

user-pic
0

連休中の5月5日だったと思うが、何かWindows Updateが来ていないかと確認したら、Windows 10バージョン1803の機能更新プログラムというのが来ていた。それで、軽い気持ちで導入を開始したら、なんとこれが重いこと、3時間もかかってしまったが、どうもこれは年に2回のWindows 10のメジャー・アップデートで、もしかしてSpring Updateというものかもしれない。というのは、これがめでたく終わると、いろいろい追加された新機能を紹介する画面になったからである。しかし私はずぼらなので、今までのメジャー・アップデートで、どんな機能が追加になったか確認したこともなく、今まで通り使い続けているだけである。

それと、5月9日は第2水曜日で、定例のWindows Updateの日だったが、Windows Updateフリークの私には珍しくうっかりWindows Updateを忘れていて、なんかパソコンが不安定だよね、大丈夫かな、と能天気に?不安がっていた。すると、勝手にWindows Updateが進んでいたようで、あるタイミングでWindows Updateが終わっているらしいので、パソコンを立ち上げなおしたら、事なきをえた。

いやあ、次々にいろいろなUpdateが来るもんだと怪訝になっていたところ、5月10日の夜、パソコンを立ち上げたまましばらく席を離れて戻ったところ、見慣れないスクリーンが表示されていた。なんだろう、なにかの不具合かと思ってよく見ると、なんとLenovoのBIOS Updateだった。へえ、いまどきはBIOS Updateもこうやってくるのかと思いつつOKをクリックしたところ、しばらく何かのプログラムをダウンロードしたあと、急にパソコンがリブートされ、立ち上がるとBIOS Updateの導入画面になった。

その導入画面の赤い帯の部分を眺めていてふと、もしこれが手の込んだウイルスなら、パソコンのファームウェアは完全に破壊されてお陀仏だよね、と変な不安がよぎった。だとしても、毒食らわばBIOSまでの諺?に従い、開き直って導入したら、無事立ち上がって安心した。しかし、立ち上げ時のLenovoの表示図面が更新され、赤みがかっていて無気味なのが気になる。BIOSなんて縁の下の力持ちなんだから、もっと地味な表示にすべきである。

ともかく、こうして私のパソコンライフは命尽きるまで続くことだろう。しかし、いくら私が長生きしても、AIはこういうシーンではあまり活躍しないと思われる。なぜなら、システムの低レベル(ハードウェア寄りという意味であって、簡単だとか、誰もが理解できるという意味ではない)のところは泥臭くて、人間が思考錯誤しながら作業するよりないと思われるからである。

Windows 10の運用の続きであるが、Windows 10のメモリ・リソース管理について、ちょっと思いついたことがあるので、メモしておく。

まず、Windows 10に限らず、オペレーティングシステム全般にそうだと思うが、アプリケーションを起動したり閉じたりを繰り返すと、メモリ・リークなどのためにより、次第に使えるメモリ・リソースが減ってきて、使えるメモリ・リソースがある閾値以下になると、アプリケーションが起動できなくなる、プリンタ・ドライバなどのデバイス・ドライバが見えなくなる、ネットワークとの通信ができなくなる、などのよく分からない障害が発生する。人間でもコンピュータでも、リソース不足は、一見原因不明の障害をもたらすのである。

一方、何がメモリ・リソースをよく消耗するかというと、私の見るところ、AbemaTVなどのネットワーク・ストリーム系のウェブ・システムあたりそうだったりする。他には、シンプルながらエラトステネスの篩のようなプログラムもそうだったりする。すなわち、それはa_num = new int[AR];によってメモリをアロケートするのだが、実はこの場合、125MBものメモリをアロケートすることになる。もちろんプログラムを終了するときはdelete[] a_num; で礼儀正しくメモリを解放するのだが、なぜか開放されきらなかったりする。

それで当然に私にも対策はある。それは、Windows 10付属のリソース・モニタで利用可能なメモリの残量を確認する。ところが、どうみてもリソース不足のはずなのに、リソース・モニタに反映されていなかったりすることがある。それで私は、リソース・モニタは参考にしつつ、それはあくまでも参考として、勘によりどうもリソースがまずそうだと想定したら、システムを再起動している。

ところが、システムを再起動して不思議なことが起きたりする。それは、リソース・モニタでみると、システムが使用しているメモリがかなり増えていたりする。あるときは、システムが2747MBも使用していた。スタンバイが756MB、空きが458MBである。再起動直後で、何も明示的にアプリケーションを起動していないのに、である。明示的に、というのはバックグラウンドでいろいろアプリケーションが走っているからである。

ともかく、その状態でEmbacadero C++ Builder XE4でコンパイルさせてみると、案の定、LME288エラーが起きた。

そこでさらに再起動させると、今度はシステムが使用するメモリが2202MBに減って、スタンバイが1416MB、空きが326MBとなった。この状態では、LME288エラーが起きなくて、万歳である。

それにしても、再起動した直後は普通、常に同じ状態のはずなのに、なぜメモリの状態が毎回異なるのだろうか。もしかして、リソース消耗が激しい状態では仮想記憶システムがディスクにスワップ領域を確保して、それは一度再起動しただけでは開放されないのだろうか、とかいろいろ考えてみるのだが、インサイドWindows 10のような事情に全く疎く、なにひとつ分からないのであった。

エラトステネスの篩で、以下のように、プログラム実行時間について言及した。

このような方針で作成したプログラムのソースコードは以下のとおりである。すると、10億までの素数の計算が、わずか38秒程度で完了した。

しかし、考えてみると、エラトステネスの篩といっても、バイト列のビットにエラトステネスの篩のアルゴリズムに従い1を立てる処理と、その位置に従い、実際の素数をファイルとして書き出す処理とに分かれる。その各々にどれだけ時間がかかるか測定しようとしたら最早、コードを書き込むしかなく、その方法をウェブで調べたので、実装してみた。また、折角なので、多少コードも改良してみた。下記の赤字で示したのが、処理時間計測用のコードである。

要するに、時間計測用の変数をclock_tとして宣言しておき、clock()で、1/1000秒単位でパソコンの内部時間を取得して変数に代入し、後でそれらの変数の差を取ることで、時間差が分かる。好適には、差を1000で割ることで、時間単位を秒に変換する。

すると、ある計測では、プログラム開始から篩計算の終了までが12.984秒で、プログラムの終了までが33.718秒だった。ある計測では、といちいちことわる理由は、試行ごとに時間が多少変わるからである。なぜなら、マルチタスクなので、同時に他のプログラムが走っているからである。

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

 unsigned int MASK[32];
 int IsPrime(int n);

int main(void)        {
 FILE *fo;
 int i,j,k,n, p_max;
 int a,b;
 int N, LIM, AR;
 int *a_num, *p_num;
 clock_t t_begin, t_end1, t_end2;

// MASKの用意
 t_begin = clock();
 
 if (( fo = fopen("sieve_n.txt","wt")) == NULL )
   return -1;
 
 n = 1;
 i = 0;
 do { MASK[i++] = n;
      n <<= 1;
 } while ( i < 32 );
 
N = 1000000000;
AR = N / 32 + 1;
LIM = int(sqrt((double) N)) + 1;

// LIM以下の素数の配列を用意

 n = int((double) LIM*1.02/(log(LIM) - 1.0)) + 1;
 p_num = new int[n];
 j = 0;
 for ( i = 2; i <= LIM; i++)   {
	if (IsPrime(i) == 1 )
	  p_num[j++] = i;
 }
 p_max = j;
// 素数配列完了

 a_num = new int[AR];
 
// 初期化

 for ( i = 0; i < AR; i++ )
  a_num[i] = 0;


//篩ビットのセット
 for ( i = 0; i < p_max; i++ )    {
   j = 2;
     while (( k = p_num[i]*j) <= N )    {
     	a = k/32;
    	b = k%32;
	    a_num[a] = a_num[a] | MASK[b];
	    j++;
    }
 }
 
 t_end1 = clock();
 
k = 0;
 for ( i = 0; i < AR; i++ )                {
   if ( (a_num[i] & 1) == 1 ) k++;
   else { fprintf(fo,"%d\n",k++); }
   for ( j = 1; j < 32; j++ )  {
	  a_num[i] >>= 1;
	if ((a_num[i] & 1) == 1 )      {
	    k++;
	} else { fprintf(fo,"%d\n",k++);           }
	 if ( k > N-1 ) goto EXT;
   }
 }
EXT:

delete[] a_num;
delete[] p_num; 

t_end2 = clock();

printf("篩計算まで、%6.3f秒\n",(double) (t_end1 - t_begin)/1000);
printf("素数出力完了まで、%6.3f秒\n",(double) (t_end2 - t_begin)/1000);

return 1;
}

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

}

HEROZ

user-pic
0

数日前、新聞をみていたら、ポナンザの開発者の山本さんが参加している、HEROZという会社が、東証マザーズに上場したのだという。以下、改めて日経の記事を引用させて戴く。

人工知能(AI)開発のHEROZ(ヒーローズ)が20日、東証マザーズに上場した。投資家の期待は大きく、終日買い気配で初値はつかなかった。将棋ゲームで名高い同社は、NEC出身の2人が2009年に創業。メンバーには将棋有段者のエンジニアも多い。将棋AIの開発で培った技術や知見を惜しみなく産業応用に注ぎ込む。
「AIの可能性は無限大だ」。ヒーローズの林隆弘最高経営責任者(CEO)は、20日の上場セレモニー後の会見でこう語った。ヒーローズは、将棋ソフト「将棋ウォーズ」などで知られる。これまで実績を持つB2C(消費者向け)サービスに加え、B2B(産業向け)の事業の柱に据えて、AIを建設業や金融業などに活用する。

私は、ここまで読んで、目を瞠った。というのは、CEOの林隆弘さんといえば、私の知る限り、ガチのアマ強豪なのである。そんな将棋の本格強豪なのに、CEOが務まるなんて、なんというパワーなんだろう。

それにしても、将棋AIの開発で培った技術や知見を活用する、というのがいいよね。一将棋ファンとしてうれしくなる。

「僕らは血へどを吐くような実装をしているんですよね」。リードエンジニアの山本一成氏は笑う。東大将棋部出身の山本氏が開発した将棋AI「Ponanza(ポナンザ)」は、17年に名人に勝利した。

確かに、強い将棋AIの開発には、とてつもないハードウェア・リソースとアルゴリズムの工夫が必要であることをポナンザで言及した山本氏の著書を読んで知った。そのスキルとノウハウは、他の分野での問題解決にも役立ちそうである。

山本氏は、DeepZenGoの開発にも関与されているそうであるが、その開発が終了したことで、その開発グループもまた、何かソリューションビジネスを始めるのだろうか。それとも、アカデミックな世界に還るのだろうか。

ともかく、囲碁や将棋は、AIの威力を世間に見せつける最高のshow caseであったことは間違いない。

ただ、国立情報学研究所教授で東大合格を目指すAI「東ロボくん」の開発で有名な新井紀子氏によれば、自然言語が関与するとAIはかなり限界があって、シンギュラリティなんてありえないという。私は、解ける問題はすさまじく解き、そうでない問題は周囲をくるくる回って核心になかなか至らないのがAIだと思っているが、このあたり、ブレークスルーがあるのかどうか、行く末を眺めたい。

カラーリング・スキームで適当な十六進数のカラーリング・スキームを適用してみて、だいたいもういいかなとも思ったが、ちょっと乱数を入れたカラーリング・スキームも作ってみたくて、それで最初はカラーリング・スキーム関数を呼び出す毎に乱数が入るようにしてみたところ、霧吹きのような図になって少なくとも私の趣味にはあわないので棄却した。

それで考え直して、少なくともプログラムが呼ばれたとき乱数に基づきカラーリング・スキームの値を決めて、それを通すようにした。それでも確かに毎回実行する度に色が変わるので当初の趣旨は活かしているとはいえる。ただ、いろいろコードを弄っているうちに奇妙奇天烈なコードになってしまったが、それはとりあえず最後に示す。

それと、ともかくこのプログラムを何度も実行してみて、比較的自分の好みにあった5個の画像を以下に示してみる。







//---------------------------------------------------------------------------

#include <vcl.h>
#include <stdlib.h>
#include <time.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
TColor ColorN(int n);
int ColorW(int n);
int RND,RND1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button2Click(TObject *Sender)
{
double x,y,tx,ty,x0,y0;
// double a = -0.541;
// double b = 0.5112;
double a = -0.15652;
double b = 1.03225;
int n,A,B;
int px, py;

srand(time(NULL));
RND = rand()%3;
RND1 = rand()%0xf35c;

 for ( A = -300; A < 300; A++ )      {
   for ( B = -300; B < 300; B++ )  {
	 x = (double) A/500.0; y = (double) B/500.0;
	 n = 0;
     while ( x*x + y*y < 4.0 )   {
      tx = x*x - y*y + a;
      ty = 2.0*x*y + b;
      x = tx;
      y = ty;
      n++;
      if ( n > 10000 ) break;
    }
    px = 320 + A; py = 320 + B;
    if ( n > 10000 ) Canvas->Pixels[px][py] = clBlack;
	else Canvas->Pixels[px][py] = ColorW(n);

  }
 }	
  
}

TColor ColorN(int n)   {
 int i;
 i = (n%7);
 
 switch(i)           {
  case 0:
   return clRed;
  case 1:
   return clYellow;
  case 2:
   return clGreen;
  case 3:
   return clLime;
  case 4:
   return clPurple;
  default:
   return clSkyBlue;
 }
 return clWhite;
}

int ColorW(int n)   {
 int i,j;
 j = RND;
 i = (n%9);
 
 switch(i)         {
  case 0:
   return 0x00dc143c + RND * RND1;
  case 1:
   return 0x00cd5c5c + RND * RND1;
  case 2:
   return 0x00f7d358 + RND * RND1;
  case 3:
   return 0x00fa8072 + RND * RND1;
  case 4:
   return 0x00ff6347 + RND * RND1;
  case 5:
   return 0x00d2691e + RND * RND1;
  case 6:
   return 0x00b18904 + RND * RND1;
  case 7:
   return 0x00ff69b4 + RND * RND1;
  default:
   return 0x006090ef + RND * RND1;
 }
  return 0x00ffc0cb;
}
//---------------------------------------------------------------------------

最近の私の愛読書の一つである、黒川信重「ラマヌジャン ζの衝撃」現代数学社を眺めていたら、ラマヌジャンの素数公式の業績が言及されていた。実はその素数公式には決定的な間違いがあって、ハーディはそのことを以て、ラマヌジャンが複素関数論さえ知らないと失望したそうなのだが、間違いとはいっても、リーマンの素数公式の第一項までは正しいのであり、寧ろ独学でそこまで導いただけでも凄いことだと黒川氏は仰られている。ラマヌジャンのファンの私は、そうだそうだと応援したくなる。きちんとした数学教育のバックグランドがなかっただけなのである。

そのリーマンの素数公式には、ラマヌジャンが気づかなかった、ゼータ関数の虚零点の全体に亘ってある式を足す、という式が含まれている。リーマンはその虚零点の実部がすべて1/2だと予想した。それが名にし負う、世紀の難問リーマン予想なのだった。よって、リーマン予想が解けると、素数公式の厳密解が得られることになる。

そのラマヌジャンの式にもリーマンの素数公式にも、メビウス関数μ(m)というのが出てくる。私は不勉強でこれが何か分からないので、「ラマヌジャン ζの衝撃」を注意深く読むと、メビウス関数μ(m)とは下記のように、ゼータ関数ζ(s)の逆数を級数展開したときの係数にあたるものらしい。

このとき、ゼータ関数ζ(s)のオイラー積の公式を用いると、以下とおりとなる。すなわち、自然数mを素因数分解して、
mが1または、異なる偶数個の素数の積ならμ(m)は、+1
mが、異なる奇数個の素数の積ならμ(m)は、-1
その他の場合、すなわち任意の素数の二乗で割り切れるなら、μ(m)は、0

そこでμ(m)を計算するプログラムを書いてみた。下の方にリストしてみる。 また、その結果の一部を示すと以下のとおりである。

μ(1) = 1
μ(2) = -1
μ(3) = -1
μ(4) = 0
μ(5) = -1
μ(6) = 1
μ(7) = -1
μ(8) = 0
μ(9) = 0
μ(10) = 1
μ(11) = -1
・・・
μ(9989) = 1
μ(9990) = 0
μ(9991) = 1
μ(9992) = 0
μ(9993) = 1
μ(9994) = -1
μ(9995) = 1
μ(9996) = 0
μ(9997) = 1
μ(9998) = 1
μ(9999) = 0

結果は、1と-1と0しか出てこないので単調のようであるが、なんと以下の式が成立するそうである。



  • この式は、アダマールと、ド・ラ・ヴァレ・プーサンが独立に証明した素数定理と同等な結果であるという意外に深い式なのだった。

ということで、可能な範囲でメビウス関数について考察してみたい。

#include <stdio.h>

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

int main(void)       {
 int i,n, P_max;


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

 P_max = i;
 
 for ( i = 1; i < 10000; i++ )
  printf("μ(%d) = %d\n",i,Moebius(i));
  
  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 ((m%Prime[k]) == 0 )  {
	m /= Prime[k];
	if ((m%Prime[k]) == 0 ) return 0;
	s++;
  }
   k++;	
 }
 if ((s%2) == 0 ) return 1;
 return -1;
}

月別 アーカイブ

このアーカイブについて

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

次のカテゴリはスピリチュアルです。

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