パンの木を植えて

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

素数判定から素因数分解へ - Ex7.7(b)(c)

\[ %%% 黒板太字 %%% \newcommand{\A}{\mathbb{A}} %アフィン空間 \newcommand{\C}{\mathbb{C}} %複素数 \newcommand{\F}{\mathbb{F}} %有限体 \newcommand{\N}{\mathbb{N}} %自然数 \newcommand{\Q}{\mathbb{Q}} %有理数 \newcommand{\R}{\mathbb{R}} %実数 \newcommand{\Z}{\mathbb{Z}} %整数 %%% 2項演算 %%% \newcommand{\f}[2]{ \frac{#1}{#2} } \]

本題に入る前に少しおしらせ.

無料で読める本を一冊追加しました.数論と代数の本です.


さて今日もシルヴァーマン『はじめての数論』の演習問題をやっていきます.

今回のテーマも素因数分解です.

\[ %%% 黒板太字 %%% \newcommand{\R}{\mathbb{R}} \newcommand{\C}{\mathbb{C}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\F}{\mathbb{F}} \newcommand{\N}{\mathbb{N}} %%% カリグラフィー %%% \newcommand{\calf}{\mathcal{F} } \newcommand{\calg}{\mathcal{G} } %%% 引数を取るもの %%% \newcommand{\f}[2]{ \frac{#1}{#2} } \newcommand{\im}{\operatorname{im} } \]

問題文(共通)

この練習問題ではあなたは,正の整数 $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$ の間にあるすべての整数について,完全な素因数分解を与えよ.

問bの回答

素因数分解アルゴリズムを考える

難しくなってきました.

まず,前回のEx7.7(a)を解いたときに私が書いたアルゴリズムって,素数判定しかしないアルゴリズムだったんですよね.つまり実行させると

  1. 与えられた数 $n$ は素数である,と返すか

  2. あるいは $n=ab$ を満たす $a,b \geq 2$ を返すか

のどちらかの挙動をするアルゴリズムでした.

これって,因数を返してはいるんですけど素数とは限りませんから,素因数分解ではなかったわけです.

今回求められているのは明確に素因数分解です.素数判定ではないです.だから素数判定を使って,素因数分解アルゴリズムを書くところから始めないといけないですね.


いろいろと弄っているうちに,次のアルゴリズムを思いつきました.

このアルゴリズムは,再帰的な考え方をしています.$N$ の約数 $d$ が見つかったなら,あとは $N/d$ の素因数分解をすればいいという考え方です.再帰文ではありませんけど.

この手続きが有限時間内に停止することは明らかだと思います.$d$ が無限に大きくなるわけはないので.

このアルゴリズムがどのように動作するかは,次の $n=135$ のときの例を見ていただければわかりやすいかと.

つまり,常にwhileループ終了時に

$$ n = N \prod_{p \in A} p $$

が成り立っています.

最終的に $N=1$ となるので,停止時には $n$ の素因数分解が得られるという寸法です.

$A$ に追加されるのが常に素数であり,合成数が追加されないことを疑問に思うかもしれませんが,それは小さいものから順に $d$ で割っているからです.もしも合成数 $d$ が $A$ に追加されたとすると,それ以前に $d$ の約数で割ったときに $N$ が割れていなかったことになって矛盾します.


後は実装していきましょう.

いつものように python で書いていきます.

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

試しにモンスター群の位数

$$ n = 808017424794512875886459904961710757005754368000000000 $$

を代入してみましたが, $0.2$ 秒でちゃんと分解してくれました.

ここで [] というのは空な配列のことです..append で配列の最後尾に要素を追加しています.

ただし,このままだと出力は

[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 7, 7, 7, 7, 7, 7, 11, 11, 13, 13, 13, 17, 19, 23, 29, 31, 41, 47, 59, 71]

というように,重複した素因数がそのままダブって表示されてしまいます.読みにくいですね.


素因数の指数まで出力したい

というわけで改良していきましょう.

問題文を見返してみると $2 \times r$ 行列として保持することが推奨されているっぽいですね.従おうかな?

でもPython で行列を扱うのって結構面倒なんですね.調べてみると NumPy というのがあるそうですが,今回は別に行列である必要はなくて,素因子のキーに対してその指数を格納するようなデータがあればいいだけだと思うので,躊躇われますね.

行列の積とか計算したいわけじゃないしなあ.

そういうわけで,何回割り切れたかを数える変数を新たに用意することにしました.

こんな感じです.

# 素因数分解
## 入力
n = 
## 手続き
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
print(A)

コードの中にある count という変数が「何回おなじ $d$ で割り切れたか」を記憶しています.$d$ で1回以上割り切れた場合には配列に count も格納され,$d$ が変更されたときには直後に count も初期化されます.

これに先ほどのように,モンスター群の位数を代入してみましょう.出力はこうなります.

[2, 46, 3, 20, 5, 9, 7, 6, 11, 2, 13, 3, 17, 1, 19, 1, 23, 1, 29, 1, 31, 1, 41, 1, 47, 1, 59, 1, 71]

いい感じです.

でも最後の素因数 $71$ の指数である $1$ がきちんと出力されていませんね.これは $N=1$ になると while ループを抜けてしまって,count を配列に加えるという操作ができないからですね.

そこで,試し割りする数 $d$ が $N$ より大きくなってしまったら,そのときの count の値を配列に滑り込みで追加するように仕様変更します.

# 素因数分解
## 入力
n = 
## 手続き
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(A)

これでちゃんと「素因数,指数,素因数,指数,・・・」という順で行儀よく出力が並ぶようになりました.


指数を指数らしく出力する

コンピュータで扱う分にはこれでいいのですが,出力を人間が読むことを考えると,もう少し綺麗に出力して欲しいものです.

そこで指数部分に ^ を連結して指数っぽく出力するように,次のように書き直しました.

# 素因数分解
## 入力
n = 
## 手続き
### 素因数分解
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)
### 結果をみやすく出力する
for i in range(0, len(A),2):
    if A[i+1] >= 2:
        print(str(A[i]) +  "^" + str(A[i+1]),end=' ')
    else:
        print(str(A[i]),end=' ') 

str() というのは整数として格納されている値を文字列に変換するための命令です.

print の後に end= という謎のオプションがついていますが,これは出力の最後に改行ではなく空白を入れなさいという命令です.

あとfor文のところに range(0,len(A),2) とあるのは,$0$ から配列 $A$ の長さまで,2ずつ動かしなさいという命令です.

これをモンスター群の位数に対して実行すると

2^46 3^20 5^9 7^6 11^2 13^3 17 19 23 29 31 41 47 59 71

という出力が返ってきます.

だいぶ体裁が整ってきました.

なんとなくやっていただけなんですが,本を見たところこれが問c そのものだったようです.

(c) $n$ の約数をきれいなフォーマットで出力するようなサブルーチンをつくれ.指数部分は指数らしく出力できると望ましいが,それが無理なら,素因数分解の出力は,たとえば $n = 75400 = 2^ 2 \cdot 5 \cdot 7^ 3 \cdot 11$ なら 2^2 * 5 * 7^3 * 11 のように出力せよ.(出力を読みやすくするために,指数部が1であれば出力しないようにせよ)


本題 100万からの素因数分解

さていよいよ本題に戻ってきました.

あなたのプログラムを使って,$100$万$0000$ から $100$万$0030$ の間にあるすべての整数について,完全な素因数分解を与えよ.

この問題をやっていきましょう.

素数をリストアップすることでスピードアップするという件は,今回はもうだいぶ長くなってるので次回に回すことにします.

出力結果はこんな感じでした.

Output exceeds the size limit. Open the full output data in a text editor
1000000 = 2^6 5^6 
1000001 = 101 9901 
1000002 = 2 3 166667 
1000003 = 1000003 
1000004 = 2^2 53^2 89 
1000005 = 3 5 163 409 
1000006 = 2 7 71429 
1000007 = 29 34483 
1000008 = 2^3 3^2 17 19 43 
1000009 = 293 3413 
1000010 = 2 5 11 9091 
1000011 = 3 333337 
1000012 = 2^2 13 19231 
1000013 = 7 373 383 
1000014 = 2 3 166669 
1000015 = 5 200003 
1000016 = 2^4 62501 
1000017 = 3^2 23 4831 
1000018 = 2 500009 
1000019 = 47 21277 
1000020 = 2^2 3 5 7 2381 
1000021 = 11 90911 
1000022 = 2 107 4673 
1000023 = 3 333341 
1000024 = 2^3 125003 
...
1000027 = 7 19 73 103 
1000028 = 2^2 250007 
1000029 = 3 31 10753 
1000030 = 2 5 100003 

サイズが大きすぎて全部は表示できないと.そうですか.意外と器が小さいのね.

あと $100$万$0003$ って素数だったんだね.

使ったコードはこんな感じ.

# 素因数分解
## 入力

## 手続き
for n in range(1000000,1000031):
    ### 素因数分解
    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(str(n) + " =", end = ' ')
    for i in range(0, len(A),2):
        if A[i+1] >= 2:
            print(str(A[i]) +  "^" + str(A[i+1]),end=' ')
        else:
            print(str(A[i]),end=' ') 
    print()  

最後に書かれてる空っぽの print は改行を出力するために置いてあります.


感想とまとめ

今回は素因数分解のコードをpythonで書きました.

まだ問題の途中ですが,記事が長くなってきたのでこの辺で切ります.