タイトルにある Feynman point ちうのは、円周率の一部に 9 が連続する箇所があり、それが何故か Feynman point と呼ばれているとの話であります。
ref. ファインマン・ポイント - Wikipedia
https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A1%E3%82%A4%E3%83%B3%E3%83%9E%E3%83%B3%E3%83%BB%E3%83%9D%E3%82%A4%E3%83%B3%E3%83%88
円周率の小数点以下 762桁から "999999" が現れるとの事。
そういや、CASIO の Topics ページに、以下の記事がありました。
ref. 小学3年生の小原さん 関数電卓でPythonを動かす「シン・電卓アート」を制作 - CASIO Topics
https://www.casio.co.jp/topics/article/2024/K-055/
python で script を組んで、円周率 1000桁を憶えようとしている小学3年生の記事であります。
時代はここまで進んでおりました。マイッタ。
過去に、100円ショップの電卓をいじって興じている幼児をみて「電卓で遊ぶ幼児は理学者の夢を見るか」というネタを書いた憶えがあるのですが、現実はもっと進んでいたという次第。
25歳の頃「円周率を憶えるゾ !」と頑張ってはみたものの、16桁程度で満足してしまったバカチンの当方としては、オツムの体操と称して、円周率を求める script を使って遊んでみたのですが、その辺りの話を書いておきますヨ。
円周率を計算する upython script は以下の様であります。(実は再掲)
def machin(n):これを使うと、(原理的には)所望の桁数分、円周率を計算をしてくれるのですが、実際には丸め誤差があるので、余裕をみて計算桁を指定するのがよいです。
sm=0
term_1=10**n//5
term_2=10**n//239
flg=0
for j in range(n):
dv=1+2*j
if flg==0:
sm+=term_1//dv*4
sm-=term_2//dv
term_1=term_1//5//5
term_2=term_2//239//239
else:
sm-=term_1//dv*4
sm+=term_2//dv
term_1=term_1//5//5
term_2=term_2//239//239
flg=1-flg
return (sm*4)
例えば 500桁を計算するならば、510桁とか。
で、fx-CG50 の upython でコレを実行すると、およそ 508までは計算できました ... Feynman point の 762桁には届かない。ウーム。
800桁程度まで計算できるならば、Feynman Point の部分を取り出すため、つぎの追加をする所なのですが、
print(str(machin(800))[760:769])CASIO python では無理ゲーでした、残念。
しかし、諦めが悪いので Python Extra などで試してみましたら、ナント ! 計算できてしまいますネ ! コリャーエエ。
調子に乗って、XCAS でも試してみましたら、これでも計算できます。いや、不可解。
メモリ管理などの方法が違う、という感じなのかも知れません。
この 「Feynman point を表示する」という設問については、C.BASIC で実現できそうな気がするのですが、そのために、C.BASIC のマニュアルを眺めていたら、その機能の多さに圧倒されてしまい、勉強が進んでおりません。
折を見て、書けたらと思います。
ス・キ・魔 :
fx-9750GIII の CASIO python では、machin(254) まで計算してくれましたが、メモリの少ない 9750GIII では致し方のない所であります。
と・こ・ろ・が ! XCAS では実行できますネ。ただ、Clock UP の仕掛けがないので、40秒くらいの時間が掛かります。
str(machin(800))[760:769]
"349999998"
akatuki様
返信削除fx－CG100ですと n=688 までは計算でて、出力が 872146844 となりました。
689以上になると Memory Error となります。
整数要素のリストは1000要素は問題ないことは確認済み。
するとスタック不足が疑われるので、計算方法を工夫すればなんとかなりそうですね。
akatuki様
削除こんなん出来ました。
fx-CG100で、990桁から1000桁までを出力できました。
>>> from PiCalc10 inport *
From digit#990 to 1000: 09216420198
このページを見ると、一致するので正しそう：
https://ryo.blue/archive/%e5%86%86%e5%91%a8%e7%8e%87-1000%e6%a1%81/
スクリプトは以下：
def n_th_digit_of_pi(n):
def s(j, n):
s = 0.0
for k in range(n + 1):
r = 8 * k + j
s = (s + pow(16, n - k, r) / r) % 1.0
for k in range(n + 1, n + 10):
s = (s + 16**(n - k) / (8 * k + j)) % 1.0
return s
# Get digit #n using BBP equasion
x = (4 * s(1, n) - 2 * s(4, n) - s(5, n) - s(6, n)) % 1.0
return int(x * 10) # Return an approximation as a dev\cimal value
def print_pi_range_optimized(start, end):
print("From digit#"+str(start)+" to "+str(end)+": ", end="")
# streaming-ish calculation by saving memory
q, r, t, k, m, x = 1, 0, 1, 1, 3, 3
count = 0
while count <= end:
if 4 * q + r - t < m * t:
count += 1
if start <= count <= end:
print(m, end="")
q, r, t, k, m, x = (10 * q, 10 * (r - m * t), t, k,
(10 * (3 * q + r)) // t - 10 * m, x
)
else:
q, r, t, k, m, x = (
q * k, (2 * q + r) * x, t * x, k + 1,
(q * (7 * k + 2) + r * x) // (t * x), x + 2
)
print()
# Calculate
print_pi_range_optimized(990, 1000)
=====
def s(j, n)では、
sum_{k=0}^{n} (16^(n-k) mod (8k+j)) / (8k+j) により
巨大な16^nを直接扱わず、pow(16, n-k, mod) でメモリを節約
def n_th_digit_of_pi(n)では、
BBP公式によるn桁目の16進数算出を行って、
10進数への変換精度を確保するため少し前から計算します。
で、10進数としての近似値を返すようにします。
def print_pi_range_optimized(start, end)では、
実際には、10進数で特定の桁を直接出すのは数学的に複雑なため、
1000桁程度の「整数演算」をメモリ最小限で行う手法を適用します。
数値が得られたらストリーミング的にメモリを使わずに計算する感じです。
一番最後の行で、表示開始の桁数と終了の桁数を即値で指定すると、色々計算できます。
akatuki様
削除fx-CG100を使って、PiCalc10.py で、最後の行を以下のようにして、実行してみました。
print_pi_range_optimized(762, 772)
すると、 Shell画面で以下のようになりました：
>>> from PiCalc10 import*
From digit#762 to 772: 49999998372
763桁目から9が6連続になるようです。
akatuki様
返信削除インデックス (count) は0から始まるから、最後の関数の start と end は、実際の桁数より1少ない値にしないと、ですね。
akatuki様
返信削除このスクリプトを応用すると、円周率の多桁出力ができてしまいそうです。Shell画面への出力のスタックは300行までOKなので、1行10桁とすると3000桁まで、1行20桁だと6000桁まで出力できる計算です。
面白そうです。
akatuki様、
返信削除3つ目以降の私のコメントを削除しただけませんか。連投でご迷惑をおかけして、すみません。
やす (Krtyski) 様、早々のお越し、多謝であります !
返信削除> fx－CG100ですと n=688 までは計算でて、出力が 872146844 となりました。
> 689以上になると Memory Error となります。
fx-CG50 より計算桁数が多いでありますネ。
python で扱う work memory (RAM) が多いのかも ?
> 整数要素のリストは1000要素は問題ないことは確認済み。
> するとスタック不足が疑われるので、計算方法を工夫すればなんとかなりそうですね。
仰るように、多桁整数を複数保持しているので、memory がアップアップしているのだと、ですネ。
> こんなん出来ました。
泉アツノさんじゃないですか !?
それは冗談でありますが、圧巻であります !
BBP計算って、知りませんでした。最新の成果でありますネ。不勉強の誹りを ...
> 763桁目から9が6連続になるようです。
> インデックス (count) は0から始まるから、最後の関数の start と end は、実際の桁数より1少ない値にしないと、ですね。
最初の「3」(index : 0) は小数点数ではなく、indexが小数点以下の桁と一致する筈ですから、763目桁が「小数点以下 762桁目」となる様であります。
> このスクリプトを応用すると、円周率の多桁出力ができてしまいそうです。Shell画面への出力のスタックは300行までOKなので、1行10桁とすると3000桁まで、1行20桁だと6000桁まで出力できる計算です。
素晴らしいであります !
多桁出力もそうでありますが、BBP計算、親分の blog に掲載して公開するのが断然よい成果であります !
( 既に作業中であるとは思いますが )
akatuki様、
削除泉アツノをご存知とは...
只者では、ござらんな
ところで、fx-CG100のshell画面出力のスタックは200行で、CG50と同じでした。上で300行と書いたのは誤りでした。訂正致します。
akatuki様、
返信削除ご指摘のように、円周率の最初の3を桁数に入れてしまっていたので、小数点以下の桁数で処理をするように修正しました。
これまで count としていた桁数変数を decimal_count に変更し、これは小数点以下の桁数にしました。そして decimal_count = 0 の時は最初の '3' が現れる時なので、表示もカウントアップもしない (つまり、何もしないでループを次に回す) といった手抜き改造をしたら、うまくゆきました。
BBP公式を使った計算があることを知って、始めてトライしました。以前、ここで円周率の多桁表示で盛り上がった時は、まさかカシオ電卓で Python が走るなんて思いもしなかったので、今回は Casio Python のスクリプトを追加できたので、楽しかったです！
ありがとうございます(^^)/
スクリプトファイルを以下からダウンロードできるようにしました：
https://egadget2.web.fc2.com/python/PiCalc10/PiCalc10.html
電卓では日本ご表示ができないので、キャプションを無理矢理英語にしているので、分かりにくいです(ゴメンナサイ)。