e-Gadget - プログラム関数電卓　(旧) Casio 関数電卓の素因数分解

(旧) Casio 関数電卓の素因数分解

追記 2021/5/28

本ブログでは、素因数分解プログラムについて取り上げています。

　・ fx-5800P で素因数分解
　・ fx-5800P で素因数分解再び
　・ fx-5800P 素因数分解 - バグ修正と表示変更
　・ fx-9860GII への移植 - 素因数分解
　・ VBAで素因数分解
　・ fx-5800P 素因数分解 - 高速化 [2020/01/07 追記]
　・ fx-5800P 素因数分解 - さらなる高速化 [2020/12/25 追記]
　・高速素因数分解を Casio Basic から Casio Python へ移植 [2021/05/24 追記]
　・fx-CG50 / fx-9750GIII - C.Basic で素因数分解15桁対応＆高速化 [2021/05/28 追記]

これらの記事の発端は、カシオのスタンダード関数電卓 fx-995ES に因数分解機能が搭載されたことです。素因数分解機能が内蔵されたのなら、プログラム電卓で作ってみようと思ったわけです。但し、自作プログラムは、逆立ちしたって内蔵機能よりも速く計算できません。これは、チョット気になっていました。

そこで、今回は、素因数分解のアルゴリズムの話題です。

実は、すけっぴぃ様から、素因数分解アルゴリズムの効率化に関するコメント（ここ）を頂いており、すぐには役に立つ情報を返信できなかったのですが、心の隅に引っかかっていたこともあって、この話題を少し掘り下げてみます。

[2020/12/25 追記]
すけっぴぃ様からのコメントを実際に実現できたのが、上に挙げた 2020/01/07 追記の記事に相当します。このアルゴリズムを利用して素因数検索効率をさらに向上させ、fx-5800P のメモリ容量ギリギリまで使うことでさらに高速化したのが、上で挙げた 2020/12/25 追記の記事になります。

以下は、fx-5800P 用高速素因数分解のアルゴリズム以前にアレコレと考えていた時の話です。

私の手持ちの電卓で素因数分解ができるのは、最初から機能が内蔵されている fx-995ES と fx-JP900、それに作ったプログラムが走る fx-5800P と fx-9860GII の４機種。

順に、fx-995ES、fx-JP900、fx-5800P、fx-9860GII

最近、fx-JP900 を入手し、素因数分解機能を試したところ、fx-995ES よりも大幅に計算が速くなっていて、さらに得られる素因数の桁数が増えています。

幾つかの素因数分解結果を、fx-995ES、fx-JP900、fx-5800P のプログラム、カシオの高精度計算サイトKe!san で行った結果の一覧表を再掲載します。← Casio fx-JP900 （その３）から抜粋

整数	fx-955ES	fx-JP900	fx-5800P	Keisan
1,234,567	127x(9721)	127x9,721	127x9721	127x9721
98,765,432	2³x37x(333667) (1.7秒)	2³x37x333,667 (0.4秒)	2³x37x333667 (27秒)	2³x37x333667
9,516,208,473	3²x17²x(3,658,673)	3²x17²x(3,658,673)	3²x17²x3658673	3²x17²x3658673
123,456,789	3²x(13717421)	3²x(13,717,421)	3²x3607x3803	3²x3607x3803

因数分解できない場合は ( ) 付きで表示される仕様で、そこは fx-955ES も fx-JP900 も同じ。
但し、fx-995ES は素因数が４桁以上で ( ) 付きになっていましたが、fx-JP900 では素因数の桁数制限が大幅に緩和され、上の例ですと６桁まではOK。

また、fx-JP900 の取扱説明書１７ページを見ると、素因数が 1,018,081 以上の素因数を持つ時は計算エラーになると書かれていますが、上の１つめ、２つめ、３つめの例では、( ) 付きで結果が表示されています。これ以上計算できないが、たまたまそれが素因数だったと解釈すれば良さそうです。( ) 付きの場合は、それが正しいかどうかの保証が無いと言うことでしょう。

４つめの例では、( ) 付きの結果は、さらに因数分解できるが、これ以上計算できないことを示しています。

[2015/06/18 追記]
コメント欄での sentaro様とのやりとりから、求める素因数の桁数を制限すること、そしてCPUの演算速度の簡単な比較から、２と３以上の奇数で順次割り算する単純なロジックても、それぞれの実行時間をほぼ説明できそうだ、と今のところの結論です。(追記終わり)

[2020/01/24 追記]
上記のすけっぴぃさんのご提案にあるように、一旦3で割ればその後は3の倍数でない奇数で割り算すると効率があがります。これを拡張して、割る数として奇数、かつ 3でない、かつ 5でない、かつ 7でない整数を順次求めて行くアルゴリズムで3～4倍の高速化ができました。⇒ fx-5800P 素因数分解 - 高速化

素数は整数論の一分野で研究されていて、素数が無限個存在することは、紀元前３００年頃に証明されています。しかし、未だに素数を求める公式が見つかっていないし、素数の密度を求める式も見つかっていません。ケースバーケースの近似的な式があるだけです。それだけ神秘的な数が素数と言えるのですが、だからこそ素数は暗号通信の要として、実用上極めて重要なものになっています。応用工学で重要な素数だからこそ、数学教育において素数は重要だと考え、カシオは素因数分解機能を内蔵したのかも知れません。

ところで、遅いCPUを搭載した fx-5800P で走らせた素因数分解プログラムは、処理が遅いので、あまり凝ったアルゴリズムを実装できません。それに比べて内蔵機能の計算が桁違いに速いことが、気になっていました。

そんなとき、fx-JP900 の評価をしていて、正しく計算する条件として求める素数の桁数に制限があることを改めて考えてみました。桁数の制限をかけるアルゴリズムがあって、それを使えば Casio Basic でも速いアルゴリズムを実装できるかもしれないと思いました。そこで、具体的なプログラムにするには、まだ不完全ですが、取りあえず書いてみます。

fx-995ES では３桁以内の素因数に分解して表示することができるのですが、正しい計算が保証できるときの素因数が３桁と言うのは、どういうことか？　例えば、何十個かの素数の表をメモリ内に持っていて、それを使って計算すれば、かなり計算量が節約できるかも知れません。３桁つまり1～999の範囲にある素数は168個で、最大の素数は 997 です。168個程度のデータならメモリに入れておくのは現実的です。与えられた自然数をこの168個の素数で次々と割り算してゆけば、２と奇数で順次割り算するよりも、遙かに計算量が少なそうです。

ある整数の素因数分解を行う時、素因数はその整数の平方根以下になることを使い、さらに２と奇数を順に割り算して素因数を求める計算を、fx-5800P や fx-9860GII 用に書いたプログラムで行っています。このとき、３桁の素数で最大は 997 なので、997² = 994009 以下の素因数分解は、メモリに保存された168個の素数で順割り算すれば、かなり計算量は節約できます。fx-5800P や fx-9860GIi で作ったプログラムでは、２と奇数で割り算するので、499個の数で割り算しています。つまり、計算量は 168÷499 = 0.3367 つまり約33%、1/3 にに抑えられます。
Casio Basic での変数参照や配列参照のアクセスに比べて、関数計算する際のテーブル参照は遙かに速いと考えられるので、これで内蔵機能が速いことが説明できそうです。この方法を Casio Basic のプログラムに反映する際は、算術演算に 1～1.5ms 程度かかり、配列変数や行列へのアクセスに20ms程度かかること、そしてこの差を考慮して、本当に速くなるのかどうかを検討する必要があります。

一方 fx-JP900 は、1,018,081 以上の素因数を持っている整数ではエラーになる仕様です。ちなみに 1,018,081 は素数でなくて、これを素因数分解すると 1009² です。fx-5800P のプログラムで計算してみると、1,018,081 未満で最大の素数は、1,018,057 だと分かります。1,000,000 (100万)までの素数は 78,498個あるので、8万個近くある素数の表を不揮発メモリに持っておくのは、関数電卓としてはあまり現実的ではないでしょう。1,018,081 以上でエラーになることから、ひょっとしてこれを実装していることも考えられますが、素因数分解だけのために原価を押し上げるメモリを使うのは疑わしいわけです。

そこで、メモリに保存された素数テーブルを使わず、完全ではないものの、素数を計算して、それで順次割り算してゆく方法は、うまくすると計算量を減らして、高速化できるのかも知れません。或いは、メモリ上の素数テーブル参照と計算の併用も現実的な折衷案かも知れません。

素数の計算については、以下の式が非常に効率よく素数を計算するものとして知られています。

これはオイラーが見つけた、素数 p を求める式です。全ての素数をもれなく見つけることはできません。
ただ、この式の面白いのは、x が 0 から 39 の時に得られる数が全て素数だと言う点にあります。但し、x=39 の時得られる素数 1601 以下には、この式で得られない素数が沢山あります。1601 以下の全ての素数を算出できるわけではありません。

さらに、x が 0 から 60 の時は、x = 40, 41, 42, 44, 49, 50, 57 の７個の x 以外で、 p は素数になります。素数を算出する効率は、x が 0 から 60 で計算した61個の値のうち88.5% が素数になります。繰り返しますが、これで得られる3071 以下には、この式の結果以外に多くの素数があります。

x	p	x	p	x	p
0	41	13	223	26	743
1	43	14	251	27	797
2	47	15	281	28	853
3	53	16	313	29	911
4	61	17	347	30	971
5	71	18	383	31	1033
6	83	19	421	32	1097
7	97	20	461	33	1163
8	113	21	503	34	1231
9	131	22	547	35	1301
10	151	23	593	36	1373
11	173	24	641	37	1447
12	197	25	691	38	1523
				39	1601
				･･･	･･･
				60	3701
				1008	1017113

ちなみに、x = 1009 の時 p = 1,019,131、x = 1008 の時 p = 1,017,113 となります。つまり、x = 1008 で得られる1,017,113 がこの式でエラーにならない最大素因数 1,018,057 に最も近いものだと分かります。この式では、100万までの素数の47.5% が求められることが調べられていて、効率は半分以下に落ちます。しかし、現在 fx-5800P や fx-9860GII に実装しているプログラムのロジック「２と奇数で約 50万回近く割り算する」よりも、1000回程度の割り算をする方が、計算量が 1/500程度になります。仮に３分かかっていた計算は１秒以下になりそうです。

与えられた整数に対して、先ず先に、このオイラーの式で得られる計算値から素数を選んで、それで順次割り算して、残った余りについて、２と奇数で順次割り算して計算すると言うアルゴリズムが考えられます。オイラーの計算値から少ない計算量で素数を見つけられれば、テーブルと計算の併用で、高速化の可能性があります。x をどの範囲まで使うのか、38 までとするか、60 とするか、1008 までとするか、その中間のどこかにするか、も実際の計算量とコマンドの処理速度から最適点を検討する必要があります。仮に、x = 1000 まで使うなら、1000の47.5% に相当する 475 個の素数で順次割り算するので、500 / 500,000,000 = 0.001% となり、他の計算量と計算時間を低く抑えられれば、大幅な高速化ができるかも知れません。

プログラムの実装は、そのうちやってみようと思います。もし先にプログラムを作られた方は、是非コメント欄で発表してください。お待ちしております。

ところで、余計な話になりますが、

The Asahi Shinbun Blobe - 数学という力 (残念ながら削除されたようです) というサイトでは、素数と円周率の関係、素数と宇宙の真理の関係について解説されていて、なかなか面白いです。素数は現代数学の花形の１なのですね。

これは、リーマンと言う数学者が名付けたゼータ関数というものです。一番右の項は、オイラー積という総積計算で、素数 p がしっかり出てきます。

この式で、s=2 の時は、

となって、素数 p のオイラー積と円周率 π が結びつくことが発見された歴史的な式をゼータ関数で表したものです。素数と円周率の関係式をより一般化したゼータ関数で表現することで、より深く調べることができるというわけです。上の記事によれば、この式を巡った素数の熱い研究が進められているようですね。簡単に素因数分解できる公式が存在するのかどうか？まだよく分かっていません。そのうち見つかる筈と言う人もいれば、存在しないかも知れないと言う人もいます。あまり簡単に計算できると、インターネットなどで使われている暗号が簡単に破られることに繋がるので、大問題です。

脱線すると、このゼータ関数で s=-1 の時は、

なんて、なってしまいます。自然数を無限に足すと、-1/12 になる、なんとも受け入れられない結果です。実は、ゼータ関数は複素数の世界のものなので、一見あり得ない結果に見えるわけです。

カシオの高精度計算サイト Ke!san のココでも取り上げられています。
ここでは、明確に書かれていませんが、s=-1 で、実数の世界を複素数の世界に拡張（解析接続）しています。実はゼータ関数といっても、色々な派生バージョンがあって、s>1 でしか成り立たないゼータ関数では s=-1 は定義域の外側、言い換えれば s=-1 としてはいけません。そこで、このゼータ関数を拡張して(解析接続して) s=1 以外で成り立つようにできます。この時のζ関数は、拡張前のゼータ関数とは別のものです。そして拡張した別のζ関数では、s=-1 としても良く、-1/12 に収束します。殆ど詐欺というか、紛らわしい話ですね。

1+2+3+4+･･･は実数の世界と誰でも思うので、自然数の無限和が負の数である -1/12 に収束するなどと言うのは、実数の世界では間違いです。でも複素数の世界では収束するので、ゼータ関数は物理の計算で重宝されているわけです。

素数は大きくなると、まばらになるのか？その密度はどうなっているのか？　これもはっきりと分からない問題でしたが、2014年には素数は極端な偏りがなく万遍なく分布することが発見されました。→ こちら

素数の性質の研究は、非常にホットな分野ですね。まぁそれだけ分からないからこそ、暗号に使われるわけです。
深淵な研究は数学者に任せるとして、取りあえずプログラム電卓で素因数分解の高速化が出来るかも知れないということで、一旦区切ることにします。

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