HSPポータル
サイトマップ お問い合わせ


HSPTV!掲示板


未解決 解決 停止 削除要請

2007
0404
匿名希望莫大な少数の扱い9解決


匿名希望

リンク

2007/4/4(Wed) 23:05:14|NO.6982

またもや失礼。
今度は super π のような円周率を求めるソフトを
作成しようと思っているのですが、莫大な桁数の
少数を保管できる方法がわかりません。
よって、ここに質問させていただきます。



この記事に返信する


匿名希望

リンク

2007/4/4(Wed) 23:08:29|NO.6983

追伸、

扱う桁数は少なくても70000桁以上、
循環しない少数でお願いします。
図々しいようですが、よろしくお願いいたします。



GENKI

リンク

2007/4/4(Wed) 23:20:39|NO.6984

…以下、思い付きで言っているだけなので軽く流してください。orz
テキストデータとしてsdimで… なんて、やっぱりだめですか?(^ ^;
数値として保持するならある範囲の桁ごとに分けて複数の変数に保持するとかでしょうか…?
普通はどうやってるんでしょうね。



93

リンク

2007/4/5(Thu) 03:09:19|NO.6986

おもしろそうだったので作りました。
仕様で小数点はまだ入れてません。
あと、小さな値(差)同士でないバグる状態なのでそこの改善をしてみてください。

結構はまって考えてました。おもしろいっす。


#module // // 桁を保存する擬似変数 // // p1 = 計算する数(文字列) // p2 = 計算する数(文字列) // p3 = 桁数 // #deffunc ExpCacle str p1, str p2, int p3 // 文字列を実数として考える sdim vtr, 64, 2 vtr(0) = p1 vtr(1) = p2 dim ptr, 2 // .(小数点)の位置を後ろから数えて repeat 2 ptr(cnt) = instr(vtr(cnt), 0, ".") : ptr(cnt) = (strlen(vtr(cnt))-ptr(cnt)) loop // とりあえず.(小数点)を消す repeat 2 // 25.020だったら 25を取得 front = strmid(vtr(cnt), 0, strlen(vtr(cnt))-ptr(cnt)) // 25.020だったら 020を取得 back = strmid(vtr(cnt), strlen(vtr(cnt))-ptr(cnt)+1, strlen(vtr(cnt))) // あわせる 25020 vtr(cnt) = front+back loop // 遠い方の回数を保存 if (ptr(0) > ptr(1)) { ev = ptr(0) } else { ev = ptr(1) } // 文字列のを実数になおして(回数分-1)の10を掛ける repeat 2 : i = cnt repeat ev-ptr(i) vtr(i) = vtr(i)+"0" loop loop // 桁数分の数値を予約 sdim exp, 2*p3 // 数の予約 dim ptr, 2 dim ans, 2 // 数を代入(strからint) repeat 2 ptr(cnt) = int(vtr(cnt)) loop // 繰り返し、割った数の余りを計算 repeat p3 : n = cnt // 倍数 ans(0) = ptr(0)/ptr(1) // 余り ans(1) = ptr(0)\ptr(1) // 桁n = v1 poke exp, n, str(ans(0)) // 桁の判断 if (ans(0) == 0)|(ans(1) < ptr(0)) { // もし倍数が 0 か、余りが前より小さいなら次に持ち越す ptr(0) = int(""+ans(1)+"0") } else { // もし倍数が 0 以外なら終わり break } loop // 小数点をいれる // ..考え中 return // // 擬似変数から桁を取得 // // p1 = はじめの桁 // p2 = 取得する桁数 // #defcfunc Expinfo int p1, int p2 // 桁数分の文字列を予約 sdim p, p2-p1 // 桁を取得 if (p1 == 0)&(p2 == 0) { p = strmid(exp, 0, strlen(exp)) } else { p = strmid(exp, p1, p2) } return p #global // ここは関係なし sdim buf, 70000, 2 pos 0, 0 : mesbox buf(0), ginfo(12), 200 pos 0, 220 : mesbox buf(1), ginfo(12), 200 // ここから dialog "7万桁の計算をします。10秒くらいまて!" // 計算がexp型にはいる ExpCacle str(55.1), str(8.0116), 70000 // exp型の中身を問い合わせる objprm 0, Expinfo() // objprm 1, (55.1/8.0116)



93

リンク

2007/4/5(Thu) 03:34:16|NO.6987

DOUBLE型と最後の桁が違うのは四捨五入されてるからです。(多分

あと70000桁の計算遅いように感じますが・・・objprmが遅いだけです。
mesに置き換えたら0.5〜1秒程度で終わります。

PCも早くなったもんだ。



りさ

リンク

2007/4/5(Thu) 21:22:02|NO.7015

こんにちわ。

>円周率を求めるソフトを作成しようと思っているのですが、
>莫大な桁数の少数を保管できる方法がわかりません。

GENKIさんが仰るように、dim または sdim を使用するのが良いかと思います。

具体的な方法は、円周率を求める方法によって変わってくるのではないでしょうか。
ただ、アークタンジェント法にしろ、算術幾何平均法にしろ、
無限に続く円周率を求めるのはHSPには向かないかと思います。

しかし、もしそのアルゴリズムを理解なさっていて、
円周率の近似値を求めるのが目的でないのであれば、近いことはできそうです。


一番簡単な方法を例に挙げます。

HSPは、標準では double型が限界で、円周率の近似値を sin関数で求めると、


mes " PI = 3.14159265\n" repeat 5 x = 3. repeat ( cnt ) x += sin( x ) loop mes str( cnt + 1 ) + "回目 : " + x loop stop

3回目で限界に達します。
これは sin(x) にも double型 にも原因があるためです。

sin関数は初等関数なので、実質上は無限回連続微分が可能なはずです。
そこで、sin関数をTaylor展開によって、置き換えてみます。

sin(x) = ... + { ( -1 )~( n - 1 ) } * { ( x~( 2n - 1 ) )/( ( 2n - 1 )! ) } + ...

( チルダ(~) は累乗を意味し、{ } は大括弧を意味します)

そうすると、
[A] = ( -1 )~( n - 1 )
[B] = ( x~( 2n - 1 ) )
[C] = ( ( 2n - 1 )! )

に分解できるので、それぞれを計算します。

n (int型・正の整数)を求める項数、t (double型)を合計とすると、

repeat ( n ) _cnt = cnt + 1 lp = ( 2*_cnt - 1 ) // [A] sign = ( _cnt\2 )*2 - 1 // [B] vp = 1.0 repeat ( lp ) vp *= x loop // [C] vw = 0 repeat ( lp ) vw += lp lp -- loop // sin(x) t += ( vp / vw * sign ) loop

で求まるはずです。

これを見ると、実際に無限大の実数が必要なのは、
vp *= x で“実数と整数の乗算”
t += ( vp / vw * sign ) で“実数と整数の除算”と“実数と実数の加減算”
の4点だけで済むことになります。

勿論、この方法でも double型の限界があります。
なので、double型もsdim型(str)に置き換えてみます。

無限長の型は、sdim命令で実現します。
これを「XLong型」とします。

以下は、
vp *= x で使うための、XLong × Int のサンプルモジュールです。


#define ctype XLongOut( %1, %2 ) ( %1 + "." + %2 ) #module XLONG_ // ■■■■■■■■■■■■■■■■■■/* by Lisa*/■ // ■ 無限長 XLong型 のサンプル // ■■■■■■■■■■■■■■■■■■■■■■■■■ #deffunc XLongDeclare var p1, var p2, str p3, str p4 sdim p1, strlen( p3 ) sdim p2, strlen( p4 ) p1 = p3 p2 = p4 return #deffunc XLongMultipleXI var p1, var p2, int p3 cry = 0 // 繰上げ用変数 repeat 2 if ( cnt ) { // 計算済みの小数部を、元の小数部変数に戻す sdim p2, strlen( val ) p2 = val // 受け取った整数部変数をコピー len = strlen( p1 ) sdim val, len val = p1 } else { // 受け取った小数部変数をコピー len = strlen( p2 ) sdim val, len val = p2 } // 乗算のループ repeat ( len ) // 後尾から先頭へインデックスを移動する idx = len - cnt - 1 // 読み出し、数値に置き換え、乗算をして、繰り上がりを足す num = ( peek( val, idx ) - '0' ) * p3 + cry // 一桁目を文字コードに戻して、書き換える poke val, idx, ( ( num \ 10 ) + '0' ) // 二桁目を繰り上げ変数へ保存する cry = int( num / 10 ) loop loop // 整数部での繰り上がりがあるか if ( cry ) :inc = str( cry ) :else :inc = "" // 計算済みの整数部を、元の整数部変数に戻す sdim p1, ( strlen( val ) + strlen( inc ) ) p1 = inc + val return #global /*// XLong型を使用してみる //*/ // 321.823546879012345678900001 を XLong型に代入 XLongDeclare i, d, "821", "323546879012345678900001" mes mes "元になる値\n" + XLongOut( i, d ) // XLong変数に 2 を掛ける XLongMultipleXI i, d, 2 mes mes "計算した値\n" + XLongOut( i, d ) stop

しかし、ご存知かと思いますが、円周率の世界では
このように筆算アルゴリズムで処理する方法はタブーとされています。
(処理が桁数の二乗になりますので)

これ以上詳しくは、調べてみないと分かりませんが、
このサンプルが何かのヒントになればと思い、長々と書かせて貰いました。

長文・乱文、本当に申し訳ありません。



りさ

リンク

2007/4/5(Thu) 21:31:19|NO.7016

間違ってたら、ごめんなさいっ(><;



匿名希望

リンク

2007/4/6(Fri) 19:24:59|NO.7065

皆さん、回答ありがとうございます。
りささん、アルゴリズムは ガウス・ルジャンドル4次の収束を使用します。
ソースコードまで書いていただいたおかげで解決したので、解決のしるしを入れておきます。
この質問に回答してくださった方すべてにお礼申し上げます。



匿名希望

リンク

2007/4/6(Fri) 19:32:21|NO.7066

追伸

りささんがHSPに向かないと申しておりましたが、そのとうりだと思います。
しかし、HSPで製作する事で、HSPの可能性を示そうと考えたから製作に走ったのです。
HSPでできない事はないのです。私もそう思っております。
皆さん本当にありがとうございます。



314

リンク

2007/5/14(Mon) 21:56:29|NO.8255

アルゴリズムにもよりますが7万桁のような少ない桁数の場合、
算術幾何平均の方法で円周率を算出するよりアークタンジェント法の方が
早く求まる可能性があります。

普通、乗算や開平の演算時間は桁数の2乗に比例して長くなるのですが、
AGMやFFT(高速フーリエ変換)といったアルゴリズムを使うことで演算時間を短縮し
計算スピードを早くさせることができます。(算術幾何平均の場合のみ)
すなわち、ギネスに挑戦するほどの多長倍の桁数を求めたいならFFTによるアルゴリズムの影響が
大きく、アークタンジェント法よりよっぽど早く求まるのですが
7万桁程度ならアークタンジェント法でも問題ないかと思われます。

このスクリプトは浮動小数の配列変数を使って円周率πを任意の桁数もとめるプログラムです。
参考程度にしてください・・・(むっちゃ見にくいけど)

screen 0,200,100
kt=100
mes "何桁求める?"
input kt,50,20,6
button goto "OK牧場!",*bokujo
stop
;ここで桁数入力
*bokujo
clrobj
await 1
color 0,0,0:boxf
kt=(kt+1)/9+1
ktt=kt*30.11-29
e=0.0:zxx=1
qa=0.0
;何回ループするかを算出
dimtype a,vartype("double"),kt
a=1.0
repeat ktt
vo=cnt
c=1.0*ktt-cnt
dc=c*2+1.0
e=0.0
repeat kt
a.cnt=a.cnt*c
e+=a.cnt
qa=int(e/dc)
a.cnt=1.0*qa
e=e-dc*qa
e=e*1000000000.0
loop
a+=1.0
loop
;グレゴリー級数のオイラー変換によるアークタンジェント法で任意の桁数のπを計算
c=0:qa=0
repeat kt
c=kt-cnt
a.c=a.c*2+qa
qa=int(a.c/1000000000):a.c-=1000000000.0*qa
loop
;繰り上がり
a+=2.0
sdim e,9,kt
repeat kt
if a.cnt<100000000:e.cnt+="0":if a.cnt<10000000:e.cnt+="0":if a.cnt<1000000:e.cnt+="0":if a.cnt<100000:e.cnt+="0":if a.cnt<10000:e.cnt+="0":if a.cnt<1000:e.cnt+="0":if a.cnt<100:e.cnt+="0":if a.cnt<10:e.cnt+="0"
e.cnt+=str(int (a.cnt))
loop
;文字列に変換
dim ktt,1:ktt=5
sdim aa,(kt-1)*9.5+10
aa="π=3."
repeat kt-1,1
if cnt\10=1:aa+="\n":ktt+=2
if cnt<kt-1:aa+=e.cnt:ktt+=9
if cnt=kt-1:aa+=strmid(e.cnt,0,7):ktt+=9
loop
bsave "円周率.txt",aa,ktt
end


もしガウス・ルジャンドルの4次収束の公式より遅かったらスイマセン・・・
ちなみにFFTは僕分かりませんが、カラツバ乗算による乗算時間短縮アルゴリズム
ならわかります。
xを配列一つあたりの丸め込み上限値として
2log10x桁の数
ax+b

cx+d
の乗算アルゴリズムを考えます。
普通に計算すれば4回の乗算a*c,a*d,b*c,c*dを計算して
acx^2+(ad+cb)x+bd,箸垢襪里、
カラツバ乗算は,鯤儼舛靴
acx^2+{-(a-b)(c-d)+ac+bd}x+bdとして
結果的に3回の乗算ac,bd,(a-b)(c-d)で計算することでlog10x桁の乗算1回分が消えることになります。
これをさらに分割して帰納的に用いることでN^2に比例した演算量が、N^(log23)≒N^1.58に減ることがわかります。
しかしこのカラツバ乗算はどんなに分割してもFFTより早くすることは出来ません。



ONION software Copyright 1997-2021(c) All rights reserved.