複天一流:どんな手を使ってでも問題解決を図るブログ

宮本武蔵の五輪書の教えに従い、どんな手を使ってでも問題解決を図るブログです(特に、科学、数学、工学の問題についてですが)

彗星観測やっと成功

19世紀の初めに発見された70年周期の楕円軌道の彗星、それが今回北半球の空に浮かんでいるPons-Brooks彗星である。 これまで2回挑戦したが、いずれも彗星の位置を確認することができなかった。最初は雲がかかり、2回目は彗星の位置がわからなかった。

しかし、牡羊座と三角座のだいたいの位置がわかってきたこともあり、今回は成功の可能性がだいぶ高まったと感じていた。なにより、雨上がりの青空が終日広がり観測の好機であった。とはいえ、地平線に近い彗星は日に日に夜の空から消えつつある。しかも、この時期は日没がどんどん遅くなっている。もはや6時半でも夕空になってしまっている。これは冬が終わり春が近づいているという意味だから、普通の人には嬉しい変化である。しかし、日没の地平線に浮かぶ彗星を狙うアマチュア天文家には、恨めしい状況なのである。今日失敗すれば「徳俵状態」になるだろう。

木星から牡羊座、そして三角座の先端はあっという間に見つかった。あとは少しずつずらしながら彗星を探すだけである。広角レンズでなんとなくそれらしいものが写ったが、はっきりしない。はっきりさせようと思って、露出時間を長くすると昼間のような青空の写真になってしまう。夕方の地平線間際の天体観測は、残光がとても邪魔である。

しかたないので、しばらく待ってから、今度は望遠レンズで狙ってみることにした。視野が狭くなるので天体の位置がはっきりわかっていないとなかなか天体をレンズの中に入れることができないのだが、局所的に光をたくさん集めたほうが淡い天体は補足しやすいと考え、望遠レンズ使用に踏み切った。

何度か試し撮影をした後、ついにその美しい姿がカメラに映ったのであった!やはり5等星ならばデジタル写真には綺麗に映る。

画像処理にチャレンジ(part 1):彗星観測の失敗と無謀な挑戦

Pons-Brook彗星の観測失敗

今年は大きめの彗星が2つほどやってくると噂に聞いた。実は一つ目はすでに観測可能であり、ネットを検索すると綺麗な緑青色のコマをもった彗星の画像が溢れている。

www.theguardian.com

(Arundel城の頭上に光るこの写真はおそらく合成された偽写真だと思うが)こんな感じで見えたら最高である!ということで、遅まきながら私も観測に参加することにした。ちなみに、この彗星は19世紀の初めに発見されたが、最初がフランス、その70年後が英国の天文学者による観測だったという。つまり70年ほどの周期で太陽の周りを周回する楕円軌道の彗星である。

彗星の観測といえば、吉田氏のHPである。光度グラフをみると5月1日あたりが最大光度のようだが、関東周辺の場合、現在は夕方になると地平線の上方にわずかに見える位置にある。おそらく4月に入ると夜間は地平線に沈み、昼間の空に彗星の位置が移動するだろうから、(北半球からは)観測できなくなるであろう。ということで、「今しかない」のである。

さっそく観測してみたが、現在5等星(最大光度の場合でも4.5等星)なので肉眼でみるのはきついが、観測写真なら十分捕捉可能な明るさである。しかし、現在の彗星存在地点が馴染みの薄い場所であり探しにくいという問題点が(私の場合は)ある。

ということで、広角レンズで、少し広めに夜空を切り取ってみた。 今、もっとも目立つのが木星である(写真の中央からやや下めに位置する明るい点)。その少し上にプレアデス星団(すばる)が見える。すばるの左斜め上にあるV字型の星座が牡牛座でアルデバランの赤い点が比較的よく目立つ。昨年は、アルデバランのすぐ横にZTF2023E3彗星が位置した時期があり、撮影に成功した。アルデバランの赤(朱色)と彗星の緑青のコントラストが見事であった。そのさらに右上にずーっと視線を動かすとオリオン座が右下を覗き込むように傾いて光っている。

しかし、今回のPons-Brooks彗星は「逆側」の領域にある。逆という意味は、木星とスバルをむすぶ垂直線に対して右側という意味である。右側に目を写し、上方を仰ぎ見ると比較的明るめの恒星カペラがみえる(ぎょしゃ座)。写真の中央右上にぎょしゃ座の5角形は丸ごと写っている。カペラの下にあるのがペルセウス座で、星がぐちゃぐちゃっと写っている領域である。その手にメデューサの首ならぬアルゴル変光星が握られている。カペラとアルゴルを結ぶ直線の先に2つの星が縦に並んでいるのが見える。これが牡羊座のハマル(上)とシェラタン(下)である。写真中央部の地平線近くにのっぽの針葉樹が見えるが、そのすぐ左上にある、縦に並んだ2つの恒星である。この星の右上に、鋭角三角形が下向きに位置している。三角座である。彗星は、牡羊座と三角座の中間に位置しているはずである。...が、残念ながらこの写真には写っていない。

拡大写真

彗星は動き回る天体なので、毎回毎回その位置ぎめが大変な作業である。星座によく習熟しておかないと、さっぱりの場合もある(今回がそれである)。いつも思うのが、星のパターンを自動認識してくれるソフトウェアがないかな、という願望である。それよりなにより、この写真に写っている星はいったい幾つあるのだろうか?明るい星はどのくらいで、暗い星は何個あるだろうか?そういう分析こそ、画像処理が得意ではないのか?

....と思ったのがきっかけである。デジタル画像処理について興味が湧いたので、さっそく自分でできる範囲までやってみようと思い立ったのである(無謀?)

画像処理の教科書を納屋から引っ張り出す

遥か以前、画像処理をやってみようと思い立った時期があっった。あの頃は、変光星に興味があり、定期的に観測して、その光度がどのように時間変化するのか興味があったのだ。フリーウェアがいくつか見つかったのだが、そのほとんどがWindows向けのものであり、LinuxからmacOSに軸足を移しつつあった私は、その試みを諦めたのである。吉田氏の彗星光度分析にあるような画像解析をやりたかったのだった。

そこで、自分でなんとかならないかと考え、いろいろな書籍を購入してみたのだが、そのほとんどが(やはり)windows向けの書籍であり、内容を読んでもさっぱりわからず、再び挫折してしまったのである。

この時購入したいくつかの書籍が埃をかぶって納屋に眠っているのを思い出し、もしかすると今なら多少はわかるようになっているかも、と根拠のない楽観に基づき、再再度挑戦することにしたのである。

doomとnanosaursの経験が生きているかも

この間チャレンジしたdoomnanosaursというゲームのコンパイルをやっている最中に、アプリケーションのコードがどのように書かれているのか、多少見聞きする機会があった。特にlinuxxdoomでは、C言語のプログラムのバグを手で修正したので、その構造や書き方について(少しではあるが)その慣習を知ることができた。

昔買った教科書のひとつが「画像処理プログラミング」というものである。以前はこのプログラムの意味がまったくわからず、コンパイルも通らず途方に暮れたのだが、今見ると、どうしてうまくいかなかったのかすぐにわかった。関数のコードリストはそのままコンパイルできないこともあるし、includeファイルの設定や、main()の記述などもどこかでやらねばならないのである。

doomの場合には、mainはほとんど形式的な利用になっていて、コードを見ると、すぐにD_DoomMain()を呼び出して実行的な処理はそちらでやるようになっていた(main関数はi_main.cにある)。

int main( int argc,  char **argv ) 
{ 
    myargc = argc; 
    myargv = argv; 
 
    D_DoomMain (); 

    return 0;
} 

さて、教科書に掲載されたプログラムであるが、まずはmain関数が載っているプログラムを探さないといけない。ずっとページを追いかけていくと、第1章の終わりに近いところにあった。ただ、ヘッダーファイルのリストが教科書のどこにも載っていないことがわかった。これではコンパイルが通らないわけである。

よくよく調べると、この教科書のコードはすべてソフトバンクのサーバーで公開されているそうで、まずはそちらをダウンロードしてみることにした。すると、そこにヘッダーファイルが含まれていたのであった。しかし、これだけではまだ不十分である。ダウンロードパッケージにはmakefileが含まれていなかったので、自分で書く必要があった。ファイルの依存関係を見出すのがなかなか骨折りであったが、mainを含むコードを読んで、必要な関数が含まれているコードを探し、それを個別にオブジェクトファイルにコンパイルできるかどうかチェックし、最後にリンクしてみたら、うまくいった!

できあがったコードを走らせると、小さなBMPファイルが生成されていた。教科書で解説していた市松模様の画像であった!

灰色の四角形に見えるが、拡大すると市松模様になっているBMP画像ファイル

これで、やっとこの教科書の読み方がわかった。画像処理の技術はひとつひとつが関数として分割して書かれており、見栄えが悪いヘッダーファイルは教科書では省略されてしまっている。コードの書き方だけを学んでね、という書き方である。

実際にコンパイルしてコンピュータで動かすには、ヘッダーファイルを書き足し、必要なメイン関数を書いて、関数ファイルとリンクさせる必要がある。主な作業はダウンロードパッケージに含まれたファイルを借用すれば省けるが、ファイルの依存性だけは自分でやらないとだめである。コードの中身を読み、教科書を読みながら、自分で分析してmakefileを書く。これらを突破し最後までは辿り着けば、実行ファイルが生成されて、BMP画像が作られるのである。

教科書をまねる

どっちみち自分でコードの構造を理解し、分析しないと動かないのであるから、思い切って同じ内容のプログラムを自分で書き直してみることにした。

まずは、画像データの先頭に組み込まれているデータ情報の読み込みをやってみることにした。画像ファイルをプログラムで読み込めなければ、なにも始まらない。

この教科書はC言語で画像処理を実装している。画像形式としては圧縮しないタイプのBMPが選ばれている。jpegとかpngとかをやりたいところであるが、これら圧縮された画像データの処理はかなり複雑であり、ライブラリを介してやらないと大変骨が折れる(ライブラリを使っても大変だと思うが)。画像処理プログラミングの最初の一歩としては(現在はあまり利用されなくなったとはいえ)BMPが適しているのかもしれない。BMP利用が減少した最大の理由がその大きなファイルサイズであったが、近年の技術革新により、BMPサイズの大きさがあまり気にならなくなってきたというのもある。Windows3.1やOS/2が利用されていた「大昔」においては4MBのフロッピーディスクに画像を何枚しまえるかどうかが問題だったが、今では4TBの外付けHDDや256GBのSDカードが1万円から2万円程度で買えるようになった。もはや画像はBMPのままでも問題ないのかもしれない(ネット利用でなければjpegは必要ないかも)。

BMPファイルを用意する

まずはBMPファイルを用意する必要がある。教科書のプログラムを使って生成した(上の)市松模様のファイルがその一つである。もう一つは、macOS Montereyのバージョンを示す画像をpngスクリーンショットした後、imagemagickをつかってBMPに変換したものである。

まずは、imagemagickをつかって、用意した画像のプロパティを確認しておこう。

% identify out.bmp

out.bmp BMP3 100x100 100x100+0+0 8-bit sRGB 30054B 0.010u 0:00.000

いろいろな情報が表示されたが、わかるのは「形式がBMP3、大きさが100x100ピクセル、ファイルサイズは30054byte」ということであろう。

pngからbmpに変換したものについても同じようにやってみると、

%identify monterey.bmp

monterey.bmp BMP 654x400 654x400+0+0 8-bit sRGB 789826B 0.010u 0:00.003

となった。ファイル形式は(BMP3ではなく)BMPとある。いろいろなバージョンがBMPにもあるのであろう。「大きさは654x400ピクセル、そしてファイルサイズが789826byte」であった。

その他の情報については、BMPについて詳しくなれば次第にわかるようになるであろう。

ファイルを読み込むコード

BMPはバイナリファイルであるから、まずはC言語でバイナリファイルを開いてから読み込むことになる。

K&Rを久しぶりに読み直してみると、「ファイルはポインタで定義される」とある。バッファの位置や、その現在の文字位置、加えて読み書きどちらなのかといった情報、ファイルの終端(EOF)が来たのか来てないのか、エラーのあるなしなど、ファイルに関する様々なプロパティが構造体の形でまとめられているのだという。そしてその情報の先頭が置かれているアドレスがファイルポインタである。したがって、まずはファイルポインタを宣言し、そのインスタンスを生成する。ちなみに、このファイル構造体の定義は、stdio.hに書かれているという(見たことはないが)。

ということで、C言語でファイルの読み出しを行う第一歩は次のようになる。

#include <stdio.h>

int main(){
    FILE *fp;
    return 0;
}

この状態ですでにコンパイルでき、実行できる(なにも起きないが)。

% gcc openFile.c

次に、BMPデータファイルの情報を、このファイルポインタ構造体に結びつける。このために必要なのが、stdioライブラリに収められているfopenである。ファイル名とファイルアクセスモード(読み込みか書き込みかなど)を指定して実行する。

FILE *fp;
fp = fopen("out.bmp","rb");

rは読み込み(read)、bはバイナリ(binary)形式、を意味する。BMPはバイナリファイルなのでbが必要である。

「ファイルオープン」つまりファイルとファイルポインタの結びつきが失敗すると、ポインタfpにはNULLが結びつく。それ以外の場合は成功と見做される。K&Rでは、ファイルオープンの成功失敗の検査も込みした例文が書き下されているので、それを採用してみた。

#include <stdio.h>

int main(){
    FILE *fp;
    if( (fp = fopen("out.dat","rb")) == NULL){
        printf("File not found. Failsure.\n");
        return 1;
    }else{
        printf("File found. Success.\n");
   }
    return 0;
}

最後にファイルの中身にアクセスしてみる。といっても最初2バイトのみとする。これは、BMPデータファイルの最初の2バイトには、識別のための"B"と"M"が必ず書き込まれているからである(そういう仕様になっていて「マジックナンバー(魔法数)」と呼ばれる)。

ASCII文字コード

C言語は米国で開発されたプログラミング言語であるため、扱う文字は米国系のアルファベットと関連する記号である。ドイツ語だったり、北欧系、東欧系の諸国で利用されている文字や記号は採用されてない。したがって、128通りの区別があれば、米国系の英文に限ってならコンピュータで表現できる(らしい)。ということで、文字を表す表現としては27=128の情報をもつASCIIコードが採用された。たとえば、Aは整数65、Bは整数66....によって識別する。ちなみに、「1」は49, 「2」は50、...といった具合に「数字」についても「文字」として認識され、ASCIIコードが割り振られている。

BMPの最初の2バイトに書き込まれているBMという2つの文字は、ASCIIでは66と77に相当する(16進数では0x42と0x4dである)。

現代のコンピューターアーキテクチャーでは、1byte=8bitであるから、28=256文字を処理することができる。ASCIIコードの最大値は127=27-1、つまり7ビットであるから、大幅に「余り」が出てしまう。この余の部分に各国の言語特有の記号を埋め込んだ新しいASCIIコードというものが制定されてはいるようだが、今考えるべきはオリジナルのASCIIコード7ビット分、つまり0.5バイトの情報量である。しかし、0.5バイトというのは中途半端な数字である。ファイルに文字を書き込む時、実際はどのようにやっているのか疑問に思うのである。

putc()関数

ファイルポインタfpで定義されたファイルの中へ一文字つまり1バイトのASCIIデータを書き込む関数はputc()である。

似たような関数にfputc()があるため、どちらを使ってよいか初心者は非常に困惑する。色々読んで調べてみたが、結局はどちらでもよいということになった(少なくとも初心者目線では)。教科書ではfputcが採用されているが、あえてputcを使って書き換えてみることにした。どちらでも同じ、というならば、ちゃんと動くだろう。

次のようなコードを書いて、ファイルに1文字書き込んでみた。果たして0.5byteになるのか、それとも1byteになるのか?

#include <stdio.h>

int main(){
  int c;
  FILE *fp;
  fp = fopen("sample.dat","wb");
  c = putc('B',fp);
  fclose(fp);
  return 0;
}

wbというのはバイナリ形式での(新規ファイルへの)書き込みモード、という意味である。書き込む先のファイルは「sample.dat」と名付けた。1文字を表すときは一重引用符で囲む。つまり'B'とする。二重引用符で囲む、つまり"B"とすると、C言語では文字列扱いになってしまうので、終端記号'\0'が自動的に付加されて2バイトになってしまう。putcが書き込めるのは、あくまで「一文字」だけであり、文字列は別の関数を使わねばならない。

さて出来上がったファイルを見ると、

% ls -l sample.dat

-rw-r--r--  1 user  group  1  3 26 12:50 sample.dat

% xxd sample.dat

00000000: 42                                       B

となっていて、1バイト、文字Bひとつだけが書き込まれたファイルになっている。では2文字書き込んだらどうなるか?

#include <stdio.h>

int main(){
  int c;
  FILE *fp;
  fp = fopen("sample.dat","wb");
  c = putc('B',fp);
  c = putc('M',fp);
  fclose(fp);
  return 0;
}

結果は、

% xxd sample.dat

00000000: 424d                                     BM

% ls -l sample.dat

-rw-r--r--  1 user  group  2  3 26 12:53 sample.dat

たしかに2バイトである。

この実験により明らかになったのは、ASCII文字コードは7bit仕様(128文字)であるが、現代のコンピュータアーキテクチャでは、8bit(256文字)=1byteとして扱っているということである。かなり無駄をしているように見えるが、無駄が無駄にならないくらい技術が進歩したのであろう。

ちなみに、「無駄」な情報といえば、生体アミノ酸を思い出す。生体アミノ酸は20種類あるそうだが、塩基配列RNAの遺伝子情報)では4進数表現(G,C,A,U)の3ビット(例えばGCCはアラニンを表す)で表現されているそうである。43=64なので、20種類のアミノ酸を表すには「多すぎる」のだが、この冗長さが進化論的に有利に働いたのであろう。これは、コロナウイルスが変異しやすい理由などと若干関係しているのかもしれないが、よくしらないのでこれ以上の深入はやめておこう.....

さて、putc()はint型なので、返値は整数で与えられる。前回確認したようにC言語の整数表現は4byte、つまり32ビット表現がデフォルトである。putcが返す整数値は、書き込んだASCIIコードになっているらしいことが実験をしてみるとわかる。しかし、ASCIIコードは0.5byteで十分だから、4byteの整数表現では「ゆるゆる」ではないだろうか? 

まずは、文字がASCIIコード(つまり0-127の整数)と等価であることを確認しておこう。putc('B',fp)の代わりに、BのASCIIコードである66を使っても同じ結果になる。つまり、

c = putc(66, fp);

c = putc('B',fp);

は同じ結果となるのである。K&Rにも

int putc( int, FILE *)

と書いてある。返値cを印字させると66と出力される。予想通り、書き込んだASCIIコードが返値になっているのである。

ここで実験である。ASCIIコードをはるかに超える大きな数、たとえば578を使ってファイルに書き込んでみたらどうなるだろうか?つまり、

c = putc(578, fp);

である。コンパイルはうまく通るので、これは間違いではないらしい。書き込んだファイルを覗いてみると、Bが書き込まれており、さらにcを出力すると66となっているではないか。つまり578は66なのである(?)どういうことかというと、現代のマシンではASCIIコードが8ビット(つまり256種)に拡張されており、256を超える大きな数は「周期的に」再定義されているらしいのである。つまり、

578 % 256 = 66

として元に戻ってくるように見えるのだ。ちなみに、putc(834,fp)もputc(21314,fp)も、Bを書き込む命令putc(66,fp)と等価であることを確認した。どれも余りが66となる整数である。

BMPの最初の2byteのデータを読み込む

ということで、bmp画像ファイルの最初の2バイトに書き込んであるという「魔法数」を読み取り、それがBMになっていることを確認してみよう。プログラムする前に、先ほど使ったxxdでBMPファイルをダンプしてみた。

% xxd out.bmp

00000000: 424d 6675 0000 0000 0000 3600 0000 2800  BMfu......6...(.
00000010: 0000 6400 0000 6400 0000 0100 1800 0000  ..d...d.........
00000020: 0000 0000 0000 2c01 0000 2c01 0000 0000  ......,...,.....
00000030: 0000 0000 0000 0000 00ff ffff 0000 00ff  ................
00000040: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
00000050: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
00000060: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
00000070: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
00000080: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
00000090: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
000000a0: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
000000b0: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
000000c0: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
000000d0: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
000000e0: 00ff ffff 0000 00ff ffff 0000 00ff ffff  ................
000000f0: 0000 00ff ffff 0000 00ff ffff 0000 00ff  ................
00000100: ffff 0000 00ff ffff 0000 00ff ffff 0000  ................
(以下省略)

確かに最初の2byteがBMになっている。しばらくはヘッダー情報が続くが、やがて0000と00ffの繰り返しになっていて、これが市松模様に対応していることがわかる。

では読み取ってみよう。

#include <stdio.h>

int main(){
  int c;
  FILE *fp;
  if( (fp = fopen("out.bmp","rb")) == NULL) {
     printf("File not found.\n");
     return 1;
  }
  printf("File found.\n");
  c = getc(fp);
  printf("%c %d %x\n", c, c, c);
  c = getc(fp);
  printf("%c %d %x\n", c, c, c);
  fclose(fp);
  return 0;
}

コンパイルして実行すると、次の結果を得る。

% a.out 

File found.
B 66 42
M 77 4d

out.bmpBMP画像データであることが確認できた!

BMPファイルの初めの14byteには、このような識別データがいくつか書き込まれているとBMPの仕様書には書いてある。2つ目のデータは4byteで構成されるファイルサイズのデータである。3つ目と4つ目は2byteずつ確保されているそうだが、零値が充てあられる規約となっているという。5つ目が4byteの情報で、画像データがどこから始まるか書いてあるそうである。

これらのデータは、BMという2バイトの識別データを読み込んだ要領と同様に読み取っていけばよい。

まずは、画像ファイルのヘッダー部分の読み取りができるようになった。この要領で、画像データ本体も読み取ってみたいと思うが、それは次回のテーマである。

東大数学2024問題6 (part 7): 数値的に問題6を解く

前回のあらすじ

前回までに素数のデータベースの取り込み、エラトステネスの篩を使った素数生成、両者の比較検証、などが済み、いよいよ問題に直接チャレンジできるところまできた。今回でけりを付ける。

コードを分ける

エラトステネスの篩の部分(q6_eratos.py)、関数の定義(q6_func.py)、そして数学の問題を解く部分(q6_main.py)、の3つにプログラムファイルを分けることにする。

まずは数学の問題を解く部分(q6_main.py)。

import q6_eratos as q6e
import q6_func as q6f

Nmax  = 21000000
prime_table = [ True for i in range(Nmax)]
prime_table[0:2] = [False, False, True]

prime_table = q6e.eratos(prime_table,Nmax)

# Negative-integer check
IMAX=101
print("\n f max",q6f.f(-IMAX))
for i in range(IMAX):
    im = -i
    fm = q6f.f(im)
    if fm > 0:
        print(im,fm,prime_table[fm])
# Positive-integer check
IMAX=271
print("\n f max",q6f.f(IMAX))
for i in range(IMAX):
    fi = q6f.f(i)
    if prime_table[fi] :
        print(i,fi,prime_table[fi])

「篩」の表の中に含む整数は21,000,000まで(2100万個)とした。このときの最大の素数は20,999,999と算出されたが、これがちゃんと素数であることはThe PrimePagesによって個別に確認した。この設定に対して、篩が見つけた素数の数は1,329,943個。表の中に含まれる整数の1/10程度のオーダーになっているようだ(これは色々と試しに計算してみて、経験的に会得した推測である)。

負の整数の場合は$n=-101$までを試しているが、実は問題で与えられた関数は$x\rightarrow\infty$で$f(x)\rightarrow x^3$のように振る舞うので、関数$f(x)$自体が素数、つまり正数となるのは限られた$x<0$の値だけであることが(やってみると)わかる。数値計算の結果がその様子を抉り出すようにプログラミングされている。

正の整数については$n=271$まで試している。これは$f(n)<2100万$となるように決めた値である(もちろん、試行錯誤で)。

次は、エラトステネスの篩の部分のコード(q6_eratos.py)である。

def eratos(prime_table, Nmax):
    Pmax = 1
    MaxPrime = 2
    PrimeMax = 2
    while MaxPrime < Nmax:
        n = MaxPrime
        while n <= Nmax/MaxPrime:
            prime_table[MaxPrime * n] = False
            n += 1
        while True:
            MaxPrime += 1
            if MaxPrime > Nmax -1:
                break
            if prime_table[MaxPrime] == True:
                Pmax += 1
                PrimeMax = MaxPrime
                break

    print("Number of primes found: ",Pmax)
    print("Largest prime found: ",PrimeMax)

    return prime_table

global変数とlocal変数の違いをきちんと理解しておかないと、意外に関数というのは使うのが難しい。今回は、素数表とその最大整数値をglobal変数とし、外から(引数として)読み取ることにした。返すのは、関数内で処理した素数表であるが、これ自体はlocalなので、q6_main.pyの中で「値渡し」している。関数内ではglobal変数を参照することはできるが、「いじる(変化させる)」ことはできない。

最後に、問題文で与えられた関数のコード(q6_func.py)である。

def f(x):
    return x*(x**2 + 10*x + 20)

多項式の計算は、本来はホーナーのアルゴリズムというのを使ってやるべきなのだが、計算するのは因数分解できる3次式なので、直感的にわかりやすい表式で済ませた。

これで以上である。実行してみると次のようになった。

% python3 q6_main.py

Number of primes found:  1329943
Largest prime found:  20999999

 f max 1134331
-3 3 True
-4 16 False
-5 25 False
-6 24 False
-7 7 True

 f max 20642341
1 31 True

前半部分が負の整数$n<0$を試した結果である。抜き出された結果は、関数の値が正値をとる場合のみである。予想通り、$n=-3,-4,-5,-6,-7$の5つの場合のみに限定されている。$n\rightarrow -\infty$で$f(n)\rightarrow n^3 < 0$であるから、どんなに大きな$n$を試しても無駄なことになる。計算では$n=-101$までを試しているが、$n=-10$まで調べればもう十分である。ちなみに、$n=-101$のとき、$f(-101)=-930,311$であった。

$n<0$に対し関数値$f(n)$が正となる、これら5つの場合のうち、$f(n)$が素数となるのは$n=-3,-7$の2つだけであった。Trueというのは、エラトステネスの篩で見つかった素数の表の中に含まれている素数である、という意味である。

後半部分は$n>0$の場合であり、こちらは$n=1$の場合のみが適合した。試した整数の最大値は$n=271$で、このとき$f(271)=20,642,341$であった。これは見つけた最大の素数20,999,999に肉薄する値である。$n$のオーダーが100程度であっても、3次関数の値自体は非常に大きくなってしまうのは、数値計算的には頭の痛い内容である。

とはいえ、この計算に要した時間はわずか5.8秒弱である。もうちょっと頑張って、$n=1000$くらいまで試してみてもいいかもしれない。

ということで、問題6の前半はおしまいである。後半に関しては丁寧に分析していくだけであろうから、数値計算の観点からはあまり面白そうに感じられない。来年の受験シーズンまではまだ時間あるので、ゆっくりやっていくことにする(つまり、気が向いたら戻ってくるが、しばらくはやらない予定、ということ)。

東大数学2024問題6 (part 6): データベースとの比較

前回のあらすじ

エラトステネスの篩のコードが書き終わった。最初の30個の素数(2,3,....,113)は正確に計算することができたが、果たしてたまたまなのか、それともちゃんと機能しているのか、あらかじめ公開されている正しい素数のデータベースと比較してみたい。

The PrimePagesのデータベースとの比較

The PrimePagesからダウンロードした最初の1万個の素数を(自分で作ったエラトステネスの古いのコードが)再現できるかどうか確認してみる。

データを読み込む部分と、エラトステネスの篩で計算する部分を融合させる。データベースは、10個の素数をまとめて1行で表現しているので、10個ずつinput()関数で読み込み、プログラムが用意した「配列」(Pythonではlistと呼ぶ)に収納していく形を採用する。

Nt = 10000
Nx = 10
Ny = Nt/Nx

prime_data = [ [0 for i in range(Nx)] for j in range(Ny)]

初期値は0とした。ネストの最も内側にあるデータが連続配置されており、アクセスするには

np = int(prime_data[j][i])

とする。連続した部分は、右側の添字で制御するという形式である。

npの位置にあるprime_tableの値の真偽を調べ、もし真であればnpは正しく計算されデータベースと同じ結果になっていることになる。

print(prime_table[np])

もちろん、np素数でない場合も検査する必要がある。その場合にはFalseが出ればOKである。

ということで、検査部分は次のような感じになった。

ip=0
for nj in range(Ny):
    for ni in range(Nx):
        np = int(prime_data[nj][ni])
        for i in range(ip,np+1):
            if i != np and prime_table[i] == True:
                print(" Error A: ", nj, ni, i, ip, np)
                sys.exit()
            if i == np and prime_table[i] == False:
                print(" Error B: ", nj, ni, i, ip, np)
                sys.exit()
        print(" {0:7d}  {1:7d}".format(np,i), end='')
        ip = np + 1
    print()
print()
print('Success: There is no inconsistency\n')
print("Prime-number search completed.")

計算結果

10000.txtについては問題なくSuccessが出た。

次に100000.txtについても試してみたが、やはりSuccessとなった。

ここまでやってみて、たぶん大丈夫ではないか、という自信が出てきた。 これでデータベースが尽きる500,000あるいは1,000,000くらいまでの素数を利用して、東大の問題を解いてみたいと思う。

東大数学2024問題6 (part 5): エラトステネスの篩のコーディング

前回のあらすじ

ちょっと前のブログ記事の最後に英語を使って「付録」を書いてみたら、海外からのアクセスが来るようになって意外に感じている。東大の入試問題が世界からも注目されているのは、ちょっとした驚きであった。ということで、今回は英語で書く分量を増やしてみて、どんな反応が返ってくるかみてみようと思う。

  • Summary of the previous article: numerical investigations were made on the largest integer that a given computer language can handle, such as the C language, FORTRAN and Python. The largest integer in FORTRAN (based on the 128-bit expression) was already impressive, but the one handled by Python was even more impressive: it has NO limitation in principle. The largest number that Python can handle can go as far as the hardware of a PC/mac configuration can allow it. Amazing.

  • Using pointers, one can investigate what type of representation is employed for the integer expression in computers. However, it is widely believed that the concept of the pointer is intrinsic to the C language, so that it seems at first glance that the other languages are not suitable for this approach.

  • In the last article, I nevertheless managed to find ways to enable this investigation with the pointers even in Fortan95 and Python. A special note is that Python employs representations of integers with 256-bit (32byte) format by default, which is already gigantic. Therefore, storing even a single integer requires 32 bytes, four times greater than an ordinary integer number in C. In this sense, Python is a modern "kid" living in the world where plenty of resources are available for users with reasonable prices.

  • My original plan was to write a code in Python, which is based on the algorythm called the Sieve of Eratosthnes, in order to solve the Maths problem in the entrance examination 2024 of University of Tokyo. After those playful technical attempts and detours, I finally come to a point that I can show you the code.

A code for the Sieve of Eratosthenes

I try to be as pedantic as possible to the original procedure of the Sieve of Eratosthenes.

  1. Namely, a table containing integer numbers, say, up to 10000. Then, 0 and 1 are slashed in the table because they are not prime number by definition.
  2. Next, the first prime number "2" is chosen and its multipliers (including 2 itself) are all marked in the table with a slash. In this process, all even integers except 2 are excluded from the list of prime numbers.
  3. The smallest integer surviving the previous processes is the second prime number, which turns out to be "3". The multipliers of 3 are again slashed in the table.
  4. The then surviving smallest integer (that is, 5) is the third prime number and its multipliers are .... these processes are repeated until the smallest integer surviving the previous processes exceeds the largest integer originally set in the table.
  5. In the end, we have a list of prime numbers included in the table.

Our attempt is to write a code in our own way based on the original algorithm, but this means that there can be errors in the code. To verify what we obtained with the code are correct, we crosscheck the database provided by the PrimePages, the online database created mainly by mathematicians in the University of Tennessee. The PrimePages provides a table containing the first 1000 small primes, so that we use the data for the verification.

I ran the code several times and learned that it is necessary to set the maximum integer in the initial table needs to exceed Nmax=110000 if we demand the first 10000 smallest prime numbers. So, we begin with writing the code as,

Nmax  = 110000
Npmax = 10000

We define a python list named prime_table as the table of integers up to Nmax.

prime_table = [True for i in range(Nmax)]
prime_table[0:2] = [0,0,1]

The list is for a logic variable and the index of prime_table corresponds to the integers in the table. That is, if prime_table[n]=True, then the integer n is one of the prime numbers. If the list returns a false value, then the number does not belong to the prime. The first two elements of the list are set "False" (0) in the initial setting, because n=0,1 are not prime numbers. On the other hand, prime_table[2]=1 (True) is set because n=2 is a prime number.

In the begining, MaxPime=2 is set and its multipliers MaxPrime * n are to be "slashed", that is, set be "False" in the list prime_table. Due to a commutability of a product for two integers, the multipliers can be larger than and equal to MaxPrime. Pmax records the total number of the prime numbers found so far.

Pmax = 1
MaxPrime = 2
while MaxPrime < Nmax:
    n = MaxPrime
    while n <= Nmax / MaxPrime:
        prime_table[MaxPrime * n] = 0
        n += 1
    while True:
        MaxPrime += 1
        if MaxPrime > Nmax-1:
            print("End of Search at tmp MaxPrime =", MaxPrime)
            break
        if prime_table[MaxPrime] == True:
            Pmax += 1
            break

In the end of the programme, a search for the next prime number, that is, MaxPrime, is attempted by referring to the data_table having the True value.

This is the contents of the Sieve of Eratosthenes.

To show the prime numbers we found with this algorithm, the additional section is placed right after the above code.

for i in range(Nmax):
   if data_table[i] == True:
      print(i)

Test run

As a test, we set a pair of small numbers for Nmax and Npmax and run the code.

Nmax = 120
Npmax = 30

The result is as the following.

$python3 eratos.py 

       2       3       5       7      11      13      17      19      23      29
      31      37      41      43      47      53      59      61      67      71
      73      79      83      89      97     101     103     107     109     113

The first 30 prime numbers are detected for Nmax less than and equal to 120. These prime numbers are verified with the database downloaded from the PrimePages! So, it seems the code works properly. But we need to repeat this check for larger numbers, too. That would be the topic of the next article.

東大数学2024問題6 (part 4): ポインタを使って整数表現を調べる

前回のあらすじ

素数を生成する「エラトステネスの篩」アルゴリズムに従って、自分でコードを書いてみようと思ったわけだが、どの程度まで大きな整数が数値的に扱えるのかを確認する必要があると感じ、C言語FORTRANについて調べてみたのが前回である。

C言語についてはlimits.hというヘッダーファイルがシステムに置かれていて、それを参照するだけでよかったのだが、Fortran95に関してはそのようなファイルはなく、計算を用いて確認することになった。128bit表現までが可能で、垓の桁を遥かに凌駕するまで計算することができることがわかった。

一方、Pythonについては「制限がない」という噂をつかみ、今回はその確認を行う予定である。

しかし、その前に、ポインタを使って整数表現が何ビットで表されているか調べる方法を思いついたので、まずはそちらについてメモっておきたい。

C言語におけるポインタの利用

C言語はOS基部やコンピュータシステムそのものにアクセスできる言語であり、本来は科学技術計算というよりは、システム構築のために利用されてきた。以前USB接続機器の制御プログラムを書いたことがあるが、もちろんCで書いた。FORTRANでは不可能だろう(しかし、PythonならRaspberry PIでハードウェア制御することができる!)。

そのため、メモリ管理についてはどんな素人プログラマーでも(理論上は)制御することができる。しかし、その代償として、素人が迂闊なプログラムを書くことで、システム全体の脆弱性を高めてしまったりする危険性もあるから注意が必要だ。

今回は、メモリリークやネットワークハックなどの問題とは無縁のコード内容であるため、それほど緊張する必要はないが、メモリ管理の基本に触れる内容であり、ある意味C言語の「一番危ないところ」の初歩を理解する上では重要な内容である。

まず普通の整数型変数を定義する。この時、配列として定義すると、メモリ上に連続して配置してくれるので計算し易くなる。

limits.hを利用した前回の分析で、普通の整数型は32bitで表現されることがわかった。情報関係でよく利用される単位byteは、8ビット)で表される最大整数28=256の情報量を指しているようである。 \begin{equation} \text{1 byte} \leftrightarrow \text{8 bit} \end{equation}

したがって、C言語の32ビット整数表現は4byteの情報量をもつことになるが、これは28*4=232=4294967296種類の数を取り扱えるという意味になる。

整数の32bit表現をFORTRANでは"integer*4"と宣言するが、これは「4バイトの整数表現という」意味でもあったのだ。いまさらながら、納得である。

さて、では実際に確かめてみよう。

#include <stdio.h>

int main(){
  int qint[5];
  unsigned int uint[5];

  printf("%p %p\n",&qint, &qint[1]);
  printf("%p %p\n",&uint, &uint[1]);
}

qintは4byteつまり32ビット表現の整数である。配列とはメモリ中の連続した情報のことである。したがって、整数型の配列であるqintの先頭アドレスと2つ目のアドレスの間には32ビット(4バイト)の間隔がメモリ中にできていると予想される。実行してみると、

0x16f86f6e4 0x16f86f6e8
0x16f86f6d0 0x16f86f6d4

となった。すべての数は16進数表現である(数字の先頭に0xがつくと16進数という意味になる)。最初の列は配列先頭のアドレスで、2つ目の列が配列2つ目のアドレスである。したがって、その差が整数表現の大きさに対応する。ちなみに、1行目は32ビット表現の場合、2行目が符号なしの32ビット表現の場合である。

最初の行の数字を使って、アドレス間の間隔を計算してみる。数字が長いので下2桁のみを引用する。最初のアドレスから2つ目のアドレスまでを「指折り数えると」、e4,e5,e6,e7,e8となるから、2つのアドレスの感覚は「4」だけ空いている。これは4byteの差、つまり4x8=32ビットの間隔が空いていることに相当する。言い換えれば、一つの整数を格納するのに32ビット必要であるという意味である。

次に2行目のアドレスも見てみる。d0, d1, d2, d3, d4であるからやっぱり4byteの間隔である。これは「符号なし」の場合であるが、符号ありとなしの場合は、先頭ビットを符号に当てるかどうかの違いであるから、一つの数字を表現するビット数(バイト数)は共通になるのは当然である。

16進数の計算は慣れてないと面倒だし、不要に難しく見える。じゃあ10進数でやってみたらいいじゃないか、という意見が飛んできそうであるが、これは不可能ではない。printfの識別子%pを%dに変えるだけの話である。10進数の方が読みやすいと思う人はこれでいいと思う。ただ、gccコンパイルしようとするとwarningが出てしまう。どうもコンピュータ関係の人は、10進数が嫌いのようである(もしかすると、理論上では明確な違いがあるのかもしれないが、素人に近い私にとってはその機微については予想すらもできないのである)。

1840707300 1840707304
0x16db6f6e4 0x16db6f6e8

(同じ内容を上の行では10進数で、下の行では16進数で表したもの。どちらも4バイトの間隔がある。)

配列のアドレスは、厳密には「ポインタ」ではない。ポインタの宣言の仕方は

int *pntr;

である。ポインタというのは「メモリの塊」の先頭部分のアドレス情報のことであるから、実質、配列の先頭のアドレスだと解釈しても(一般のプログラマにとっては)実害はないだろう。OSやカーネルの開発をやっているプログラマからしてみれば、ポインタの実装については配列と違うかもしれないから、区別して理解するのは重要だろう。実際、ポインタが示すアドレス情報自体は、配列とは無関係のアドレスに格納されている。上のプログラムでポインタを使って、整数の情報を引き出す場合には、たぶん次のようにやるのだと思う。

#include <stdio.h>

int main(){
  int qint[5];
  int *pntr;
  
  printf("%p %p\n",&qint, &qint[1]);

  printf("%p\n",pntr);
  ++pntr;
  printf("%p\n",pntr);

  pntr = qint;
  printf("%p %p \n", pntr,qint);

  ++pntr;
  printf("%p %p \n", pntr,qint);

}

実行結果は次のようになる。

0x16ef4f6e4 0x16ef4f6e8
0x100ebc060
0x100ebc064
0x16ef4f6e4 0x16ef4f6e4 
0x16ef4f6e8 0x16ef4f6e4 

最初の行は、これまでのプログラムと同じ内容である。整数型の配列の先頭と2つめのアドレスである。その間隔は4バイト、つまり32ビットである。

次に、ポインタ $*pntr$が定義され、そのアドレスをそんまま印字させている。結果は0x100ebc060となり、整数型の配列とは無関係の、新しいアドレスが表示されている。これが宣言当初にポインタが割り当てられたアドレスであるが、どこも向いてない無意味なアドレスである。ただし、この無意味なポインタを一つ「横に」ずらすと($++pntr$)、整数一つ分だけ「ずれる」ことが確認でき、そのずれ幅は4バイト(32ビット)になっていることがわかる。

最後に、ポインタを配列の先頭にセットする。ポインタが示す「方向」は配列qintの先頭であり、「意味のある方向」を向いたことになる。ポインタを印字すると、配列の先頭アドレスに一致していることがわかる(0x16ef4f6e4)。

ポインタを一つ「横にずらす」と、qint[1]のアドレスにポインタの位置はずれる(0x16ef4f6e8)。そのずれ幅は、またもや4バイト(32ビット)である。

Fortran95における「ポインタ」の利用

C言語と同じようなことをFortran95でやってみたいと思っても、なかなか難しい。というのは、Fortranの本質は、システムの詳細にあまりこだわらなくても、簡単に数値計算が実施できることだからである。したがって、意図的にメモリ構造やそのアドレス情報といったものにアクセスできないように言語が設計されている。

しかし、時にはメモリ構造を念頭に置きながらコーディングした方が、効率的なプログラムを書ける場合がある。その場合、Fortran95でも、どのように変数がメモリ空間に配置されているか知りたいケースも出てくるだろう。

いろいろ調べてみると、「非正規な」関数ではあるが、loc()という関数が利用できることがわかった。C言語のようにポインタ操作ができるわけではないが、配列や変数のアドレスを印字する関数だというから、C言語の&に相当するものである。

Fortranの場合は、4byte表現, 8byte表現, 16byte表現を一気に確認してみたい。使ったコンパイラgnu gfortranであるが、非正規の内部関数loc()を含んでいた。

program pointer_address
  implicit none
  integer :: ad1, ad2
  integer, dimension(5) :: quadInt
  integer*8, dimension(5) :: octInt
  integer*16, dimension(5) :: hexInt

  ad1 = loc(quadInt(1))
  ad2 = loc(quadInt(2))
  print '(z10 z10 i4)',ad1, ad2, ad2-ad1
  print *

  ad1 = loc(octInt(1))
  ad2 = loc(octInt(2))
  print '(z10 z10 i4)',ad1, ad2, ad2-ad1
  print *

  ad1 = loc(hexInt(1))
  ad2 = loc(hexInt(2))
  print '(z10 z10 i4)',ad1, ad2, ad2-ad1
  print *
end program pointer_address

Fortran95で16進数表現の表示をする場合は、フォーマットで指定する(記号zをprint文で利用)。

計算結果は次のような感じになる。

  6DCFF6C0  6DCFF6C4   4

  6DCFF6D8  6DCFF6E0   8

  6DCFF700  6DCFF710  16

配列の先頭と2つ目の間が、メモリ空間中でどのくらい空いているか計算したものである。予想通り、4byte, 8byte, 16byteとなり、整数の表現が32bit, 64bit, 128bitになっていることが確認できた。

Pythonの場合

Pythonにもloc()に似た関数が存在し、id()という。同じようなことをやってみよう。ただし、Pythonでは整数はデフォルトで無制限に扱えるため、64bit表現とか32bit表現とか、配列の宣言で規定する必要がない。これがPythonが人気があり、他の言語に比べて優れている点であろう。とにかく、変数型の宣言をしないとプログラムというのは、こんなにも簡単になるのである!

i1=[1,2]
ad1=id(i1[0])
ad2=id(i1[1])
print(hex(ad1),hex(ad2),int(ad2)-int(ad1))

16進数表示するための関数hex()を利用して印字しているので、結果は次のようになった。

0x101738e48 0x101738e68 32

アドレスの引き算(差)は0x20であるから、これを10進数に変換すると$2\cdot 16^1+0\cdot 16^0=32$となる。

つまり、Pythonでは最初から32byte表現、つまり256 bit表現を採用しているらしいのである。

ある意味、これはメモリ容量をかなり逼迫するような設計である。たった一つの整数を保存するだけでも256ビットを取られてしまうのである。明らかに、pythonは現代の計算機環境に根ざした「資産が潤沢なプログラマ向け」の言語である。

東大数学2024問題6 (part 3): エラトステネスの篩(ふるい)と数値計算で扱える最大の整数

前回のあらすじ

前回のブログでは、The PrimePagesで公開されている素数のデータベースをダウンロードし、これを読み込む部分を(pythonを用いて)書きあげるところまで完了した。

今回は、いよいよエラトステネスの篩と呼ばれる、古来より知られる有名なアルゴリズムを用いて、最初の10000個(1万個)の素数を計算してみようと思う。(計算結果は、The PrimePagesからダウンロードしたデータベースと比較し、正しく計算できたか確認する計画。)

エラトステネスの篩

高校数学でやったという記憶はないが、最近の高校の教科書では紹介されているようである。 うろ覚えながら、その手順は次のような感じとなる。

  1. 0は除く。
  2. 1は除く。
  3. 素数を見つけたい範囲を決める(たとえば0から100までの中に含まれる素数、といった具合)
  4. 2が最初の素数
  5. 最初の素数(2)の倍数を、素数を見つけたい範囲から消す(2,4,6,8,...100、つまり100以下の偶数を消す)。
  6. 4.のプロセスで消去されなかった最初の数が、次の素数となる(3)。
  7. 2番目の素数(3)の倍数を、素数を見つけたい範囲から消す(3,6,9,.....,99)。
  8. ここまでで生き延びた最小の整数を3番目の素数とする(5)。
  9. 3-5のプロセスを繰り返し、最大の素数が100を超えたらアルゴリズムを終了する(97の次の素数が101となるので、最大の素数は97として計算終了)。

3-5のプロセスがメインの繰り返しジョブとなり、8の終了判定が重要である。

このアルゴリズムをより効率よく書き換えて、数値的に計算が楽になる有名なコードもおそらく存在するであろう。しかし、今回はまず自分でプログラムを組んでみることにする。

pythonによる実装

pythonはどの程度まで大きな整数を扱うことができるのであろうか?FORTRANやCではその最大値が決まっていて、gccの場合はlimits.hで調べることができるという。

www.gnu.org

実際にコードを書いて確かめてみた。

#include <stdio.h>
#include <limits.h>

int main(){
  printf("Min for signed int: %d\n", INT_MIN);
  printf("Max for signed int: %d\n", INT_MAX);
  printf("Max for unsigned int: %u\n", UINT_MAX);
}

実行結果は次のようになった。

Min for signed int: -2147483648
Max for signed int:  2147483647
Max for unsigned int: 4294967295

素数の計算では負数は出てこないから、unsigned intを用いることができる。となると、扱える最大の整数は4,294,967,295、つまり43億程度まで扱えるということである(裏金の帳簿管理プログラムにはギリギリ使える値だが、国家予算の帳簿プログラムには使えないかも)。この値は32ビットの最大値に対応しており、bcで確かみるとその通りである。

$ echo '2^32' | bc -ql

> 4294967296

符号付きだと1ビット持って行かれてしまうので、31ビットとなり、上の値の半分となる。

$ echo '2^31' | bc -ql 

> 2147483648

ここで実験をやってみる。最大の整数に対して、あえて1を足してみるのである。水をなみなみと注いだグラスに、一滴水を垂らすと何が起きるであろうか?たぶん、これはマシンによって結果は異なると思うが。

#include <stdio.h>
#include <limits.h>

int main(){
  int n;
  unsigned int m;

  n = INT_MAX;
  m = UINT_MAX;

  printf("Out of bounds: %d %d\n",n,n+1);
  printf("Out of bounds: %u %u\n",m,m+1);
}

ここでのポイントは、符号なし整数の識別子は%uと表す点である(今日まで知らなかった....)。結果は次のとおり。

Out of bounds: 2147483647 -2147483648
Out of bounds: 4294967295 0

どうやら、丸い地球を平らだと思って探検に出かけ、水平線を越えてみたら出発点に戻ってしまった、という感じである。

FORTRANでは整数を使った数値計算はあまりやらないのだが、たぶんC言語と同じような状況になっているはずである。特にgfortranはgccの付加体として定義されているからなおのことである。(もしかするとgfortranなら#include <limits.h>を受け付けるかもしれないと思って、やってみたがダメであった。)

いろいろとFortranについての資料を探してみたが、整数型についての説明はあまりなかった。これは当然で、科学技術計算では浮動点小数計算が主流であり、あまり整数を使った計算はやらないからだ。とはいえ、最近では素数を利用した暗号技術などが重要になってきたので、Fortranで色々計算している人も多少はいるのではないかと想像する。また、デジタル情報の処理は基本的に整数でできるから、工学系の人で整数を扱うFortranプログラミングをしている人もいるのかもしれない。とはいえ、調べた限りでは、手持ちのgfortranの整数がどの程度の範囲まで扱えるのか調べきることができなかった。

こういう時は、自分で試してみるしかない。

Fortran95で整数の最大値を確かめる

たとえば、整数の4ビット表現の場合、ビットが4つあるので \begin{equation} [a][b][c][d] \end{equation} となる。ビット表現は2進数なので、a,b,c,dは1か0の値しかとれない。たとえば、 $[1][0][0][1]$といったような感じになる。これを10進数に変換する方法は(高校数学でも解説があるように) \begin{equation} 1\cdot 2^3 + 0\cdot 2^2 + 0\cdot 2^1 + 1\cdot 2^0 \end{equation} つまり、 $8+0+0+1 =9$である。したがって、4ビットの一般的な表現を10進数に変換する公式は \begin{equation} [a][b][c][d] = a\cdot 2^3 + b\cdot 2^2 + c\cdot 2^1 + d\cdot 2^0 \end{equation} である。

4ビットの場合は、2の3乗、つまり(二進数の)3桁までが登場する。したがって、4ビット表現で表せる最大の整数は $8+4+2+1=15$となる。これは$2^4-1$とも書けるが、よくよく考えると2の等比級数(n次まで)になっているから当然である。 \begin{equation} S_n = 1+2+2^2+\cdots +2^n = \frac{1-2^{n+1}}{1-2} \end{equation} しかし、この「当然」はプログラミングにおいては「当然」ではない。なぜなら、もし4ビット表現が採用されていた場合、$2^4$を計算した瞬間にオーバーフローが発生し、そこから1を引いたとしても正しい値に戻ってきてくれるかどうかは保証されないからである。したがって、上の等比級数の公式を「鵜呑み状態」で利用するのは危険である。より安全なのは、等比級数の定義式、つまり$2^3$までの等比数列文字通り足しあげるアプローチである。

数値実験のためのプログラムは次のようなものである。

program max_int_limit
  implicit none
  integer :: n
  integer  :: r,sum
  !integer*4 :: r,sum
  !integer*8 :: r,sum
  !integer*16 :: r,sum

  sum = 0
  r = 1
  do n=1,35
     sum = sum + r
     print *,n,r,sum, sum + 1
     r = r*2
  enddo
  
end program max_int_limit

integerの宣言にはさまざまなタイプがあるようだ。デフォルトのもの、$*4, *8, *16$などがある。C言語にある「符号なし」の宣言は見つからなかった。

プログラム中の変数の意味は、$n$が次数を表し(ビット数に相当)、$r$が等比数列、そして$sum$が等比級数である。

プログラム後半のループでは、$n$を一つずつ増加させながら、整数が正しく計算され表現できているかを確かめていく。破綻した表現が出た場所が、Fortran95が扱う整数の限界(最大値)となる。最大値がどこかを確かめれば、何ビット表現になっているかわかる、という狙いである。

さて、gfortranでコンパイルし、実行してみる。デフォルト、4, 8, *16の順番でやってみる。

デフォルトのinteger :: r, sumの場合

まずは$r$と$sum$を、デフォルトのintegerとして宣言して計算した場合を見てみる。結果は次のようになった。 (出力は、$n, r_n, S_n, S_n+1$、つまり次数、等比数列級数級数に1を足したもの、の順番)

           1           1           1           2
           2           2           3           4
           3           4           7           8
           4           8          15          16
           5          16          31          32
           6          32          63          64
           7          64         127         128
           8         128         255         256
           9         256         511         512
          10         512        1023        1024
          11        1024        2047        2048
          12        2048        4095        4096
          13        4096        8191        8192
          14        8192       16383       16384
          15       16384       32767       32768
          16       32768       65535       65536
          17       65536      131071      131072
          18      131072      262143      262144
          19      262144      524287      524288
          20      524288     1048575     1048576
          21     1048576     2097151     2097152
          22     2097152     4194303     4194304
          23     4194304     8388607     8388608
          24     8388608    16777215    16777216
          25    16777216    33554431    33554432
          26    33554432    67108863    67108864
          27    67108864   134217727   134217728
          28   134217728   268435455   268435456
          29   268435456   536870911   536870912
          30   536870912  1073741823  1073741824
          31  1073741824  2147483647 -2147483648
          32 -2147483648          -1           0
          33           0          -1           0
          34           0          -1           0
          35           0          -1           0

$n=30$までは順調に計算は進んでいるが、$n=31$のとき、sum+1の値がおかしくなっている。つまり、gfortranのintegerは、32ビット表現になっていて、扱える最大の整数値は$2147483647 = 1+2+4+\cdots 2^{31}$というわけである。

$n=32$ではすでにオーバーフローしてしまい、$r$自体の値が変になっていて、級数(和)を計算してももはや何の意味もないことが確認できる。

integer*4 の場合

次にinteger4を試してみると、上の計算結果と同じとなった。つまり、デフォルトのintegerはinteger4と同じであり、この2つは32ビット表現になっているのである。

integer*8の場合

ではinteger*8を試してみよう。おそらく64bit表現になっているはずだから、nのループの最大値を65にしてみる。

program max_int_limit
  implicit none
  integer :: n
  integer*8 :: r,sum

....

  do n=1,65
     ...
  enddo

end program max_int_limit

結果はこうなった(前半部分省略)

          56    36028797018963968    72057594037927935    72057594037927936
          57    72057594037927936   144115188075855871   144115188075855872
          58   144115188075855872   288230376151711743   288230376151711744
          59   288230376151711744   576460752303423487   576460752303423488
          60   576460752303423488  1152921504606846975  1152921504606846976
          61  1152921504606846976  2305843009213693951  2305843009213693952
          62  2305843009213693952  4611686018427387903  4611686018427387904
          63  4611686018427387904  9223372036854775807 -9223372036854775808
          64 -9223372036854775808                   -1                    0
          65                    0                   -1                    0

予想通りである。$n=63$の$sum$までは計算できているが、そこに1を足すと変なこと(つまり負数になってしまう)が発生し、その後の計算は「もうダメ」状態である。つまり、integer*8は64bit表現である。「900京」をこえる整数が扱えるというのは、すごいことだ!

この上となると「垓」の世界である(アボガドロ数の世界へ近づくということでもある!)。integer*16でコンパイルエラーが出なければ、とんでもなく大きさな整数をgfortranは扱えるということになりそうである。果たしてどうだろうか?

integer*16にトライ!

予想ではinteger*16は128bit表現になるはずである。$2^{128}=10^x$とおいて$x$を計算すると \begin{equation} x = 128\log_{10}2 \sim 128\cdot 0.3010 = 38.53... \end{equation} つまり$10^{38}$の世界である!アボガドロ数のオーダーが$10^{23}$であったから、それを遥かに上回る「宇宙の果て」の世界に似た感じがする大きさの整数である。

実行結果は次のようになって、ちゃんと「宇宙の果て」の雰囲気を醸し出す巨大な整数が表示されたのであった! まちがいなく、128bitである!

ということで、どうやらFortran95で整数の研究をやっている数学者/情報工学者がいそうな感じに見えてきた。しかし、暗号資産とか電子通貨といったものはこれを超える桁での素数を研究しているはずだ。

pythonではどうなのか?

大きな素数の研究をするならFortran95でやった方がよいのだろうか?といった気持ちが芽生えてきたわけだが、果たして我々が利用しているPythonはどうなっているのであろうか?

ググってみると、どうやら「限界なし」というおそろしいことになっているらしいのである。

python-codes.net

この辺りは次回にもう少し詳しく見てみたいと思う。有用そうな資料としては次のブログも注目である。

note.nkmk.me

エラトステネスの篩のコードまでいけなかった....

コードはできているのだが、最大整数の扱いのところで知らなかったことが山ほど出てきてしまい、今回はそこで一杯一杯となってしまった。次回でPythonでの最大整数の扱いについて議論した後、その次でエラトステネスへいけるのではないかと思う。