fx-5800P 素因数分解 - 高速化

追記修正 2020/05/25
バグフィックス：2021/01/02
追記修正：2021/06/07

[2020/01/24] ソースコード表示を追加し、C.Basicでの実行例の説明を追加しました。
[2020/05/25] 動作可能モデルの記述を更新しました。
[2021/01/02] FACTOR-F1 の終了時に配列変数を解放しない問題を修正しました。
[2021/06/07] fx-7400GIII 用に行列をリストに書き換えて移植したプログラムを公開

今回は、fx-5800P の素因数分解プログラムの高速化の話題です。

素因数分解の一般的な方法は、"エラストテレスのふるい" と言われるものです。最初に入力した数 N を小さい順に素数で割り算を繰り返す方法です。最大の素因数は N の平方根以下の整数なので、平方根以下の整数を対象に２と３以上の奇数で割り算を行う方法で fx-5800P での素因数分解プログラムを作っています。
⇒ fx-5800P 素因数分解 - バグ修正と表示変更

さらに高速化する方法は無いものかと、色々と考えたり、試したりしていました。
　⇒ Casio 関数電卓の素因数分解



▋最速の素因数分解プログラム

最近、目からうろこの高速化を達成した記事を見つけました。
A fast prime factorizing program for Casio fx-5800P 
このトピックで、この作者は、2つのプログラムを公開されています。いずれもプログラム名は FACTOR です。
オリジナルは、上のサイトをご覧ください。

これらには、面白い工夫がなされています。その部分には手を付けず、表示やユーザーインターフェースの不具合(と私が思うだけですが．．．)を解消するために、チョット変更しました。


▍1本目のプログラムの変更版：
結果の表示を 3＾２ と乗数を表示し、List にも乗数の結果が FREQ に反映するように変更しました。
▶ FACTOR-M1 のダウンロード

  
FACTOR-M1 の結果出力 - 左: 乗数表示、- 右: [MODE] [3] で現れる List表示


▍2本目のプログラムの変更版：
オリジナルプログラムでは、結果表示画面を、複数ページで切り替えて自由に見られるようになっています。[EXE] [▲] [▶] [+] で次のページを表示、[(ー)] [▼] [◀] で前のページを表示、[EXIT] [DEL] でプログラムを終了するようになっています。計算が高速化しているだけでなく、結果表示も良くなっています。

このキーの使い方がチョット馴染まないので、[▼] で次のページ、[▲] で前のページ、[+] は使えないようにし、指定以外のキーを押した時の誤動作を抑制するように修正しました。時間をかけて計算できた結果が、誤動作で見えなくなるのは嫌ですから．．．
さらに、オリジナルはキーを軽く押しても応答せず長押しが必要、つまり応答がとても悪いので、キー入力待ちを最小のループにして応答を十分高速にしました。軽くキーを1回押すだけで、チョット待ちますが必ず画面が変わります。
またオリジナルでは結果表示一覧で素因数として 1 が表示されますが、1 は素数ではないので、素因数として 1 が表示されないように修正しました。
▶ FACTOR-F1 のダウンロード [2021/01/02 バグフィックス]
▶ ソースコードはこのページの一番下に掲載

   
FACTOR-F1 の結果表示 - 全結果を画面切り替えで確認、左: 2/1ページ、右: 2/2ページ


▍比較のための以前作った PRIME DECOMP：
▶ PRIME DECOMP のダウンロード

  


上記でダウンロードしたZIPファイルには、それぞれ CCL ファイルと TEXT ファイルが含まれています。CCL ファイルは CcLinker を使って fx-5800P に転送できます。或いは、下のリンクからテキストファイルを参考にしてください。



▋計算時間の比較

先ずは、fx-5800P でどのくらい高速化されたかの結果を示します。

	PRIME DECOMP ソースコード	FACTOR-M1 ソースコード	FACTOR-F1 ソース (メインルーチン) ソース (サブルーチン)
123,456,789 = 32 x 3607 x 3803	170 秒	60 秒 3 倍高速化	42 秒 4 倍高速化
6,666,666,667 = 19 x 1627 x 215659	77 秒	27 秒 3 倍高速化	20 秒 4 倍高速化
7,849,516,203 = 32 x 9811 x 88897	458 秒	165 秒 3 倍高速化	111 秒 4 倍高速化

１本目の FACTOR-M1 で３倍高速化、さらに２本目の FACTOR-F1 は４倍高速化されていることが分かります。

素因数分解は、与えられた数を小さい数から順に割ってゆく時、その操作の回数を減らせば、高速化に繋がります。
PRIME DECOMP では入力した数に、先ず平方根をとって一気に探査範囲を狭め、２と３以上の奇数で小さい方から順に割り算してゆく作戦です。

一方で、理想的なのは、小さい素数から素数だけで順に割り算してゆくことです。それには素数リストが必要ですが、それがあれば苦労しません。素数を算出する計算式などありません。

さて、この作者の工夫は、高い確率で素数を見つけるだけでなく、割り算する候補 (例えば、2 や 3 だけでなく、一旦割り算で使った素数の倍数) を効果的にふるい分ける手法にあります。そして、最初の FACTOR-M1 よりも 次の FACTOR-F1 の方が、素数を見つける確率が高くなっているので、さらに高速化されています。



▋最速プログラムでの工夫 [2019/07/31]

読者のまつ様から、高速プログラムの考え方をご説明頂きました。私はこれに大変納得しましたので、それを掲載致します。以前の記述は撤回させて頂きます。

素因数分解プログラムで通常行われている割り算では，「2,および,3以上の奇数」を割る数としています。
これですと，ご存知のように，例えば3で割り切った後に9でも割るという無駄が出てしまいます。
FACTOR-M1は，「2, 3, 3より大きい3の倍数を除く奇数」を割る数としており，3で割ったあと9や15で割ることはないようにしています。

3より大きい3の倍数を除く奇数について考えてみます。

まず，3より大きい3の倍数でない整数は次の(1),(2)のどちらかです。
　3m+1 (1)
　3m+2 (2)
　(m>=1)

次に(1),(2)が奇数となる式をそれぞれ求めます。

(1)は (2m+1)+m と変形できます。2m+1 は奇数ですから，(1)が奇数であるためには m が偶数である必要があります。そこで m=2n (n>=1) とおけば，(1)は 6n+1 となります。

(2)は 2(m+1)+m と変形できます。2(m+1) は偶数ですから，(2)が奇数であるためには m が奇数である必要があります。そこで m=2n+1 (n>=1) とおけば，(2)は 6n+5 となります。

5とこれらの式を並べ，さらに，隣り合う式の差も書き加えると次のようになります。(n=1とします)

式の並び	隣り合う式の差
5	2
6n + 1	4
6n + 5	2
6(n+1) + 1	4
6(n+1) + 5	2
6(n+2) + 1	4
6(n+2) + 5	2
: :	: :

このような訳で，2, 4, 2, 4, ... と加算していると思われます。
この繰り返しは，6n+1 および 6n+5 の式からも分かるように，2と3の最小公倍数6を周期としています。

FACTOR-F1は，これを拡張して，割る数として
　奇数，かつ，3の倍数でない，かつ，5の倍数でない，かつ，7の倍数でない整数
を順番に求めていく方法をとっていると思われます。
2,3,5,7の最小公倍数は210ですから，隣り合う割る数の差(2,4,6,8など)の並びは210が1周期です。
FACTOR-F1の場合は，13+210n〜13+210(n+1)-1 [n>=0]の範囲で，2,3,5,7のいずれの倍数でもない整数を並べて，隣り合う整数の差をリストにしていると思われます。

ちなみに，13〜222の210個の整数のうち，
　奇数，かつ，3の倍数でない，かつ，5の倍数でない，かつ，7の倍数でない整数
の個数は，FACTOR-F1の「While 1」の下にある
　B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1
のパターンの行数と同じ48個です。

48個となる理由は次の通りです。

まず，2,3,5,7の倍数に関わるそれぞれの個数を求めます。

2の倍数の個数=210/2=105
3の倍数の個数=210/3= 70
5の倍数の個数=210/5= 42
7の倍数の個数=210/7= 30
　
2と3の公倍数の個数= 6の倍数の個数=210/ 6=35
2と5の公倍数の個数=10の倍数の個数=210/10=21
2と7の公倍数の個数=14の倍数の個数=210/14=15
3と5の公倍数の個数=15の倍数の個数=210/15=14
3と7の公倍数の個数=21の倍数の個数=210/21=10
5と7の公倍数の個数=35の倍数の個数=210/35= 6
　
2,3,5の公倍数の個数= 30の倍数の個数=210/ 30=7
2,3,7の公倍数の個数= 42の倍数の個数=210/ 42=5
2,5,7の公倍数の個数= 70の倍数の個数=210/ 70=3
3,5,7の公倍数の個数=105の倍数の個数=210/105=2
　
2,3,5,7の公倍数210の個数=210/210=1

重複を考慮すると，
2,3,5,7のいずれかの倍数の個数
=(105+70+42+30)-(35+21+15+14+10+6)+(7+5+3+2)-1
=162

従って，2,3,5,7のうちどの倍数でもない整数の個数は，
　210-162=48
で 48個となります。


▋素因数分解の正しい動作を検証

この作者は自分のプログラムを検証するテストプログラムまで作って、ソースコード (C++) を公開しています。fx-5800P 用なので、素因数分解する数は最大１０桁と決まっています。この条件下でアルゴリズムの正しさを検証しています。

検証プログラムのソースコードが公開されています。そこで Visual Studio 2019 Community でビルドしました。


▶ FACTOR-F1 のテストプログラム Factor.exe のダウンロード [2019/08/28 リンクを修正]

コマンドプロンプトで、Factor.exe のあるデイレクトリに移動し、そこで

 Factor 1000 50000 4

と入力してエンターキーで実行すると、1000 から 50000まで 4 刻みで入力を変化させて FACTOR-F1 を実行した結果を一気に連続して自動実行してくれます。そして、結果が素数がどうか判定し、素数でないものが出てくるとエラーを出し、素数であれば実行を継続します。正常終了すれば、正しく素因数分解が実行されたことになります。


正常終了しているので、1000 から 50000 までの素因数分解は問題ないことが検証されました。


fx-5800P 用のプログラムなので、入力値は 1 ～ 9,999,999,999 (10桁) の範囲なので、完璧にテストするには、

 Factor 1 9999999999 1

とすれば良いのですが、試しに私のPCだと、１晩で150億回分の検査が済みました。 10桁に相当する (1000億 - 1) 個全部のテストは、夜のみ終夜運転で計算させるとして１週間程度必要になりそうです。


1 ～ 9,999,999,999 までのテスト中

テストが正常終了すれば、このアルゴリズムが正しいことが検証されます。但し、他のパラメータ設定でも正しい組み合わせは有りそうです。



▋グラフ関数電卓用に移植 [2020/05/25 更新]

FACTOR-F1 をグラフ関数電卓用に移植しました。素因数分解の計算部分は変更の余地がありませんが、結果表示は７行全部を使うように変更しました。

▶ グラフ関数電卓用 FactorG のダウンロード [2021/01/02 バグフィックス]
　 バグフィックス：終了時に行列解放しない点を修正しました。

 
FactorG の結果表示

ここで作成したプログラムを CFX-9860G 以降の全てのグラフ関数電卓モデルで動作するように、CFX-9850Gシリーズと fx-7400GIII 用には若干の手直しを行いました。これら以外のモデルでは、上記のダウンロードがそのまま適用できます。 

▶ 実際に動作確認したモデル (若干の手直し移植を含む)
　・CFX-9850G / CFX-9850GC PLUS: 移植のために手直ししたプログラムのダウンロード
　・fx-9860G / fx-9860G Slim / fx-9860GIIシリーズ : 上記ダウンロードをそのまま適用
　・fx-CGシリーズ: 上記ダウンロードをそのまま適用
　・fx-9750GIII / fx-9860GIII:  上記ダウンロードをそのまま適用
　・fx-7400GIII: 移植のために手直ししたプログラムのダウンロード / 転送方法 [2021/06/07 追記]

なお、上の画面は、fx-CG50 にインストールした アドインCasio Basic  - C.Basic for CG で実行した画面です。C.Basic は純正Casio Basic のほぼ上位互換なので、グラフ関数電卓用 FactorG.g3m がそのまま実行できます。C.Basic は純正Casio Basic プログラムが動作し、なおかつ非常に高速に動作するので、上記の例だと一瞬で結果が表示されます。
⇒ アドインCasio Basic - トップページ



▋Factor-F1 (fx-5800P用) のソースコード

FACTOR-F1 – fx-5800P code
0→DimZ↵
22→DimZ↵
"NUMBER"?→F↵
If F<1 Or F≧1x₁₀10↵
Then "NUMBER□MUST□BE□□≧1□ And <1x₁₀10"↵
Stop↵
IfEnd↵
If F≠Int(F)↵
Then "NUMBER□MUST□BE□□AN□INTEGER"↵
Stop↵
IfEnd↵
For 1→E To 22↵
0→Z[E]↵
Next↵
0→Z[1]↵
0→Z[12]↵
0→E↵
F→A↵
↵
Int(√(A))→C↵
2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
3→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
5→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
7→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
11→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
↵
While 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+8→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+8→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+10→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+10→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
WhileEnd↵
↵
Lbl 1↵
If A>1↵
Then Isz E↵
A→Z[E]↵
1→A↵
1→Z[E+11]↵
IfEnd↵
Int(E÷3)→D↵
E-3ｘD>0⇒Isz D↵
1→C↵
↵
Lbl 2↵
Cls↵
Locate 1,1,F↵
Locate 12,1,C↵
Locate 13,1,":"↵
Locate 14,1,D↵
3ｘ(C-1)+1→B↵
Locate 1,2,Z[B]↵
Locate 11,2,"^("↵
Locate 13,2,Z[B+11]↵
Locate 16,2,")"↵
If B+1≦E↵
Then Locate 1,3,Z[B+1]↵
Locate 11,3,"^("↵
Locate 13,3,Z[B+12]↵
Locate 16,3,")"↵
IfEnd↵
If B+2≦E↵
Then Locate 1,4,Z[B+2]↵
Locate 11,4,"^("↵
Locate 13,4,Z[B+13]↵
Locate 16,4,")"↵
IfEnd↵
↵
Do↵
Getkey→K↵
LpWhile K=0↵
↵
0→M↵
If K=34 Or K=73↵
Then 1→M↵
Else If K=86 Or K=85 Or K=47↵
Then 2→M↵
Else If K=83 Or K=84 Or K=57↵
Then 3→M↵
Else Goto 2↵
IfEnd:IfEnd↵
IfEnd↵
↵
If M=1↵
Then Cls↵
"DONE"↵
Else If M=2↵
Then Isz C↵
C>D⇒1→C↵
Goto 2↵
Else If M=3↵
Then C-1→C↵
C<1⇒D→C↵
Goto 2↵
IfEnd:IfEnd↵
IfEnd↵
↵
0→DimZ↵


WFSUB – fx-5800P code
E+1→E↵
B→Z[E]↵
Do↵
D→A↵
Z[E+11]+1→Z[E+11]↵
A÷B→D↵
LpWhile Frac(D)=0↵
Int(√(A))→C↵
Return↵

---

今回は、Universal Casio Forum への投稿プログラムの解説記事になってしまいましたが、私としては十分楽しませてもらったので、記事にして残そうと思いました。



応援クリックをお願いします。励みになるので．．．

にほんブログ村

人気ブログランキング

FC2ブログランキングへ

 



keywords: プログラム電卓、Casio Basic、fx-5800P、素因数分解

---

リンク集　|　ブログ内マップ