ｍの平方根の連分数展開なのですがうまく動きません。

質問

subroutine keisan(m)
integer a,a0,n,i,m
real t,X,Y,Z
t=sqrt(real(m))
a0=int(t)
write(*,*) 'a0=',a0
a=a0
X=1
Y=0
Z=1
do 10 n=1,20
if (Z==0.0) exit
X=X*Z
Y=Z*(-Y+real(a)*Z)
Z=X*X*real(m)-(real(a)*Z-Y)**2
do 100 i=1,X
if(real(i)==X) then
X=X
else
if(((mod(int(X),i))==0).AND.((mod(int(Y),i))==0).AND.((mod(int(Z),i))==0)) then
X=X/real(i)
Y=Y/real(i)
Z=Z/real(i)
else
X=X
endif
endif
100 continue
a=int((X*t+Y)/Z)
write(*,*) 'a',n,'=',a
10 continue
end subroutine
サブルーチンのみ載せてますけど・・・・どうなのでしょう？？

回答No.1

a-kuma · ベストアンサー

なんで、こんな難しくなるの？

一番単純な形だと、

　　　　　　　　１
√ｍ ＝ ａ1 ＋ ──
　　　　　　　　Ｘ

で、ａ1 は、√ｍの整数部。で、この Ｘ は、√ｍ－ａ1 の逆数で、
それを又、この形にしてゆく、の繰返しでしょ。

　　　subroutine keisan(m)

　　　integer i,m,a;
　　　double x;

　　　x = sqrt(real(m))
　　　a = int(x);

　　　write(*,*) 'a=',a

　　　do 10 i = 1, 20
　　　　if (x .eq. a) return;
　　　　x = 1.0 / (x - a)
　　　　a = int(x)
　　　　write(*,*) 'a=',a
　10　continue

　　　return
　　　end

＃ fortran は、ずいぶん久しぶりなので、文法間違いがあったらごめんなさい

お礼

無事にプログラム完成致しました。ありがとうございました♪

補足

ありがとうございます・・・でも、「誤差なく」値を出せということだったんでこんなに複雑になったんです。
ヒントとして“Ｘ＝（Ｘ＋Ｙ√ｍ）／Ｚのカタチで書けることが帰納法で証明出来る”というのがあったんで・・・・・