円周率

　　円周率 匿名希望 2006/12/11 12:05:41 
　　├数列の和の計算の仕方が 白石　和夫 2006/12/11 13:28:18 
　　├ライプニッツは、収束が極めて遅いです。 SECOND 2006/12/11 14:19:18 
　　│├!オマケ。 SECOND 2006/12/12 02:08:36 
　　│└!オマケ２収束の速いマチンの公式 SECOND 2006/12/12 20:43:13 
　　│　└最強公式 河川屋 2006/12/14 00:01:22 
　　│　　└びっくりの速さ。 SECOND 2006/12/14 02:08:34 
　　│　　　└無理数 河川屋 2006/12/15 12:36:10 
　　│　　　　└できれば、２０００桁のサンプル・プログラ... SECOND 2006/12/15 16:38:36 
　　│　　　　　├これで動くには動きますが.... 河川屋 2006/12/15 20:09:25 
　　│　　　　　│└十進BASICで動くものを、お待ちします。 SECOND 2006/12/16 00:19:32 
　　│　　　　　│　└嫌です. 河川屋 2006/12/23 20:13:43 
　　│　　　　　└SUBAVE(A(),B(),C(),KETA) 河川屋 2006/12/15 20:11:01 
　　├式をどう見るかでいろいろなプログラムがで... 山中和義 2006/12/12 09:30:57 
　　│├「辛さ」ですか！？わたしは老眼の辛さだけ... SECOND 2006/12/12 21:35:30 
　　│└MicrosoftAccessのOLEも加えて頂けますか。... SECOND 2006/12/13 01:51:19 
　　└Re:円周率 河川屋 2006/12/23 19:50:10

　　円周率 匿名希望 2006/12/11 12:05:41  ツリーへ

円周率 返事を書く

匿名希望 2006/12/11 12:05:41

ライプニッツの公式による円周率を近似するのに
第ｎ項まで計算して近似値を与える関数を作るにはどのようにすればいいのですか？
教えてください

　　├数列の和の計算の仕方が 白石　和夫 2006/12/11 13:28:18  ツリーへ

Re: 円周率 返事を書く

白石　和夫 2006/12/11 13:28:18

数列の和の計算の仕方が
http://hp.vector.co.jp/authors/VA008683/tutorial/section2.htm
にあります。
また，関数定義は
http://hp.vector.co.jp/authors/VA008683/tutorial/section4.htm#4.1
にあるので，参考にしてください。
要は，級数の第k項をkの式で書けるかどうかの問題でしょう。

　　├ライプニッツは、収束が極めて遅いです。 SECOND 2006/12/11 14:19:18  ツリーへ

Re: 円周率 返事を書く

SECOND 2006/12/11 14:19:18

ライプニッツは、収束が極めて遅いです。
ubasic の結果ですが、ご参考。

list
100 'ライプニッツの級数　π/4= Σk=1~n (-1)^(k-1)/(2k-1)
110 word 540
120 N=10000000
130 Sum=0
140 for K=1 to N
150 Sum=Sum+(-1)^(K-1)/(2*K-1)
160 'print 4*Sum;
170 next K
180 print 4*Sum
190 end
OK
run
長変数のワード数は 540 (内部計算のワード数は 540)
3.1415925535897932388
OK

　　│├!オマケ。 SECOND 2006/12/12 02:08:36  ツリーへ

Re: ライプニッツは、収束が極めて遅いです。 返事を書く

SECOND 2006/12/12 02:08:36

!オマケ。
!ライプニッツの級数 π/4= Σk=1~n (-1)^(k-1)/(2k-1)

!LET n= 10000000
!PRINT 4*ライプニッツ(n)
!END
!n=10000000 の結果 3.14159255358802　　　 ← 十進BASIC
!n=10000000 の結果 3.1415925535897932388　← ubasic

!----ステップ・モニター
FOR k=1 TO 10
PRINT "k=";k;"項"; TAB(15);(-1)^(k-1);"/";(2*k-1); TAB(30);4*ライプニッツ(k)
NEXT K

END

!----
EXTERNAL FUNCTION ライプニッツ(n)
LET sum= 0
FOR k=1 TO n
LET sum= sum + (-1)^(k-1)/(2*k-1)
NEXT k
LET ライプニッツ= sum
END FUNCTION

　　│└!オマケ２収束の速いマチンの公式 SECOND 2006/12/12 20:43:13  ツリーへ

Re: ライプニッツは、収束が極めて遅いです。 返事を書く

SECOND 2006/12/12 20:43:13

!オマケ２収束の速いマチンの公式
!----
!Pai/4= 4*arctan(1/5) -arctan(1/239)

!　　 = 4*Σk=1~∞ { (-1)^(k-1)*(1/　5)^(2k-1)/(2k-1) }
!　　　　-Σk=1~∞ { (-1)^(k-1)*(1/239)^(2k-1)/(2k-1) }

OPTION ARITHMETIC DECIMAL_HIGH !10進1000桁

LET sum=0
FOR k=1 TO 10 !maxim=212
LET sum=sum+( 4/5^(2*k-1) -1/239^(2*k-1) ) *(-1)^(k-1) /(2*k-1)
PRINT 4*sum
NEXT K
PRINT TAB(2);4*sum-PI
PRINT TAB(3);".00… の所までが誤差無し。"

!CALL 十進1000桁π !必要なら実行。
END

!----
EXTERNAL SUB 十進1000桁π
OPTION ARITHMETIC DECIMAL_HIGH !10進1000桁
LET w$=STR$(PI)
PRINT
PRINT "十進BASIC 内蔵の 1000 桁のπ"
PRINT left$(w$,2)
FOR m=3 TO 1001 STEP 64
FOR n=m TO m+63 STEP 8
PRINT mid$(w$,n,8);" ";
NEXT N
PRINT
NEXT M
END SUB

　　│　└最強公式 河川屋 2006/12/14 00:01:22  ツリーへ

Re: !オマケ２収束の速いマチンの公式 返事を書く

河川屋 2006/12/14 00:01:22

最強公式

収束回数が少なくていいのは、ガウスールジャンドル公式ですね。
ループを回る毎に正しい桁が２倍になっていくから、000桁求めるのにループは１０回くらい。
ただし、式の性質上、最後の何桁かは正確でないかも。

!' 円周率の計算　　ガウス－ルジャンドル公式
OPTION ARITHMETIC DECIMAL_HIGH !10進1000桁
LET A=1
LET B=1/SQR(2)
LET T=0.25
LET X=1
LET PI0=0 !１回前のπの近似値
LET GOSA=1
LET COUNT=0
DO WHILE GOSA>1e-990 ! ABS(A-B)>1e-990 !後ろのほうはπ計算をループ外に出した場合。
LET Y=A
LET A=(A+B)/2
LET B=SQR(B*Y)
LET T=T-(Y-A)^2/X
LET X=X/2
LET PAI=((A+B)/2)^2/T !ループの外でも可だが、収束を見るため中に記述。
LET GOSA=ABS(PAI-PI0) !GOSA=ABS(PAI-PI) !ここでPIを使うのは一般論として反則技。
LET PI0=PAI
LET COUNT=COUNT+1
PRINT USING "##### -%.###^^^^^":COUNT,PAI-PI !ループ回数と誤差を表示
LOOP
PRINT PAI !途中で折り返すルーチンにする必要があるが省略。
END

なお、サンプルプログラムの
SUUGAKUB\PIBAS　にも収束計算例が載っていますが、コレは
内接する正方形、８角形、１６角形....
の辺長を順に求める方法だから、ヴィエトの公式の系譜です。

　　│　　└びっくりの速さ。 SECOND 2006/12/14 02:08:34  ツリーへ

Re: 最強公式 返事を書く

SECOND 2006/12/14 02:08:34

びっくりの速さ。
でも、開平計算が、級数に含まれているのが、残念です。

例えば、２０００桁のπを求めようとすると、どうしても
演算レジスターなど、BASIC 側からのサービスが受けられ
なくなりますので、
プログラムで桁の展開をする事になりますね、その場合、開平
計算で、時折発生する無理数を、どうすればよいでしょう。

　　│　　　└無理数 河川屋 2006/12/15 12:36:10  ツリーへ

Re: びっくりの速さ。 返事を書く

河川屋 2006/12/15 12:36:10

無理数

＞開平計算で、時折発生する無理数を、どうすればよいでしょう。
2000桁の答が欲しいなら、√の答も2000桁正しければいいわけであり、
NEWTON法(サンプルプログラムの中にあります)を多倍長用に展開すれば
４則演算のみで計算はできます。
（開平法で解くのは面倒だし計算量が多い。）

あと、2000桁あればいいと書きましたが、
計算誤差対策のため、全計算を少し多めの桁(2050桁とか。)で計算しておき、
最後の数桁は使わないようにしないと、
正しく2000桁までは求められないです。

　　│　　　　└できれば、２０００桁のサンプル・プログラ... SECOND 2006/12/15 16:38:36  ツリーへ

Re: 無理数 返事を書く

SECOND 2006/12/15 16:38:36

できれば、２０００桁のサンプル・プログラムを御願いできますか。
すみません、ご面倒がけて・・

　　│　　　　　├これで動くには動きますが.... 河川屋 2006/12/15 20:09:25  ツリーへ

Re: できれば、２０００桁のサンプル・プログラ... 返事を書く

河川屋 2006/12/15 20:09:25

これで動くには動きますが....
あくまで動くだけ。
配列１個が１桁で覚える方式なので非常に遅いです。
多倍長計算もインチキ(とりあえず用を足せるだけ)です。
おまけに、QuickBasicです。

DEFINT A-Z
DECLARE SUB SQRT (A(), B(), C(), KETA) 'A=SQRT(B) Cは初期値
DECLARE SUB DIV (A(), B(), C(), KETA) 'A=B/C
DECLARE SUB MUL (A(), B(), C(), KETA) 'A=B*C
DECLARE SUB MOV (A(), B(), KETA) 'A=B
DECLARE SUB AVE (A(), B(), C(), KETA) 'A=(B+C)/2
DECLARE SUB SUBT (A(), B(), C(), KETA) 'A=B-C
DECLARE SUB CHK (A(), B(), FOREVER, KETA) 'Aとの一致桁数をFOREVERに返す
DECLARE SUB PRT (A(), KETA)

INPUT "KETA="; KETA
DIM A(2010), B(2010), C(2010), D(2010), E(2010), T(2010)
DIM X(2010), Y(2010), PI(2010), Z(2010)
B(0) = 2
PI(0) = 1: PI(1) = 4: 'PI=√2の近似値
X(0) = 1
K2 = KETA * 1.02 + 2 '桁は少し余分に。
CALL SQRT(A(), B(), PI(), K2)
CALL MOV(B(), A(), K2)
CALL DIV(B(), X(), A(), K2) ' B=1/SQR(2)
FOR I = 0 TO KETA: A(I) = 0: NEXT
A(0) = 1 ' A=1
T(1) = 2: T(2) = 5 ' T=0.25
FOR I = 0 TO K2: X(I) = 0: NEXT
X(0) = 1 ' X=1
FOREVERR = -1
WHILE FOREVERR
CALL MOV(Y(), A(), K2) ' Y=A
CALL AVE(C(), A(), B(), K2)
CALL MOV(A(), C(), K2) ' A=(A+B)/2
CALL MOV(C(), B(), K2)
CALL MUL(PI(), B(), Y(), K2)
CALL SQRT(B(), PI(), A(), K2) ' B=SQR(B*Y)
CALL SUBT(C(), Y(), A(), K2)
CALL MOV(PI(), C(), K2)
CALL MUL(D(), C(), PI(), K2)
CALL DIV(C(), D(), X(), K2)
CALL SUBT(D(), T(), C(), K2)
CALL MOV(T(), D(), K2) ' T=T-(Y-A)^2/X
CALL AVE(D(), X(), Z(), K2)
CALL MOV(X(), D(), K2) ' X=X/2
CALL CHK(A(), B(), STAT, K2)
IF STAT >= KETA THEN FOREVERR = 0 ELSE FOREVERR = STAT
WEND
CALL AVE(D(), A(), B(), K2)
CALL MOV(E(), D(), K2)
CALL MUL(C(), D(), E(), K2)
CALL DIV(PI(), C(), T(), K2)
PRINT "PI=3.";
FOR I = 1 TO KETA
PRINT HEX$(PI(I));
IF I MOD 5 = 0 THEN PRINT " ";
IF I MOD 50 = 0 THEN PRINT : PRINT " ";
NEXT
END

　　│　　　　　│└十進BASICで動くものを、お待ちします。 SECOND 2006/12/16 00:19:32  ツリーへ

Re: これで動くには動きますが.... 返事を書く

SECOND 2006/12/16 00:19:32

十進BASIC で動くものを、お待ちします。
いそがなくでも、ゆっくりで・・
有難うございました。

　　│　　　　　│　└嫌です. 河川屋 2006/12/23 20:13:43  ツリーへ

Re: 十進BASICで動くものを、お待ちします。 返事を書く

河川屋 2006/12/23 20:13:43

嫌です.

>十進BASIC で動くものを、お待ちします。
いえ、そこまでやる気がないから、QuickBasicのままなのです。
理由の1
あれは相当大昔に組んだもので、プロちゃんが組んだマチンの公式
による2000桁計算に速度でベタ負けします。
10進basicは、整数型も固定小数点型も無いから、
整数でいいところを実数計算するため、もっと悲惨になるでしょう。
※実は、JIS_BASIC規格に固定小数点型というのもあるが、実装必須でないため(?)
今のところサポートされていない。

理由の2
私はハンドル名のとおり河川屋なので、科学技術計算(土木技術計算?)
onlyであって、倍精度実数(約15～17桁)あればほとんど用が足りるので、
多倍長にはあまり興味がない。
かわりに、有効桁の制限がある中で、どうやれば誤差が累積しないか
というのをいれておきます。

理由の3
当初の記述
>例えば、２０００桁のπを求めようとすると、
>プログラムで桁の展開をする事になりますね、
>その場合、開平計算で、時折発生する無理数を、どうすればよいでしょう。
というのを、
・無理数をどうすれば良いか(適当に打ちきればええねん)
次には
・4則演算の多倍長計算なら自力で組めるが√は組めない

という意味にとったため。
どちらの意味でも、２０００桁のπのプログラムをそっくり
提示する必要は無いと感じ、あとからそっくり提示してくれ、
というのは反則に思えるのです。

　　│　　　　　└SUBAVE(A(),B(),C(),KETA) 河川屋 2006/12/15 20:11:01  ツリーへ

Re: できれば、２０００桁のサンプル・プログラ... 返事を書く

河川屋 2006/12/15 20:11:01

SUB AVE (A(), B(), C(), KETA)
AMARI = 0
FOR I = 0 TO KETA
J = I
T = (AMARI * 10 + B(I) + C(I))
AMARI = T MOD 2
A(I) = T \ 2
WHILE A(J) >= 10
A(J - 1) = A(J - 1) + A(J) \ 10
A(J) = A(J) MOD 10
J = J - 1
WEND
NEXT
END SUB

SUB CHK (A(), B(), FOREVER, KETA)
FOREVER = KETA
FOR I = 1 TO KETA
IF A(I) <> B(I) THEN FOREVER = I: EXIT SUB
NEXT
END SUB

SUB DIV (A(), B(), C(), KETA)
DIM X(2010), BB(2010), CC(2010)
CALL MOV(BB(), B(), KETA)
CALL MOV(CC(), C(), KETA)
FOR I = 0 TO KETA
FOR J = 0 TO 10
CALL SUBT(X(), BB(), CC(), KETA)
IF X(0) < 0 THEN
A(I) = J
EXIT FOR
ELSE
CALL MOV(BB(), X(), KETA)
END IF
NEXT
FOR J = KETA TO I + 1 STEP -1
CC(J) = CC(J - 1)
NEXT
CC(I) = 0
NEXT
END SUB

SUB MOV (A(), B(), KETA)
FOR I = 0 TO KETA: A(I) = B(I): NEXT
END SUB

SUB MUL (A(), B(), C(), KETA)
FOR I = 0 TO KETA: A(I) = 0: NEXT
FOR I = 0 TO KETA
FOR J = 0 TO KETA - I
K = I + J
A(K) = A(K) + B(I) * C(J)
IF A(K) >= 10 THEN A(K - 1) = A(K - 1) + A(K) \ 10: A(K) = A(K) MOD 10
NEXT
NEXT
FOR K = KETA TO 0 STEP -1
IF A(K) >= 10 THEN A(K - 1) = A(K - 1) + A(K) \ 10: A(K) = A(K) MOD 10
NEXT
END SUB

SUB PRT (A(), KETA)
FOR I = 0 TO KETA
PRINT HEX$(A(I));
NEXT
PRINT
END SUB

SUB SQRT (A(), B(), C(), KETA)
DIM X(2010), Y(2010), CC(2010)
CALL MOV(CC(), C(), KETA)
FOREVER = -1
WHILE FOREVER
CALL DIV(X(), B(), CC(), KETA + 1)
CALL AVE(Y(), X(), CC(), KETA + 1)
CALL CHK(X(), Y(), STAT, KETA + 1)
CALL MOV(CC(), Y(), KETA + 1)
IF STAT >= KETA THEN FOREVER = 0 ELSE FOREVER = STAT
WEND
CALL MOV(A(), Y(), KETA)
END SUB

SUB SUBT (A(), B(), C(), KETA)
AMARI = 0
FOR I = 0 TO KETA
A(I) = B(I) - C(I)
IF A(I) < 0 THEN
J = I - 1
WHILE A(J + 1) < 0 AND J >= 0
A(J) = A(J) - 1: A(J + 1) = A(J + 1) + 10
J = J - 1
WEND
END IF
NEXT
END SUB

　　├式をどう見るかでいろいろなプログラムがで... 山中和義 2006/12/12 09:30:57  ツリーへ

Re: 円周率 返事を書く

山中和義 2006/12/12 09:30:57

式をどう見るかでいろいろなプログラムができます。

!ケース１
!ライプニッツの公式　π/4=ATN(1)=1-1/3+1/5-1/7+1/9-1/11+ …

DEF A(n)=(-1)^(n-1)/(2*n-1) !第n項 n=1,2,3,…

LET t=TIME !計測開始
LET S=0
FOR i=1 TO 1000000 !収束が遅い
LET S=S+A(i)
NEXT i
PRINT S*4 !S=π/4
PRINT TIME-t !計測終了、表示

!ケース２
!π/4=ATN(1)=(1)/1+(-1)/3+(1)/5+(-1)/7+(1)/9+(-1)/11+ …

DEF Den(n)=(2*n-1) !第n項の分母
LET t=TIME
LET S=0
LET D=1 !分子（符号を含む）
FOR i=1 TO 1000000 !収束が遅い
LET S=S+D/Den(i)
LET D=-D !反転
NEXT i
PRINT S*4 !S=π/4
PRINT TIME-t

!分母の奇数をnの式（第n項）で書けない場合
LET t=TIME
LET S=0
LET D=1 !符号
FOR i=1 TO 2000000 STEP 2 !収束が遅い
LET S=S+D/i
LET D=-D !反転
NEXT i
PRINT S*4 !S=π/4
PRINT TIME-t

!ケース３
!π/4=ATN(1)=(1/1-1/3)+(1/5-1/7)+(1/9-1/11)+ …
LET t=TIME
LET S=0
LET i=1
DO WHILE i<=2000000 !収束が遅い
LET S=S+1/i !カッコ内の前の値
LET i=i+2
LET S=S-1/i !後の値
LET i=i+2
LOOP
PRINT S*4 !S=π/4
PRINT TIME-t

!精度　ケース２
LET t=TIME
LET S=0
LET D=1 !符号
LET i=1
DO WHILE Den(i)<10000000 !分母が大きければ０に近くなる
LET S=S+D/Den(i)
LET D=-D !反転
LET i=i+1
LOOP
PRINT "第";i-1;"項まで",S*4 !S=π/4
PRINT TIME-t

END

アイディアを出すなどはプログラミングの「楽しみ」で、
メモリ、実行時間、読みやすさ、言語仕様などの制約があり、これは「辛さ」です。

正解だけではなくいろいろ考えてみるとおもしろいです。

　　│├「辛さ」ですか！？わたしは老眼の辛さだけ... SECOND 2006/12/12 21:35:30  ツリーへ

Re: 式をどう見るかでいろいろなプログラムがで... 返事を書く

SECOND 2006/12/12 21:35:30

「辛さ」ですか！？わたしは老眼の辛さだけですが・・
もっと書いてください。おかげで、十進BASIC の魅力の
沢山のことが、覚えられました。有難うございます。

　　│└MicrosoftAccessのOLEも加えて頂けますか。... SECOND 2006/12/13 01:51:19  ツリーへ

Re: 式をどう見るかでいろいろなプログラムがで... 返事を書く

SECOND 2006/12/13 01:51:19

Microsoft Access の OLE も加えて頂けますか。全て重宝しております。

　　└Re:円周率 河川屋 2006/12/23 19:50:10  ツリーへ

Re: 円周率 返事を書く

河川屋 2006/12/23 19:50:10

Re:円周率

>ライプニッツの公式による円周率を近似するのに
>第ｎ項まで計算して近似値を与える関数を作るには
以下を回答例とします。

要点は以下のとおり。
1.関数(Function)で定義する。
2.第ｎ項までを足す。nはメインルーチンから引数で与える。
※4*(1-1/3+1/5) まで計算した場合、第3項までと数えている。
ここまでが質問に対する最低必要な要件。
以下、オマケの要点。
3.データ型の範囲内で、なるべく正確な答を返す。
このためには、n項から遡って足していく。
※正確な値とは、有効桁無限大の計算機でn項まで足した値と定義する。
※※2進モード(IEEE754倍精度)と10進15桁モードを比較すると
数学的にたまたま整数解となる場合は10進15桁が正確であるが、
循環小数が混じると、
4則演算の場合2進モードのほうが正確で、
超越関数を使うと10進15桁モードのほうが正確となるもよう。
※※※第1兆項あたりまで計算しない限り、あまり意味はないかも。
4.第n項めの値を計算するには、(-1)^(n-1)/(2*n-1)でも
数学では正しいけれど、一般にべき計算は遅いのでif文で代用した。
(整数演算は組込関数内部で別ルーチンを使っていればこの限りでは
ないが、念のためこうしている。)

DECLARE EXTERNAL FUNCTION LEIBNIZ
LET N=2 !第n項まで
LET PAI=LEIBNIZ(N)
PRINT PAI
END

EXTERNAL FUNCTION LEIBNIZ(N)
LET T=0
FOR I=N TO 1 STEP -1
IF mod(I,2)=1 THEN
LET T=T+1/(2*I-1)
ELSE
LET T=T-1/(2*I-1)
END IF
NEXT I
LET LEIBNIZ=4*T
END FUNCTION

インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。

円周率	返事を書く
匿名希望 2006/12/11 12:05:41
ライプニッツの公式による円周率を近似するのに第ｎ項まで計算して近似値を与える関数を作るにはどのようにすればいいのですか？教えてください

Re: 円周率	返事を書く
白石　和夫 2006/12/11 13:28:18
数列の和の計算の仕方が http://hp.vector.co.jp/authors/VA008683/tutorial/section2.htm にあります。また，関数定義は http://hp.vector.co.jp/authors/VA008683/tutorial/section4.htm#4.1 にあるので，参考にしてください。要は，級数の第k項をkの式で書けるかどうかの問題でしょう。

Re: 円周率	返事を書く
SECOND 2006/12/11 14:19:18
ライプニッツは、収束が極めて遅いです。 ubasic の結果ですが、ご参考。 list 100 'ライプニッツの級数　π/4= Σk=1~n (-1)^(k-1)/(2k-1) 110 word 540 120 N=10000000 130 Sum=0 140 for K=1 to N 150 Sum=Sum+(-1)^(K-1)/(2K-1) 160 'print 4Sum; 170 next K 180 print 4*Sum 190 end OK run 長変数のワード数は 540 (内部計算のワード数は 540) 3.1415925535897932388 OK

Re: ライプニッツは、収束が極めて遅いです。	返事を書く
SECOND 2006/12/12 02:08:36
!オマケ。 !ライプニッツの級数 π/4= Σk=1~n (-1)^(k-1)/(2k-1) !LET n= 10000000 !PRINT 4ライプニッツ(n) !END !n=10000000 の結果 3.14159255358802　　　 ← 十進BASIC !n=10000000 の結果 3.1415925535897932388　← ubasic !----ステップ・モニター FOR k=1 TO 10 PRINT "k=";k;"項"; TAB(15);(-1)^(k-1);"/";(2k-1); TAB(30);4ライプニッツ(k) NEXT K END !---- EXTERNAL FUNCTION ライプニッツ(n) LET sum= 0 FOR k=1 TO n LET sum= sum + (-1)^(k-1)/(2k-1) NEXT k LET ライプニッツ= sum END FUNCTION

Re: ライプニッツは、収束が極めて遅いです。	返事を書く
SECOND 2006/12/12 20:43:13
!オマケ２収束の速いマチンの公式 !---- !Pai/4= 4arctan(1/5) -arctan(1/239) !　　 = 4Σk=1~∞ { (-1)^(k-1)(1/　5)^(2k-1)/(2k-1) } !　　　　-Σk=1~∞ { (-1)^(k-1)(1/239)^(2k-1)/(2k-1) } OPTION ARITHMETIC DECIMAL_HIGH !10進1000桁 LET sum=0 FOR k=1 TO 10 !maxim=212 LET sum=sum+( 4/5^(2k-1) -1/239^(2k-1) ) (-1)^(k-1) /(2k-1) PRINT 4sum NEXT K PRINT TAB(2);4sum-PI PRINT TAB(3);".00… の所までが誤差無し。" !CALL 十進1000桁π !必要なら実行。 END !---- EXTERNAL SUB 十進1000桁π OPTION ARITHMETIC DECIMAL_HIGH !10進1000桁 LET w$=STR$(PI) PRINT PRINT "十進BASIC 内蔵の 1000 桁のπ" PRINT left$(w$,2) FOR m=3 TO 1001 STEP 64 FOR n=m TO m+63 STEP 8 PRINT mid$(w$,n,8);" "; NEXT N PRINT NEXT M END SUB

Re: !オマケ２収束の速いマチンの公式	返事を書く
河川屋 2006/12/14 00:01:22
最強公式収束回数が少なくていいのは、ガウスールジャンドル公式ですね。ループを回る毎に正しい桁が２倍になっていくから、000桁求めるのにループは１０回くらい。ただし、式の性質上、最後の何桁かは正確でないかも。 !' 円周率の計算　　ガウス－ルジャンドル公式 OPTION ARITHMETIC DECIMAL_HIGH !10進1000桁 LET A=1 LET B=1/SQR(2) LET T=0.25 LET X=1 LET PI0=0 !１回前のπの近似値 LET GOSA=1 LET COUNT=0 DO WHILE GOSA>1e-990 ! ABS(A-B)>1e-990 !後ろのほうはπ計算をループ外に出した場合。 LET Y=A LET A=(A+B)/2 LET B=SQR(B*Y) LET T=T-(Y-A)^2/X LET X=X/2 LET PAI=((A+B)/2)^2/T !ループの外でも可だが、収束を見るため中に記述。 LET GOSA=ABS(PAI-PI0) !GOSA=ABS(PAI-PI) !ここでPIを使うのは一般論として反則技。 LET PI0=PAI LET COUNT=COUNT+1 PRINT USING "##### -%.###^^^^^":COUNT,PAI-PI !ループ回数と誤差を表示 LOOP PRINT PAI !途中で折り返すルーチンにする必要があるが省略。 END なお、サンプルプログラムの SUUGAKUB\PIBAS　にも収束計算例が載っていますが、コレは内接する正方形、８角形、１６角形.... の辺長を順に求める方法だから、ヴィエトの公式の系譜です。

Re: 最強公式	返事を書く
SECOND 2006/12/14 02:08:34
びっくりの速さ。でも、開平計算が、級数に含まれているのが、残念です。例えば、２０００桁のπを求めようとすると、どうしても演算レジスターなど、BASIC 側からのサービスが受けられなくなりますので、プログラムで桁の展開をする事になりますね、その場合、開平計算で、時折発生する無理数を、どうすればよいでしょう。

Re: びっくりの速さ。	返事を書く
河川屋 2006/12/15 12:36:10
無理数＞開平計算で、時折発生する無理数を、どうすればよいでしょう。 2000桁の答が欲しいなら、√の答も2000桁正しければいいわけであり、 NEWTON法(サンプルプログラムの中にあります)を多倍長用に展開すれば４則演算のみで計算はできます。（開平法で解くのは面倒だし計算量が多い。）あと、2000桁あればいいと書きましたが、計算誤差対策のため、全計算を少し多めの桁(2050桁とか。)で計算しておき、最後の数桁は使わないようにしないと、正しく2000桁までは求められないです。

Re: 無理数	返事を書く
SECOND 2006/12/15 16:38:36
できれば、２０００桁のサンプル・プログラムを御願いできますか。すみません、ご面倒がけて・・

Re: できれば、２０００桁のサンプル・プログラ...	返事を書く
河川屋 2006/12/15 20:09:25
これで動くには動きますが.... あくまで動くだけ。配列１個が１桁で覚える方式なので非常に遅いです。多倍長計算もインチキ(とりあえず用を足せるだけ)です。おまけに、QuickBasicです。 DEFINT A-Z DECLARE SUB SQRT (A(), B(), C(), KETA) 'A=SQRT(B) Cは初期値 DECLARE SUB DIV (A(), B(), C(), KETA) 'A=B/C DECLARE SUB MUL (A(), B(), C(), KETA) 'A=BC DECLARE SUB MOV (A(), B(), KETA) 'A=B DECLARE SUB AVE (A(), B(), C(), KETA) 'A=(B+C)/2 DECLARE SUB SUBT (A(), B(), C(), KETA) 'A=B-C DECLARE SUB CHK (A(), B(), FOREVER, KETA) 'Aとの一致桁数をFOREVERに返す DECLARE SUB PRT (A(), KETA) INPUT "KETA="; KETA DIM A(2010), B(2010), C(2010), D(2010), E(2010), T(2010) DIM X(2010), Y(2010), PI(2010), Z(2010) B(0) = 2 PI(0) = 1: PI(1) = 4: 'PI=√2の近似値 X(0) = 1 K2 = KETA 1.02 + 2 '桁は少し余分に。 CALL SQRT(A(), B(), PI(), K2) CALL MOV(B(), A(), K2) CALL DIV(B(), X(), A(), K2) ' B=1/SQR(2) FOR I = 0 TO KETA: A(I) = 0: NEXT A(0) = 1 ' A=1 T(1) = 2: T(2) = 5 ' T=0.25 FOR I = 0 TO K2: X(I) = 0: NEXT X(0) = 1 ' X=1 FOREVERR = -1 WHILE FOREVERR CALL MOV(Y(), A(), K2) ' Y=A CALL AVE(C(), A(), B(), K2) CALL MOV(A(), C(), K2) ' A=(A+B)/2 CALL MOV(C(), B(), K2) CALL MUL(PI(), B(), Y(), K2) CALL SQRT(B(), PI(), A(), K2) ' B=SQR(B*Y) CALL SUBT(C(), Y(), A(), K2) CALL MOV(PI(), C(), K2) CALL MUL(D(), C(), PI(), K2) CALL DIV(C(), D(), X(), K2) CALL SUBT(D(), T(), C(), K2) CALL MOV(T(), D(), K2) ' T=T-(Y-A)^2/X CALL AVE(D(), X(), Z(), K2) CALL MOV(X(), D(), K2) ' X=X/2 CALL CHK(A(), B(), STAT, K2) IF STAT >= KETA THEN FOREVERR = 0 ELSE FOREVERR = STAT WEND CALL AVE(D(), A(), B(), K2) CALL MOV(E(), D(), K2) CALL MUL(C(), D(), E(), K2) CALL DIV(PI(), C(), T(), K2) PRINT "PI=3."; FOR I = 1 TO KETA PRINT HEX$(PI(I)); IF I MOD 5 = 0 THEN PRINT " "; IF I MOD 50 = 0 THEN PRINT : PRINT " "; NEXT END

Re: これで動くには動きますが....	返事を書く
SECOND 2006/12/16 00:19:32
十進BASIC で動くものを、お待ちします。いそがなくでも、ゆっくりで・・有難うございました。

Re: 十進BASICで動くものを、お待ちします。	返事を書く
河川屋 2006/12/23 20:13:43
嫌です. >十進BASIC で動くものを、お待ちします。いえ、そこまでやる気がないから、QuickBasicのままなのです。理由の1 あれは相当大昔に組んだもので、プロちゃんが組んだマチンの公式による2000桁計算に速度でベタ負けします。 10進basicは、整数型も固定小数点型も無いから、整数でいいところを実数計算するため、もっと悲惨になるでしょう。 ※実は、JIS_BASIC規格に固定小数点型というのもあるが、実装必須でないため(?) 今のところサポートされていない。理由の2 私はハンドル名のとおり河川屋なので、科学技術計算(土木技術計算?) onlyであって、倍精度実数(約15～17桁)あればほとんど用が足りるので、多倍長にはあまり興味がない。かわりに、有効桁の制限がある中で、どうやれば誤差が累積しないかというのをいれておきます。理由の3 当初の記述 >例えば、２０００桁のπを求めようとすると、 >プログラムで桁の展開をする事になりますね、 >その場合、開平計算で、時折発生する無理数を、どうすればよいでしょう。というのを、・無理数をどうすれば良いか(適当に打ちきればええねん) 次には・4則演算の多倍長計算なら自力で組めるが√は組めないという意味にとったため。どちらの意味でも、２０００桁のπのプログラムをそっくり提示する必要は無いと感じ、あとからそっくり提示してくれ、というのは反則に思えるのです。

Re: できれば、２０００桁のサンプル・プログラ...	返事を書く
河川屋 2006/12/15 20:11:01
SUB AVE (A(), B(), C(), KETA) AMARI = 0 FOR I = 0 TO KETA J = I T = (AMARI * 10 + B(I) + C(I)) AMARI = T MOD 2 A(I) = T \ 2 WHILE A(J) >= 10 A(J - 1) = A(J - 1) + A(J) \ 10 A(J) = A(J) MOD 10 J = J - 1 WEND NEXT END SUB SUB CHK (A(), B(), FOREVER, KETA) FOREVER = KETA FOR I = 1 TO KETA IF A(I) <> B(I) THEN FOREVER = I: EXIT SUB NEXT END SUB SUB DIV (A(), B(), C(), KETA) DIM X(2010), BB(2010), CC(2010) CALL MOV(BB(), B(), KETA) CALL MOV(CC(), C(), KETA) FOR I = 0 TO KETA FOR J = 0 TO 10 CALL SUBT(X(), BB(), CC(), KETA) IF X(0) < 0 THEN A(I) = J EXIT FOR ELSE CALL MOV(BB(), X(), KETA) END IF NEXT FOR J = KETA TO I + 1 STEP -1 CC(J) = CC(J - 1) NEXT CC(I) = 0 NEXT END SUB SUB MOV (A(), B(), KETA) FOR I = 0 TO KETA: A(I) = B(I): NEXT END SUB SUB MUL (A(), B(), C(), KETA) FOR I = 0 TO KETA: A(I) = 0: NEXT FOR I = 0 TO KETA FOR J = 0 TO KETA - I K = I + J A(K) = A(K) + B(I) * C(J) IF A(K) >= 10 THEN A(K - 1) = A(K - 1) + A(K) \ 10: A(K) = A(K) MOD 10 NEXT NEXT FOR K = KETA TO 0 STEP -1 IF A(K) >= 10 THEN A(K - 1) = A(K - 1) + A(K) \ 10: A(K) = A(K) MOD 10 NEXT END SUB SUB PRT (A(), KETA) FOR I = 0 TO KETA PRINT HEX$(A(I)); NEXT PRINT END SUB SUB SQRT (A(), B(), C(), KETA) DIM X(2010), Y(2010), CC(2010) CALL MOV(CC(), C(), KETA) FOREVER = -1 WHILE FOREVER CALL DIV(X(), B(), CC(), KETA + 1) CALL AVE(Y(), X(), CC(), KETA + 1) CALL CHK(X(), Y(), STAT, KETA + 1) CALL MOV(CC(), Y(), KETA + 1) IF STAT >= KETA THEN FOREVER = 0 ELSE FOREVER = STAT WEND CALL MOV(A(), Y(), KETA) END SUB SUB SUBT (A(), B(), C(), KETA) AMARI = 0 FOR I = 0 TO KETA A(I) = B(I) - C(I) IF A(I) < 0 THEN J = I - 1 WHILE A(J + 1) < 0 AND J >= 0 A(J) = A(J) - 1: A(J + 1) = A(J + 1) + 10 J = J - 1 WEND END IF NEXT END SUB

Re: 円周率	返事を書く
山中和義 2006/12/12 09:30:57
式をどう見るかでいろいろなプログラムができます。 !ケース１ !ライプニッツの公式　π/4=ATN(1)=1-1/3+1/5-1/7+1/9-1/11+ … DEF A(n)=(-1)^(n-1)/(2n-1) !第n項 n=1,2,3,… LET t=TIME !計測開始 LET S=0 FOR i=1 TO 1000000 !収束が遅い LET S=S+A(i) NEXT i PRINT S4 !S=π/4 PRINT TIME-t !計測終了、表示 !ケース２ !π/4=ATN(1)=(1)/1+(-1)/3+(1)/5+(-1)/7+(1)/9+(-1)/11+ … DEF Den(n)=(2n-1) !第n項の分母 LET t=TIME LET S=0 LET D=1 !分子（符号を含む） FOR i=1 TO 1000000 !収束が遅い LET S=S+D/Den(i) LET D=-D !反転 NEXT i PRINT S4 !S=π/4 PRINT TIME-t !分母の奇数をnの式（第n項）で書けない場合 LET t=TIME LET S=0 LET D=1 !符号 FOR i=1 TO 2000000 STEP 2 !収束が遅い LET S=S+D/i LET D=-D !反転 NEXT i PRINT S4 !S=π/4 PRINT TIME-t !ケース３ !π/4=ATN(1)=(1/1-1/3)+(1/5-1/7)+(1/9-1/11)+ … LET t=TIME LET S=0 LET i=1 DO WHILE i<=2000000 !収束が遅い LET S=S+1/i !カッコ内の前の値 LET i=i+2 LET S=S-1/i !後の値 LET i=i+2 LOOP PRINT S4 !S=π/4 PRINT TIME-t !精度　ケース２ LET t=TIME LET S=0 LET D=1 !符号 LET i=1 DO WHILE Den(i)<10000000 !分母が大きければ０に近くなる LET S=S+D/Den(i) LET D=-D !反転 LET i=i+1 LOOP PRINT "第";i-1;"項まで",S*4 !S=π/4 PRINT TIME-t END アイディアを出すなどはプログラミングの「楽しみ」で、メモリ、実行時間、読みやすさ、言語仕様などの制約があり、これは「辛さ」です。正解だけではなくいろいろ考えてみるとおもしろいです。

Re: 式をどう見るかでいろいろなプログラムがで...	返事を書く
SECOND 2006/12/12 21:35:30
「辛さ」ですか！？わたしは老眼の辛さだけですが・・もっと書いてください。おかげで、十進BASIC の魅力の沢山のことが、覚えられました。有難うございます。

Re: 式をどう見るかでいろいろなプログラムがで...	返事を書く
SECOND 2006/12/13 01:51:19
Microsoft Access の OLE も加えて頂けますか。全て重宝しております。

Re: 円周率	返事を書く
河川屋 2006/12/23 19:50:10
Re:円周率 >ライプニッツの公式による円周率を近似するのに >第ｎ項まで計算して近似値を与える関数を作るには以下を回答例とします。要点は以下のとおり。 1.関数(Function)で定義する。 2.第ｎ項までを足す。nはメインルーチンから引数で与える。 ※4(1-1/3+1/5) まで計算した場合、第3項までと数えている。ここまでが質問に対する最低必要な要件。以下、オマケの要点。 3.データ型の範囲内で、なるべく正確な答を返す。このためには、n項から遡って足していく。 ※正確な値とは、有効桁無限大の計算機でn項まで足した値と定義する。 ※※2進モード(IEEE754倍精度)と10進15桁モードを比較すると数学的にたまたま整数解となる場合は10進15桁が正確であるが、循環小数が混じると、 4則演算の場合2進モードのほうが正確で、超越関数を使うと10進15桁モードのほうが正確となるもよう。 ※※※第1兆項あたりまで計算しない限り、あまり意味はないかも。 4.第n項めの値を計算するには、(-1)^(n-1)/(2n-1)でも数学では正しいけれど、一般にべき計算は遅いのでif文で代用した。 (整数演算は組込関数内部で別ルーチンを使っていればこの限りではないが、念のためこうしている。) DECLARE EXTERNAL FUNCTION LEIBNIZ LET N=2 !第n項まで LET PAI=LEIBNIZ(N) PRINT PAI END EXTERNAL FUNCTION LEIBNIZ(N) LET T=0 FOR I=N TO 1 STEP -1 IF mod(I,2)=1 THEN LET T=T+1/(2I-1) ELSE LET T=T-1/(2I-1) END IF NEXT I LET LEIBNIZ=4*T END FUNCTION