tribonacci_closed_form.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# License: GPLv3
# Date: 20 January 2017
# https://github.com/trizen

# Tribonacci numbers - closed form.

# See also:
#   https://oeis.org/A000073

# Formula from Wolfram|Alpha
#   https://www.wolframalpha.com/input/?i=a(0)+%3D+0,a(1)%3D0,+a(2)%3D1,+a(n)+%3D+a(n-1)+%2B+a(n-2)+%2B+a(n-3)

define m = 1/3
define p = 2/3

define a = (99 + 19*sqrt(33))
define b = (19 + 3*sqrt(33))**m
define c = (19 - 3*sqrt(33))**m
define d = (4 * 33**(p))
define e = (33*a)**(m)
define f = a**(m)
define g = (1/6)*c
define h = (1/6)*b
define i = Complex.i
define j = i*sqrt(3)
define k = (1 - j)
define l = (1 + j)

func tribonacci(n) {
    [
      (k/e - l*f/d) * (m - g*k - h*l)**n,
      (l/e - k*f/d) * (m - g*l - h*k)**n,
      (f/(2 * 33**p) - 2/(33 * a)**m) * (3/(1 + b + c))**(-n)
    ].sum
}

for n in (2 .. 20) {
    say "T(#{n}) = #{tribonacci(n)}"
}

assert_eq(tribonacci( 8).round(-20), 24)
assert_eq(tribonacci( 9).round(-20), 44)
assert_eq(tribonacci(10).round(-20), 81)
assert_eq(tribonacci(20).round(-20), 35890)

say("Tribonacci constant: ", tribonacci(1e3+1)/tribonacci(1e3))

__END__
T(2) = 1
T(3) = 1
T(4) = 2
T(5) = 4
T(6) = 7
T(7) = 13
T(8) = 24
T(9) = 44
T(10) = 81
T(11) = 149
T(12) = 274
T(13) = 504
T(14) = 927
T(15) = 1705
T(16) = 3136
T(17) = 5768
T(18) = 10609
T(19) = 19513
T(20) = 35890
Tribonacci constant: 1.83928675521416113255185256465328660042417874609759