５．Lucas列の計算方法

y_nの値を直接計算する方法を考える。

y_nについては、以下のような式が成り立つ。

(1) y2n＝2ynyn+1−ayn2
(2) y2n+1＝yn+12＋byn2
(3) y2n+2＝ayn+12＋2bynyn+1

n番目と n+1番目から2n番目、2n+1番目、2n+2番目が生成されるということで、
以前、べき乗計算のときに使った、ｎを２進数に展開して、各桁のビットに応じて処理を行う、
という手法が使えそうである。

23を例に考えてみよう。
奇数番目を生成する式は、(2)しかない。

y₂₃＝y₁₂＋y₁₁

つまり、12番目と11番目を計算しておく必要がある。

次に、11番目は奇数番目なので、生成する式は、(2)しかない。

y₁₁＝y₆＋y₅

12番目をこの６番目と５番目から生成するには、(3)の式を使えばよい。

次に、５番目は奇数番目なので、生成する式は、(2)しかない。

y₅＝y₃＋y₂

６番目をこの３番目と２番目から生成するには、(3)の式を使えばよい。

次に、３番目は奇数番目なので、生成する式は、(2)しかない。

y₃＝y₂＋y₁

２番目は、０番目と１番目から生成されるが、(1)式より、２番目と１番目からも計算できる。

以上の流れをまとめると、以下のようになる。

w0：０→(2)→１→(1)→２→(2)→５→(2)→11→(2)→23
w1：１→(3)→２→(2)→３→(3)→６→(3)→12→(3)→24

23＝(10111)₂と、式の選択を見比べると、以下のアルゴリズムでｎ番目を求められることがわかる。

1. ｎを２進展開する。y₀＝０、y₁＝１とする
2. 最上位ビットから順に見ていき、

• １なら w₀は(2)式、w₁は(3)式より求める
• ０なら w₀は(1)式、w₁は(2)式より求める

3. w₀が求める解

プログラムは以下のとおり。

 10   ' lucas sequence
 20   ' y(0)=0, y(1)=1, y(n+1)=ay(n)+by(n-1)
 30   fnYn1(N)
 40   local I,A,B,W0,W1,Y0,Y1
 50   dim W%(len(N))
 60   if N=0 then return(0)
 70   if N=1 then return(1)
 80   I=0
 90   while N>0
100     inc I:W%(I)=N@2:N=N\2
110   wend
120   W%(0)=I
130   '
140   A=1:B=1:Y0=0:Y1=1
150   for I=W%(0) to 1 step -1
160     if W%(I)=1 then W0=Y1^2+B*Y0^2:W1=A*Y1^2+2*B*Y0*Y1
170     if W%(I)=0 then W0=2*Y0*Y1-A*Y0^2:W1=Y1^2+B*Y0^2
180     Y0=W0:Y1=W1
190   next I
200   return(Y0)

０番目、１番目に対する例外処理は不要である。
また、UBASICでは、bit(i,n)により、ｎを２進展開したときの最下位からｉ番目のビットを取得できる。

これらを使って書き直すと以下のようになる。

 10   ' lucas sequence
 20   ' y(0)=0, y(1)=1, y(n+1)=ay(n)+by(n-1)
 30   fnYn2(N)
 40   local I,A,B,W0,W1,Y0,Y1
 50   dim W%(len(N))
 60   A=1:B=1:Y0=0:Y1=1
 70   for I=len(N)-1 to 0 step -1
 80     if bit(I,N)=1 then W0=Y1^2+B*Y0^2:W1=A*Y1^2+2*B*Y0*Y1
 90     if bit(I,N)=0 then W0=2*Y0*Y1-A*Y0^2:W1=Y1^2+B*Y0^2
100     Y0=W0:Y1=W1
110   next I
120   return(Y0)