Casio Python - 大きな数の計算:高速素因数分解(3)

Python Casio Python
 Casioグラフ関数電卓の Python を使ってみる
     - 大きな数の計算:高速素因数分解(3) 
目次


初版:2020/07/28
追記修正:202011/14

前の記事 - 9. Python らしい反復処理 |  次の記事 - 11. 関数呼出オーバーヘッド


10. 大きな数の計算:高速素因数分解(3)

今回は、Casio Python では、より大きな桁数の整数で精度を確保できることを紹介し、15桁入力に対応した高速素因数分解スクリプトに改造します。

前回までに作った高速素因数分解スクリプトは、オリジナルの Casio Basic プログラムから移植したものなので、入力値の桁数はオリジナルと同じように10桁が上限になっています。

Casio Basic は仕様上10進数10桁を超えると計算精度の保証がないので、入力値は最大10桁に制限されます。一方、Casio Python では10進数10桁を超えても精度が確保されます。


例えば、15桁の入力を2で除算してみます。
999,999,999,999,999,999 / 2 = 49,999,999,999,999.8
15digit 

入力15桁に対して、1桁多い16桁で出力されます。
正解は 499999999999999.5 なので、16桁の最下位に誤差が発生しています。


次に、16桁の入力を2で除算してみます。
9,999,999,999,999,999,999 / 2 = 50,000,000,000,000.0
16digit 

16桁の入力に対して、出力は1桁多い17桁になっています。
奇数を2で除算して整数になることから、少なくとも16桁の入力では、計算精度に問題があることが分かります。

続いて、1/3 の計算結果は、出力16 桁になります。
digit_one_third_shell 

円周率 pi も 16桁出力になります。
digit_pi_shell 

Casio Python での浮動小数点計算は、10進数15桁以内の入力なら、計算結果は15桁の精度が確保できるようです。

そこで、高速素因数分解の入力値の桁数を 10 桁から 15 桁に拡張することにします。


10.1 変更箇所の洗い出し

前回作った高速素因数分解 - 10桁対応版:FactorG2.py のダウンロード

15桁対応にするために、前回作った下記のスクリプトに対して、下記の2項目を変更します。

(1) 増加する素因数に合わせて表示行数を拡張する - disp() の変更
(2) 入力桁数の制限を拡張し、それに合わせてエラー表示を変更する


これら2つの項目に対応して変更箇所を以下に赤文字で示します。

fx-CG50 Pythonモード:高速素因数分解 - FactorG2.py
"""Sample script

 Exercise;
 ported from Casio Basic
 "Prime Factor 10 digits"
  FactorG.py
   ver 2.0

 by Krtyski/e-Gadget
"""

from u import *

def ckpwr():
 global a,b,c,d,e,z
 e+=1
 z[e]=b
 while 1:
  a=int(d)
  z[e+11]+=1
  d=a/b
  if frac(d):
   break
 c=int(sqrt(a))

def
disp():
 global a,e,f,z
 clear_screen()
 locate(00f3'm'0)
 locate(160': ' str(len(inp)) + ' digits'3'm'0)
 line(0,15,383,15,1,0)
 for i in range(112):  ◀ (1) 増加する素因数に対応して、表示行数を拡張する
  if i<=e:
   locate(,,z[i], 1'm'0)
   locate(10i'^('2'm'0)
   locate(12iz[i+21], 1'm'0)
   locate(15i')'2'm'0)
 locate
(00''0'm'1)


while
1:
 try:
  inp=str(eval(input('Number:')))
 except (SyntaxError,TypeError,NameError) as e:
  print(e)
  print('*must be number or\n expression')
  continue
 if '.' in inp:
  print('*must be integer')
  continue
 elif inp.isdigit():
  if len(inp)>10:    ◀ (2) 桁数制限を拡張する
   print('*must be 10 or less') ◀ (2) 拡張した桁数制限に合わせて変更
   continue
  else:
   f=int(inp)
   break
 else:
  continue

z=list(range(23))
for e in range(1,23):
 z[e]=0
e=0
a=f
c=int(sqrt(a))

prime_list=[2,3,5,7,11]
for b in prime_list:
 d=a/b
 if frac(d)==0ck_pwr()
 if b>c: break


increment=\
[2,4,2,4,6,2,6,4,
 2,4,6,6,2,6,4,2,
 6,4,6,8,4,2,4,2,
 4,8,6,4,6,2,4,6,
 2,6,6,4,2,4,6,2,
 6,4,2,4,2,10,2,10]

while 1:
 for i in increment:
  b+=id=a/b
  if frac(d)==0ck_pwr()
  if b>cbreak
 if b>c
break

if a>1:
 e+=1
 z[e]=a
 a=1
 z[e+11]=
1

disp
()



10.2 増加する素因数に合わせて表示行数を拡張する - disp() の変更

素数の積を調べてみます。

2から41までの素数の積;
2*3*5*7*11*13*17*19*23*29*31*37*41 = 304,250,263,527,210
15桁になります。

2から43までの素数の積;
2*3*5*7*11*17*19*23*29*31*37*41*43 = 21,306,211,311,576,906
17桁になります。

これらから、素数13個で15桁になることが確認でき、乗数を含めて出力には最大 13行必要になります。

前回作った disp() を示します。

def disp():
 global a,e,f,z
 clear_screen()
 locate(00f3'm'0)
 locate(160': ' str(len(inp)) + ' digits'3'm'0)
 line(0,15,383,15,1,0)
 for i in range(112):  ◀ (1) 増加する素因数に対応して、表示行数を拡張する
  if i<=e:
   locate(,,z[i], 1'm'0)
   locate(10i'^('2'm'0)
   locate(12iz[i+21], 1'm'0)
   locate(15i')'2'm'0)
 locate
(00''0'm'1)


Ver2_10lines1 

1画面あたり、素因数は11個まで表示できるので、for i in range(12): としました。
これを13行まで拡張するためには、以下のように2列にして13行まで表示できるように変更します。

13lines 

従って、反復処理の for 文は、for i in range(14): で良いですが、1行余裕をみて

for i in range(15):

としておきます。しかし、これだけでは、12行目 と 13行目 は、右側に表示されないので、4行ある locate() の引数を適切に変更してゆきます。
 
i の値が 12 未満では左側に表示、12 以上の場合は右側に表示すれば良いですね。
そこで、if i <12: / else: と場合分けをする方法もあるでしょうが、もう少しスッキリと簡潔に記述するために、locate() で指定する 桁数と行数の増分 dxdyi の値に応じて計算で求める方式を考えました。
 dx=int(i/12)*16  i<12 のとき dx=0i≧12 のとき dx=16
 dy=int(i/12)*11
  i<12 のとき dy=0i≧12 のとき dy=11

これを使って、スクリプトを以下のように変更しました。

def disp():
 global a,e,f,z
 clear_screen()
 locate(00f3'm'0)
 locate(160': ' str(len(inp)) + ' digits'3'm'0)
 line(0153831510)
 for i in range(115):  ◀ 変更部分
  if i<=e:
   dx=int(i/12)*16  ◀ 追加部分: i < 12 のとき dx=0、i ≧ 12 のとき dx=16
   dy=int(i/12)*11    ◀ 追加部分: i < 12 のとき dy=0、i ≧ 12 のとき dy=11
   locate(0+dxi-dyz[i], 1'm'0)    ◀ 第1引数に +dx を、第2引数に -dy を追加 
   locate(10+dxi-dy'^('2'm'0)    ◀ 第1引数に +dx を、第2引数に -dy を追加
   locate(12+dxi-dyz[i+21], 1'm'0) ◀ 第1引数に +dx を、第2引数に -dy を追加
   locate(15+dxi-dy')'2'm'0)    ◀ 第1引数に +dx を、第2引数に -dy を追加 
 locate
(00''0'm'1)


他にもスッキリと記述する方法があると思いますが、これで進めます。


10.3 入力桁数の制限を拡張し、それに合わせてエラー表示を変更する

前回までに作った "入力と入力チェックの記述" を以下に示し、今回変更する部分を赤文字で示します。

while 1:
 try:
  inp str(eval(input('Number:')))
 except (SyntaxErrorTypeErrorNameError) as e:
  print(e)
  print('*must be number or\n expression')
  continue
 if '.' in inp:
  print('*must be integer')
  continue
 elif inp.isdigit():
  if len(inp)>10:    ◀ (2) 桁数制限を拡張する
   print('*must be 10 or less') ◀ (2) 拡張した桁数制限に合わせて変更
   continue
  else:
   f int(inp)
   break
 else:
  continue

ここで、赤文字で追記している2箇所は、10 を 15 に変更すればOKです。

今回は、変数 digit を使って、これに 15 を代入して使うことにします。

digit=15

をスクリプトの冒頭に書いておけば、代入する数値を書き換えるだけで入力桁数の制限を簡単に変更できます。

そこで、上記2箇所は、以下のように変更します。

if ken(inp) > digit:
 print('*must be '+str(digit)+' digit\n or less')

 continue

出力文字列に エスケープシーケンス \n (バックスラッシn) による改行を入れているので、狭い画面の電卓で以下のように見やすくなります。
*must be 15
 or less



入力と入力チェック部分は、まとめると以下のようになります。

digit=15

・・・
・・・

while
1:
 try:
  inp str(eval(input('Number:')))
 except (SyntaxErrorTypeErrorNameError) as e:
  print(e)
  print('*must be number or\n expression')
  continue
 if '.' in inp:
  print('*must be integer')
  continue
 elif inp.isdigit():
  if len(inp)>digit:
   print('*must be '+str(digit)+' digit\n or less')
   continue
  else:
   f int(inp)
   break
 else:
  continue


10.4 素因数探索回数の表示を追加する

ここまで作ったスクリプトを実際に以下の2通りの入力値で走らせてみると、実行時間が大きく異なります。

NoSearch_ver3_fast_14digits NoSearch_ver3_slow_13digits 

左の14桁に比べて、右の13桁は、桁数が1つ少ないにもかかわらず、かなり時間がかかります。
素因数の桁数が大きいほど時間がかかると考えれば、素因数探索の while ループを回る回数が大きいほど時間がかかると理解できます。

そこで、実際に素因数探索の回数をカウントして、以下のように表示しようと思います。

Search_ver3_fast_14digits Search_ver3_slow_13digits 

search は探索回数です。

変数 search の初期化
スクリプト冒頭の digit=15 の直後に search=0 を追加し、変数 search を探索回数とします。

from u import *
digit=
15
search=0


 変数 search のカウントアップ
変数 search をカウントアップ (serch+=1) するのは、以下の赤文字で示した 2 箇所になります。

prime_list=[2,3,5,7,11]
for b in prime_list:
 search+=1    ◀ ここで search をカウントアップ
 d=a/b
 if frac(d)==0ck_pwr()
 if b>c: break


increment=\
[2,4,2,4,6,2,6,4,
 2,4,6,6,2,6,4,2,
 6,4,6,8,4,2,4,2,
 4,8,6,4,6,2,4,6,
 2,6,6,4,2,4,6,2,
 6,4,2,4,2,10,2,10]

while 1:
 for i in increment:
  search+=1   ◀ ここで search をカウントアップ
  b+=id=a/b
  if frac(d)==0ck_pwr()
  if b>cbreak
 if b>c
break



変数 search の表示
カウントアップした search の表示は、disp() 関数の内で、以下の赤文字で示したようにします。

def disp():
 global a,e,f,z
 clear_screen()
 locate(00f3'm'0)
 locate(160': ' str(len(inp)) + ' digits'3'm'0)
 locate(16, 11, 'search:'+str(search), 2, 'm', 0)   ◀ search の表示
 line(0,15,383,15,1,0)
 for i in range(115):
  if i<=e:
   dx=int(i/12)*16
   dy=int(i/12)*11
   locate(0+dxi-dyz[i], 1'm'0)
   locate(10+dxi-dy'^('2'm'0)
   locate(12+dxi-dyz[i+21], 1'm'0)
   locate(15+dxi-dy')'2'm'0)
 locate
(00''0'm'1)



10.5 素因数が10桁を超えるときの表示の変更 - disp() の変更

入力桁数を15桁まで拡張したので、見つかった素因数が10桁を超えることがあります。
例えば、以下のようなケースです。

11digit_prime 

4つめの素因数 58402449151 (11桁) が、乗数表示 ^(1 ) の表示と被ってしまっています。
これでは、ちょっと格好悪いので、素因数が10桁を超えた時は、超えた桁数に応じて乗数表示を右に移動して、被らないようにします。

disp() 関数内では、素因数は z[i] なので、素因数の桁数は、len(str(z[i])) で分かります。
この桁数が 10 を超えるとき、その超過分は len(str(z[i])) - 10 で求められます。これを変数 r に代入しておきます。

この r を使えば、乗数表示にかかわる locate() の第1引数に r を加えれば、被らなくなりなります。

この変更を行った disp() 関数は以下のようになります。変更部分は赤文字で示します。

def
disp():
 global a,e,f,z
 clear_screen()
 locate(00f3'm'0)
 locate(160': ' str(len(inp)) + ' digits'3'm'0)
 line(0,15,383,15,1,0)
 for i in range(115):
  if i<=e:
   r=0
   if len(str(z[i]))>10:
    r=len(str(z[i]))-10
   dx=int(i/12)*16
   dy=int(i/12)*11
   locate(0+dxi-dyz[i], 1'm'0)
   locate(10+dx+ri-dy'^('2'm'0)    ◀ 第1引数に r を加算する
   locate(12+dx+ri-dyz[i+21], 1'm'0)   ◀ 第1引数に r を加算する
   locate(15+dx+ri-dy')'2'm'0)     ◀ 第1引数に r を加算する
 locate
(00''0'm'1)



変更したスクリプトで、上記の "被ってしまった" ケースを入力すると、以下のように改善されたことが分かります。
Case1_14digit 


別のケースとして、素因数が15桁になるときは、以下のように乗数表示は17桁めに移動して、被りません。
15digit_prime 

なお、このケースでは、探索回数 search が、なんと434万回を超えています。

==========

以上で、15桁入力への拡張版が完成です。

※ 高速素因数分解 15桁入力対応版 - FactorG3.py のダウンロード

fx-CG50 Pythonモード:高速素因数分解 - FactorG3.py
"""Sample script

Exercise;
ported from Casio Basic
"Prime Factor 10 digits"
factorG.py
ver 3.0

by Krtyski/e-Gadget
"""

from u import *
digit=
15
search=0

def ckpwr():
 global a,b,c,d,e,z
 e+=1
 z[e]=b
 while 1:
  a=int(d)
  z[e+11]+=1
  d=a/b
  if frac(d):
   break
 c=int(sqrt(a))

def
disp():
 global a,e,f,z
 clear_screen()
 locate(00f3'm'0)
 locate(160': ' str(len(inp)) + ' digits'3'm'0)
 locate(16, 11, 'search:'+str(search), 2, 'm', 0)
 line(0,15,383,15,1,0)
 for i in range(115):
  if i<=e:
   dx=int(i/12)*16
   dy=int(i/12)*11
   locate(0+dxi-dyz[i], 1'm'0)
   locate(10+dxi-dy'^('2'm'0)
   locate(12+dxi-dyz[i+21], 1'm'0)
   locate(15+dxi-dy')'2'm'0)
 locate
(00''0'm'1)


while
1:
 try:
  inp=str(eval(input('Number:')))
 except (SyntaxError,TypeError,NameError) as e:
  print(e)
  print('*must be number or\n expression')
  continue
 if '.' in inp:
  print('*must be integer')
  continue
 elif inp.isdigit():
  if len(inp)>digit:
   print('*must be '+str(digit)+' digit\n or less')
   continue
  else:
   f=int(inp)
   break
 else:
  continue

z=list(range(23))
for e in range(1,23):
 z[e]=0
e=0
a=f
c=int(sqrt(a))

prime_list=[2,3,5,7,11]
for b in prime_list:
 search+=1
 d=a/b
 if frac(d)==0ck_pwr()
 if b>c: break


increment=\
[2,4,2,4,6,2,6,4,\
 2,4,6,6,2,6,4,2,\
 6,4,6,8,4,2,4,2,\
 4,8,6,4,6,2,4,6,\
 2,6,6,4,2,4,6,2,\
 6,4,2,4,2,10,2,10]

while 1:
 for i in increment:
  search+=1
  b+=id=a/b
  if frac(d)==0ck_pwr()
  if b>cbreak
 if b>c
break

if a>1:
 e+=1
 z[e]=a
 a=1
 z[e+11]=
1

disp
()






目 次

前の記事 - 9. Python らしい反復処理

次の記事 - 11. 関数呼出オーバーヘッド





応援クリックをお願いします。励みになるので...
にほんブログ村 IT技術ブログ 開発言語へ


 


keywords: fx-CG50Pythonプログラミング入門プログラミングプログラム関数電卓

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

関連記事

テーマ : プログラム関数電卓
ジャンル : コンピュータ

コメントの投稿

非公開コメント

最新記事
検索フォーム
最新コメント
カテゴリ
C# (3)
Visitors
Online Counter
現在の閲覧者数:
プロフィール

やす (Krtyski)

Author:やす (Krtyski)
since Oct 30, 2013


プログラム電卓は、プログラムを作って、使ってナンボ!

プログラム電卓を実際に使って気づいたこと、自作プログラム、電卓での Casio Basic や Casio Python プログラミングについて書いています。

なお管理人はカシオ計算機の関係者ではなく、Casio Basicが面白いと感じる1ユーザーです。


写真: 「4駆で泥んこ遊び@オックスフォード郊外」

リンク
月別アーカイブ
Sitemap

全ての記事を表示する

ブロとも申請フォーム

この人とブロともになる

QRコード
QR