From Fibonacci to Horadam

by Joseph De Kerf

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 122

Abstract¶

Details are given about the generalization of the Fibonacci sequence as investigated by A. F. Horadam. Defined functions generating these so-called Horadam sequences, respectively based on implicit recursion, programming of a recursive loop, and explicit calculation from the auxiliary equation are given. Benchmarks demonstrate that explicit calculation is much more efficient in CPU time than the calculation based on recursion.

Fibonacci Sequences¶

The Fibonacci sequence is a special case of the sequence [1]

\[\begin{align*} T_1 &= p \\ T_2 &= q \\ \rlap{\ldots}\phantom{T_n} & \\ T_n &= T_{n-2} + T_{n-1} \qquad (\text{for } n > 2) \end{align*}\]

the domain of \(p\) and \(q\) being the set of integers. For \(p = 1\) and \(q = 1\), the sequence is called the Fibonacci sequence \(F_n\). For \(p = 1\) and \(q = 3\), the sequence is known in the literature as the Lucas sequence \(L_n\). It may be shown that:

\[\begin{equation*} T_n = pF_{n-2} + qF_{n-1} \end{equation*}\]

The sequences may be generated by application of implicit recursion or by programming a recursive loop. J. Sullivan [2] has shown that, as far as the Fibonacci sequence is concerned, explicit calculation based on the golden section formula of Binet, is much more efficient in CPU time. J. De Kerf [3] has indicated how, apart from some overhead, by calculating only the first term of the golden section formula and rounding off, the efficiency may be speeded up by an extra factor of 2:

\[\begin{equation*} F_n = \frac{1}{\sqrt{5}} \left( \frac{1+\sqrt{5}}{2} \right)^{\!n} \qquad (\text{rounded off}) \end{equation*}\]

F←⌊0.5+(÷R)×(0.5×1+R←5*0.5)*⍳N

APL defined functions based on implicit recursion (FLT1), programming of a recursive loop (FLT2), and explicit calculation from \(F_{n-2}\) and \(F_{n-1}\) as described above (FLT3), are given below:

[0]   Z←P FLT1 N
[1]   →(N≤2)/0,Z←(2⌊N)⍴P
[2]   Z←P FLT1 N-1
[3]   Z←Z,+/¯2↑Z

[0]   Z←P FLT2 N
[1]   →(N≤2)/0,Z←(2⌊N)⍴P
[2]  END:→(N>⍴Z←Z,+/¯2↑Z)/END

[0]   Z←P FLT3 N;F;R
[1]   →(N≤2)/0,Z←(2⌊N)⍴P
[2]   F←⌊0.5+(÷R)×(0.5×1+R←5*0.5)*⍳N-1
[3]   Z←Z,(P[1]×¯1↓F)+P[2]×1↓F

The left argument P is a two-element vector containing the parameters \(p\) and \(q\). The right argument N is the non-negative integer scalar or one-element vector specifying the number of terms of the sequence required. The explicit result Z is a vector containing those terms. If N is 0, the empty numeric vector is returned. If N is a negative integer or not an integer, an error report is generated.

Examples:

      1 3 FLT1 10
1 3 4 7 11 18 29 47 76 123
      1 3 FLT2 10
1 3 4 7 11 18 29 47 76 123
      1 3 FLT3 10
1 3 4 7 11 18 29 47 76 123

      3 1 FLT1 10
3 1 4 5 9 14 23 37 60 97
      3 1 FLT2 10
3 1 4 5 9 14 23 37 60 97
      3 1 FLT3 10
3 1 4 5 9 14 23 37 60 97

the first of the two examples given being the Lucas sequence.

Horadam Sequences¶

A much wider generalization of the Fibonacci sequence has been thoroughly discussed by A.F. Horadam [4]:

\[\begin{align*} H_1 &= p \\ H_2 &= q \\ \rlap{\ldots}\phantom{T_n} & \\ H_n &= rH_{n-2} + sH_{n-1} \qquad\qquad (\text{for } n > 2) \end{align*}\]

Originally, the domain of the parameters \(p\ q\ r\ s\) was restricted to integers. In most current publications however, the domain is extended to the set of real numbers. Specific sequences from the literature are:

Sequence	\(p\)	\(q\)	\(r\)	\(s\)
Fibonacci (\(F_n\))	1	1	1	1
Lucas (\(L_n\))	1	3	1	1
Pell (\(P_n\))	1	2	1	2
Pell-Lucas (\(Q_n\))	1	3	1	2

\[\begin{alignat*}{2} &\text{such that:}\quad & F_n &= 1\ 1\ 2\ 3\ 5\ 8\ 13\ 21\ 34\ 55 \\ && L_n &= 1\ 3\ 4\ 7\ 11\ 18\ 29\ 47\ 76\ 123 \\ && P_n &= 1\ 2\ 5\ 12\ 29\ 70\ 169\ 408\ 985\ 2378 \\ && Q_n &= 1\ 3\ 7\ 17\ 41\ 99\ 239\ 577\ 1393\ 3363 \end{alignat*}\]

To calculate the sequence explicitly we may solve the auxiliary equation \(h^2 - sh - r = 0\), the solution depending on the signum of the discriminant \(s^2 + 4r\) [5]:

\[\begin{alignat*}{2} &\text{(a)}\quad s^2 + 4r > 0\colon\quad H_n = ca^{n-1} + db^{n-1} &\quad\text{with}\quad &\begin{aligned}[t] \Delta &= \sqrt{s^2 + 4r} \\ a &= (s + \Delta)/2 \\ b &= (s - \Delta)/2 \\ c &= (q - pb)/\Delta \\ d &= (pa - q)/\Delta \end{aligned} \\[1ex] &\text{(b)}\quad s^2 + 4r = 0\colon\quad H_n = (c + d)(s/2)^{n-2} &\quad\text{with}\quad &\begin{aligned}[t] c &= ps - q \\ d &= n(2q - ps)/2 \end{aligned} \\[1ex] &\text{(c)}\quad s^2 + 4r < 0\colon\quad H_n = (c + d)\rho^{n-1} &\quad\text{with}\quad &\begin{aligned}[t] \Delta &= \sqrt{\lvert s^2 + 4r \rvert} \\ \rho &= \pm\sqrt{-r} \qquad (+ \text{ for } s \ge 0 \text{ and } - \text{ for } s < 0) \\ \theta &= \arcsin(\Delta/2\rho) \\ c &= -(2p\rho.\sin(n - 2)\theta)/\Delta \\ d &= (2q.\sin(n - 1)\theta)/\Delta \end{aligned} \end{alignat*}\]

Defined Functions¶

APL defined functions for calculating the sequence \(H_1 \ldots H_n\), based on implicit recursion (HOR1), programming of a recursive loop (HOR2), and explicit calculation from the auxiliary equation (HOR3) are given below:

[0]   Z←P HOR1 N
[1]   →(N≤2)/0,Z←(2⌊N)⍴P
[2]   Z←P HOR1 N-1
[3]   Z←Z,+/(2↓P)×¯2↑Z

[0]   Z←P HOR2 N
[1]   →(N≤2)/0,Z←(2⌊N)⍴P
[2]  END:→(N>⍴Z←Z,+/(2↓P)×¯2↑Z)/END

[0]   Z←P HOR3 N;A;B;C;D;R;T;∆
[1]   →(N≤2)/0,Z←(2⌊N)⍴P
[2]   →((P[4]*2)>¯4×P[3])/LABDP
[3]   →((P[4]*2)<¯4×P[3])/LABDN
[4]   C←(P[1]×P[4])-P[2]
[5]   D←(P[2]-P[1]×P[4]÷2)×N←2↓⍳N
[6]   →0,Z←Z,(C≠-D)×(C+D)×(P[4]÷2)*N-2
[7]  LABDP:∆←((P[4]*2)+4×P[3])*0.5
[8]   B←P[4]-A←0.5×P[4]+∆
[9]   C←(P[2]≠P[1]×B)×(P[2]-P[1]×B)÷∆
[10]  D←(P[2]≠P[1]×A)×((P[1]×A)-P[2])÷∆
[11]  →0,Z←Z,(C×A*N)+D×B*N←1↓⍳N-1
[12] LABDN:∆←(|(P[4]*2)+4×P[3])*0.5
[13]  R←((-P[3])*0.5)×¯1*P[4]<0
[14]  T←¯1○∆÷2×R
[15]  C←-(2×P[1]×R×1○T×N←⍳N-2)÷∆
[16]  D←(2×P[2]×1○T×N+1)÷∆
[17]  Z←Z,(C≠-D)×(C+D)×R*N+1

The left argument P is a four-element vector containing the parameters \(p\ q\ r\ s\). The right argument N is the non-negative integer scalar or one-element vector specifying the number of terms of the sequence requested. The explicit result is a vector containing those terms. If N is 0, the empty numeric vector is returned. If N is a negative integer or not an integer, an error report is generated. Results are not rounded off, as the parameters may be any reals. Examples:

      1 2 1 2 HOR1 10
1 2 5 12 29 70 169 408 985 2378
      1 2 1 2 HOR2 10
1 2 5 12 29 70 169 408 985 2378
      1 2 1 2 HOR3 10
1 2 5 12 29 70 169 408 985 2378

      1 3 1 2 HOR1 10
1 3 7 17 41 99 239 577 1393 3363
      1 3 1 2 HOR2 10
1 3 7 17 41 99 239 577 1393 3363
      1 3 1 2 HOR3 10
1 3 7 17 41 99 239 577 1393 3363

The two examples given being the Pell sequence and the Pell-Lucas sequence respectively.

With HOR3 some pitfalls are cured. If \(s^2 + 4r\) is 0, the calculated value (P[4]*2)+4×P[3] may be slightly different from 0. If compared with 0, the program would jump to LABDP if the approximation is positive and to LABDN if the approximation is negative. To avoid this, in lines 2 and 3 the values of P[4]*2 and ¯4×P[3] are compared. In addition, depending on the signum of \(s^2 + 4r\):

(a) \(s^2 + 4r > 0\). If \(c = 0\) and \(\lvert b\rvert < 1\) or \(d = 0\) and \(\lvert a\rvert < 1\), the sequence \(H_n\) converges to 0. If \(c = 0\) or \(d = 0\) however, their calculated value may be slightly different from 0, such that for \(\lvert a\rvert > 1\) and \(\lvert b\rvert > 1\) respectively, the calculated sequence \(H_n\) converges to \(\pm\infty\). Correction is done in line 9 by setting C=0 when P[2]=P[1]×B and in line 10 by setting D=0 when P[2]=P[1]×A.

Examples:

p=1.24 q=-0.558 r=0.81 s=1.35
(c=0 a=1.80 b=-0.45)
p=1.24 q=0.558 r=0.81 s=-1.35
(d=0 a=0.45 b=-1.80)

(b) \(s^2 + 4r = 0\). If \(n\) is a positive integer such that \(c + d = 0\), then \(H_n = 0\) for that value of \(n\). The calculated value of \(c + d\) however, may be slightly different from 0, and if \(\lvert s/2\rvert > 1\) the value of \(H_n\) returned may be considerably different from 0. Correction is done in line 6 by setting the corresponding value of the explicit result Z to 0 when C = -D.

Examples:

p=0.6  q=0.98 r=-2.7225 s=3.3
p=-0.6 q=0.98 r=-2.7225 s=-3.3
(|s/2|=1.65 and c+d=0 for n=100)

(c) \(s^2 + 4r < 0\). If for specific values of \(n\) the sum \(c + d = 0\), then \(H_n = 0\) for those values of \(n\). The calculated value of \(c + d\) however, may be slightly different from 0, and, if \(\lvert\rho\rvert > 1\), the value of \(H_n\) returned may be considerably different from 0, the error increasing in absolute value with \(n\). Correction is done in line 17 by setting the corresponding value of Z to 0 when C = -D.

Examples:

p=2 q=3 r=-3 s=3
p=2 q=-3 r=-3 s=-3
(|ρ| = √3 and c+d=0 for n=4(6)...)

Finally, if the value of \(H_n\) is only requested for one specific value of \(n\), when applying recursion, the entire sequence \(H_1 \ldots H_n\) must be calculated. When applying explicit calculation from the auxiliary equation however, \(H_n\) must only be calculated for that specific value of \(n\). A modified version HOR4 of the defined function HOR3, that only calculates \(H_n\) for a given non-negative scalar or one-element vector is given below:

[0]   Z←P HOR4 N;A;B;C;D;R;T;∆
[1]   →(N≤2)/0,Z←P[2⌊N]
[2]   →((P[4]*2)>¯4×P[3])/LABDP
[3]   →((P[4]*2)<¯4×P[3])/LABDN
[4]   C←(P[1]×P[4])-P[2]
[5]   D←(P[2]-P[1]×P[4]÷2)×N
[6]   →0,Z←(C≠-D)×(C+D)×(P[4]÷2)*N-2
[7]  LABDP:∆←((P[4]*2)+4×P[3])*0.5
[8]   B←P[4]-A←0.5×P[4]+∆
[9]   C←(P[2]≠P[1]×B)×(P[2]-P[1]×B)÷∆
[10]  D←(P[2]≠P[1]×A)×((P[1]×A)-P[2])÷∆
[11]  →0,Z←(C×A*N-1)+D×B*N-1
[12] LABDN:∆←(|(P[4]*2)+4×P[3])*0.5
[13]  R←((-P[3])*0.5)×¯1*P[4]<0
[14]  T←¯1○∆÷2×R
[15]  C←-(2×P[1]×R×1○T×N-2)÷∆
[16]  D←(2×P[2]×1○T×N-1)÷∆
[17]  Z←(C≠-D)×(C+D)×R*N-1

Note: Rather unknown is the fact that applying recursion, as done in the defined functions HOR1 and HOR2, also gives erroneous results for the special cases described as pitfalls above. Results for those special cases are compared with the exact solutions as returned by HOR3 in the accompanying Appendix. It would be a challenge to incorporate the corrections needed such as is done in the defined functions HOR3 and HOR4 (if possible!).

Benchmarks¶

Benchmarks have been done for N = 250(250)1000, in order to compare the effectiveness in CPU times for the functions defined. Results are given below (in ms):

`N`	250	500	750	1000
`HOR1`	98.6	217.6	359.2	524.7
`HOR2`	64.6	145.2	248.0	366.8
`HOR3` (\(s^2 + 4r > 0\))	4.12	6.61	9.26	12.12
`HOR3` (\(s^2 + 4r = 0\))	3.77	5.77	7.82	10.18
`HOR3` (\(s^2 + 4r < 0\))	7.60	12.91	18.55	24.16
`HOR4` (\(s^2 + 4r > 0\))	1.99	1.99	1.99	1.99
`HOR4` (\(s^2 + 4r = 0\))	1.72	1.72	1.72	1.72
`HOR4` (\(s^2 + 4r < 0\))	2.14	2.14	2.14	2.14

Depending on N, the CPU times are from 34% to 30% lower for HOR2 than for HOR1. Depending on N and the signum of \(s^2 + 4r\), the efficiency factor in CPU time is from 9 to 36 for HOR3 versus HOR2, and from 30 to 213 for HOR4 versus HOR2. Of course, those figures depend on the hardware and software used.

Functions defined conform to the standard ISO APL [6]. Programming, debugging, illustrations and benchmarks have been done on a MicroLine Pentium S-100, with mathematical co-processor, using Dyalog APL/W Version 7.1.2, under Windows 3.11 [7]. Comparison Tolerance was set to 1E¯13, the default value being 1E¯14. Benchmarks have been done using the system function QdMONITOR.

References¶

V.E. Hoggatt: Fibonacci and Lucas Numbers. Houghton Mifflin Company, Boston, Mass. 1969
J. Sullivan: Fast Fibbing. Vector Vol.3 No.1, July 1986, pp 113–114
J. De Kerf: More on "Fast Fibbing". Vector Vol.3 No.4, April 1987, p 120
A.F. Horadam: Basic Properties of a Certain Generalized Sequence of Numbers. The Fibonacci Quarterly, Vol.3 No.3, August 1965, pp 161–176
A. Brousseau: Linear Recursion and Fibonacci Sequences. The Fibonacci Association, San Jose State College, San Jose, California, 1971
International Standard ISO 8485: 1989 (E), First Edition; 1989-11-01; Programming Languages - APL. International Organization for Standardization ISO, Geneva, Switzerland, 1989
Dyalog APL/W Reference Manual - Version 7. Dyadic Systems Ltd., Basingstoke, 1994

Appendix¶

Results for the special cases described in the text, returned by HOR1 and HOR2, are compared with the exact solutions returned by HOR3. Results shown are those for Dyalog APL/W [7].

(a) \(s^2 + 4r > 0\)

      P←1.24 ¯0.558 0.81 1.35
      R1←¯5↑P HOR1 502
      R2←¯5↑P HOR2 502
      R3←¯5↑P HOR3 502
      (¯5↑⍳502),R1,R2,[1.5]R3
498 ¯2.359343255E110 ¯2.359343255E110 ¯5.495927961E¯173
499 ¯4.246817860E110 ¯4.246817860E110  2.473167583E¯173
500 ¯7.644272148E110 ¯7.644272148E110 ¯1.112925412E¯173
501 ¯1.375968987E111 ¯1.375968987E111  5.008164355E¯174
502 ¯2.476744176E111 ¯2.476744176E111 ¯2.253673960E¯174

      P←1.24 0.558 0.81 ¯1.35
      R1←¯5↑P HOR1 502
      R2←¯5↑P HOR2 502
      R3←¯5↑P HOR3 502
      (¯5↑⍳502),R1,R2,[1.5]R3
498  2.359343255E110  2.359343255E110 5.495927961E¯173
499 ¯4.246817860E110 ¯4.246817860E110 2.473167583E¯173
500  7.644272148E110  7.644272148E110 1.112925412E¯173
501 ¯1.375968987E111 ¯1.375968987E111 5.008164355E¯174
502  2.476744176E111  2.476744176E111 2.253673960E¯174

Results returned by HOR1 and HOR2 converge to \(\pm\infty\) and are erroneous. Results returned by HOR3 converge to 0 and are correct.

(b) \(s^2 + 4r = 0\)

      P←0.6 0.98 ¯2.7225 3.3
      R1←¯5↑P HOR1 102
      R2←¯5↑P HOR2 102
      R3←¯5↑P HOR3 102
      (¯5↑⍳102),R1,R2,[1.5]R3
 98  1.511780155E19  1.511780155E19  1.511780155E19
 99  1.247218628E19  1.247218628E19  1.247218628E19
100 ¯1.005223936E9  ¯1.005223936E9   0.000000000E0
101 ¯3.395552715E19 ¯3.395552715E19 ¯3.395552715E19
102 ¯1.120532396E20 ¯1.120532396E20 ¯1.120532396E20

      P←¯0.6 0.98 ¯2.7225 ¯3.3
      R1←¯5↑P HOR1 102
      R2←¯5↑P HOR2 102
      R3←¯5↑P HOR3 102
      (¯5↑⍳102),R1,R2,[1.5]R3

 98  1.511780155E19  1.511780155E19  1.511780155E19
 99 ¯1.247218628E19 ¯1.247218628E19 ¯1.247218628E19
100 ¯1.005223936E9  ¯1.005223936E9   0.000000000E0
101  3.395552715E19  3.395552715E19  3.395552715E19
102 ¯1.120532396E20 ¯1.120532396E20 ¯1.120532396E20

Results returned by HOR1 and HOR2 for N = 100 are considerably different from 0 and are erroneous. Results returned by HOR3 for N = 100 are 0 and are correct.

(c) \(s^2 + 4r < 0\)

      P←2 3 ¯3 3
      R1←¯5↑P HOR1 252
      R2←¯5↑P HOR2 252
      R3←¯5↑P HOR3 252
      (¯5↑⍳252),R1,R2,[1.5]R3

248 ¯1.455578343E59 ¯1.455578343E59 ¯1.455578343E59
249 ¯1.455578343E59 ¯1.455578343E59 ¯1.455578343E59
250 ¯3.568119232E44 ¯3.568119232E44  0.000000000E0
251  4.366735029E59  4.366735029E59  4.366735029E59
252  1.310020509E60  1.310020509E60  1.310020509E60

      P←2 ¯3 ¯3 ¯3
      R1←¯5↑P HOR1 252
      R2←¯5↑P HOR2 252
      R3←¯5↑P HOR3 252
      (¯5↑⍳252),R1,R2,[1.5]R3

248  1.455578343E59  1.455578343E59  1.455578343E59
249 ¯1.455578343E59 ¯1.455578343E59 ¯1.455578343E59
250  3.568119232E44  3.568119232E44  0.000000000E0
251  4.366735029E59  4.366735029E59  4.366735029E59
252 ¯1.310020509E60 ¯1.310020509E60 ¯1.310020509E60

Results returned by HOR1 and HOR2 for N = 250 are totally different from 0 and are erroneous. Results returned by HOR3 for N = 250 are 0 and are correct.