素因数分解

プログラム

use strict;
my @prime = (2, 3);

while (1) {
  print 'Input number(2 - 999,999,999 or expr[+, -, *, **]): ';
  my $n = <STDIN>; chomp $n; $n =~ s/,//g;
  $n = eval $n if $n =~ /[1-9]\d* *(\+|-|\*|\*\*) *[1-9]\d*/;
  exit if $n !~ /^[1-9]\d{0,8}$/;
  print "$n";
  my %result = ();

  for (my $i = 0; $prime[$i] ** 2 <= $n; $i++) {
    unless ($n % $prime[$i]) {
      $result{$prime[$i]}++;
      $n /= $prime[$i]; $i--;
      next;
    }
    if ($i == $#prime) {
      my $p = $prime[-1] + 2;
      for (my $j = 1; $prime[$j] ** 2 <= $p; $j++) {
        unless ($p % $prime[$j]) {
          $p += 2; $j = 1;
          redo;
        }
      }
      push @prime, $p;
    }
  }

  if (%result) {
    $result{$n}++;
    my $expr = join ' x ',
               map { $result{$_} == 1 ? "$_" : "$_**$result{$_}" }
               sort { $a <=> $b } keys %result;
    print " = $expr\n\n";
  } else {
    print ": prime\n\n";
  }
}

プログラムの説明

今回は、ごく易しいプログラムです。まずは、素因数分解を取り上げます。 このプログラムを実行すると、プロンプトを出して入力を促します。入力できる桁数は9桁までで、 カンマで区切ることができます。また、式で入力することもできます (厳密なチェックはしていません)。 使える演算子は、+, -, *, ** のみです。式で入力した場合は、計算の結果が9桁を越えたり、 負の数になったり、小数点が含んではいけません。入力は、以下のようにします。

Input number(2 - 999,999,999 or expr[+, -, *, **]): 999,999,999
Input number(2 - 999,999,999 or expr[+, -, *, **]): 7 ** 5 + 1
Input number(2 - 999,999,999 or expr[+, -, *, **]): 2 * 3 * 5 * 7 - 1
Input number(2 - 999,999,999 or expr[+, -, *, **]): 2 * 3 * 5 * 7 * 11 * 13 + 1

プロンプトに対して整数または式を入力すると、素因数分解した解 (あるいは素数である旨) を表示して、再び入力待ちになります。プログラムを終了するには、(本当は何でもいいのですが) quit, end 等と入力します。ところで、最後のような式を入力した場合には、何か新しい発見があるかもしれません。 (というのも、素数の無限性に関する有名な式 x = 2 * 3 * ... * p + 1 の x は、必ず本物の素数であると長い間間違って覚えていました。もちろん素数の場合もあるのですが、p よりも大きな素数の約数を持つ場合もあることに気がつきました。それでも、p よりも大きな素数が存在することの証明になっているのですね。閑話休題。)

素因数分解をするには、入力された整数を 2, 3, 5 ... と小さな素数から順に割っていきます。 そのために、素数を格納する配列 @prime を用意します。また、分解した結果を格納するハッシュ %result も用意します。%result のキーには素数が入り、値には分解数が入ります。分解は、次のコードで行っています。

for (my $i = 0; $prime[$i] ** 2 <= $n; $i++) {
  unless ($n % $prime[$i]) {
    $result{$prime[$i]}++;
    $n /= $prime[$i]; $i--;
    next;
  }
  ...
}

ここでは、for ループを使って分解を進めていきます。分解を進めて $n が素数になったら、それで分解は終了します。ここで使っている条件式 $prime[$i] ** 2 <= $n では、素数を判定しています。 素数の判定には、「平方根以下の素数に約数がない場合は素数である」という定義を使っています。 この定義通りの条件式は、$prime[$i] <= int(sqrt $n) になりますが、両辺を2乗してあるだけで意味は同じです。

個々の分解の処理は、unless ブロックの中で行います。分解できた素因数を %result に登録して、$n を再設定します。そして、次の分解に移るのですが、素因数分解の場合には分解できた素因数で (分解できなくなるまで) 再び試みる必要があります。このような場合には、言葉の意味からも redo を使うのが適当と思いますが、ここでは制御変数を1つ戻して next を使っています。redo を使ったときのコードは、次のようになります。

for (my $i = 0; $prime[$i] ** 2 <= $n; $i++) {
  unless ($n % $prime[$i] || $prime[$i] == $n) {
    $result{$prime[$i]}++;
    $n /= $prime[$i];
    redo;
  }
  if ($#prime == $i) {
    last if $prime[$i] ** 2 > $n;
  ...
}

太字で示した最初の部分 $prime[$i] == $n は、素数を素数自身で割ることを防ぐためです。これは、最大値の素因数の分解数が2以上の場合に発生します。 2番目の部分 last if $prime[$i] ** 2 > $n は、まだ説明していませんが余分な素数を生成するのを防止するものです。太字で示した2つの部分は、for の条件式を通すことで必要がなくなります。

プログラムの2行目では、配列 @prime に素数 2 と 3 だけしか入れていません。配列 @prime を使いきってしまって素数が不足する場合は、その都度1つずつ生成します。 以下のコードで、次の素数を生成しています。

for (my $i = 0; $prime[$i] ** 2 <= $n; $i++) {
  ...
  if ($i == $#prime) {
    my $p = $prime[-1] + 2;
    for (my $j = 1; $prime[$j] ** 2 <= $p; $j++) {
      unless ($p % $prime[$j]) {
        $p += 2; $j = 1;
        redo;
      }
    }
    push @prime, $p;
  }
}

配列 @prime に残りがない ($i == $#prime) ときには、@prime の最後の要素に 2 をプラスして素数探しを開始します。ここでも、for ループを使って平方根以下の素数に約数がないか調べます。約数があって素数でないと分かったら、さらに $p に 2 を加えて素数探しを見つかるまで続行します。素数が見つかったら、配列 @prime に見つかった素数を追加して、外側の for ループに移って追加された素数を使って素因数分解を試みます。 このようにして、最終的に素因数分解が終了します。

配列 @prime は、毎回最初から素数を生成するのではなく、前に生成した素数を保持しています。 9桁の最大の素数 999,999,937 を入力すると、必要なすべての素数は生成され、 それ以降の入力では素数の生成をしません。そのときの配列 @prime の内容は、2, 3, 5 ... 31601, 31607, 31627 (要素数 3402) となります。31607 が平方根以下の最大の素数であり、31627 が for ループを終了させるための平方根を越える最初の素数になります。

素因数分解が終了したら、結果を表示します。while ループの最初に %result を空リストで初期化してあるので、その中身を見ることで素数であるか否かが判別できます。%result が空のままであれば、入力された数は素数ということになります。 分解できたときの表示用の書式は、次のコードで生成します。

my $expr = join ' x ',
           map { $result{$_} == 1 ? "$_" : "$_**$result{$_}" }
           sort { $a <=> $b } keys %result;

1行で書いてしまうと少々分かりにくいので、シュウォーツ変換風に書いてみました (単なる好みですが)。シュウォーツ変換は、下から読んでいきます。sort では、分解した素因数を数値順に並べ替えます。map は、sort で並べ替えたリストを受け取って分解数が 2 以上のときに要素を書き換えます。最後に join で連結します。以下は、プログラムの実行の様子です。

Input number(2 - 999,999,999 or expr[+, -, *, **]): 7
7: prime

Input number(2 - 999,999,999 or expr[+, -, *, **]): 77
77 = 7 x 11

Input number(2 - 999,999,999 or expr[+, -, *, **]): 777
777 = 3 x 7 x 37

Input number(2 - 999,999,999 or expr[+, -, *, **]): 7777
7777 = 7 x 11 x 101

Input number(2 - 999,999,999 or expr[+, -, *, **]): 77777
77777 = 7 x 41 x 271

Input number(2 - 999,999,999 or expr[+, -, *, **]): 777777
777777 = 3 x 7**2 x 11 x 13 x 37

引き算

sed について

2番目は、引き算のプログラムです。と言っても、Perl プログラムではなく、Perl の祖先の1つである sed スクリプトです。sed で出来ることはすべて Perl で出来てしまいますので、sed を覚えたり使ったりすることは少なくなっているでしょう。sed は Perl よりも古くからあって、 一時期熱中したことを覚えています。sed は、文字列の編集加工用に作られているので、 数値演算の機能を備えていません。sed で引き算をするには、文字列として処理する必要があります。

sed の動作を、簡単に説明しましょう。sed はフィルタと呼ばれるプログラムで、 ファイル (または標準入力) から1行ずつ読み込んで、編集加工して標準出力に出力します。 sed は、パターンスペースとホールドスペースという2つの作業領域を使います。ファイルからの1行は、 パタースペースに読み込まれます。そして、スクリプトに書かれているコマンドが適用されて、 スクリプトの末尾に達したら出力されます。出力されたら次の行をパタースペースに読み込んで、 同じことを繰り返してというふうに、ファイル全体を処理します。パターンスペースは、Perl の while(<>) {...} の $_ を思い浮かべると理解しやすいでしょう。 もう1つのホールドスペースは、パタースペースの内容を一時的に保存しておくために使用します。 これについては、後ほど触れます。

sed のコマンドは、アルファベットの1文字になっています。 今回のスクリプトで使っているコマンドを、説明しておきます。使っているコマンドは、少しだけで b, h, x, G だけです。

   [address[,address]]b [label]
b コマンドは、address にマッチしたら、指定した label に分岐します。address には、行アドレスあるいは正規表現アドレスを指定することができます。label には、ジャンプ先のラベル名を書きます。ラベル名は、先頭にコロンを書いて :label のようにします。address と label は、共に省略することができます。address を省略した場合には、無条件にラベルにジャンプします。label を省略した場合には、スクリプトの末尾にジャンプします。 そのときには、パターンスペースの内容を出力します。

   h(H), g(G), x
h(H), g(G), x の3つのコマンドは、パターンスペースとホールドスペースのやり取りに関するものです。h は、パターンスペースの内容をホールドスペースにコピーします。 そのときには、ホールドスペースの内容は破棄されます。大文字の H は、ホールドスペースの内容は破棄せずに、間に改行を挿入してホールドスペースの後ろに追加コピーします。g と G は、ホールドスペースからパターンスペースへと逆方向にコピーします。x は、パターンスペースとホールドスペースの内容を交換します。

sed の特徴は、何と言っても正規表現です。今回のスクリプトも、 ほとんどが正規表現のコードからなっています。sed の正規表現は、Perl にも受け継がれています。ここで、sed の正規表現と Perl の正規表現の違いに触れておきます。sed と Perl では、メタ文字の書き方に違いがあります。sed は初期の正規表現に基づいて作られたため、 後から追加したメタ文字には前に \ が付きます。., *, ^, $, [, ] 等の書き方は sed も Perl も同じですが、(, ), +, ?, | などは sed と Perl で書き方が異なります。sed では、\(, \), \+, \?, \| のように書く必要があります。例えば /\([1-9][0-9]\)/ と書いた場合に、sed では捕捉用の丸カッコとなりますが、Perl ではリテラルの丸カッコになります。 もう1つ、s/.../.../ の置換部での捕捉した文字列の参照の仕方に違いがあります。Perl では $1, $2, ... と書きますが、sed では正規表現部と同じように \1, \2, ... と書きます。

それでは、スクリプトを紹介することにしましょう。 スクリプトは、次のようにして実行することとします。

$ echo "4321 - 1234" | sed -f minus.sed

スクリプト

minus.sed

/[1-9][0-9]* *- *[1-9][0-9]*/!{
  s/^.*$/Input error/
  b
}
s/^[^0-9]*\([1-9][0-9]*\)[^0-9]*\([1-9][0-9]*\)[^0-9]*$/\1-\2/

h
s/-.*$//

:n2z
s/z/zzzzzzzzzz/g
s/^9/zzzzzzzzz/
s/^8/zzzzzzzz/
s/^7/zzzzzzz/
s/^6/zzzzzz/
s/^5/zzzzz/
s/^4/zzzz/
s/^3/zzz/
s/^2/zz/
s/^1/z/
s/^0//
/[0-9]/{
  s/^\(z*\)\([0-9].*\)$/\2\1/
  b n2z
}

x
/-/{
  s/^.*-//
  b n2z
}

G
s/\n/-/

s/^\(zz*\)\(z*\)-\1/\2-/

/^-$/{
  s/-/0/
  b
}

s/-$//

:z2n
s/zzzzzzzzzz/y/g
s/y\([0-9]*\)$/y0\1/
s/zzzzzzzzz/9/
s/zzzzzzzz/8/
s/zzzzzzz/7/
s/zzzzzz/6/
s/zzzzz/5/
s/zzzz/4/
s/zzz/3/
s/zz/2/
s/z/1/
s/y/z/g
/z/b z2n

スクリプトの説明

このスクリプトは、主に3つの部分から成ります。

  1. 左右の2つの整数を文字列に変換する。
  2. 正規表現を使って文字列として引き算をする。
  3. 引き算の結果の文字列を整数に変換する。

整数から文字列へまたはその逆の変換

今回のスクリプトでは最初に整数を文字列に変換して、 計算をした結果を最後に整数に変換しています。まとめて、説明することにします。使用する文字は、z にしました。次が、変換をするコードです。

整数 --> 文字列                            文字列 --> 整数
  :n2z                                       :z2n
  s/z/zzzzzzzzzz/g                           s/zzzzzzzzzz/y/g
  s/^9/zzzzzzzzz/                            s/y\([0-9]*\)$/y0\1/
  s/^8/zzzzzzzz/                             s/zzzzzzzzz/9/
  s/^7/zzzzzzz/                              s/zzzzzzzz/8/
  s/^6/zzzzzz/                               s/zzzzzzz/7/
  s/^5/zzzzz/                                s/zzzzzz/6/
  s/^4/zzzz/                                 s/zzzzz/5/
  s/^3/zzz/                                  s/zzzz/4/
  s/^2/zz/                                   s/zzz/3/
  s/^1/z/                                    s/zz/2/
  s/^0//                                     s/z/1/
  /[0-9]/{                                   s/y/z/g
    s/^\(z*\)\([0-9].*\)$/\2\1/              /z/b z2n
    b n2z
  }

整数は、ループを使って上位の位から1つずつ数字を変換していきます。 その際に、すでに変換されている上位の位の z の数を 10 倍にします。 こうすることで、10 の位の z は1の位を変換するときに 10 倍され、100 の位の z は 10 の位を変換するときに 10 倍され、さらに1の位の変換のときに 10 倍されますので、結果的に 100 倍されることになります。

文字列を整数に変換する場合は、下位の位から数字に変換していきます。10 個の z を1個の y にすべて変換すると、残っている z の文字数がその位の数字ということになります。もちろん、z が 0 のときもありますので、0 への変換は最初に行っておきます。後は、大きな数字から変換するように書いておきます。 そして、1つの位を変換し終わったら、y を z に戻してループを続けます。ループを繰り返して、z がなくなったら変換は終了です。

正規表現による引き算

整数を z 文字列に変換したら、次に正規表現を使って引き算をすることになります。z 文字列に変換された式は、zz...zz-zz...zz のように、間をハイフン (マイナス) で繋いであります。z 文字列の引き算は、左右同数の z を削除して、片側 (もしくは両側) の z を 0 にすることによって行います。 引き算をする正規表現を3つほど考えてみました。

1)  :decr
    s/z-z/-/
    /z-z/b decr

2)  s/^\(zz*\)\(z*\)-\1/\2-/   zzzzzzzz-zzzzz

3)  s/\(zz*\)-\1/-/            zzzzzzzz-zzzzz

1) の方法は、ループを使ってハイフンの両側の z を1つずつ削除していきます。2) と 3) は、一度に引き算を行ってしまいますが、マッチの仕方が少し異なります。8 - 5 を例にすると、2) は先頭の z 5個が削除されますが、3) ではハイフン前の5個が削除されます (上の例の下線部)。

上の3つの正規表現は、どれを使っても正常にスクリプトを実行できます。 しかし、計算できる数の範囲には大きな違いがあります。2) と 3) では、5000 以上で異常終了してしまうことがあります。1) では、9999 - 9999 でも正常に実行できますし、19999 - 19999 でも大丈夫です。簡素で単純な正規表現の方がよい場合も、あることの1つの例です。

では、実行速度の面ではどうでしょうか。私の遅いパソコン (Pentium II 800MHz) で、実行してみました。実行した sed は、Vine Linux 3.1 上の v4.1.2 です。2つの整数は、4321 と 1234 を使用しました。

 4321 - 43214321 - 12341234 - 4321
loop s/z-z/-/ 1992
s/^\(zz*\)\(z*\)-\1/\2-/ 1682
s/\(zz*\)-\1/-/ 12550

1) の結果はわかりやすいもので、ループの回数が多いほど、マッチ位置が後ろのものほど時間が 掛かっています。2) の結果は、4321 - 4321 が 16 秒も掛かり、4321 - 1234 よりも多く掛かっているのは意外でした。 4321 - 4321 の場合には、先頭から1度でマッチできるはずです。それに比べて 4321 - 1234 の場合には、\1 に捕捉した 4321 個の z が 1234 個になるまでバックトラックするので時間が掛かるはずです。 正規表現の動作の内部までは分かりませんので、説明ができないのが残念です。3) の結果は、極端なものになりました。4321 - 4321 と 1234 - 4321 の場合は、ほとんど時間が掛かりません。 この場合には、先頭からの1度のマッチで成功します。それに比べて 4321 - 1234 の場合には、255 秒という膨大な時間が掛かります。これは、少し時間が掛かりすぎる感じもしますが、先頭からハイフン前の z が 1234 個になるまでマッチが失敗し、その1つ1つでバックトラックが必要となるためと思われます。

    ---------------------------------------------

今回、初めて sed スクリプトを紹介しました。sed は、Perl ほど難しいことが多くはありません。sed 理解のポイントは、正規表現です。それも Perl ほど高機能でもなく拡張もされていないので、基本を学ぶことができます。機会をみて、sed を再び取り上げることにしましょう。

最後に、上の3つの正規表現を Perl で実行した結果を記しておきます。Perl の正規表現は、sed の正規表現に比較して、高速で高機能です。上の表の正規表現を Perl で実行すると、全部が瞬時に終了してしまいます。そこで、もう1ケタ増やして 54321 と 12345 を使った結果です。Perl のバージョンは、Vine Linux 3.1 上の v5.8.2 です。

$_ = 'z' x 54321 . '-' . 'z' x 54321;

s/z-z/-/ while /z-z/;

if (/^-$/)   { print "0\n"; }
elsif (/^-/) { print '-', length($_) - 1, "\n"; }
else         { print length($_) - 1, "\n"; }

 54321 - 5432154321 - 1234512345 - 54321
s/z-z/-/ while /z-z/ 82254
s/^(z+)(z*)-\1/$2-/ 0260
s/(z+)-\1/-/ 0400

(2005/12/15)

TopPage