平方根の連分数展開

TAOCPの第4.5.3項の演習問題12は, quadratic irrationalityという話題である. よく分からぬが, どうやら平方数でない整数の平方根の連分数展開をしているらしい.

そのアルゴリズムを, さっそくSchemeに書き直した.

(define (cf d n)
 (let* ((x (sqrt d)) (v 1)
   (a (inexact->exact (floor x))) (u a))
  (define (loop)
   (display (list 'v v 'a a 'u u)) (newline)
   (set! v (/ (- d (* u u)) v))
   (set! a (inexact->exact (floor (/ (+ x u) v))))
   (set! u (- (* a v) u))
   (set! n (- n 1))
   (if (> n 0) (loop) 'ok))
  (loop)))

関数cfは, dの平方根の連分数の部分商をn回計算するものである. 31の平方根は

(cf 31 13)
(v 1 a 5 u 5)
(v 6 a 1 u 1)
(v 5 a 1 u 4)
(v 3 a 3 u 5)
(v 2 a 5 u 5)
(v 3 a 3 u 4)
(v 5 a 1 u 1)
(v 6 a 1 u 5)
(v 1 a 10 u 5)
(v 6 a 1 u 1)
(v 5 a 1 u 4)
(v 3 a 3 u 5)
(v 2 a 5 u 5)

となる. このaの値の5,1,1,3,5,3,1,1,10,1,1,3,5が, √31の連分数の部分商である. √31=5+1/(1+1/(1+1/(3+1/(5+....)))).

アルゴリズムに戻って, 変数はv,a,uの3個(nは繰返し数). let*で設定するvの初期値は1, aは√dの整数部分, つまり最初の部分商, uもaであって,
v_n+1=(d-u_n²)/v_n,
a_n+1=floor((√d+u_n)/v_n+1),
u_n+1=a_n+1v_n+1-u_n
を繰り返す. aに√dがあるが, 整数部分をとるので下の方はいらない. vもuも整数のままで, 誤差が入らないから, いつまでも計算が続く.

この演習問題の次の設問は, ある回数から先は, 0<u_n<√d, 0<v_n<2√dを証明せよというのだ. 私は証明しようとも思わないが, そういうことはこの計算が周期的になるわけで, 平方根は無理数なのに, 連分数にすると繰り返すというのは面白い.

さて, このアルゴリズムは何をやっていたか調べてみる. √dをx₀, 最初の部分商をa₀とすると, (1)のように書ける. (TEXのプログラムでは, D,A,U,Vは大文字になっている.) ただし, u_-1, v₀は(0)のようだ. uの添字がvのより1だけ小さいのに注意. x=(√d+u)/vを各xの標準形とする.



(1)の添字をnだけ増やすと, 一般形x_nが(2)のように得られる. a_nを左辺に移すと, 1/x_n+1が(3)である. 上下を引っくり返せば(4)のx_n+1になるが, 分母に平方根があるので, 有理化し, 左上に√dが来るように, v_nで割る. これを(4)の最後のような標準形だと思うと, 対応する項から(5)の漸化式を得る. u_nを書いたら, 添字を1増やすとu_n+1が出来る. v_n+1も作れるが, u_nの式を入れると, vの漸化式になり, たしかにschemeで書いたプログラムが得られる.

というわけで, 今回もアルゴリズムが解明出来た.