C言語系/「デーモン君のソース探検」読書メモ/06, factor(1)

• [ Prev ]
• [ Next ]
• [ C言語系 ]

お題：素因数分解を行うfactorコマンドで、素因数分解を保証している最大値を調査せよ。

• まずfactorを使ってみる
• factor/factor.cのmain関数より、"Result too large"を追跡する。
• factor/factor.cのpr_fact()関数で素因数分解の処理を追跡する。
• 素因数分解可能な最大値の調査

まず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 <bsd.prog.mk>

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だった。

$ cat sizeofs.c
#include <stdio.h>
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関数より、まずコマンドライン引数または標準入力から数値を読み込む部分を見てみる。

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" :

/* 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()関数：

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" を読んでいた段階で次のようなコメントがあり、ヒントが書かれていたりした。

/*
 * prime - prime table
 *
 * By: Landon Curt Noll   chongo@toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
 *
 *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
 *
 * We are able to sieve 2^32-1 because this table has primes up to 65537
 * and 65537^2 > 2^32-1.
 */

今回のお題については、ここまで。

• [ Prev ]
• [ Next ]
• [ Up ]
• [ C言語系 ]