多倍長数のＬＯＧの計算

　　多倍長数のＬＯＧの計算 会社員 2007/04/30 10:54:29  (修正1回)
　　├!続き１ 会社員 2007/04/30 10:55:44 
　　├!続き２ 会社員 2007/04/30 10:58:37 
　　├!続き３ 会社員 2007/04/30 10:59:42 
　　├一般論ですが，有理数モードで計算して誤差... 白石　和夫 2007/04/30 14:02:24 
　　├Re:多倍長のＬＯＧ 河川屋 2007/05/03 22:34:03 
　　├バグの箇所がわかりました。 会社員 2007/05/04 23:34:48 
　　│├DELETED  会社員  2007/05/05 00:05:41  (削除)
　　│├DELETED  会社員  2007/05/05 00:09:15  (削除)
　　│├DELETED  会社員  2007/05/05 15:33:09  (削除)
　　│├DELETED  会社員  2007/05/05 15:34:14  (削除)
　　│└小数点以下の場合に誤動作していたところを... 会社員 2007/05/16 23:03:53 
　　│　├!続き 会社員 2007/05/16 23:04:58 
　　│　└!続き２ 会社員 2007/05/16 23:05:47 
　　└LOG（π）の結果です。 会社員 2007/05/16 23:10:37

　　多倍長数のＬＯＧの計算 会社員 2007/04/30 10:54:29  (修正1回) ツリーへ

多倍長数のＬＯＧの計算 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/04/30 10:54:29 ** この記事は1回修正されてます

多倍長数のＬＯＧのプログラムを作成中です。

８桁あたりまで正しくでてきているのですが、
それ以降が正しくでてこないようです。

もし、バグの理由がわかれば教えてください。

下記は、log５　を　１００桁求める場合の試験プログラムです。

! 多倍長のLOG計算 10桁区切り
!
! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}
!
!

OPTION ARITHMETIC DECIMAL_HIGH ! 10進 1000桁
OPTION BASE 0

! ＝＝＝＝＝＝＝＝＝初期設定テスト用初期値 100桁＝＝＝＝＝＝＝＝＝＝
!DO
! INPUT PROMPT "いくらの桁で計算しますか？":KTA_L

LET KTA_L=100 ! テスト用初期値 100桁
print "桁数=";KTA_L;"桁"
! LOOP UNTIL MOD(KTA_L,20)=0 AND KTA_L>=100
LET HYRETU=ROUND(KTA_L*1.05/10,0) ! 5％のマージンをつけおく
! ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

! 10桁区切り
PRINT "HYRETU=";HYRETU

DIM N(HYRETU),N_1(HYRETU),INV_NP1ANS(HYRETU),KAKE_ANS(HYRETU)
DIM NM1_NP1(HYRETU),NM1_NP1_2(HYRETU)
DIM LG_N(HYRETU),TLR_K(HYRETU)

! ＝＝＝＝＝＝＝＝＝初期設定テスト用初期値 LOG(5)＝＝＝＝＝＝＝＝＝＝
LET N(0)=5 !
FOR I=1 TO HYRETU
LET N(I)=0
NEXT I
! ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}

　　├!続き１ 会社員 2007/04/30 10:55:44  ツリーへ

Re: 多倍長数のＬＯＧの計算 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/04/30 10:55:44

!　続き１

LET N_1(0)=N(0)-1 !(N-1)

CALL INV_N(N,UBOUND(N),INV_NP1ANS) !１／(N+1) INV_NP1()が帰ってくる

FOR I=0 TO HYRETU
PRINT "N=";N(I)
NEXT I

FOR I=0 TO HYRETU
PRINT "N-1=";N_1(I)
NEXT I

FOR I=0 TO HYRETU
PRINT "1/(N+1)=";INV_NP1ANS(I)
NEXT I

CALL KAKEZAN(N_1,INV_NP1ANS,UBOUND(N_1),KAKE_ANS) ! (N-1)／(N+1)を計算 KAKE_ANS()が帰ってくる

FOR I=0 TO HYRETU
PRINT "(N-1)/(N+1)=";KAKE_ANS(I)
LET NM1_NP1(I)=KAKE_ANS(I)
NEXT I

CALL KAKEZAN(NM1_NP1,NM1_NP1,UBOUND(NM1_NP1),KAKE_ANS) ! (N-1)^2／(N+1)^2 答えが、KAKE_ANS()で帰ってくる

FOR I=0 TO HYRETU
LET NM1_NP1_2(I)=KAKE_ANS(I)
NEXT I

PRINT "(N-1)/(N+1)の結果"
FOR I=0 TO HYRETU
PRINT "NM1_NP1(";I;")=";NM1_NP1(I)
NEXT I

PRINT "(N-1)^2/(N+1)^2の結果"
FOR I=0 TO HYRETU
PRINT "NM1_NP1_2(";I;")=";NM1_NP1_2(I)
NEXT I

! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}

FOR S=1 TO 100000000 STEP 2

IF S=1 THEN
FOR I=0 TO HYRETU
LET TLR_K(I)=NM1_NP1(I) ! (N-1)/(N+1)を設定
NEXT I
END IF

! テーラー項の計算テーラー項の計算テーラー項の計算テーラー項の計算テーラー項の計算
! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}

IF S=>3 THEN
CALL KAKEZAN(TLR_K,NM1_NP1_2,UBOUND(TLR_K),KAKE_ANS) ! (N-1)/(N+1)^3/3 答えが、KKZ_ANS()で帰ってくる

FOR I=0 TO HYRETU
LET TLR_K(I)=KAKE_ANS(I)
NEXT I

FOR I=0 TO HYRETU
LET TLR_K(I)=TLR_K(I)*(S-2) ! S-2 倍する
NEXT I
CALL KTA_L_SYORI(TLR_K,UBOUND(TLR_K)) !桁処理
END IF

FOR I=0 TO HYRETU
IF I=<HYRETU-1 THEN
LET TLR_K(I)=INT(TLR_K(I)/S) ! Sで割る
LET TLR_K(I+1)=MOD(TLR_K(I),S)*10^10+TLR_K(I+1) ! 余りを次の桁へ渡す
END IF
IF I=HYRETU THEN
LET TLR_K(I)=INT(TLR_K(I)/S) ! Sで割る
END IF
NEXT I

CALL KTA_L_SYORI(TLR_K,UBOUND(TLR_K)) !桁処理

PRINT "＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝"
PRINT "S=";S

FOR I=0 TO HYRETU
PRINT "TLR_K(";I;")=";TLR_K(I)
NEXT I

FOR I=0 TO HYRETU
LET LG_N(I)=LG_N(I)+TLR_K(I)
NEXT I

CALL KTA_L_SYORI(LG_N,UBOUND(LG_N)) !桁処理

! PRINT "桁処理後"
! FOR I=0 TO HYRETU
! PRINT "LG_N(";I;")=";LG_N(I)
! NEXT I

IF TLR_K(HYRETU/1.05+1)=<10^9 AND TLR_K(HYRETU/1.05)=0 THEN
PRINT "TLR_K(HYRETU/1.05+1)<10^9になりました"
EXIT FOR
END IF

NEXT S

FOR I=0 TO HYRETU
LET LG_N(I)=2*LG_N(I)
NEXT I
CALL KTA_L_SYORI(LG_N,UBOUND(LG_N)) !桁処理

! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}

PRINT "---計算終了！--- LOGの結果"
FOR I=0 TO HYRETU
PRINT "LG_N(";I;")=";LG_N(I)
NEXT I

　　├!続き２ 会社員 2007/04/30 10:58:37  ツリーへ

Re: 多倍長数のＬＯＧの計算 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/04/30 10:58:37

!続き２
!●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●●
SUB INV_N(A(),LLLL,INV_NP1()) !１/(N+1)の作成
DIM Xn1(LLLL),Xn0(LLLL),KKZ_ANS(LLLL),Xn0_2(LLLL),NP1(LLLL)

! １/(N+1)の作成
FOR I=0 TO LLLL
IF I=0 THEN LET NP1(0)=A(0)+1
IF I<>0 THEN LET NP1(I)=A(I)
NEXT I

CALL KTA_L_SYORI(NP1,UBOUND(NP1))

FOR I=0 TO LLLL
PRINT "NP1(I)=";NP1(I)
NEXT I

print "逆数スタート";time$

! Xn+1 = 2 * Xn - A * Xn^2
LET Xn0(0)=INT(1/NP1(0)) ! 初期化
LET Xn0(1)=INT(1/NP1(0)*10^10)-INT(1/NP1(0))*10^10 ! 初期化

FOR I=2 TO LLLL ! ニュートン法パラメータ初期化
LET Xn0(I)=0
LET Xn1(I)=0
NEXT I

LET CNT=0
LET KTA_L_T=KTA_L
DO
LET KTA_L_T=KTA_L_T/10
LET CNT=CNT+1
LOOP UNTIL KTA_L_T<1
LET LOOOP=CNT/0.3010+1
PRINT "ループ回数=";LOOOP;"回"

FOR JJJJ=1 TO LOOOP
PRINT "■■■■JJJJ=";JJJJ

FOR I=0 TO HYRETU
LET Xn1(I)=0
NEXT I

PRINT "Xn^2の計算開始！！-----------------------------------------"
CALL KAKEZAN(Xn0,Xn0,UBOUND(Xn0),KKZ_ANS) !乗算サブルーチン Xn^2

FOR I=0 TO HYRETU !
LET Xn0_2(I)=KKZ_ANS(I)
NEXT I

PRINT "A * Xn^2の計算開始！！-----------------------------------------"
CALL KAKEZAN(NP1,Xn0_2,UBOUND(A),KKZ_ANS) !乗算サブルーチン A * Xn^2

PRINT "Xn1(I)=2*Xn0(I)-KKZ_ANS(I) 計算開始！！"
FOR I=0 TO HYRETU !
LET Xn1(I)=2*Xn0(I)-KKZ_ANS(I)
NEXT I

CALL KTA_L_SYORI(Xn1,UBOUND(Xn1))

PRINT "Xn1(I)=2*Xn0(I)-KKZ_ANS(I) 計算終了"

FOR I=0 TO LLLL !
LET Xn0(I)=Xn1(I)
NEXT I

PRINT "================================================================="
PRINT "================================================================="

FOR I=0 TO LLLL !
PRINT "■Xn1(";I;")=";Xn1(I)
NEXT I
NEXT JJJJ

　　├!続き３ 会社員 2007/04/30 10:59:42  ツリーへ

Re: 多倍長数のＬＯＧの計算 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/04/30 10:59:42

!続き３

FOR I=0 TO LLLL
IF I=0 THEN LET INV_NP1(0)=Xn1(I)
IF I<>0 THEN LET INV_NP1(I)=Xn1(I)
NEXT I

END SUB
!●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●●

! 掛算サブルーチン始まり■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
SUB KAKEZAN(JUX(),JUY(),LLL,KKZ_ANS()) ! Z=X*Y

FOR I=0 TO LLL
LET KKZ_ANS(I)=0
NEXT I

! X * Yの計算
FOR I=0 TO LLL
FOR J=0 TO LLL-I
LET KKZ_ANS(I+J)=JUX(I)*JUY(J)+KKZ_ANS(I+J)
NEXT J
NEXT I

CALL KTA_L_SYORI(KKZ_ANS,LLL)

END SUB

SUB KTA_L_SYORI(BBB(),LLL) !桁処理

DO ! 桁処理１ --負記号の削除
LET MNSF=0
IF BBB(1)<0 THEN
LET BBB(0)=BBB(0)+INT(BBB(1)/10^10)-1
LET BBB(1)=BBB(1)-10^10*INT(BBB(1)/10^10)+10^10
END IF
FOR I=1 TO LLL
IF BBB(I)<0 THEN
LET MNSF=1
LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10)-1
LET BBB(I)=BBB(I)-10^10*INT(BBB(I)/10^10)+10^10
END IF
NEXT I
LOOP UNTIL MNSF=0

DO ! 桁処理２ --桁上げ処理
LET OBF=0
FOR I=LLL TO 1 STEP -1
LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10)
LET BBB(I)=BBB(I)-10^10*INT(BBB(I)/10^10)
IF BBB(I)>10^10 THEN
LET OBF=1
END IF
NEXT I
LOOP UNTIL OBF=0
END SUB
! 掛算サブルーチン終わり■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■

END

　　├一般論ですが，有理数モードで計算して誤差... 白石　和夫 2007/04/30 14:02:24  ツリーへ

Re: 多倍長数のＬＯＧの計算 返事を書くノートメニュー

白石　和夫 <fbdfvqwhki> 2007/04/30 14:02:24

一般論ですが，有理数モードで計算して誤差があるようなら近似式の見直し，そうでないなら，異符号数の加算あるいは同符号数の減算などによる桁落ちの可能性が高いと思います。
有理数モードでの計算結果はPRINT USINGを利用して1000桁程度までは小数に直せます。
http://hp.vector.co.jp/authors/VA008683/Aprox.htm

　　├Re:多倍長のＬＯＧ 河川屋 2007/05/03 22:34:03  ツリーへ

Re: 多倍長数のＬＯＧの計算 返事を書くノートメニュー

河川屋 <evcdumfptp> 2007/05/03 22:34:03

Re:多倍長のＬＯＧ

＞８桁あたりまで正しくでてきているのですが、
＞それ以降が正しくでてこないようです。
＞もし、バグの理由がわかれば教えてください。
う～ん、目の前にあるのは自作の多倍長ルーチンのバグですが、
それがなくなったとしても、まだダメです。
というのは、LOG(1000)あたりでさえ、級数の収束が非常に遅くなってしまうという致命傷が生じているため。
あと、1000桁モードを使用していながら、
そのうち10桁しか使っていないというのは、大いなるムダでしょう。
（メモリーも、計算量も。）
1000桁モードなら、配列1つあたり400桁でもいけるでしょう。
ただし、これらは、1000桁以上の答がほしい場合であり、
1000桁弱の答でいいなら、配列を使わないで生で計算してしまえばいいです。
たぶん、以下で正しい範囲は995～997桁くらいだと思います。

! 1000桁モードLOG計算
OPTION ARITHMETIC DECIMAL_HIGH

LET X=5 !LOG(X)を求める
LET KETA=1000 !求めたい桁数
LET EPS=1/10^KETA/10 !収束判定値
CALL EE(E) !自然対数の底

! 級数の収束を加速するためLOG(X)整数部の概略値を計算
LET GETA=0
DO WHILE X>1.65 !1.65はSQRT(E)の近似値。
LET GETA=GETA+1
LET X=X/E
LOOP
DO WHILE X<0.6 !0.6は1/SQRT(E)の近似値
LET GETA=GETA-1
LET X=X*E
LOOP

! LOG(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}
LET XX=(X-1)/(X+1)
LET XX2=XX
LET LG=0
LET LG0=XX2
LET I=1
DO WHILE ABS(LG0)>EPS ! このループを回る回数は、1000桁計算時最悪で800回強。
LET LG=LG+LG0
LET XX2=XX2*XX*XX
LET I=I+2
LET LG0=XX2/I
LOOP
LET LG=GETA+2*LG
PRINT LG
STOP

SUB EE(E) !自然対数の底の1000桁計算
DIM F(450)
! E=1+1/1!+1/2!+1/3!+.....を、桁落ち対策のため尻から足す
LET E=0
LET F(1)=1
FOR I=2 TO 450
LET F(I)=F(I-1)*I
NEXT I
FOR I=450 TO 1 STEP -1
LET E=E+1/F(I)
NEXT I
LET E=E+1
END SUB
END

　　├バグの箇所がわかりました。 会社員 2007/05/04 23:34:48  ツリーへ

Re: 多倍長数のＬＯＧの計算 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/05/04 23:34:48

バグの箇所がわかりました。
テーラー項を求める際の整数で割る手順に
誤りが見つかりました。

　　│├DELETED  会社員  2007/05/05 00:05:41  (削除) ツリーへ

Re: バグの箇所がわかりました。 返事を書く

会社員 <qiekrrepwd> 2007/05/05 00:05:41 ** この記事は削除されました

　　│├DELETED  会社員  2007/05/05 00:09:15  (削除) ツリーへ

Re: バグの箇所がわかりました。 返事を書く

会社員 <qiekrrepwd> 2007/05/05 00:09:15 ** この記事は削除されました

　　│├DELETED  会社員  2007/05/05 15:33:09  (削除) ツリーへ

Re: バグの箇所がわかりました。 返事を書く

会社員 <qiekrrepwd> 2007/05/05 15:33:09 ** この記事は削除されました

　　│├DELETED  会社員  2007/05/05 15:34:14  (削除) ツリーへ

Re: バグの箇所がわかりました。 返事を書く

会社員 <qiekrrepwd> 2007/05/05 15:34:14 ** この記事は削除されました

　　│└小数点以下の場合に誤動作していたところを... 会社員 2007/05/16 23:03:53  ツリーへ

Re: バグの箇所がわかりました。 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/05/16 23:03:53

小数点以下の場合に誤動作していたところを修正です

!多倍長のLOG計算 10桁区切り
!log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}

OPTION ARITHMETIC DECIMAL_HIGH ! 10進 1000桁
OPTION BASE 0

! ＝＝＝＝＝＝＝＝＝初期設定＝＝＝＝＝＝＝＝＝＝＝＝
DO
INPUT PROMPT "いくらの桁まで計算しますか？":KTA
print "桁数=";KTA;"桁"
LOOP UNTIL MOD(KTA,20)=0 AND KTA>=100
LET HYRE2=ROUND(KTA*1.15/10,0) ! 15％のマージンをつけおく
! 10桁区切り
PRINT "HYRE2=";HYRE2

DIM N(HYRE2)

! テスト用の初期値設定
PRINT "試験用 LOG(5)の計算"
LET N(0)=5! テスト用初期値LOG(5)
FOR I=1 TO HYRE2
LET N(I)=0
NEXT I
! ＝＝＝＝＝＝＝＝＝初期設定＝＝＝＝＝＝＝＝＝＝＝＝

DIM N_1(HYRE2),INV_NP1ANS(HYRE2),KAKE_ANS(HYRE2)
DIM NM1_NP1(HYRE2),NM1_NP1_2(HYRE2)
DIM LG_N(HYRE2),TYL_K(HYRE2)

! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}
PRINT "計算開始 ";DATE$;"　";TIME$
CALL TYLR_PRE ! 事前準備

FOR S=1 TO 100000000 STEP 2
IF S=1 THEN
FOR I=0 TO HYRE2
LET TYL_K(I)=NM1_NP1(I) ! 初項に (N-1)/(N+1)を設定
NEXT I
END IF

IF S=>3 THEN
CALL KAKEZAN(TYL_K,NM1_NP1_2,UBOUND(TYL_K),KAKE_ANS) ! (N-1)/(N+1)^3/3 答えが、KKZ_ANS()で帰ってくる
FOR I=0 TO HYRE2
LET TYL_K(I)=KAKE_ANS(I)*(S-2) ! (N-1)/(N+1)^2*(S-2) 倍する
NEXT I
CALL KTA_SYORI(TYL_K,UBOUND(TYL_K)) !桁処理
END IF

FOR I=0 TO HYRE2
IF I=<HYRE2-1 THEN
LET AMARI=MOD(TYL_K(I),S)
LET TYL_K(I)=INT(TYL_K(I)/S) ! Sで割る
LET TYL_K(I+1)=AMARI*10^10+TYL_K(I+1) ! 余りを次の桁へ渡す
END IF

IF I=HYRE2 THEN
LET TYL_K(I)=INT(TYL_K(I)/S) ! Sで割る
END IF
NEXT I
CALL KTA_SYORI(TYL_K,UBOUND(TYL_K)) !桁処理

FOR I=0 TO HYRE2
LET LG_N(I)=LG_N(I)+TYL_K(I)
NEXT I

CALL KTA_SYORI(LG_N,UBOUND(LG_N)) !桁処理

IF TYL_K(HYRE2)=<10^9 AND TYL_K(HYRE2-1)=0 THEN
PRINT "終了！ ";DATE$;"　";TIME$
EXIT FOR
END IF

NEXT S

　　│　├!続き 会社員 2007/05/16 23:04:58  ツリーへ

Re: 小数点以下の場合に誤動作していたところを... 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/05/16 23:04:58

!続き
FOR I=0 TO HYRE2
LET LG_N(I)=2*LG_N(I)
NEXT I
CALL KTA_SYORI(LG_N,UBOUND(LG_N)) !桁処理

PRINT ""
FOR I=0 TO HYRE2
IF I=0 THEN
PRINT STR$(LG_N(I));"."
END IF
IF I>0 THEN
PRINT REPEAT$("0",10-LEN(STR$(LG_N(I)))) ;STR$(LG_N(I));" ";
END IF

IF I>0 AND I/5=INT(I/5) THEN PRINT "：";STR$(I*10)
NEXT I

!==========サブルーチン============= 事前準備
SUB TYLR_PRE

FOR I=0 TO HYR2
IF I=0 THEN LET N_1(0)=N(I)-1 !(N-1)
IF i>0 THEN LET N_1(I)=N(I)
NEXT I

CALL INV_N(N,UBOUND(N),INV_NP1ANS) !１／(N+1)INV_NP1()が帰ってくる
CALL KAKEZAN(N_1,INV_NP1ANS,UBOUND(N_1),KAKE_ANS) ! (N-1)／(N+1)を計算 KAKE_ANS()が帰ってくる

FOR I=0 TO HYRE2
LET NM1_NP1(I)=KAKE_ANS(I)
NEXT I

CALL KAKEZAN(NM1_NP1,NM1_NP1,UBOUND(NM1_NP1),KAKE_ANS) ! (N-1)／(N+1)^2を計算 KAKE_ANS()が帰ってくる

FOR I=0 TO HYRE2
LET NM1_NP1_2(I)=KAKE_ANS(I)
NEXT I

end SUB

!●●●●1/(N+1)のルーチン●●●●
SUB INV_N(A(),LLLL,INV_NP1()) !１/(N+1)の作成
DIM Xn1(LLLL),Xn0(LLLL),KKZ_ANS(LLLL),Xn0_2(LLLL),NP1(LLLL)

! １/(N+1)の作成
FOR I=0 TO LLLL
IF I=0 THEN LET NP1(0)=A(0)+1
IF I<>0 THEN LET NP1(I)=A(I)
NEXT I
CALL KTA_SYORI(NP1,UBOUND(NP1))

! Xn+1 = 2 * Xn - A * Xn^2
LET Xn0(0)=INT(1/NP1(0)) ! 初期化
LET Xn0(1)=INT(1/NP1(0)*10^10)-INT(1/NP1(0))*10^10 ! 初期化

FOR I=2 TO LLLL ! ニュートン法パラメータ初期化
LET Xn0(I)=0
LET Xn1(I)=0
NEXT I

LET CNT=0
LET KTA_T=KTA
DO
LET KTA_T=KTA_T/10
LET CNT=CNT+1
LOOP UNTIL KTA_T<1
LET LOOOP=CNT/0.3010+2
PRINT "ループ回数=";ROUND(LOOOP,2);"回"

FOR JJJJ=1 TO LOOOP
FOR I=0 TO HYRE2
LET Xn1(I)=0
NEXT I
CALL KAKEZAN(Xn0,Xn0,UBOUND(Xn0),KKZ_ANS) !乗算サブルーチン Xn^2

FOR I=0 TO HYRE2
LET Xn0_2(I)=KKZ_ANS(I)
NEXT I

　　│　└!続き２ 会社員 2007/05/16 23:05:47  ツリーへ

Re: 小数点以下の場合に誤動作していたところを... 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/05/16 23:05:47

!続き２
! A * Xn^2の計算開始！
CALL KAKEZAN(NP1,Xn0_2,UBOUND(A),KKZ_ANS) !乗算サブルーチン A * Xn^2

! Xn1(I)=2*Xn0(I)-KKZ_ANS(I) 計算開始！
FOR I=0 TO HYRE2 !
LET Xn1(I)=2*Xn0(I)-KKZ_ANS(I)
NEXT I
CALL KTA_SYORI(Xn1,UBOUND(Xn1))

FOR I=0 TO LLLL !
LET Xn0(I)=Xn1(I)
NEXT I
NEXT JJJJ

FOR I=0 TO LLLL
IF I=0 THEN LET INV_NP1(0)=Xn1(I)
IF I<>0 THEN LET INV_NP1(I)=Xn1(I)
NEXT I
END SUB

! 掛算サブルーチン始まり■■■■
SUB KAKEZAN(JUX(),JUY(),LLL,KKZ_ANS())
FOR I=0 TO LLL
LET KKZ_ANS(I)=0
NEXT I
! X * Yの計算
FOR I=0 TO LLL
FOR J=0 TO LLL-I
LET KKZ_ANS(I+J)=JUX(I)*JUY(J)+KKZ_ANS(I+J)
NEXT J
NEXT I
CALL KTA_SYORI(KKZ_ANS,LLL)

END SUB

SUB KTA_SYORI(BBB(),LLL) !桁処理
DO ! 桁処理１ --負記号の削除
LET MNSF=0
IF BBB(1)<0 THEN
LET BBB(0)=BBB(0)+INT(BBB(1)/10^10)-1
LET BBB(1)=BBB(1)-10^10*INT(BBB(1)/10^10)+10^10
END IF
FOR I=1 TO LLL
IF BBB(I)<0 THEN
LET MNSF=1
LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10)-1
LET BBB(I)=BBB(I)-10^10*INT(BBB(I)/10^10)+10^10
END IF
NEXT I
LOOP UNTIL MNSF=0

DO ! 桁処理２ --桁上げ処理
LET OBF=0
FOR I=LLL TO 1 STEP -1
LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10)
LET BBB(I)=BBB(I)-10^10*INT(BBB(I)/10^10)
IF BBB(I)>10^10 THEN
LET OBF=1
END IF
NEXT I
LOOP UNTIL OBF=0
END SUB

END

　　└LOG（π）の結果です。 会社員 2007/05/16 23:10:37  ツリーへ

Re: 多倍長数のＬＯＧの計算 返事を書くノートメニュー

会社員 <qiekrrepwd> 2007/05/16 23:10:37

LOG（π）　の結果です。

1.
1447298858 4940017414 3427351353 0587116472 9481291531 ：50
1571513623 0714721377 6988482607 9783623270 2754897077 ：100
0200981222 8697989159 0482055279 2345658727 9081078810 ：150
2868252763 9391426634 5902902484 7733588699 3778920311 ：200
9630824756 7940119160 2821722737 9888126563 1780498236 ：250
9731331069 5003600064 4054872638 8022327009 6433504959 ：300
5118150662 3725246834 3391269896 5797514047 7703857799 ：350
5399825842 5660228485 0148136217 9159252505 6707638686 ：400
0280763456 8897505123 3436078143 9914144264 2959671289 ：450
7781136526 4523450410 5900716081 8570824981 1881831868 ：500
9767284592 8110257656 8751724223 3833718927 3043288217 ：550
3486510427 6153237516 1028392221 3401436967 1758561644 ：600
2473718780 5060466920 5628337731 0133621627 4515898752 ：650
0151299654 5465739691 5282523916 9585245379 3594601400 ：700
3799565196 6603653800 0112659858 5001297656 9906074466 ：750
7455472671 0450849506 6855874339 0774251341 5924126523 ：800
1777178491 7799588095 7678805102 9644475090 1508911403 ：850
2780807683 3733793894 9488075152 8900918753 6376608670 ：900
7435833345 1081392325 3557406768 4327431198 0496339997 ：950
6180304622 1286361595 8598364047 5800986179 9938264629 ：1000

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

多倍長数のＬＯＧの計算	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/04/30 10:54:29 ** この記事は1回修正されてます
多倍長数のＬＯＧのプログラムを作成中です。８桁あたりまで正しくでてきているのですが、それ以降が正しくでてこないようです。もし、バグの理由がわかれば教えてください。下記は、log５　を　１００桁求める場合の試験プログラムです。 ! 多倍長のLOG計算 10桁区切り ! ! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・} ! ! OPTION ARITHMETIC DECIMAL_HIGH ! 10進 1000桁 OPTION BASE 0 ! ＝＝＝＝＝＝＝＝＝初期設定テスト用初期値 100桁＝＝＝＝＝＝＝＝＝＝ !DO ! INPUT PROMPT "いくらの桁で計算しますか？":KTA_L LET KTA_L=100 ! テスト用初期値 100桁 print "桁数=";KTA_L;"桁" ! LOOP UNTIL MOD(KTA_L,20)=0 AND KTA_L>=100 LET HYRETU=ROUND(KTA_L*1.05/10,0) ! 5％のマージンをつけおく ! ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ! 10桁区切り PRINT "HYRETU=";HYRETU DIM N(HYRETU),N_1(HYRETU),INV_NP1ANS(HYRETU),KAKE_ANS(HYRETU) DIM NM1_NP1(HYRETU),NM1_NP1_2(HYRETU) DIM LG_N(HYRETU),TLR_K(HYRETU) ! ＝＝＝＝＝＝＝＝＝初期設定テスト用初期値 LOG(5)＝＝＝＝＝＝＝＝＝＝ LET N(0)=5 ! FOR I=1 TO HYRETU LET N(I)=0 NEXT I ! ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・}

Re: 多倍長数のＬＯＧの計算	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/04/30 10:55:44
!　続き１ LET N_1(0)=N(0)-1 !(N-1) CALL INV_N(N,UBOUND(N),INV_NP1ANS) !１／(N+1) INV_NP1()が帰ってくる FOR I=0 TO HYRETU PRINT "N=";N(I) NEXT I FOR I=0 TO HYRETU PRINT "N-1=";N_1(I) NEXT I FOR I=0 TO HYRETU PRINT "1/(N+1)=";INV_NP1ANS(I) NEXT I CALL KAKEZAN(N_1,INV_NP1ANS,UBOUND(N_1),KAKE_ANS) ! (N-1)／(N+1)を計算 KAKE_ANS()が帰ってくる FOR I=0 TO HYRETU PRINT "(N-1)/(N+1)=";KAKE_ANS(I) LET NM1_NP1(I)=KAKE_ANS(I) NEXT I CALL KAKEZAN(NM1_NP1,NM1_NP1,UBOUND(NM1_NP1),KAKE_ANS) ! (N-1)^2／(N+1)^2 答えが、KAKE_ANS()で帰ってくる FOR I=0 TO HYRETU LET NM1_NP1_2(I)=KAKE_ANS(I) NEXT I PRINT "(N-1)/(N+1)の結果" FOR I=0 TO HYRETU PRINT "NM1_NP1(";I;")=";NM1_NP1(I) NEXT I PRINT "(N-1)^2/(N+1)^2の結果" FOR I=0 TO HYRETU PRINT "NM1_NP1_2(";I;")=";NM1_NP1_2(I) NEXT I ! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・} FOR S=1 TO 100000000 STEP 2 IF S=1 THEN FOR I=0 TO HYRETU LET TLR_K(I)=NM1_NP1(I) ! (N-1)/(N+1)を設定 NEXT I END IF ! テーラー項の計算テーラー項の計算テーラー項の計算テーラー項の計算テーラー項の計算 ! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・} IF S=>3 THEN CALL KAKEZAN(TLR_K,NM1_NP1_2,UBOUND(TLR_K),KAKE_ANS) ! (N-1)/(N+1)^3/3 答えが、KKZ_ANS()で帰ってくる FOR I=0 TO HYRETU LET TLR_K(I)=KAKE_ANS(I) NEXT I FOR I=0 TO HYRETU LET TLR_K(I)=TLR_K(I)(S-2) ! S-2 倍する NEXT I CALL KTA_L_SYORI(TLR_K,UBOUND(TLR_K)) !桁処理 END IF FOR I=0 TO HYRETU IF I=<HYRETU-1 THEN LET TLR_K(I)=INT(TLR_K(I)/S) ! Sで割る LET TLR_K(I+1)=MOD(TLR_K(I),S)10^10+TLR_K(I+1) ! 余りを次の桁へ渡す END IF IF I=HYRETU THEN LET TLR_K(I)=INT(TLR_K(I)/S) ! Sで割る END IF NEXT I CALL KTA_L_SYORI(TLR_K,UBOUND(TLR_K)) !桁処理 PRINT "＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝" PRINT "S=";S FOR I=0 TO HYRETU PRINT "TLR_K(";I;")=";TLR_K(I) NEXT I FOR I=0 TO HYRETU LET LG_N(I)=LG_N(I)+TLR_K(I) NEXT I CALL KTA_L_SYORI(LG_N,UBOUND(LG_N)) !桁処理 ! PRINT "桁処理後" ! FOR I=0 TO HYRETU ! PRINT "LG_N(";I;")=";LG_N(I) ! NEXT I IF TLR_K(HYRETU/1.05+1)=<10^9 AND TLR_K(HYRETU/1.05)=0 THEN PRINT "TLR_K(HYRETU/1.05+1)<10^9になりました" EXIT FOR END IF NEXT S FOR I=0 TO HYRETU LET LG_N(I)=2*LG_N(I) NEXT I CALL KTA_L_SYORI(LG_N,UBOUND(LG_N)) !桁処理 ! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・} PRINT "---計算終了！--- LOGの結果" FOR I=0 TO HYRETU PRINT "LG_N(";I;")=";LG_N(I) NEXT I

Re: 多倍長数のＬＯＧの計算	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/04/30 10:58:37
!続き２ !●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●● SUB INV_N(A(),LLLL,INV_NP1()) !１/(N+1)の作成 DIM Xn1(LLLL),Xn0(LLLL),KKZ_ANS(LLLL),Xn0_2(LLLL),NP1(LLLL) ! １/(N+1)の作成 FOR I=0 TO LLLL IF I=0 THEN LET NP1(0)=A(0)+1 IF I<>0 THEN LET NP1(I)=A(I) NEXT I CALL KTA_L_SYORI(NP1,UBOUND(NP1)) FOR I=0 TO LLLL PRINT "NP1(I)=";NP1(I) NEXT I print "逆数スタート";time$ ! Xn+1 = 2 * Xn - A * Xn^2 LET Xn0(0)=INT(1/NP1(0)) ! 初期化 LET Xn0(1)=INT(1/NP1(0)10^10)-INT(1/NP1(0))10^10 ! 初期化 FOR I=2 TO LLLL ! ニュートン法パラメータ初期化 LET Xn0(I)=0 LET Xn1(I)=0 NEXT I LET CNT=0 LET KTA_L_T=KTA_L DO LET KTA_L_T=KTA_L_T/10 LET CNT=CNT+1 LOOP UNTIL KTA_L_T<1 LET LOOOP=CNT/0.3010+1 PRINT "ループ回数=";LOOOP;"回" FOR JJJJ=1 TO LOOOP PRINT "■■■■JJJJ=";JJJJ FOR I=0 TO HYRETU LET Xn1(I)=0 NEXT I PRINT "Xn^2の計算開始！！-----------------------------------------" CALL KAKEZAN(Xn0,Xn0,UBOUND(Xn0),KKZ_ANS) !乗算サブルーチン Xn^2 FOR I=0 TO HYRETU ! LET Xn0_2(I)=KKZ_ANS(I) NEXT I PRINT "A * Xn^2の計算開始！！-----------------------------------------" CALL KAKEZAN(NP1,Xn0_2,UBOUND(A),KKZ_ANS) !乗算サブルーチン A * Xn^2 PRINT "Xn1(I)=2Xn0(I)-KKZ_ANS(I) 計算開始！！" FOR I=0 TO HYRETU ! LET Xn1(I)=2Xn0(I)-KKZ_ANS(I) NEXT I CALL KTA_L_SYORI(Xn1,UBOUND(Xn1)) PRINT "Xn1(I)=2*Xn0(I)-KKZ_ANS(I) 計算終了" FOR I=0 TO LLLL ! LET Xn0(I)=Xn1(I) NEXT I PRINT "=================================================================" PRINT "=================================================================" FOR I=0 TO LLLL ! PRINT "■Xn1(";I;")=";Xn1(I) NEXT I NEXT JJJJ

Re: 多倍長数のＬＯＧの計算	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/04/30 10:59:42
!続き３ FOR I=0 TO LLLL IF I=0 THEN LET INV_NP1(0)=Xn1(I) IF I<>0 THEN LET INV_NP1(I)=Xn1(I) NEXT I END SUB !●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●1/(N+1)のルーチン●●●●● ! 掛算サブルーチン始まり■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ SUB KAKEZAN(JUX(),JUY(),LLL,KKZ_ANS()) ! Z=XY FOR I=0 TO LLL LET KKZ_ANS(I)=0 NEXT I ! X Yの計算 FOR I=0 TO LLL FOR J=0 TO LLL-I LET KKZ_ANS(I+J)=JUX(I)JUY(J)+KKZ_ANS(I+J) NEXT J NEXT I CALL KTA_L_SYORI(KKZ_ANS,LLL) END SUB SUB KTA_L_SYORI(BBB(),LLL) !桁処理 DO ! 桁処理１ --負記号の削除 LET MNSF=0 IF BBB(1)<0 THEN LET BBB(0)=BBB(0)+INT(BBB(1)/10^10)-1 LET BBB(1)=BBB(1)-10^10INT(BBB(1)/10^10)+10^10 END IF FOR I=1 TO LLL IF BBB(I)<0 THEN LET MNSF=1 LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10)-1 LET BBB(I)=BBB(I)-10^10INT(BBB(I)/10^10)+10^10 END IF NEXT I LOOP UNTIL MNSF=0 DO ! 桁処理２ --桁上げ処理 LET OBF=0 FOR I=LLL TO 1 STEP -1 LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10) LET BBB(I)=BBB(I)-10^10INT(BBB(I)/10^10) IF BBB(I)>10^10 THEN LET OBF=1 END IF NEXT I LOOP UNTIL OBF=0 END SUB ! 掛算サブルーチン終わり■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ END

Re: 多倍長数のＬＯＧの計算	返事を書くノートメニュー
白石　和夫 <fbdfvqwhki> 2007/04/30 14:02:24
一般論ですが，有理数モードで計算して誤差があるようなら近似式の見直し，そうでないなら，異符号数の加算あるいは同符号数の減算などによる桁落ちの可能性が高いと思います。有理数モードでの計算結果はPRINT USINGを利用して1000桁程度までは小数に直せます。 http://hp.vector.co.jp/authors/VA008683/Aprox.htm

Re: 多倍長数のＬＯＧの計算	返事を書くノートメニュー
河川屋 <evcdumfptp> 2007/05/03 22:34:03
Re:多倍長のＬＯＧ＞８桁あたりまで正しくでてきているのですが、＞それ以降が正しくでてこないようです。＞もし、バグの理由がわかれば教えてください。う～ん、目の前にあるのは自作の多倍長ルーチンのバグですが、それがなくなったとしても、まだダメです。というのは、LOG(1000)あたりでさえ、級数の収束が非常に遅くなってしまうという致命傷が生じているため。あと、1000桁モードを使用していながら、そのうち10桁しか使っていないというのは、大いなるムダでしょう。（メモリーも、計算量も。） 1000桁モードなら、配列1つあたり400桁でもいけるでしょう。ただし、これらは、1000桁以上の答がほしい場合であり、 1000桁弱の答でいいなら、配列を使わないで生で計算してしまえばいいです。たぶん、以下で正しい範囲は995～997桁くらいだと思います。 ! 1000桁モードLOG計算 OPTION ARITHMETIC DECIMAL_HIGH LET X=5 !LOG(X)を求める LET KETA=1000 !求めたい桁数 LET EPS=1/10^KETA/10 !収束判定値 CALL EE(E) !自然対数の底 ! 級数の収束を加速するためLOG(X)整数部の概略値を計算 LET GETA=0 DO WHILE X>1.65 !1.65はSQRT(E)の近似値。 LET GETA=GETA+1 LET X=X/E LOOP DO WHILE X<0.6 !0.6は1/SQRT(E)の近似値 LET GETA=GETA-1 LET X=XE LOOP ! LOG(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・} LET XX=(X-1)/(X+1) LET XX2=XX LET LG=0 LET LG0=XX2 LET I=1 DO WHILE ABS(LG0)>EPS ! このループを回る回数は、1000桁計算時最悪で800回強。 LET LG=LG+LG0 LET XX2=XX2XXXX LET I=I+2 LET LG0=XX2/I LOOP LET LG=GETA+2LG PRINT LG STOP SUB EE(E) !自然対数の底の1000桁計算 DIM F(450) ! E=1+1/1!+1/2!+1/3!+.....を、桁落ち対策のため尻から足す LET E=0 LET F(1)=1 FOR I=2 TO 450 LET F(I)=F(I-1)*I NEXT I FOR I=450 TO 1 STEP -1 LET E=E+1/F(I) NEXT I LET E=E+1 END SUB END

Re: 多倍長数のＬＯＧの計算	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/05/04 23:34:48
バグの箇所がわかりました。テーラー項を求める際の整数で割る手順に誤りが見つかりました。

Re: バグの箇所がわかりました。	返事を書く
会社員 <qiekrrepwd> 2007/05/05 00:05:41 ** この記事は削除されました

Re: バグの箇所がわかりました。	返事を書く
会社員 <qiekrrepwd> 2007/05/05 00:09:15 ** この記事は削除されました

Re: バグの箇所がわかりました。	返事を書く
会社員 <qiekrrepwd> 2007/05/05 15:33:09 ** この記事は削除されました

Re: バグの箇所がわかりました。	返事を書く
会社員 <qiekrrepwd> 2007/05/05 15:34:14 ** この記事は削除されました

Re: バグの箇所がわかりました。	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/05/16 23:03:53
小数点以下の場合に誤動作していたところを修正です !多倍長のLOG計算 10桁区切り !log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・} OPTION ARITHMETIC DECIMAL_HIGH ! 10進 1000桁 OPTION BASE 0 ! ＝＝＝＝＝＝＝＝＝初期設定＝＝＝＝＝＝＝＝＝＝＝＝ DO INPUT PROMPT "いくらの桁まで計算しますか？":KTA print "桁数=";KTA;"桁" LOOP UNTIL MOD(KTA,20)=0 AND KTA>=100 LET HYRE2=ROUND(KTA1.15/10,0) ! 15％のマージンをつけおく ! 10桁区切り PRINT "HYRE2=";HYRE2 DIM N(HYRE2) ! テスト用の初期値設定 PRINT "試験用 LOG(5)の計算" LET N(0)=5! テスト用初期値LOG(5) FOR I=1 TO HYRE2 LET N(I)=0 NEXT I ! ＝＝＝＝＝＝＝＝＝初期設定＝＝＝＝＝＝＝＝＝＝＝＝ DIM N_1(HYRE2),INV_NP1ANS(HYRE2),KAKE_ANS(HYRE2) DIM NM1_NP1(HYRE2),NM1_NP1_2(HYRE2) DIM LG_N(HYRE2),TYL_K(HYRE2) ! log(N)= ２{ (N-1)/(N+1) + (N-1)/(N+1)^3/3 + (N-1)/(N+1)^5/5 +・・・} PRINT "計算開始 ";DATE$;"　";TIME$ CALL TYLR_PRE ! 事前準備 FOR S=1 TO 100000000 STEP 2 IF S=1 THEN FOR I=0 TO HYRE2 LET TYL_K(I)=NM1_NP1(I) ! 初項に (N-1)/(N+1)を設定 NEXT I END IF IF S=>3 THEN CALL KAKEZAN(TYL_K,NM1_NP1_2,UBOUND(TYL_K),KAKE_ANS) ! (N-1)/(N+1)^3/3 答えが、KKZ_ANS()で帰ってくる FOR I=0 TO HYRE2 LET TYL_K(I)=KAKE_ANS(I)(S-2) ! (N-1)/(N+1)^2(S-2) 倍する NEXT I CALL KTA_SYORI(TYL_K,UBOUND(TYL_K)) !桁処理 END IF FOR I=0 TO HYRE2 IF I=<HYRE2-1 THEN LET AMARI=MOD(TYL_K(I),S) LET TYL_K(I)=INT(TYL_K(I)/S) ! Sで割る LET TYL_K(I+1)=AMARI10^10+TYL_K(I+1) ! 余りを次の桁へ渡す END IF IF I=HYRE2 THEN LET TYL_K(I)=INT(TYL_K(I)/S) ! Sで割る END IF NEXT I CALL KTA_SYORI(TYL_K,UBOUND(TYL_K)) !桁処理 FOR I=0 TO HYRE2 LET LG_N(I)=LG_N(I)+TYL_K(I) NEXT I CALL KTA_SYORI(LG_N,UBOUND(LG_N)) !桁処理 IF TYL_K(HYRE2)=<10^9 AND TYL_K(HYRE2-1)=0 THEN PRINT "終了！ ";DATE$;"　";TIME$ EXIT FOR END IF NEXT S

Re: 小数点以下の場合に誤動作していたところを...	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/05/16 23:04:58
!続き FOR I=0 TO HYRE2 LET LG_N(I)=2LG_N(I) NEXT I CALL KTA_SYORI(LG_N,UBOUND(LG_N)) !桁処理 PRINT "" FOR I=0 TO HYRE2 IF I=0 THEN PRINT STR$(LG_N(I));"." END IF IF I>0 THEN PRINT REPEAT$("0",10-LEN(STR$(LG_N(I)))) ;STR$(LG_N(I));" "; END IF IF I>0 AND I/5=INT(I/5) THEN PRINT "：";STR$(I10) NEXT I !==========サブルーチン============= 事前準備 SUB TYLR_PRE FOR I=0 TO HYR2 IF I=0 THEN LET N_1(0)=N(I)-1 !(N-1) IF i>0 THEN LET N_1(I)=N(I) NEXT I CALL INV_N(N,UBOUND(N),INV_NP1ANS) !１／(N+1)INV_NP1()が帰ってくる CALL KAKEZAN(N_1,INV_NP1ANS,UBOUND(N_1),KAKE_ANS) ! (N-1)／(N+1)を計算 KAKE_ANS()が帰ってくる FOR I=0 TO HYRE2 LET NM1_NP1(I)=KAKE_ANS(I) NEXT I CALL KAKEZAN(NM1_NP1,NM1_NP1,UBOUND(NM1_NP1),KAKE_ANS) ! (N-1)／(N+1)^2を計算 KAKE_ANS()が帰ってくる FOR I=0 TO HYRE2 LET NM1_NP1_2(I)=KAKE_ANS(I) NEXT I end SUB !●●●●1/(N+1)のルーチン●●●● SUB INV_N(A(),LLLL,INV_NP1()) !１/(N+1)の作成 DIM Xn1(LLLL),Xn0(LLLL),KKZ_ANS(LLLL),Xn0_2(LLLL),NP1(LLLL) ! １/(N+1)の作成 FOR I=0 TO LLLL IF I=0 THEN LET NP1(0)=A(0)+1 IF I<>0 THEN LET NP1(I)=A(I) NEXT I CALL KTA_SYORI(NP1,UBOUND(NP1)) ! Xn+1 = 2 * Xn - A * Xn^2 LET Xn0(0)=INT(1/NP1(0)) ! 初期化 LET Xn0(1)=INT(1/NP1(0)10^10)-INT(1/NP1(0))10^10 ! 初期化 FOR I=2 TO LLLL ! ニュートン法パラメータ初期化 LET Xn0(I)=0 LET Xn1(I)=0 NEXT I LET CNT=0 LET KTA_T=KTA DO LET KTA_T=KTA_T/10 LET CNT=CNT+1 LOOP UNTIL KTA_T<1 LET LOOOP=CNT/0.3010+2 PRINT "ループ回数=";ROUND(LOOOP,2);"回" FOR JJJJ=1 TO LOOOP FOR I=0 TO HYRE2 LET Xn1(I)=0 NEXT I CALL KAKEZAN(Xn0,Xn0,UBOUND(Xn0),KKZ_ANS) !乗算サブルーチン Xn^2 FOR I=0 TO HYRE2 LET Xn0_2(I)=KKZ_ANS(I) NEXT I

Re: 小数点以下の場合に誤動作していたところを...	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/05/16 23:05:47
!続き２ ! A * Xn^2の計算開始！ CALL KAKEZAN(NP1,Xn0_2,UBOUND(A),KKZ_ANS) !乗算サブルーチン A * Xn^2 ! Xn1(I)=2Xn0(I)-KKZ_ANS(I) 計算開始！ FOR I=0 TO HYRE2 ! LET Xn1(I)=2Xn0(I)-KKZ_ANS(I) NEXT I CALL KTA_SYORI(Xn1,UBOUND(Xn1)) FOR I=0 TO LLLL ! LET Xn0(I)=Xn1(I) NEXT I NEXT JJJJ FOR I=0 TO LLLL IF I=0 THEN LET INV_NP1(0)=Xn1(I) IF I<>0 THEN LET INV_NP1(I)=Xn1(I) NEXT I END SUB ! 掛算サブルーチン始まり■■■■ SUB KAKEZAN(JUX(),JUY(),LLL,KKZ_ANS()) FOR I=0 TO LLL LET KKZ_ANS(I)=0 NEXT I ! X * Yの計算 FOR I=0 TO LLL FOR J=0 TO LLL-I LET KKZ_ANS(I+J)=JUX(I)JUY(J)+KKZ_ANS(I+J) NEXT J NEXT I CALL KTA_SYORI(KKZ_ANS,LLL) END SUB SUB KTA_SYORI(BBB(),LLL) !桁処理 DO ! 桁処理１ --負記号の削除 LET MNSF=0 IF BBB(1)<0 THEN LET BBB(0)=BBB(0)+INT(BBB(1)/10^10)-1 LET BBB(1)=BBB(1)-10^10INT(BBB(1)/10^10)+10^10 END IF FOR I=1 TO LLL IF BBB(I)<0 THEN LET MNSF=1 LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10)-1 LET BBB(I)=BBB(I)-10^10INT(BBB(I)/10^10)+10^10 END IF NEXT I LOOP UNTIL MNSF=0 DO ! 桁処理２ --桁上げ処理 LET OBF=0 FOR I=LLL TO 1 STEP -1 LET BBB(I-1)=BBB(I-1)+INT(BBB(I)/10^10) LET BBB(I)=BBB(I)-10^10INT(BBB(I)/10^10) IF BBB(I)>10^10 THEN LET OBF=1 END IF NEXT I LOOP UNTIL OBF=0 END SUB END

Re: 多倍長数のＬＯＧの計算	返事を書くノートメニュー
会社員 <qiekrrepwd> 2007/05/16 23:10:37
LOG（π）　の結果です。 1. 1447298858 4940017414 3427351353 0587116472 9481291531 ：50 1571513623 0714721377 6988482607 9783623270 2754897077 ：100 0200981222 8697989159 0482055279 2345658727 9081078810 ：150 2868252763 9391426634 5902902484 7733588699 3778920311 ：200 9630824756 7940119160 2821722737 9888126563 1780498236 ：250 9731331069 5003600064 4054872638 8022327009 6433504959 ：300 5118150662 3725246834 3391269896 5797514047 7703857799 ：350 5399825842 5660228485 0148136217 9159252505 6707638686 ：400 0280763456 8897505123 3436078143 9914144264 2959671289 ：450 7781136526 4523450410 5900716081 8570824981 1881831868 ：500 9767284592 8110257656 8751724223 3833718927 3043288217 ：550 3486510427 6153237516 1028392221 3401436967 1758561644 ：600 2473718780 5060466920 5628337731 0133621627 4515898752 ：650 0151299654 5465739691 5282523916 9585245379 3594601400 ：700 3799565196 6603653800 0112659858 5001297656 9906074466 ：750 7455472671 0450849506 6855874339 0774251341 5924126523 ：800 1777178491 7799588095 7678805102 9644475090 1508911403 ：850 2780807683 3733793894 9488075152 8900918753 6376608670 ：900 7435833345 1081392325 3557406768 4327431198 0496339997 ：950 6180304622 1286361595 8598364047 5800986179 9938264629 ：1000