#navi_header|C言語系| お題：素因数分解を行うfactorコマンドで、素因数分解を保証している最大値を調査せよ。 #more|| #outline|| ---- * まずfactorを使ってみる $ factor 5 bash: factor: command not found 「デーモン君のソース探検」によれば、"/usr/games"にPATHをのばせば使えるとのこと。 ・・・ところが、インストールセットで選択し忘れたのか、"/usr/games"以下にfactorコマンドが存在しない。 # find /usr/games /usr/games /usr/games/hide ソースコードは "/usr/src/games/factor/" ディレクトリにあった。 # cat Makefile # $NetBSD: Makefile,v 1.7 1997/10/22 04:43:26 lukem Exp $ # @(#)Makefile 8.1 (Berkeley) 5/31/93 PROG= factor SRCS= factor.c pr_tbl.c CPPFLAGS+=-I${.CURDIR}/../primes MAN= factor.6 MLINKS+=factor.6 primes.6 .PATH: ${.CURDIR}/../primes .include Makefileを見てみると、 "/usr/src/games/primes/pr_tbl.c" も使っている。ここでmakeを実行してみたが、システムビルド用の一時的なtool群が生成されていなかった為失敗。 ソースがあれば何とかなると思い、"/usr/src/games/{factor|primes}" を一般ユーザの適当なディレクトリにコピー。 $ cp -rp factor ~/lang.c/daemons_src_adv/ $ cp -rp primes ~/lang.c/daemons_src_adv/ で、試しにコピーした先でprimesをmakeしてみたらそのままコンパイル出来てしまった。 $ cd ~/lang.c/daemons_src_adv/primes $ make cc -O2 -Werror -c pattern.c cc -O2 -Werror -c pr_tbl.c cc -O2 -Werror -c primes.c cc -o primes pattern.o pr_tbl.o primes.o -lm $ ./primes 3 10 3 5 7 同様に、コピーした先でfactorコマンドもコンパイルに成功した。 $ make cc -O2 -Werror -I/home/msakamoto/lang.c/daemons_src_adv/factor/../primes \ -c factor.c cc -O2 -Werror -I/home/msakamoto/lang.c/daemons_src_adv/factor/../primes \ -c /home/msakamoto/lang.c/daemons_src_adv/factor/../primes/pr_tbl.c cc -o factor factor.o pr_tbl.o groff -Tascii -mtty-char -mandoc factor.6 > factor.cat6.tmp && \ mv factor.cat6.tmp factor.cat6 $ ./factor 1000 1000: 2 2 2 5 5 5 準備が出来たところで、「デーモン君のソース探検」に出てくるフェルマーの素数 f(n) = 2^2^n + 1 を試してみる。 f(1) = 2^2^1 + 1 = 5 f(2) = 2^2^2 + 1 = 17 f(3) = 2^2^3 + 1 = 257 f(4) = 2^2^4 + 1 = 65537 f(5) = 2^2^5 + 1 = 4294967297 → 5, 17は明らかに素数なので、f(3), f(4), f(5)を試してみる。 $ ./factor 257 257: 257 $ ./factor 65537 65537: 65537 $ ./factor 4294967297 factor: 4294967297: Result too large 「デーモン君のソース探検」によれば、f(5)については2の32乗 + 1、つまり32bitで表現出来る値を超えている。このため32bitマシンでは計算しきれない。かといって、64bitマシンで計算出来るかというとその保証は無いらしい。「デーモン君のソース探検」ではなぜ64bitマシンで計算出来る保証がないか、という点が詳しく解析されている。こちらも、それに沿ってソースを読んでいく。 なお手元には64bitマシンが無いため、実際に動かすことが出来ない。また、手元のNetBSD1.6は32bit Pen4上のVMware内で動かしているが、longのサイズは4バイト = 32bitだった。 #pre||> $ cat sizeofs.c #include int main() { printf("sizeof(long) = %d\n", sizeof(long)); return 0; } $ cc -Wall -o sizeofs sizeofs.c $ ./sizeofs sizeof(long) = 4 ||< * factor/factor.cのmain関数より、"Result too large"を追跡する。 factor/factor.cのmain関数より、まずコマンドライン引数または標準入力から数値を読み込む部分を見てみる。 #code|c|> int main(argc, argv) int argc; char *argv[]; { ubig val; int ch; char *p, buf[100]; /* > max number of digits. */ /* ... */ if (argc == 0) for (;;) { if (fgets(buf, sizeof(buf), stdin) == NULL) { /* ... */ } /* ... *p が終端に設定される */ val = strtoul(buf, &p, 10); /* start = buf, end = p, base = 10 */ /* ... */ } else /* ... */ val = strtoul(argv[0], &p, 10); /* ... */ } ||< ここで、次のエラーメッセージからstrtoul()がERANGEを返していることが分かる。 factor: 4294967297: Result too large grepしてみる： $ cd /usr/include $ grep -r "Result too" * sys/errno.h:#define ERANGE 34 /* Result too large */ よって入力値はunsinged long, 2の32乗 - 1までと分かる。 また、valが ubig val; と定義されている。ubigは "primes/primes.h" で "unsinged long" がtypedefされたものである。 "primes/primes.h" : #code|c|> /* ubig is the type that holds a large unsigned value */ typedef unsigned long ubig; /* must be >=32 bit unsigned value */ #define BIG ULONG_MAX /* largest value will sieve */ /* bytes in sieve table (must be > 3*5*7*11) */ #define TABSIZE 256*1024 ||< ということで、試しに境界値付近を調べてみる。 ex1) : 2 ^ 32 - 1 = 4294967295 $ ./factor 4294967295 4294967295: 3 5 17 257 65537 ex2) : 2 ^ 32 = 4294967296 $ ./factor 4294967296 factor: 4294967296: Result too large ん～～、2の32bitは4294967296ではあるけれど、"0"が一つ分あるため、表現出来るのは 2^32 - 1 まで。2^32はstrtoul()がERANGEになってしまう。 とはいえ、64bitマシンの場合は恐らく 2^32+1 も問題なくstrtoul()で変換出来るはず。 * factor/factor.cのpr_fact()関数で素因数分解の処理を追跡する。 とりあえず "Result too large" の調査はここまでで充分なので、続いて実際に素因数分解をしている箇所を見てみる。 "factor/factor.c"のpr_fact()関数： #code|c|> void pr_fact(val) ubig val; /* Factor this value. */ { const ubig *fact; /* The factor found. */ /* ... */ /* Factor value. */ (void)printf("%lu:", val); for (fact = &prime[0]; val > 1; ++fact) { /* Look for the smallest factor. */ do { if (val % (long)*fact == 0) break; } while (++fact <= pr_limit); /* Watch for primes larger than the table. */ if (fact > pr_limit) { (void)printf(" %lu", val); break; } /* Divide factor out until none are left. */ do { (void)printf(" %lu", *fact); val /= (long)*fact; } while ((val % (long)*fact) == 0); /* Let the user know we're doing something. */ (void)fflush(stdout); } (void)putchar('\n'); } ||< ここで "prime"という配列や "pr_limit" というのが出てくる。これは factor.c の冒頭でextern宣言されており、実体は "primes/pr_tbl.c" で定義されている。(factor/Makefile中で "../primes/" ディレクトリ指定とともに、ソースターゲットに指定されている) まずprime[]は、 const ubig prme[] = { ... } として、2^16 + 1 = 65537までの素数が入ったテーブルになっている。 続いて pr_limit は、 /* pr_limit - largest prime in the prime table */ const ubig *pr_limit = &prime[(sizeof(prime)/sizeof(prime[0]))-1]; として、つまり配列の終端要素のポインタを指すようになっている。ここまで分かれば、factor.cのpr_fact()のコードは読める。 まず、次のforの初期化部分で fact に primeテーブルの先頭の素数を指すように初期化する。つまり一番小さい素数からtryする。 for (fact = &prime[0]; val > 1; ++fact) { 続いて次のループで、小さい素数から順に、0で割り切れるか試していく。割り切れるものがあれば、一旦ループを抜ける。 do { if (val % (long)*fact == 0) break; } while (++fact <= pr_limit); もしこの時点で pr_limit を超えてしまっていれば、用意した素数テーブルだけでは処理出来なかったものとして、元の数値を表示して終了する。 /* Watch for primes larger than the table. */ if (fact > pr_limit) { (void)printf(" %lu", val); break; } pr_limitの範囲内なら、見つかった素数を使って割りきれなくなるまで分解していく。 /* Divide factor out until none are left. */ do { (void)printf(" %lu", *fact); val /= (long)*fact; } while ((val % (long)*fact) == 0); 割り切れなくなったら、またforループの先頭まで戻り、より大きい素数で割りきれるかtryしていく。 * 素因数分解可能な最大値の調査 「factorで保証している、素因数分解可能な範囲はどこまでか？」という本当のテーマが、ようやくここで「デーモン君のソース探検」に出てくる。ひとまず64bitマシンと想定し、strtoul()でも以下の値を変換出来るものと仮定しておく。(long=64bit) 4294967296 (= 2^32) 4294967297 (= 2^32 + 1) まずprimeテーブルが65537以下の素数を定義しているということから、65537以下の素数の合成数までは素因数分解できると考えられる。 続いて、65537以下の素数で割り切ることが出来ない場合、表示上は素数として表示する。しかしそれが素数であることが保証されるのは 65537 x 65537 の範囲までとなる。なぜなら、65537より大きな素数が因数にあったとしても、primeテーブルの定義上、そこまで計算できないからである。 割り切れない場合に素数であるのが65537 x 65537以下であることから、 1. 1回目に割り切れた数が65537以下の素数 2. 2回目は割り切れなかったが、65537 x 65537以下 の場合も素因数分解できたことになる。「デーモン君のソース探検」で "f(5) = 2^32 + 1" に対して % factor 4294967297 4294967297: 641 6700417 となり、「素因分解出来た」に対して「たまたま出来ただけ」という風に書かれているのはこのケースに該当する。 実のところ、"primes/pr_tbl.c" を読んでいた段階で次のようなコメントがあり、ヒントが書かれていたりした。 #pre||> /* * prime - prime table * * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo * * chongo /\oo/\ * * We are able to sieve 2^32-1 because this table has primes up to 65537 * and 65537^2 > 2^32-1. */ ||< 今回のお題については、ここまで。 #navi_footer|C言語系|