新しく発言する  EXIT  インデックスへ

連分数の計算


  連分数の計算 山中和義 2008/04/06 13:13:14 

  連分数の計算 山中和義 2008/04/06 13:13:14   ツリーへ
連分数の計算  返事を書く  ノートメニュー
山中和義 <drdlxujciw> 2008/04/06 13:13:14
!連分数(continued fraction)展開

LET N=16 !桁数

DIM L(0 TO N) !連分数のリスト

INPUT PROMPT "分子a、分母bを入力してください。": a,b !有理数
CALL contfrac(a,b,n, L,k)
CALL dispCF(L,k)
CALL dispCF2(L,k)

LET a=PI !円周率π
LET b=1
CALL contfrac(a,b,n, L,k)
CALL dispCF(L,k)

LET a=EXP(1) !自然対数の底e
LET b=1
CALL contfrac(a,b,n, L,k)
CALL dispCF(L,k)

LET a=SQR(2) !平方根
LET b=1
CALL contfrac(a,b,n, L,k)
CALL dispCF(L,k)


SUB dispCF(cf(),k) !連分数a+1/(b+1/(c+...))を、[a;b,c,...]形式で表示する
PRINT a;"/";b;"=[";
FOR i=0 TO k-1
IF i>1 THEN PRINT ",";
IF i=1 THEN PRINT ";";
PRINT cf(i);
NEXT i
PRINT "]"
END SUB
SUB dispCF2(cf(),k) !a+1/(b+1/(c+...))形式で表示する
PRINT a;"/";b;"=";
FOR i=0 TO k-2
PRINT cf(i);"+1/(";
NEXT i
PRINT cf(k-1);REPEAT$(")",k-1)
END SUB

SUB contfrac(a,b,n, cf(),k) !有理数a/bの連分数cf()を求める
LET k=0 !個数

LET x=a !作業変数へ
LET y=b

DO WHILE y<>0 AND k<n !ユークリッドの互除法
LET cf(k)=INT(x/y) !商
LET k=k+1

LET t=MOD(x,y) !余り
LET x=y
LET y=t
LOOP
END SUB




!連分数[a0;a1,....,an]=pn/qnの値を求める

!行列計算による計算

!有理数a/bを連分数[a0,a1, ... ,an]に展開すると
!┌ a ┐=┌ a0 1 ┐┌ a1 1 ┐┌ a2 1 ┐...┌ an 1 ┐┌ 1 ┐
!└ b ┘ └ 1 0 ┘└ 1 0 ┘└ 1 0 ┘  └ 1 0 ┘└ 0 ┘
DIM m(2,2)

LET m(1,1)=L(0)
LET m(1,2)=1
LET m(2,1)=1
LET m(2,2)=0

DIM r(2,2) !求める分子、分母
MAT r=IDN
FOR i=1 TO N
MAT r=r*m
LET m(1,1)=L(i)
NEXT i
MAT PRINT r



!漸化式による計算

DIM p(0 TO N),q(0 TO N) !求める分子、分母

LET p(0)=1
LET q(0)=0

LET p(1)=L(0)
LET q(1)=1

FOR i=2 TO N
LET p(i)=L(i-1)*p(i-1)+p(i-2)
LET q(i)=L(i-1)*q(i-1)+q(i-2)
NEXT i

PRINT L(N-1)*p(N-2)+p(N-3);"/";L(N-1)*q(N-2)+q(N-3)
PRINT L(N)*p(N-1)+p(N-2);"/";L(N)*q(N-1)+q(N-2)




!連分数[a0,a1,....,an]の近似値を求める

LET t=L(N)
FOR i=N-1 TO 0 STEP -1
IF t<>0 THEN LET t=1/t
LET t=L(i)+t
NEXT i
PRINT t



END

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