Skip to content

The Generalized Tribonacci Sequence

Vector 17:2, page 131

mathematics and statistics performance

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

Abstract

In a previous paper we treated the Horadam Sequences, the sequences based on the second order linear recursion algorithm. In this paper we treat the Tribonacci Sequences, those based on the third order linear recursion algorithm. APL defined functions for generating these sequences are given. Benchmarks demonstrate that, just as for the Horadam Sequences, explicit calculation from the auxiliary equation is much more efficient in CPU time than the calculation based on recursion.

Introduction

In a previous paper (1) we treated the Generalized Fibonacci Sequences or Horadam Sequences (2), the sequences based on the second order linear recursion algorithm:

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

It was shown how explicit calculation of these sequences from the auxiliary equation is much more efficient in cpu time than generating the sequences from the recursion algorithm. In this paper we treat the Generalized Tribonacci Sequences, the sequences generated by the third order linear recursion algorithm.

Tribonacci Sequences

Originally, the Tribonacci Sequence was defined as the sequence generated by the third order linear recursion algorithm (3)(4):

\[\begin{align*} T_1 &= 1 \\ T_2 &= 1 \\ T_3 &= 2 \\ T_n &= T_{n-3} + T_{n-2} + T_{n-1} \qquad \text{(for } n > 3\text{)} \end{align*}\]

such that \(T_n = 1\ 1\ 2\ 4\ 7\ 13\ 24\ 44\ 81\ 149 \ldots\)

It may be shown that the sums of the ascending diagonals of the Wong and Maddockx triangle form the Tribonacci Sequence (5)(6).

It may look strange that \(\mathrm{T}_3\) is set to 2 instead of 1, but the underlying idea is the same as for the Fibonacci Sequence. The Fibonacci Sequence starts with the initial values \(F_1 = 1\) and \(F_2 = 1\). If we set \(F_0 = 0\) and \(F_1 = 1\), we get \(F_2 = F_0 + F_1 = 1\). In the same way, if we set \(T_0 = 0\) and \(T_1 = T_2 = 1\) we get \(T_3 = T_0 + T_1 + T_2 = 0 + 1 + 1 = 2\).

In most current literature (7) (8), the recursion algorithm as defined above is generalized to the form:

\[\begin{align*} T_1 &= f \\ T_2 &= g \\ T_3 &= h \\ T_n &= aT_{n-3} + bT_{n-2} + cT_{n-1} \quad \text{(for } n > 3\text{)} \end{align*}\]

where the parameters \(f\ g\ h\) and \(a\ b\ c\) are any reals. We assume however \(a \ne 0\), as for \(a = 0\) the sequence \(T_n\) becomes independent of the initial value \(T_1 = f\).

The question arises whether, such as for the Horadam Sequences, direct calculation of the Generalized Tribonacci Sequences from their auxiliary equation is much more efficient in cpu time than generating the sequences from the recursion algorithm.

The Auxiliary Equation

The auxiliary or characteristic equation associated with the third order linear recursion algorithm given above is the third degree or cubic equation:

\[\begin{equation*} t^3 - ct^2 - bt - a = 0 \end{equation*}\]

By linear transformation, this equation may be reduced to the form whose coefficient of the second degree term is zero:

\[\begin{equation*} t = s + \frac{c}{3} \end{equation*}\]
\[\begin{alignat*}{2} \text{such that}\;\; s^3 + 3ps + 2q = 0 \qquad & \text{with}\;\; && p = -\bigl(3b + c^2\bigr)/9 \\ & \text{and}\;\; && q = -\bigl(27a + 9bc + 2c^3\bigr)/54 \end{alignat*}\]

The algorithm for calculating the roots of the equation depends on the signum of the discriminant \(\Delta\):

\[\begin{equation*} \Delta = p^3 + q^2 \end{equation*}\]
  1. If \(\Delta < 0\), the equation has three different real roots:

    \[\begin{alignat*}{2} s_1 &= 2.r.\cos\frac{\varphi}{3} &\qquad& \text{with } r = \sqrt{-p} \\[1ex] s_2 &= 2.r.\cos\left(\frac{\varphi}{3} - 120^\circ\right) && \text{and } \varphi = \arccos\left(\frac{q}{pr}\right) \\[1ex] s_3 &= 2.r.\cos\left(\frac{\varphi}{3} + 120^\circ\right) \end{alignat*}\]
  2. If \(\Delta = 0\) and \(pq \ne 0\), the equation has a simple and a double root:

    \[\begin{equation*} s_1 = -2.q^{1/3} \qquad \text{and} \qquad s_2 = s_3 = q^{1/3} \end{equation*}\]
  3. If \(\Delta = 0\) and \(pq = 0\), the equation has the triple root 0:

    \[\begin{equation*} s_1 = s_2 = s_3 = 0 \end{equation*}\]
  4. If \(\Delta > 0\), the equation has a real root and two conjugate complex roots:

    \[\begin{alignat*}{3} s_1 &= u + v &\qquad s_2 &= s_r + i.s_i &\qquad s_3 &= s_r - i.s_i \\ && s_r &= -(u + v)/2 & s_i &= \sqrt{3}\bigl(\lvert u - v\rvert\bigr)/2 \end{alignat*}\]
    \[\begin{equation*} \text{with } u = \bigl(-q + \sqrt{\Delta}\bigr)^{1/3} \qquad \text{and } v = -\bigl(q + \sqrt{\Delta}\bigr)^{1/3} \end{equation*}\]

In fact, the case \(\Delta = 0\) may be treated by the same formula as for \(\Delta < 0\) and by the same formula as for \(\Delta > 0\). Treating the case separately however, is much simpler and more accurate. Finally it should be noticed that, as we assumed \(a \ne 0\) and \(a = t_1 \times t_2 \times t_3\), the roots \(t\) are always different from 0.

More details about the solution of the third degree or cubic equation may be found in (9). An APL defined function CUBIC, based on the algorithm described above but treating the case \(\Delta = 0\) with the same procedure as for \(\Delta > 0\) has been published in (10). This means that if \(\Delta = 0\), the value of \(s_i\) may be slightly different from 0, and the function CUBIC returns two conjugate complex roots instead of the double or triple root.

Explicit Calculation of \(T_n\)

Depending on the signum of the discriminant \(\Delta\) of the s-form of the auxiliary or characteristic equation, we have:

  1. If \(\Delta < 0\): \(\;T_n = xt_1^n + yt_2^n + zt_3^n\;\) with

    \[\begin{alignat*}{4} & xt_1 && {}+ yt_2 && {}+ zt_3 && {}= f \\ & xt_1^2 && {}+ yt_2^2 && {}+ zt_3^2 && {}= g \\ & xt_1^3 && {}+ yt_2^3 && {}+ zt_3^3 && {}= h \end{alignat*}\]
  2. If \(\Delta = 0\) and \(pq \ne 0\): \(\;T_n = xt_1^n + (y + nz)t_2^n\;\) with

    \[\begin{alignat*}{8} & xt_1 && {}+ (y + z)t_2 && {}= f &\qquad& \text{or} &\qquad& xt_1 && {}+ yt_2 && {}+ zt_2 && {}= f \\ & xt_1^2 && {}+ (y + 2z)t_2^2 && {}= g &&&& xt_1^2 && {}+ yt_2^2 && {}+ 2zt_2^2 && {}= g \\ & xt_1^3 && {}+ (y + 3z)t_2^3 && {}= h &&&& xt_1^3 && {}+ yt_2^3 && {}+ 3zt_2^3 && {}= h \end{alignat*}\]
  3. If \(\Delta = 0\) and \(pq = 0\): \(\;T_n = \bigl(x + ny + n^2z\bigr)t_1^n\;\) with

    \[\begin{alignat*}{5} & (x + y + z)t_1 && {}= f &\qquad& \text{or} &\qquad& x + y + z && {}= f/t_1 \\ & (x + 2y + 4z)t_1^2 && {}= g &&&& x + 2y + 4z && {}= g/t_1^2 \\ & (x + 3y + 9z)t_1^3 && {}= h &&&& x + 3y + 9z && {}= h/t_1^3 \end{alignat*}\]
  4. If \(\Delta > 0\): \(\;T_n = xt_1^n + (y\cos n\theta + z\sin n\theta)\rho^n\;\) with

    \[\begin{alignat*}{8} & xt_1 && {}+ (y\cos\theta && {}+ z\sin\theta)\rho && {}= f &\quad& \text{or} &\quad& xt_1 && {}+ y\rho\cos\theta + z\rho\sin\theta && {}= f \\ & xt_1^2 && {}+ (y\cos 2\theta && {}+ z\sin 2\theta)\rho^2 && {}= g &&&& xt_1^2 && {}+ y\rho^2\cos 2\theta + z\rho^2\sin 2\theta && {}= g \\ & xt_1^3 && {}+ (y\cos 3\theta && {}+ z\sin 3\theta)\rho^3 && {}= h &&&& xt_1^3 && {}+ y\rho^3\cos 3\theta + z\rho^3\sin 3\theta && {}= h \end{alignat*}\]
\[\begin{alignat*}{2} &\text{and}\quad & \rho &= \pm\sqrt{t_r^2 + t_i^2} \qquad (+ \text{ for } t_r \ge 0 \text{ and } - \text{ for } t_r < 0) \\ && \theta &= \arcsin(t_i/\rho) \end{alignat*}\]

Defined Functions

APL defined functions for calculating the sequence \(T_1 \ldots T_n\), based on implicit recursion (TRIB1), programming of a recursive loop (TRIB2), and explicit calculation from the auxiliary or characteristic equation (TRIB3) are given below:

    ∇ Z←L TRIB1 N
[1]   →(N≤3)/0,Z←(N⌊3)⍴L
[2]   Z←L TRIB1 N-1
[3]   Z←Z,+/(3↓L)ׯ3↑Z
    ∇

    ∇ Z←L TRIB2 N
[1]   →(N≤3)/0,Z←(N⌊3)⍴L
[2]  END:→(N>⍴Z←Z,+/(3↓L)ׯ3↑Z)/END
    ∇

    ∇ Z←L TRIB3 N;P;Q;∆;PHI;T;XYZ;UV;TR;TI;RHO;TETA
[1]   →(N≤3)/0,Z←(N⌊3)⍴L
[2]   P←-((3×L[5])+L[6]*2)÷9
[3]   Q←-((27×L[4])+L[6]×(9×L[5])+2×L[6]*2)÷54
[4]   →(0 1=×∆←(Q*2)+P*3)/∆Z,∆P
[5]   PHI←¯2○Q÷P×R←(-P)*0.5
[6]   T←(L[6]÷3)+2×R×2○(PHI+0,¯1 1×○2)÷3
[7]   XYZ←(3↑L)⌹⍉T∘.*⍳3
[8]   XYZ←(1E¯15<|XYZ)×XYZ
[9]   →0,Z←Z,(XYZ[1]×T[1]*N)+(XYZ[2]×T[2]*N)+XYZ[3]×T[3]*N←3↓⍳N
[10] ∆Z:→(Q=0)/QZ
[11]  T←(L[6]÷3)+(¯2 1 1)×Q*÷3
[12]  XYZ←(3↑L)⌹((⍳3)∘.*0 0 1)×⍉T∘.*⍳3
[13]  XYZ←(1E¯15<|XYZ)×XYZ
[14]  →0,Z←Z,(XYZ[1]×T[1]*N)+(XYZ[2]+N×XYZ[3])×T[2]*N←3↓⍳N
[15] QZ:T←L[6]÷3
[16]  XYZ←((3↑L)÷T*⍳3)⌹(⍳3)∘.*0 1 2
[17]  XYZ←(1E¯15<|XYZ)×XYZ
[18]  →0,Z←Z,(XYZ[1]+N×XYZ[2]+N×XYZ[3])×T*N←3↓⍳N
[19] ∆P:UV←(×UV)×(|UV←-Q+(¯1 1)×∆*0.5)*÷3
[20]  T←(L[6]÷3)++/UV
[21]  TR←(L[6]÷3)-0.5×+/UV
[22]  TI←(0.75*0.5)×|-/UV
[23]  RHO←(TR÷|TR)×((TR*2)+TI*2)*0.5
[24]  TETA←¯1○TI÷RHO
[25]  XYZ←(3↑L)⌹⍉(T*⍳3),[1](2 1∘.○TETA×⍳3)×2 3⍴RHO*⍳3
[26]  XYZ←(1E¯15<|XYZ)×XYZ
[27]  Z←Z,(XYZ[1]×T*N)+((XYZ[2]×2○TETA×N)+XYZ[3]×1○TETA×N)×RHO*N←3↓⍳N
    ∇

      L←1 1 2 1 1 1
      L TRIB1 10
1 1 2 4 7 13 24 44 81 149
      L TRIB2 10
1 1 2 4 7 13 24 44 81 149
      L TRIB3 10
1 1 2 4 7 13 24 44 81 149

The left argument L is a six-element vector containing the initial values \(f\ g\ h\) and the parameters \(a\ b\ c\). 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 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. Results are not rounded off, as the initial values and parameters may be any reals.

TRIB3 is a literal translation in APL of the algorithm for the calculation of the roots of the auxiliary or characteristic equation, and of the algorithm for the calculation of the Tribonacci numbers from those roots, as given in the two previous sections. Coefficients \(x\ y\ z\) are set to zero when their calculated absolute value is smaller than an arbitrarily chosen tolerance. On an experimental basis this tolerance has been set to 1E¯15. Calculation of the terms of the sequence from those coefficients may be simplified by using the outer product, but at considerable cost of performance in cpu time.

Standard ISO APL requires that for the dyadic function power A*B, when the left argument A is a negative number and the right argument B is not an integer, a domain error report is generated, this to allow complex arithmetic to be a consistent extension. As some implementations conform to this requirement, we have coded line 19 of the function TRIB3 :

[19] ∆P:UV←(×UV)×(|UV←-Q+(¯1 1)×∆*0.5)*÷3

instead of

[19] ∆P:UV←(-Q+(¯1 1)×∆*0.5)*÷3

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

    ∇ Z←L TRIB4 N;P;Q;∆;PHI;T;XYZ;UV;TR;TI;RHO;TETA
[1]   →(N≤3)/0,Z←L[N⌊3]
[2]   P←-((3×L[5])+L[6]*2)÷9
[3]   Q←-((27×L[4])+L[6]×(9×L[5])+2×L[6]*2)÷54
[4]   →(0 1=×∆←(Q*2)+P*3)/∆Z,∆P
[5]   PHI←¯2○Q÷P×R←(-P)*0.5
[6]   T←(L[6]÷3)+2×R×2○(PHI+0,¯1 1×○2)÷3
[7]   XYZ←(3↑L)⌹⍉T∘.*⍳3
[8]   XYZ←(1E¯15<|XYZ)×XYZ
[9]   →0,Z←+/XYZ×T*N
[10] ∆Z:→(Q=0)/QZ
[11]  T←(L[6]÷3)+(¯2 1 1)×Q*÷3
[12]  XYZ←(3↑L)⌹((⍳3)∘.*0 0 1)×⍉T∘.*⍳3
[13]  XYZ←(1E¯15<|XYZ)×XYZ
[14]  →0,Z←+/XYZ×(1 1,N)×T*N
[15] QZ:T←L[6]÷3
[16]  XYZ←((3↑L)÷T*⍳3)⌹(⍳3)∘.*0 1 2
[17]  XYZ←(1E¯15<|XYZ)×XYZ
[18]  →0,Z←(+/XYZ×1,N,N*2)×T*N
[19] ∆P:UV←(×UV)×(|UV←-Q+(¯1 1)×∆*0.5)*÷3
[20]  T←(L[6]÷3)++/UV
[21]  TR←(L[6]÷3)-0.5×+/UV
[22]  TI←(0.75*0.5)×|-/UV
[23]  RHO←(×TR)×((TR*2)+TI*2)*0.5
[24]  TETA←¯1○TI÷RHO
[25]  XYZ←(3↑L)⌹⍉(T*⍳3),[1](2 1∘.○TETA×⍳3)×2 3⍴RHO*⍳3
[26]  XYZ←(1E¯15<|XYZ)×XYZ
[27]  Z←+/XYZ×(1,2 1○TETA×N)×(T,2⍴RHO)*N
    ∇

      L←1 1 2 1 1 1
      L TRIB4 10
149

      (⊂L) TRIB4¨⍳10
1 1 2 4 7 13 24 44 81 149

The right argument is rounded off. If this rounded off value is less than or equal to 0, an error report is generated. Calculation of \(T_n\) from the coefficients \(x\ y\ z\) has been simplified by using reduction, the impact on performance in cpu time being negligible.

Illustrations

Illustrations are given, using TRIB3, for all possible combinations of the signum of the discriminant \(\Delta\) and of the signum of \(p\) and \(q\):

\(a\) \(b\) \(c\) \(\Delta\) \(p\) \(q\)
2 1 ¯2 - - -
¯1 ¯1 3 - - 0
¯2 1 2 - - +
2 ¯5 4 0 - -
¯3 5 ¯1 0 - +
1 ¯3 3 0 0 0
1 1 1 + - -
¯1 ¯1 2 + - +
2 ¯3 3 + 0 -
¯1 ¯3 3 + 0 +
2 ¯2 1 + + -
4 ¯6 3 + + 0
¯2 ¯2 ¯2 + + +
      L←1 2 3 2 1 ¯2
      L TRIB3 10
1 2 3 ¯2 11 ¯18 43 ¯82 171 ¯338

      L←1 2 3 ¯1 ¯1 3
      L TRIB3 10
1 2 3 6 13 30 71 170 409 986

      L←1 2 3 ¯2 1 2
      L TRIB3 10
1 2 3 6 11 22 43 86 171 342

      L←1 2 3 2 ¯5 4
      L TRIB3 10
1 2 3 4 5 6 7 8 9 10

      L←1 2 3 ¯3 5 ¯1
      L TRIB3 10
1 2 3 4 5 6 7 8 9 10

      L←1 2 3 1 ¯3 3
      L TRIB3 10
1 2 3 4 5 6 7 8 9 10

      L←1 2 3 1 1 1
      L TRIB3 10
1 2 3 6 11 20 37 68 125 230

      L←1 2 3 ¯1 ¯1 2
      L TRIB3 10
1 2 3 3 1 ¯4 ¯12 ¯21 ¯26 ¯19

      L←1 2 3 2 ¯3 3
      L TRIB3 10
1 2 3 5 10 21 43 86 171 341

      L←1 2 3 ¯1 ¯3 3
      L TRIB3 10
1 2 3 2 ¯5 ¯24 ¯59 ¯100 ¯99 62

      L←1 2 3 2 ¯2 1
      L TRIB3 10
1 2 3 1 ¯1 3 7 ¯1 ¯9 7

      L←1 2 3 4 ¯6 3
      L TRIB3 10
1 2 3 1 ¯7 ¯15 1 65 129 1

      L←1 2 3 ¯2 ¯2 ¯2
      L TRIB3 10
1 2 3 ¯12 14 ¯10 16 ¯40 68 ¯88

Benchmarks

Benchmarks have been done for \(\mathrm{N} = 250(250)1000\), in order to compare the effectiveness in cpu time for the functions defined.

Results are given below (in ms).

N 250 500 750 1000
TRIB1 97 219 365 564
TRIB2 68 155 262 390
TRIB3 \((\Delta < 0)\) 7.1 11.2 15.3 19.5
TRIB3 \((\Delta = 0 \text{ and } pq \ne 0)\) 6.3 9.2 12.0 14.8
TRIB3 \((\Delta = 0 \text{ and } p{=}q{=}0)\) 4.8 6.6 8.5 10.3
TRIB3 \((\Delta > 0)\) 9.0 14.3 19.6 24.9
TRIB4 3.8 3.8 3.8 3.8

On average, cpu times are 30% lower for TRIB2 than for TRIB1. Depending on N and on the signum of the discriminant \(\Delta\) the efficiency factor in cpu time is from 8 to 38 for TRIB3 versus TRIB2. Figures are comparable to those reported for the Generalized Fibonacci or Horadam Sequences (1).

CPU times for TRIB4 are based on a weighted average, as the impact of the signum of the discriminant \(\Delta\) is negligible. The efficiency factor in cpu time is from 18 to 105 for TRIB4 versus TRIB2, from 2.4 to 2.7 versus TRIB3. Efficiency factors are smaller than those reported for the Generalized Fibonacci or Horadam Numbers (1), due to the overhead while computing the discriminant, but are still considerable.

Functions defined conform to the standard ISO APL (11). Benchmarks have been done on a MicroLine Pentium S-100, with mathematical co-processor, using Dyalog APL/W Version 8.0.8 under Windows 95 (12). Benchmarks have been done using the defined function ⎕MONITOR.

References

  1. J. De Kerf; From Fibonacci to Horadam; Vector Vol. 15 No.3, January 1999, pp. 122-130. See also: C.A. Reiter; Exact Horadam Numbers with a Chebyshevish Accent; Vector Vol.16 No.1, July 1999, pp. 122-131.
  2. A.F. Horadam; Basic Properties of a Certain Generalized Sequence of Numbers; The Fibonacci Quarterly, Vol. 3, No. 3, 1965, pp. 161-176.
  3. M. Feinberg; Fibonacci-Tribonacci; The Fibonacci Quarterly, Vol. 1, No. 1, 1963, pp. 71-74.
  4. A. Scott, T. Delaney, and V.E. Hoggatt; The Tribonacci Sequence; The Fibonacci Quarterly, Vol. 15, No. 3, 1977, pp. 193-200.
  5. C.K. Wong and T.W. Maddockx; A Generalized Pascal’s Triangle; The Fibonacci Quarterly, Vol. 13, No. 2, 1975, pp. 134-136.
  6. J. De Kerf; Some Other Triangles and APL; Les Nouvelles d’APL, No. 30, Juillet 1999, pp. 71-99.
  7. D. Jarden; Third order recurring sequences; In "Recurring Sequences" — A Collection of Papers; Reveon Lematematika, Jerusalem, Israel, 1966, pp. 86-89.
  8. S. Rabinowitz; Algorithmic Manipulation of Third-Order Linear Recurrences; The Fibonacci Quarterly, Vol. 34, No. 5, 1996, pp. 447-464.
  9. J.V. Uspensky; Theory of Equations — Chapter V: Cubic and Biquadratic Equations; McGraw-Hill Book Company, New York, New York, 1948; pp. 82-98.
  10. U. Grenander; Mathematical Experiments on the Computer; Pure and Applied Mathematics Series; Academic Press, New York, New York, 1982; pp. 238-239.
  11. International Standard ISO 8485 : 1989 (E), First Edition; 1989—11—01; Programming Languages — APL; International Organization for Standardization ISO, Geneva, Switzerland, 1989.
  12. Dyalog APL/W Reference Manual — Version 8; Dyadic Systems Limited, Basingstoke, Hampshire, United Kingdom, 1996.