πの多桁 山中和義 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での結果２です、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
│ └!付録２、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)は正の整数４桁（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 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 ------------------------------------ !階乗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 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 !０は計算しない 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()) !０収束を確認する 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 収束を１５回としていますが、収束判定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 |
|
├かなり速いです。これは、皮肉ではありませ... 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 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桁モードで普通に四則演算すればｏｋですって。 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での結果２です、2000桁、18.62秒、B... SECOND 2007/05/30 23:11:30 ツリーへ
| Re: 計算方法の見直し |
返事を書く ノートメニュー |
| SECOND <cszcthjjdj> 2007/05/30 23:11:30 | |
|
500MHzでの結果２です、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での結果２です、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 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での結果２です、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)の計算 フェーズ１ 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の逆数の計算 フェーズ２ 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)の計算 フェーズ３ 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))の計算 フェーズ４ 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 |
|
│ │└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) ) を直接、テーラー展開しても、この級数は出てきません。 次の２つのマクローリン展開を、使用します。 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 ・・・} |
|
│ └!付録２、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 | |
|
!付録２、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 自身の再帰呼び出しとして機能し、 ! 全てのｘは、最終的に、1<=x<=3 の範囲になってから、 ! DO~LOOP の部分で級数計算されます。 |
|
│ └!参考、log.libの内部関数形式 SECOND 2007/06/03 05:40:29 (修正6回) ツリーへ
| Re: !付録２、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 ← 全てのｘ が この範囲になってから、以下の級数計算。 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 !定数１ !!!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 効率よく計算するポイントは、 ・無駄のないループ回数 ・ループ内は軽量な計算 のプログラムにする。 数学的には ・収束しやすい ・計算が少ない ・総計算回数が見当できる、収束判定ができる 初期値と漸化式を選ぶ。 |
|
└整数べきの計算 山中和義 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 !ビットが１なら 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) !２進数にする 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 |
|