素因数分解アルゴリズムを賢くする - Ex7.7(b)#2

今回も『はじめての数論』の演習問題をやっていきます．

作者:Joseph H. Silverman
丸善出版

演習7.7(b) です．(c)は前回の記事でついでにやってしまったので，残るは(b)のみです．

問題文(共通)

この練習問題ではあなたは，正の整数 $n$ を素数の積に分解するプログラムを書くことを求められている．（もし $n=0$ であれば，無限ループに入らずにエラーを返すことを忘れずに！）$n$ の素因数分解を表現する便利な方法は $2 \times r$ 行列として表すことである．つまり，

$$ n = p_1^{k_1} p_2^{k_2} \cdots p_r^{k_r} $$

のとき，$n$ の分解を行列

$$ \begin{pmatrix} p_1 & p_2 & \cdots & p_r \\ k_1 & k_2 & \cdots & k_r \end{pmatrix} $$

で表す．（もしあなたのコンピュータが動的な記憶割付を許さないときは，実行前にどれだけ多くの因数を許すか決めておく必要がある．）

問bの問題文

最初の100個（あるいはそれ以上）の素数を記憶させておいて，大きな素因数の探索以前に，まずこれらの素数をすべて $n$ から取り除くようプログラムを変更せよ．考え得る大きな因数 $d$ を試す際，$d$ が偶数，そして $3,5$ の倍数であるかのチェックをいとわなければ，スピードアップできる．効率を上げるのに，$m$ が素数であれば $2$ と $\sqrt{m}$ の間にある数では割り切れないことを利用することもできる．あなたのプログラムを使って，$100$万$0000$ から $100$万$0030$ の間にあるすべての整数について，完全な素因数分解を与えよ．

回答

100万からの連続する31個の整数の素因数分解は前回で終わっています．

残っているのは，アルゴリズムの計算時間の改善ですね．

前回は，いろいろとコードを加えてごちゃごちゃしたコードを得ましたが，根幹部分はこの疑似コードでした．

このアルゴリズムの計算時間は，$d$ が $N$ と一致するまで $d \leftarrow d+1$ し続けるので $O(n)$ です．

実際に whileループの繰り返し回数を確かめてみましょう．

N = n
d = 2
A = []
count = 0
while N >= 2:
    if N % d == 0:
        if count == 0:
            A.append(d)
        N = N // d
        count = count + 1
    else:
        if count >= 1:
            A.append(count)
        d = d+1
        count = 0
    if d > N:
        A.append(count)
print(d)

このように whileループが終わった後に print(d) と指示すれば，whileループを少なくとも何回繰り返したかが分かります．

やってみると，確かに $n$ 回繰り返してしまっていることがわかります．

前回の最終形のアルゴリズムに手を加えて修正しようかと思いましたが，どうにもうまく行きません．そこで，疑似コードからやり直すことにします．

f:id:seasawher:20220318170221j:plain

whileループの繰り返し条件ですが，最初は $N \geq 2$ だったんですけどそのままズバリ $d^ 2 \leq N$ に変更しましょう．

コードとしてはこんな感じ．

# 素因数分解
## 入力
n = 
## 手続き
N = n
d = 2
A = []
count = 0
while d * d <= N:
    r = N % d
    if r==0:
        A.append(d)
        N = N // d
    else:
        d = d+1
A.append(N)
print(A)

これで繰り返し回数は $O(\sqrt{n})$ になります．やったね．

後は，このアルゴリズムを前回の記事の要領で改善していくだけです．

結果，完成したものがこちら．

# 素因数分解 - O(√n)
## 入力
n = 
## 手続き
N = n
d = 2
A = {}
while d*d <= N :
    r = N % d
    if r==0:
        if d in A:
            A[d] = A[d] + 1
        else:
            A[d] = 1
        N = N // d
    else:
        d = d+1
#### AにNを追加
if N in A:
    A[N] = A[N]+1
else:
    A[N] = 1
## 結果を見栄えよく出力する
print(str(n) + " =", end = ' ')
primes = list(A)
for p in primes:
    if A[p] >= 2:
        print(str(p) +  "^" + str(A[p]),end=' ')
    else:
        print(p,end=' ')

今回のコードのポイントは，配列ではなく辞書を使っていることです．A={} と書いてありますが，これは空な辞書を宣言しています．

前回のコードでは素因数とその指数を配列で管理していて，偶数番目に素因数を，奇数番目に指数を格納していました．

でもそれがバグが発生しまくる原因になっていました．異なるものを同じリストにごちゃまぜにするのがいけないのです．

今回は辞書を使うことで，素因数をキーにしてその指数を呼び出すことに成功しました．かなりいじりやすいコードになったと思います．

if d in A: と書いてあるのは，辞書 A に d をキーとする値が格納されているか確認して分岐しています．直接 A[d] = A[d] + 1 と書くと，d をキーとする値がなかったときにエラーになってしまうようです．

いつものようにモンスター群の位数も分解してみました．OK．

感想とまとめ

$O(\sqrt{n})$ 時間の素因数分解アルゴリズムを実装しました．

素数をリストにあらかじめ蓄えておくとか，エラストテネスの篩みたいにサイズ $n$ の表を作って，そこから合成数を順に割って消していくとかすれば $O(\sqrt{n}/ \log n)$ まで改善されると思いますが，メモリ使用量が多くなっちゃうからな……．機会があればやってみたいですけど．

最近数論ブログでさえなくなって，プログラミングを勉強するブログのようになってしまっていて申し訳ないです．でもメインは数学ですよ．

次回以降は，pythonではなくてjuliaを使ってみようかと思っています．理由は，なんとなくイマドキだからです．新しいものには触れてみたいじゃないですか？

juliaはまだ比較的新しい言語なのでライブラリが充実していないそうですが，学習のためなんだから関係ありません．

パンの木を植えて

主として数学の話をするブログ

素因数分解アルゴリズムを賢くする - Ex7.7(b)#2

問題文(共通)

問bの問題文

回答

感想とまとめ