新しく発言する EXIT インデックスへ
πの多桁

  πの多桁 山中和義 2007/05/25 09:39:22  (修正1回)
  e、n!の多桁 山中和義 2007/05/25 09:45:16  (修正1回)
   └符号付きの計算を可能にする(一部) 山中和義 2007/05/26 20:27:43  (修正2回)
    ├!平方根の計算 山中和義 2007/05/26 20:29:13  (修正2回)
    ├LOG(5)の多桁 山中和義 2007/05/27 19:17:58  (修正4回)
    ├かなり速いです。これは、皮肉ではありませ... SECOND 2007/05/29 13:09:42  (修正3回)
    │└!LOG(5)、BASE=1E250の結果。 SECOND 2007/05/30 07:01:24 
    │ └計算方法の見直し 山中和義 2007/05/30 07:25:38  (修正2回)
    │  ├logの収束は遅いので,log(5)=2*log(sqr(5))... 白石 和夫 2007/05/30 07:35:17 
    │  ├500MHzでの結果です SECOND 2007/05/30 07:53:09  (修正1回)
    │  │└Re:500MHzでの結果です 河川屋 2007/05/30 22:06:59 
    │  │ └アルゴリズムが、無限桁であることに、魅力... SECOND 2007/05/30 23:43:16  (修正1回)
    │  ├500MHzでの結果2です、2000桁、18.62秒、B... SECOND 2007/05/30 23:11:30 
    │  │├走行条件 SECOND 2007/05/30 23:12:27 
    │  │└ビンゴです。 会社員 2007/06/02 21:33:31 
    │  │ └数値が合っていて、よかったですね。 SECOND 2007/06/03 00:32:57 
    │  ├LOG(5)やLOG(π)の計算 山中和義 2007/05/31 14:11:19 
    │  │├サンプル(サブルーチン部分は省略) 山中和義 2007/05/31 14:15:19  (修正1回)
    │  ││└続き 山中和義 2007/05/31 14:16:17 
    │  │└log(π)2000桁 会社員 2007/06/02 21:59:44 
    │  └付録、級数log(x)=2*{((x-1)/(x+1))^1/1+((... SECOND 2007/06/02 22:51:33 
    │   └!付録2、Library\log.libの使い方。(1000桁... SECOND 2007/06/03 03:53:42 
    │    └!参考、log.libの内部関数形式 SECOND 2007/06/03 05:40:29  (修正6回)
    ├多桁の逆数 山中和義 2007/05/30 20:58:13 
    └整数べきの計算 山中和義 2007/06/07 21:33:58 

  πの多桁 山中和義 2007/05/25 09:39:22  (修正1回) ツリーへ

πの多桁 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/25 09:39:22 ** この記事は1回修正されてます
!多桁の計算

!桁の数字の列を格納する配列a()を考える。
!上位桁      下位桁
!a(0),a(1),…,a(L2-1),a(L2) の構造で表すことができる。
!言い換えると、n進数表記の各桁がa()となる。
!ここでは、10000進数とする。 したがって、a(k)は正の整数4桁(0〜9999)。
!例. 1234567は123*10000+4567だから、a(0)=123、a(1)=4567。

!ladd,lsub,lmul,ldivは正の多桁の計算である。
!負の値は、符号や補数を導入して、正負の計算は、場合分けが必要である。

OPTION BASE 0
LET BASE=10000 !基数
SUB ladd(a(),b(), c()) !多桁+多桁
LET cy=0 !桁上がり
FOR i=L2 TO 0 STEP -1 !下の桁から
LET d=a(i)+b(i)+cy !同じ桁で
IF d<BASE THEN
LET c(i)=d
LET cy=0
ELSE
LET c(i)=d-BASE
LET cy=1 !上の桁へ
END IF
NEXT i
IF cy>0 THEN
PRINT "オーバーフロー"
STOP
END IF
END SUB
SUB lsub(a(),b(), c()) !多桁−多桁 ※A>B
LET brrw=0 !借り
FOR i=L2 TO 0 STEP -1 !下の桁から
LET d=a(i)-b(i)-brrw !同じ桁で
IF d>=0 THEN
LET c(i)=d
LET brrw=0
ELSE
LET c(i)=d+BASE
LET brrw=1 !上の桁から借りる
END IF
NEXT i
IF brrw>0 THEN
PRINT "A>Bではありません。"
STOP
END IF
END SUB
SUB lmul(a(),b, c()) !多桁×正の整数(0〜32767)
LET cy=0
FOR i=L2 TO 0 STEP -1 !下の桁から
LET d=a(i)*b+cy
LET c(i)=MOD(d,BASE) !この桁
LET cy=INT(d/BASE) !桁上がり
NEXT i
IF cy>0 THEN
PRINT "オーバーフロー"
STOP
END IF
END SUB
SUB ldiv(a(),b, c()) !多桁÷正の整数(0〜32767)
LET r=0 !余り
FOR i=0 TO L2 !上の桁から
LET d=a(i)+r
LET c(i)=INT(d/b) !商はこの桁
LET r=MOD(d,b)*BASE !余りを下の桁へ
NEXT i
END SUB


LET L=1000 !求める桁数
LET L1=INT(L/4)+1 !配列の大きさ
LET L2=L1+1
LET N=L/1.39794+1 !計算する項 ※10^Lより小さくなる項(収束)

DIM s(L2),w(L2),v(L2),q(L2)
FOR k=0 TO L2
LET s(k)=0
LET w(k)=0
LET v(k)=0
LET q(k)=0
NEXT k

LET w(0)=16*5 !マチン(J.Machin)の公式
LET v(0)=4*239
FOR k=1 TO N !第n項 (16/5^(2n-1)-4/239^(2n-1))/(2n-1) 符号はnが奇数なら正
CALL ldiv(w,25, w)
CALL ldiv(v,239, v)
CALL ldiv(v,239, v)
CALL lsub(w,v, q)
CALL ldiv(q,2*k-1, q)
IF MOD(k,2)<>0 THEN CALL ladd(s,q, s) ELSE CALL lsub(s,q, s)
NEXT k

PRINT s(0);"." !表示する
FOR i=1 TO L1-1
PRINT USING " %%%%":s(i);
IF MOD(i,16)=0 THEN PRINT
NEXT i
PRINT

END

  e、n!の多桁 山中和義 2007/05/25 09:45:16  (修正1回) ツリーへ

Re: πの多桁 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/25 09:45:16 ** この記事は1回修正されてます
e、n!の多桁


!eの計算

LET L=1000 !求める桁数
LET L1=INT(L/4)+1 !配列の大きさ
LET L2=L1+1
LET N=451 !計算する項 ※10^Lより小さくなる項

DIM s(L2)
FOR k=0 TO L2
LET s(k)=0
NEXT k

LET s(0)=1 !1.0
FOR k=N TO 1 STEP -1 !テイラー展開 ((…(((1/n+1)/(n-1)+1)/(n-2)+ … +1)/2+1)/1+1
CALL ldiv(s,k, s)
LET s(0)=s(0)+1
NEXT k

PRINT s(0);"." !表示する
FOR i=1 TO L1-1
PRINT USING " %%%%":s(i);
IF MOD(i,16)=0 THEN PRINT
NEXT i
PRINT



------------------------------------


!階乗n!の計算

LET L=64 !求める桁数
LET L1=INT(L/4)+1
LET L2=L1+1

DIM s(L2)
FOR k=0 TO L2
LET s(k)=0
NEXT k

LET s(L2)=1
FOR k=1 TO 50 !1〜50まで
CALL lmul(s,k, s)

PRINT USING "##":k;
PRINT "!=";
FOR i=0 TO L2 !表示する
PRINT USING "%%%%":s(i);
NEXT i
PRINT
NEXT k

   └符号付きの計算を可能にする(一部) 山中和義 2007/05/26 20:27:43  (修正2回) ツリーへ

Re: e、n!の多桁 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/26 20:27:43 ** この記事は2回修正されてます
符号付きの計算を可能にする(一部)

前出の手続きladd,lsub,lmul,ldivは、正の多桁の計算である。
小数点は、固定とする。位置はプログラムに依存する。

負の値は、符号や補数を導入して、正負の計算は、場合分けが必要である。


SUB lcopy(a(), b()) !コピー
FOR i=0 TO L2
LET b(i)=a(i) !同じ桁で
NEXT i
END SUB
FUNCTION lcomp(a(),b()) !比較 A>Bなら1、A=Bなら0、A<Bなら-1
LET lcomp=0
FOR i=0 TO L2 !上の桁から
IF a(i)>b(i) THEN
LET lcomp=1
EXIT FUNCTION
ELSEIF a(i)<b(i) THEN
LET lcomp=-1
EXIT FUNCTION
END IF
NEXT i
END FUNCTION
SUB llmul(a(),b(), c()) !多桁×多桁
FOR i=0 TO L2*2
LET c(i)=0
NEXT i

FOR j=L2 TO 0 STEP -1 !乗数:下の桁から
LET cy=0
IF b(j)>0 THEN !0は計算しない
FOR i=L2 TO 0 STEP -1 !被乗数:下の桁から
LET d=a(i)*b(j)+cy + c(i+j) !累積 ※筆算参照
LET c(i+j)=MOD(d,BASE) !この桁
LET cy=INT(d/BASE) !桁上がり
NEXT i
IF j=0 THEN
IF cy>0 THEN
PRINT "オーバーフロー!!!"
STOP
END IF
ELSE
LET c(i+j)=cy !j-1
END IF
END IF
NEXT j
END SUB
SUB slsub(a(),aSgn,b(),bSgn, c(),cSgn) !符号付多桁の減算
IF aSgn>0 THEN
IF bSgn>0 THEN !++
IF lcomp(a,b)>0 THEN !大きい数から小さい数を引く
CALL lsub(a,b, c)
LET cSgn=1
ELSE
CALL lsub(b,a, c)
LET cSgn=-1
END IF
ELSE !+−
CALL ladd(a,b, c)
LET cSgn=1
END IF
ELSE
IF bSgn>0 THEN !−+
CALL ladd(a,b, c)
LET cSgn=-1
ELSE !−−
IF lcomp(a,b)>0 THEN !大きい数から小さい数を引く
CALL lsub(a,b, c)
LET cSgn=-1
ELSE
CALL lsub(b,a, c)
LET cSgn=1
END IF
END IF
END IF
END SUB
FUNCTION lEps(a()) !0収束を確認する
FOR l=0 TO L2
IF a(l)<>0 THEN EXIT FOR
NEXT l
LET lEps=l !桁位置
END FUNCTION

    ├!平方根の計算 山中和義 2007/05/26 20:29:13  (修正2回) ツリーへ

Re: 符号付きの計算を可能にする(一部) 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/26 20:29:13 ** この記事は2回修正されてます
!平方根の計算

LET L=1000 !求める桁数
LET L1=INT(L/4)+1
LET L2=L1+1
LET N=15

LET L3=L2*2
DIM s(L3),x(L2),q(L2),v(L2)
FOR k=0 TO L3
LET s(k)=0
NEXT k
FOR k=0 TO L2
LET x(k)=0
LET q(k)=0
LET v(k)=0
NEXT k

LET AA=2 !SQR(A)
LET q(0)=1 !定数 1.0
LET xSgn=1
LET x(1)=5000 !x0=0.5
FOR k=1 TO N !ニュートン法
CALL llmul(x,x, s) !Xn+1=Xn-Xn*(A*Xn^2-1)/2 根の逆数(1/SQR(A))を求める
CALL lmul(s,AA, v)
CALL slsub(v,1,q,1, v,vSgn)
CALL llmul(x,v, s)
LET vSgn=xSgn*vSgn
CALL ldiv(s,2, v)
CALL slsub(x,xSgn,v,vSgn, x,xSgn)
NEXT k
CALL lmul(x,AA, s) !SQR(A)=Xn*A

PRINT s(0);"." !表示する
FOR i=1 TO L1-1
PRINT USING " %%%%":s(i);
IF MOD(i,16)=0 THEN PRINT
NEXT i
PRINT


収束を15回としていますが、収束判定lcompで見極めることができます。
サブルーチンslsubを使っていますが、使わなくても計算可能です。

「LOG(5)やLOG(π)の計算」を参照してください。

    ├LOG(5)の多桁 山中和義 2007/05/27 19:17:58  (修正4回) ツリーへ

Re: 符号付きの計算を可能にする(一部) 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/27 19:17:58 ** この記事は4回修正されてます
LOG(5)の多桁

1000桁を求めるのに、2900回程度の繰り返しでやっと収束します。


LET t0=TIME

!LOGの計算

LET L=1000 !求める桁数
LET L1=INT(L/4)+1
LET L2=L1+1

LET L3=L2*2
DIM s(L2),v(L2),x(L2),x2(L3),t(L3)
FOR k=0 TO L3
LET x2(k)=0
LET t(k)=0
NEXT k
FOR k=0 TO L2
LET x(k)=0
LET s(k)=0
LET v(k)=0
NEXT k

LET A=5 !LOG(A)=2*((A-1)/(A+1))^(2*k-1)/(2*k-1)
LET x(0)=A-1 !x=(n-1)/(n+1)
CALL ldiv(x, A+1, x)
LET k=1
DO
CALL ldiv(x,2*k-1, v) !s=s+x/(2*k-1)
CALL ladd(s,v, s)

CALL lmul(x,(A-1)*(A-1), x) !x=x*x^2 次へ
CALL ldiv(x,(A+1)*(A+1), x)
LET k=k+1
LOOP WHILE lEps(v)<L2
CALL lmul(s,2, s) !s=2*s


PRINT "繰り返し回数=";k
PRINT "計算時間=";TIME-t0


PRINT s(0);"." !表示する
FOR i=1 TO L1-1
PRINT USING " %%%%":s(i);
IF MOD(i,16)=0 THEN PRINT
NEXT i
PRINT

    ├かなり速いです。これは、皮肉ではありませ... SECOND 2007/05/29 13:09:42  (修正3回) ツリーへ

Re: 符号付きの計算を可能にする(一部) 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/05/29 13:09:42 ** この記事は3回修正されてます
かなり速いです。これは、皮肉ではありません。多倍長計算が
速いです。確かにLOG(5)は、収束が遅いですが、他の原因です。
基数の BASE を、もっと大きく移動できれば、絶対速くなる。

    │└!LOG(5)、BASE=1E250の結果。 SECOND 2007/05/30 07:01:24  ツリーへ

Re: かなり速いです。これは、皮肉ではありませ... 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/05/30 07:01:24
!LOG(5)、BASE=1E250 の結果。

実行時間= 20.82 (Pentium3-500MHz)、 1.5GHz~なら多分7~8秒
繰り返し回数= 2831

1 .
6094379124 3410037460 0759333226 1876395256 0135426851
7721912647 8914741789 8770765776 4630133878 0931796107
9996630302 1715562899 7240052293 2467619963 3616617463
7057275521 7963749718 3245653492 8562023415 2505727015
5193600879 7773897256 8819354071 2766154731 2218095279

4852129282 1358059722 5676722852 8724046158 9448178364
6713286739 9842463775 9593189423 8439343534 5105097505
4454194740 5013659870 8786738321 3130572972 0406594853
8383872366 2753876545 5627181615 1165993091 5243207364
9116778639 0067587258 5778766391 5838368239 5042548795

6239484031 0019826971 1749099374 1498480957 6210169110
1437886240 3354321512 7231257345 8846155978 7291980886
5706840200 6659984447 2699973217 6811865176 0122020297
8408109019 6475266997 6506892065 8978913381 5717162072
2777455977 0534320387 7874968293 7536113380 0946764048

3302850021 7477487970 8007143465 6163598970 4125890133
7641524017 0058885987 3494848778 7971054314 1783202014
5306209061 5788790671 5981251674 9904143417 2122036944
7000926855 9375978549 2118961786 4328758768 4832205835
5961996791 3166965093 2338853158 5398982370 9854314630

実行時間= 20.82
繰り返し回数= 2831



!-------------Condition
!LOG(5)の多桁

!1000桁を求めるのに、2900回程度の繰り返しでやっと収束します。

OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
LET base桁=250
LET BASE=10^base桁 !基数

LET t0=TIME

!LOGの計算

LET L=1000 !求める桁数
LET L1=INT(L/base桁)+1
LET L2=L1+1

 (
  )同じ
 (

PRINT s(0);"." !表示する
FOR i=1 TO L1-1
LET w$=USING$(REPEAT$("%",base桁),s(i))
FOR j=1 TO base桁 STEP 10
PRINT " ";mid$(w$,j,10);
IF MOD(j,50)=41 THEN PRINT
NEXT j
PRINT
NEXT i

PRINT "実行時間=";TIME-t0
PRINT "繰り返し回数=";k

END

    │ └計算方法の見直し 山中和義 2007/05/30 07:25:38  (修正2回) ツリーへ

Re: !LOG(5)、BASE=1E250の結果。 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/30 07:25:38 ** この記事は2回修正されてます
計算方法の見直し

xのべき乗を求める計算に多桁×多桁の乗算を使っていました。
これは現状、o(n^2)の計算量になります。
しかも、ループ内ですから大変な時間を費やすことに、、、


漸化式部分の抜粋


LET A=5 !LOG(A)=2*((A-1)/(A+1))^(2*k-1)/(2*k-1)
LET x(0)=A-1 !x=(n-1)/(n+1)
CALL ldiv(x, A+1, x)
LET k=1
DO
CALL ldiv(x,2*k-1, v) !s=s+x/(2*k-1)
CALL ladd(s,v, s)

CALL lmul(x,(A-1)*(A-1), x) !x=x*x^2 次へ
CALL ldiv(x,(A+1)*(A+1), x)
LET k=k+1
LOOP WHILE lEps(v)<L2
CALL lmul(s,2, s) !s=2*s

    │  ├logの収束は遅いので,log(5)=2*log(sqr(5))... 白石 和夫 2007/05/30 07:35:17  ツリーへ

Re: 計算方法の見直し 返事を書く ノートメニュー
白石 和夫 <fbdfvqwhki> 2007/05/30 07:35:17
logの収束は遅いので,log(5)=2*log(sqr(5))のように変形して計算すると速くなると思います。
この手法はサンプルプログラムのlibrary\log.libで利用しています。

    │  ├500MHzでの結果です SECOND 2007/05/30 07:53:09  (修正1回) ツリーへ

Re: 計算方法の見直し 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/05/30 07:53:09 ** この記事は1回修正されてます
500MHz での結果です

1 .
6094379124 3410037460 0759333226 1876395256 0135426851
7721912647 8914741789 8770765776 4630133878 0931796107
9996630302 1715562899 7240052293 2467619963 3616617463
7057275521 7963749718 3245653492 8562023415 2505727015
5193600879 7773897256 8819354071 2766154731 2218095279

4852129282 1358059722 5676722852 8724046158 9448178364
6713286739 9842463775 9593189423 8439343534 5105097505
4454194740 5013659870 8786738321 3130572972 0406594853
8383872366 2753876545 5627181615 1165993091 5243207364
9116778639 0067587258 5778766391 5838368239 5042548795

6239484031 0019826971 1749099374 1498480957 6210169110
1437886240 3354321512 7231257345 8846155978 7291980886
5706840200 6659984447 2699973217 6811865176 0122020297
8408109019 6475266997 6506892065 8978913381 5717162072
2777455977 0534320387 7874968293 7536113380 0946764048

3302850021 7477487970 8007143465 6163598970 4125890133
7641524017 0058885987 3494848778 7971054314 1783202014
5306209061 5788790671 5981251674 9904143417 2122036944
7000926855 9375978549 2118961786 4328758768 4832205835
5961996791 3166965093 2338853158 5398982370 9854314631

実行時間= 7.41
繰り返し回数= 3540

    │  │└Re:500MHzでの結果です 河川屋 2007/05/30 22:06:59  ツリーへ

Re: 500MHzでの結果です 返事を書く ノートメニュー
河川屋 <evcdumfptp> 2007/05/30 22:06:59
Re:500MHz での結果です

う〜ん、1000桁でいいなら、別に自力で多倍長ルーチンを組まなくても、
1000桁モードで普通に四則演算すればokですって。
4つ下で書いたプログラムでlog(5)を計算すると、

1 .
6094379124 3410037460 0759333226 1876395256 0135426851
(途中略)
5961996791 3166965093 2338853158 5398982370 985431463X

で、1桁足りないだけ。(整数部を入れてちょうど1000桁。ただし、これ以上は増やせないのと、1000桁めが合っているのは運がいいから?)
実行時間= 0.24秒(athlon3500+ 2.2Ghz)
繰り返し回数(概算)= eを求めるのに450回、 logに700回。
※log(5)でなく、log(5/e)=log(0.67667....)を計算することで収束回数を稼げば2800回もまわさなくて大丈夫。

    │  │ └アルゴリズムが、無限桁であることに、魅力... SECOND 2007/05/30 23:43:16  (修正1回) ツリーへ

Re: Re:500MHzでの結果です 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/05/30 23:43:16 ** この記事は1回修正されてます
アルゴリズムが、無限桁であることに、魅力を感じます。

    │  ├500MHzでの結果2です、2000桁、18.62秒、B... SECOND 2007/05/30 23:11:30  ツリーへ

Re: 計算方法の見直し 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/05/30 23:11:30
500MHzでの結果2です、2000桁、18.62 秒、BASE=1E400

繰り返し回数= 6805
計算時間= 18.62
1 .
6094379124 3410037460 0759333226 1876395256 0135426851
7721912647 8914741789 8770765776 4630133878 0931796107
9996630302 1715562899 7240052293 2467619963 3616617463
7057275521 7963749718 3245653492 8562023415 2505727015
5193600879 7773897256 8819354071 2766154731 2218095279
4852129282 1358059722 5676722852 8724046158 9448178364
6713286739 9842463775 9593189423 8439343534 5105097505
4454194740 5013659870 8786738321 3130572972 0406594853

8383872366 2753876545 5627181615 1165993091 5243207364
9116778639 0067587258 5778766391 5838368239 5042548795
6239484031 0019826971 1749099374 1498480957 6210169110
1437886240 3354321512 7231257345 8846155978 7291980886
5706840200 6659984447 2699973217 6811865176 0122020297
8408109019 6475266997 6506892065 8978913381 5717162072
2777455977 0534320387 7874968293 7536113380 0946764048
3302850021 7477487970 8007143465 6163598970 4125890133

7641524017 0058885987 3494848778 7971054314 1783202014
5306209061 5788790671 5981251674 9904143417 2122036944
7000926855 9375978549 2118961786 4328758768 4832205835
5961996791 3166965093 2338853158 5398982370 9854314631
1007588294 8144753148 9617891122 0148595518 5183682437
9862838066 8942003786 2023016071 9794961521 3758199640
2203356064 0972728415 6381544998 0030196573 4680232678
6239683225 0990910096 4168099970 5575316720 6826150961

5222060579 6683304425 1673080143 7890921605 5068000803
2623194994 6324147661 5174477406 1636230434 6578989656
6101224551 8259802869 6601164051 2920144258 0200794467
4735886317 6325743436 9029415038 7326580344 8006253200
2011778286 0532755204 7501919332 1327014575 8652626139
9058791660 5172464900 2200687293 4154461681 8570048196
5440585469 0522067457 4253474675 1241226962 8744080817
5207310687 9891678805 0333059079 6271009021 5533184482

7511281894 5570078255 1782812162 3902433343 3998226351
2169477366 3131235177 5823518432 6168935155 4549267357
3673891698 7974018447 3743860584 1090793408 9707442617
4633972186 9769918286 7616219433 8066238892 8991523440
1499298514 3570748060 1883820172 9033841623 6856427421
5722892362 7028751742 6794559457 0538459660 2740716554
0790275978 6273677533 9496854869 9971669202 3361423270
8222547970 2670595413 2203032071 4312534353 5367510480

    │  │├走行条件 SECOND 2007/05/30 23:12:27  ツリーへ

Re: 500MHzでの結果2です、2000桁、18.62秒、B... 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/05/30 23:12:27
走行条件

!------------- BASE を 1E400、2000桁 でも、やってみた。

OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
LET base桁=400
LET BASE=10^base桁 !基数

LET t0=TIME

!LOGの計算

LET L=2000 !求める桁数
LET L1=INT(L/base桁)+1
LET L2=L1+1

 ( 「LOG(5)の多桁」のスレッドから。(修正日付が出ないので、2007/5/30以降のコピー)
  )

PRINT "繰り返し回数=";k
PRINT "計算時間=";TIME-t0

PRINT s(0);"." !表示する
FOR i=1 TO L1-1
LET w$=USING$(REPEAT$("%",base桁),s(i))
FOR j=1 TO base桁 STEP 10
PRINT " ";mid$(w$,j,10);
IF MOD(j,50)=41 THEN PRINT
NEXT j
PRINT
NEXT i

!-------「多桁の計算」のスレッドから。

SUB ladd(a(),b(), c()) !多桁+多桁
 (
  )
SUB ldiv(a(),b, c()) !多桁÷正の整数(0〜32767)
 (
  )
END SUB

!-------「符号付きの計算を可能にする(一部)」のスレッドの全部。
 (
  )

END

    │  │└ビンゴです。 会社員 2007/06/02 21:33:31  ツリーへ

Re: 500MHzでの結果2です、2000桁、18.62秒、B... 返事を書く ノートメニュー
会社員 <qiekrrepwd> 2007/06/02 21:33:31
ビンゴです。
log5

自作プログラムにて
2000桁

1.
6094379124 3410037460 0759333226 1876395256 0135426851 : 50
7721912647 8914741789 8770765776 4630133878 0931796107 : 100
9996630302 1715562899 7240052293 2467619963 3616617463 : 150
7057275521 7963749718 3245653492 8562023415 2505727015 : 200
5193600879 7773897256 8819354071 2766154731 2218095279 : 250
4852129282 1358059722 5676722852 8724046158 9448178364 : 300
6713286739 9842463775 9593189423 8439343534 5105097505 : 350
4454194740 5013659870 8786738321 3130572972 0406594853 : 400
8383872366 2753876545 5627181615 1165993091 5243207364 : 450
9116778639 0067587258 5778766391 5838368239 5042548795 : 500
6239484031 0019826971 1749099374 1498480957 6210169110 : 550
1437886240 3354321512 7231257345 8846155978 7291980886 : 600
5706840200 6659984447 2699973217 6811865176 0122020297 : 650
8408109019 6475266997 6506892065 8978913381 5717162072 : 700
2777455977 0534320387 7874968293 7536113380 0946764048 : 750
3302850021 7477487970 8007143465 6163598970 4125890133 : 800
7641524017 0058885987 3494848778 7971054314 1783202014 : 850
5306209061 5788790671 5981251674 9904143417 2122036944 : 900
7000926855 9375978549 2118961786 4328758768 4832205835 : 950
5961996791 3166965093 2338853158 5398982370 9854314631 : 1000
1007588294 8144753148 9617891122 0148595518 5183682437 : 1050
9862838066 8942003786 2023016071 9794961521 3758199640 : 1100
2203356064 0972728415 6381544998 0030196573 4680232678 : 1150
6239683225 0990910096 4168099970 5575316720 6826150961 : 1200
5222060579 6683304425 1673080143 7890921605 5068000803 : 1250
2623194994 6324147661 5174477406 1636230434 6578989656 : 1300
6101224551 8259802869 6601164051 2920144258 0200794467 : 1350
4735886317 6325743436 9029415038 7326580344 8006253200 : 1400
2011778286 0532755204 7501919332 1327014575 8652626139 : 1450
9058791660 5172464900 2200687293 4154461681 8570048196 : 1500
5440585469 0522067457 4253474675 1241226962 8744080817 : 1550
5207310687 9891678805 0333059079 6271009021 5533184482 : 1600
7511281894 5570078255 1782812162 3902433343 3998226351 : 1650
2169477366 3131235177 5823518432 6168935155 4549267357 : 1700
3673891698 7974018447 3743860584 1090793408 9707442617 : 1750
4633972186 9769918286 7616219433 8066238892 8991523440 : 1800
1499298514 3570748060 1883820172 9033841623 6856427421 : 1850
5722892362 7028751742 6794559457 0538459660 2740716554 : 1900
0790275978 6273677533 9496854869 9971669202 3361423270 : 1950
8222547970 2670595413 2203032071 4312534353 5367510480 : 2000

    │  │ └数値が合っていて、よかったですね。 SECOND 2007/06/03 00:32:57  ツリーへ

Re: ビンゴです。 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/06/03 00:32:57
数値が合っていて、よかったですね。
ホッとしました。実は、よく見ていませんでした。
掲示されていた級数は、無断借用されています。ご容赦を。

    │  ├LOG(5)やLOG(π)の計算 山中和義 2007/05/31 14:11:19  ツリーへ

Re: 計算方法の見直し 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/31 14:11:19
LOG(5)やLOG(π)の計算

LOG(A)=2*LOG(SQR(A))より(サンプルプログラムのlibrary\log.lib参照)
・平方根 SQR(A)
・逆数 (N-1)/(N+1)を求めるため
・LOG
・2倍
の手順で算出すると、
LOG(5)の計算の繰り返し回数は、1200回ぐらいです。

LOG(π)で、900回ぐらいです。

    │  │├サンプル(サブルーチン部分は省略) 山中和義 2007/05/31 14:15:19  (修正1回) ツリーへ

Re: LOG(5)やLOG(π)の計算 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/31 14:15:19 ** この記事は1回修正されてます
サンプル(サブルーチン部分は省略)


LET L=1000 !求める桁数
LET L1=INT(L/4)+1
LET L2=L1+1


DIM c1(L2)
FOR k=0 TO L2
LET c1(k)=0
NEXT k
LET c1(0)=1 !定数 1.0


!LOGの計算 LOG(A)=2*((A-1)/(A+1))^(2*k-1)/(2*k-1)

DIM A(L2)
FOR k=0 TO L2
LET A(k)=0
NEXT k


!A=πとする

LET t0=TIME

LET N=L/1.39794+1 !計算する項 ※10^Lより小さくなる項(収束)

LET L3=L2*2
DIM s(L3)
FOR k=0 TO L2
LET s(k)=0
NEXT k
DIM w(L2),v(L2),q(L2)
FOR k=0 TO L2
LET w(k)=0
LET v(k)=0
LET q(k)=0
NEXT k

LET w(0)=16*5 !マチン(J.Machin)の公式
LET v(0)=4*239
FOR k=1 TO N !第n項 (16/5^(2n-1)-4/239^(2n-1))/(2n-1) 符号はnが奇数なら正
CALL ldiv(w,25, w)
CALL ldiv(v,239, v)
CALL ldiv(v,239, v)
CALL lsub(w,v, q)
CALL ldiv(q,2*k-1, q)
IF MOD(k,2)<>0 THEN CALL ladd(s,q, s) ELSE CALL lsub(s,q, s)
NEXT k

PRINT "繰り返し回数=";N
PRINT "計算時間=";TIME-t0

CALL lcopy(s, A)
!-----------------------------------


!B=SQR(A)の計算 フェーズ1

LET t0=TIME

LET L3=L2*2
DIM t(L3)
FOR k=0 TO L3
LET s(k)=0
LET t(k)=0
NEXT k
DIM x(L2)
FOR k=0 TO L2
LET x(k)=0
LET v(k)=0
NEXT k

LET x(1)=5000 !x0=0.5
LET k=0 !回数の確認のみ
DO !ニュートン法
CALL llmul(x,x, s) !!Xn+1=Xn+Xn*(1-A*Xn^2)/2 根の逆数(1/SQR(A))を求める
CALL llmul(s,A, t)
CALL lsub(c1,t, v)
CALL llmul(x,v, s)
CALL ldiv(s,2, v)
CALL ladd(x,v, x)
LET k=k+1
LOOP WHILE lEps(v)<L2 !Xn*(1-A*Xn^2)/2→0
CALL llmul(x,A, s) !SQR(A)=Xn*A


PRINT "繰り返し回数=";k
PRINT "計算時間=";TIME-t0

DIM B(L2)
CALL lcopy(s, B)
!-----------------------------------


!B+1の逆数の計算 フェーズ2

LET t0=TIME

FOR k=0 TO L3
LET s(k)=0
NEXT k
FOR k=0 TO L2
LET x(k)=0
LET v(k)=0
NEXT k

CALL ladd(B,c1, A) !A=B+1
LET x(1)=1000 !0.1 |x0|<2/A
LET k=0 !回数の確認のみ
DO !ニュートン法
CALL llmul(x,A, s) !Xn+1=Xn*(2-A*Xn)
CALL lsub(c1,s, v)
LET cEps=lEps(v) !右辺=Xn*(1+(1-A*Xn))と変形
LET v(0)=v(0)+1 !call ladd(c1,v, v)
CALL llmul(x,v, s)
CALL lcopy(s, x)
LET k=k+1
LOOP WHILE cEps<L2 !1-A*Xn→0


PRINT "繰り返し回数=";k
PRINT "実行時間=";TIME-t0

DIM Bp1i(L2)
CALL lcopy(s, Bp1i)
!-----------------------------------

    │  ││└続き 山中和義 2007/05/31 14:16:17  ツリーへ

Re: サンプル(サブルーチン部分は省略) 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/31 14:16:17
続き



!LOG(B)の計算 フェーズ3

LET t0=TIME

DIM x2(L3)
FOR k=0 TO L3
LET x2(k)=0
LET t(k)=0
NEXT k
FOR k=0 TO L2
LET x(k)=0
LET s(k)=0
LET v(k)=0
NEXT k

CALL lsub(B,c1, v) !x=(n-1)/(n+1)
CALL llmul(v,Bp1i, t)
CALL lcopy(t, x)
CALL llmul(x,x, x2)
LET k=1
DO
CALL ldiv(x,2*k-1, v) !s=s+x/(2*k-1)
CALL ladd(s,v, s)

CALL llmul(x,x2, t) !x=x*x^2 次へ
CALL lcopy(t, x)
LET k=k+1
LOOP WHILE lEps(v)<L2
CALL lmul(s,2, s) !s=2*s


PRINT "繰り返し回数=";k
PRINT "計算時間=";TIME-t0
!-----------------------------------


!LOG(A)=2*LOG(SQR(A))の計算 フェーズ4

LET t0=TIME

CALL lmul(s,2, s) !s=2*s

PRINT "計算時間=";TIME-t0
!-----------------------------------



PRINT s(0);"." !結果を表示する
FOR i=1 TO L1-1
PRINT USING " %%%%":s(i);
IF MOD(i,16)=0 THEN PRINT
NEXT i
PRINT

    │  │└log(π)2000桁 会社員 2007/06/02 21:59:44  ツリーへ

Re: LOG(5)やLOG(π)の計算 返事を書く ノートメニュー
会社員 <qiekrrepwd> 2007/06/02 21:59:44
log(π) 2000桁
自作プログラムにて
計算時間 11分27秒

1.
1447298858 4940017414 3427351353 0587116472 9481291531 1571513623 0714721377 6988482607 9783623270 2754897077
0200981222 8697989159 0482055279 2345658727 9081078810 2868252763 9391426634 5902902484 7733588699 3778920311
9630824756 7940119160 2821722737 9888126563 1780498236 9731331069 5003600064 4054872638 8022327009 6433504959
5118150662 3725246834 3391269896 5797514047 7703857799 5399825842 5660228485 0148136217 9159252505 6707638686
0280763456 8897505123 3436078143 9914144264 2959671289 7781136526 4523450410 5900716081 8570824981 1881831868
9767284592 8110257656 8751724223 3833718927 3043288217 3486510427 6153237516 1028392221 3401436967 1758561644
2473718780 5060466920 5628337731 0133621627 4515898752 0151299654 5465739691 5282523916 9585245379 3594601400
3799565196 6603653800 0112659858 5001297656 9906074466 7455472671 0450849506 6855874339 0774251341 5924126523
1777178491 7799588095 7678805102 9644475090 1508911403 2780807683 3733793894 9488075152 8900918753 6376608670
7435833345 1081392325 3557406768 4327431198 0496339997 6180304622 1286361595 8598364047 5800986179 9938264629 :1000桁
2776462759 4848489641 4107483132 5934620536 3507304605 5030768215 4944441547 7888455953 5228440047 8509182172
5591517990 0785243523 8371128671 3234290556 6964492585 5826231188 2422324466 1476739136 1533394142 6453460088
1979155478 9677575298 7830759323 0499751706 7853706663 1522213475 1026417324 9189065342 5737305183 5228316776
8773114429 4436810899 7522287634 5549099334 6925398102 8398378467 6950799719 6516300838 6496663274 2238867613
9294411237 9606529081 4635455024 1519364336 8404005225 6155756180 5368045961 3160686367 2262971268 4805551803
8239624057 9831384339 5588248355 6816617339 0181955089 2466778204 2898879384 6230819535 0708252369 9065543916
0296765653 4950948710 2686726405 0363448899 5781395484 0804697878 6037235600 3103351889 0166410542 2451404008
2148002607 1893924502 0777856356 9881069323 3664357379 4810929277 8193626598 0614204270 0943982983 6473376792
2501305495 4459753800 3764761751 9082652294 8577288283 4937991341 8698964043 4834570915 5046062991 2859614271
4322563776 9979432888 9523074041 4635294661 1331364188 4192574888 1893207965 7199144493 9402534883 2282628130 :2000桁

    │  └付録、級数log(x)=2*{((x-1)/(x+1))^1/1+((... SECOND 2007/06/02 22:51:33  ツリーへ

Re: 計算方法の見直し 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/06/02 22:51:33
付録、級数 log(x)= 2*{ ((x-1)/(x+1))^1 /1 +((x-1)/(x+1))^3 /3 +((x-1)/(x+1))^5 /5 + …}
は、収束が速いですが、何故かを調べてみる。

log(x) = 2*log( x^(1/2) ) と照合すると、この級数は、log( sqr(x) ) と等価ですが、

log( x^(1/2) ) を直接、テーラー展開しても、この級数は出てきません。

次の2つのマクローリン展開を、使用します。

log(1+x) = x/1-x^2/2+x^3/3-x^4/4+x^5/5 ・・・
log(1-x) = -x/1-x^2/2-x^3/3-x^4/4-x^5/5 ・・・

log(1+x) - log(1-x) =
log( (1+x)/(1-x) ) = 2*{ x/1+x^3/3+x^5/5 ・・・} ここで、高速化しています。

変数を、調整する。

 z=(1+x)/(1-x) → z(1-x)=1+x → (z-1)/(z+1)=x

log(z)= 2*{ ((z-1)/(z+1))^1 /1 +((z-1)/(z+1))^3 /3 +((z-1)/(z+1))^5 /5 ・・・}

    │   └!付録2、Library\log.libの使い方。(1000桁... SECOND 2007/06/03 03:53:42  ツリーへ

Re: 付録、級数log(x)=2*{((x-1)/(x+1))^1/1+((... 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/06/03 03:53:42
!付録2、Library\log.lib の使い方。(1000桁モードで)

!--------Main 側

DECLARE EXTERNAL FUNCTION LOG ! ←この行が不可欠です、LOG()関数の再定義。
OPTION ARITHMETIC DECIMAL_HIGH !1000桁モード

PRINT LOG(5)

END
merge "log.lib"
!------------------


!--------Library 側( Library\log.lib )

! 1000桁モードで利用する対数関数
!EXTERNAL FUNCTION LOG(x)
!OPTION ARITHMETIC DECIMAL_HIGH ! ← 脱落しているので、追加。
!IF x<=0 THEN
! CAUSE EXCEPTION 3004
!ELSEIF x<1 THEN
! LET LOG=-log(1/x)
!ELSEIF x>3 THEN
! LET LOG=2*log(SQR(x))
!ELSE ! 1<=x<=3
! LET h=(x-1)/(x+1) ! 0<=h<=0.5
! LET t=0
! LET n=1
! LET k=h
! LET h2=H^2
! DO
! LET t=t+k/n
! LET n=n+2
! LET k=k*h2
! LOOP UNTIL k<=eps(0)
! LET LOG=2*t
!END IF
!END FUNCTION

!※1000桁モードでは、本来のLOG()関数は、使えませんので
! ここでの、LOG() は、FUNCTION 自身の再帰呼び出しとして機能し、
! 全てのxは、最終的に、1<=x<=3 の範囲になってから、
! DO~LOOP の部分で級数計算されます。

    │    └!参考、log.libの内部関数形式 SECOND 2007/06/03 05:40:29  (修正6回) ツリーへ

Re: !付録2、Library\log.libの使い方。(1000桁... 返事を書く ノートメニュー
SECOND <cszcthjjdj> 2007/06/03 05:40:29 ** この記事は6回修正されてます
!参考、log.lib の内部関数形式

OPTION ARITHMETIC DECIMAL_HIGH
!DECLARE FUNCTION LOG ! 関数名を log() にする時、この行は、省略不可。

PRINT llog(5)
PRINT llog(5000)
PRINT llog(0.5)


! 1000桁モードで利用する対数関数

FUNCTION llog(x)
IF x<=0 THEN
CAUSE EXCEPTION 3004
ELSEIF x<1 THEN
LET llog=-llog(1/x) !再帰function
ELSEIF x>3 THEN
LET llog=2*llog(SQR(x)) !再帰function
ELSE ! 1<=x<=3 ← 全てのx が この範囲になってから、以下の級数計算。
LET h=(x-1)/(x+1)
LET h2=h^2
LET k=1
LET sum=0
DO
LET sum=sum+h/k
LET h=h*h2
LET k=k+2
LOOP UNTIL h<=EPS(0)
LET llog=2*sum
END IF
END FUNCTION

END

!----------------
! llog(256)の場合

! x=256 →再帰 llog(16) with SQR(x )
! x'=16 →再帰 llog( 4) with SQR(x')
! x"= 4 →再帰 llog( 2) with SQR(x")
! x"'=2 →log(2)の級数計算

! x"'のリターン引数llog= 計算値
! x"のリターン引数llog= 2* x"'のリターン引数llog
! x'のリターン引数llog= 2* x"のリターン引数llog
! x のリターン引数llog= 2* x'のリターン引数llog

!----------------
! llog(1/256)の場合

! x0=1/256 →再帰 llog(256) with 1/x0

! x=256 →再帰 llog(16) with SQR(x )
! x'=16 →再帰 llog( 4) with SQR(x')
! x"= 4 →再帰 llog( 2) with SQR(x")
! x"'=2 →log(2)の級数計算

! x"'のリターン引数llog= 計算値
! x"のリターン引数llog= 2* x"'のリターン引数llog
! x'のリターン引数llog= 2* x"のリターン引数llog
! x のリターン引数llog= 2* x'のリターン引数llog

! x0 のリターン引数llog= -1* xのリターン引数llog

    ├多桁の逆数 山中和義 2007/05/30 20:58:13  ツリーへ

Re: 符号付きの計算を可能にする(一部) 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/05/30 20:58:13
多桁の逆数


LET L=1000 !求める桁数
LET L1=INT(L/4)+1
LET L2=L1+1

LET L3=L2*2
DIM s(L3),x(L2),v(L2),A(L2),c1(L2)
FOR k=0 TO L3
LET s(k)=0
NEXT k
FOR k=0 TO L2
LET x(k)=0
LET v(k)=0
LET A(k)=0
LET c1(k)=0
NEXT k

LET c1(0)=1 !定数1
!!!LET AA=7
LET A(0)=7 !A=7.0
LET x(1)=1000 !0.1 |x0|<2/A
LET k=0 !回数の確認のみ
DO !ニュートン法
!!!CALL lmul(x,AA, s) !Xn+1=Xn*(2-A*Xn)
CALL llmul(x,A, s) !Xn+1=Xn*(2-A*Xn)
CALL lsub(c1,s, v)
LET cEps=lEps(v) !右辺=Xn*(1+(1-A*Xn))と変形
LET v(0)=v(0)+1 !call ladd(c1,v, v)
CALL llmul(x,v, s)
CALL lcopy(s, x)
LET k=k+1
LOOP WHILE cEps<L2 !1-A*Xn→0


PRINT "繰り返し回数=";k
PRINT "実行時間=";TIME-t0


PRINT s(0);"." !表示する
FOR i=1 TO L1-1
PRINT USING " %%%%":s(i);
IF MOD(i,16)=0 THEN PRINT
NEXT i
PRINT




効率よく計算するポイントは、
・無駄のないループ回数
・ループ内は軽量な計算
のプログラムにする。

数学的には
・収束しやすい
・計算が少ない
・総計算回数が見当できる、収束判定ができる
初期値と漸化式を選ぶ。

    └整数べきの計算 山中和義 2007/06/07 21:33:58  ツリーへ

Re: 符号付きの計算を可能にする(一部) 返事を書く ノートメニュー
山中和義 <drdlxujciw> 2007/06/07 21:33:58
整数べきの計算


LET t0=TIME


LET L=1000 !求める桁数
LET L1=INT(L/4)+1 !※
LET L2=L1+1


LET L3=L2*2
DIM x2(L3),s(L3),t(L3)
FOR k=0 TO L2
LET x2(k)=0
LET s(k)=0
LET t(k)=0
NEXT k
DIM x(L2)
FOR k=0 TO L2
LET x(k)=0
NEXT k


LET x(1)=8000 !0.8
LET n=12 !べき乗
LET s(0)=1 !1.0

LET k=0 !回数の確認のみ
DO UNTIL n=0
IF MOD(n,2)=1 THEN !ビットが1なら
CALL llmul(s,x, t)
CALL lcopy(t, s)
END IF
CALL llmul(x,x, x2)
CALL lcopy(x2, x)

LET n=INT(n/2) !2進数にする

LET k=k+1
LOOP


PRINT "繰り返し回数=";k
PRINT "計算時間=";TIME-t0


FOR i=0 TO L1-1 !表示する
PRINT USING " %%%%":s(i);
IF MOD(i,16)=0 THEN PRINT
NEXT i
PRINT


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