素因数の数を数える

素因数分解は、数論の中でも重要かつ難しい問題のひとつとして知られています。 素数で割っていくナイーブな方法から楕円曲線法や数対ふるい法までいくつかのアルゴリズムが考案されていますが、素因数分解そのものではなく、素因数の個数を求めたいだけならば、簡単かつそれなりに高速なアルゴリズムがあります。

このアルゴリズムは、エラトステネスの篩によく似た形をしています。 次のプログラムは、配列 factorNum の各要素に1から nNum - 1 の素因数の個数を格納します。 例えば、 factorNum[360] には6が格納されます( 360 = 2 * 2 * 2 * 3 * 3 * 5 )。

1:  void prime_factor_num(int *factorNum, int nNum){
2:    for (int i = 1; i < nNum; i++)
3:      factorNum[i] = 1;
4:
5:    for (int i = 2; i * i < nNum; i++)
6:      if (factorNum[i] == 1)
7:        for (int j = 2; i * j < nNum; j++)
8:          factorNum[i * j] = factorNum[i] + factorNum[j];
9:  }

このプログラムについて解説を加えていきましょう。

2行目と3行目では、配列の初期化を行っています。 全ての正の整数は少なくとも1つの素因数を持つ(1の素因数分解は1と定義されている)ので、1を代入しておきます。 0はそもそも素因数分解自体が定義されませんので、このプログラムでは不定の値を返すことにします。

5行目から8行目まで、特に8行目がこのプログラムの本質です。 このアルゴリズムは、「(a * b の素因数の個数) = (a の素因数の個数) + (b の素因数の個数)」という定理に基づいて計算を繰り返します。 また6行目のif文は、エラトステネスの篩の場合はbool型の値を使って素数かどうか判定するのが一般的でしたが、このアルゴリズムでは素数かどうかの判定に「素因数が1個かどうか」という判定を用いています。 i = 2 のときのループを見てみましょう。

添字 0123456789101112...
-111212131213...

5行目から始まるforループをi = 2の時まで実行完了したときの、factorNumの値をトレースしてみました。 まだ計算されていないfactorNum[9]を除いては、正しい答えが得られています。

さて、ここで疑問に思う方がいるかもしれません。
この時点では3の素因数が1個だということはわかっていないのだから、factorNum[3]を使ってもfactorNum[2 * 3]は求められないのではないか、と。

実は、その通りなのです。i = 2のループでは、factorNum[6]の値は偶然正しい値になってはいますが、6の素因数の個数だと言うことはできません。
しかしi = 3のループのとき、factorNum[6]の値はもう一度計算され上書きされます。このときにはfactorNum[2]の値は確定しているので、factorNum[3 * 2] = factorNum[3] + factorNum[2] で上書きされて、ちゃんと意味の通った値が格納されるのです。
例をひとつだけ挙げましたが、同様の論理は他の要素についても成り立ちます。

このアルゴリズムの計算量は、明らかにエラトステネスの篩と同じになります。最適化を行えばエラトステネスの篩はO(n)の計算量になるそうですが、率直に実装してもO(n3/2)で抑えられますね。

なぜ今日はこんなトピックを持ち出したのかといいますと、UVaでこれを使うと高速に解ける問題を見つけたからでした。 884番の Factorial Factors です。 ぜひお試しくださいませ。

担当: 田山 (これから新年会)