## Tetrahedral Numbers

### September 13, 2011

Our first function calculates triangular numbers:

```(define (tri n)   (let loop ((n n) (s 0))     (if (zero? n) s       (loop (- n 1) (+ s n)))))```

Alternately, you may recognize Gauss’ formula for the sum of the first n integers:

`(define (tri n) (* n (+ n 1) 1/2))`

Then we can calculate tetrahedral numbers as the sum of triangular numbers:

```(define (tet n)   (let loop ((n n) (s 0))     (if (zero? n) s       (loop (- n 1) (+ s (tri n))))))```

We won’t show the derivation, but the recurrence has a closed form:

`(define (tet n) (* n (+ n 1) (+ n 2) 1/6))`

Then binary search can calculate the desired result:

```(define (prob-18 n)   (let loop ((lo 1) (hi 2))     (if (< (tet hi) n)         (loop hi (* hi 2))         (let loop ((lo lo) (hi hi))           (let* ((mid (quotient (+ lo hi) 2))                  (tet-mid (tet mid)))             (cond ((< tet-mid n) (loop mid hi))                   ((< n tet-mid) (loop lo mid))                   (else mid)))))))```

```> (prob-18 169179692512835000) 1005000```

We’ll let the mathematicians in the group figure out how this works:

```> (floor (expt (* 6 169179692512835000) 1/3)) 1005000.0```

You can run the program at http://programmingpraxis.codepad.org/6ZT1RVxa. If you’re interested in the math behind what we have done, see Sloane’s A000292.

Pages: 1 2

### 15 Responses to “Tetrahedral Numbers”

1. Tante Hedwig said

import Data.List

tria = tail \$ scanl (+) 0 [1..]
tetra = 1:zipWith (+) tetra (tail tria)

main = do
let (Just n) = findIndex (==169179692512835000) tetra
putStrLn \$ show n

```import Data.List

triangular :: [Integer]
triangular = scanl1 (+) [1..]

tetrahedral :: [Integer]
tetrahedral = scanl1 (+) triangular

main :: IO ()
main = print . maybe 0 succ \$
findIndex (== 169179692512835000) tetrahedral
```
3. […] today’s Programming Praxis, our goal is to find the base of the three-sided pyramid that has […]

4. Graham said

For my Python solution, I wrote a linear search as well as a binary one:

def tetra(n):
return n * (n + 1) * (n + 2) / 6

def prob18_linear(target):
n = 1
while tetra(n) < target:
n += 1
return n

def bin_search(target, func):
lo, hi = 1, 10
while func(hi) < t:
hi *= 10
mid = (lo + hi) / 2
fmid = func(mid)
while fmid != target:
if fmid < target:
lo, mid = mid, (mid + hi) / 2
else:
hi, mid = mid, (lo + mid) / 2
fmid = func(mid)
return mid

def prob18_log(t):
return bin_search(t, tetra)

5. Graham said

Sorry for the formatting error, I forgot a second quotation mark. Trying again:

```def tetra(n):
return n * (n + 1) * (n + 2) / 6

def prob18_linear(target):
n = 1
while tetra(n) < target:
n += 1
return n

def bin_search(target, func):
lo, hi = 1, 10
while func(hi) < t:
hi *= 10
mid = (lo + hi) / 2
fmid = func(mid)
while fmid != target:
if fmid < target:
lo, mid = mid, (mid + hi) / 2
else:
hi, mid = mid, (lo + mid) / 2
fmid = func(mid)
return mid

def prob18_log(t):
return bin_search(t, tetra)
```
6. Lautaro Pecile said
```def triangular():
t = 0
count = 1
while True:
t += count
yield t
count += 1

def tetrahedral():
tr = triangular()
a = 0
b = tr.next()
while True:
a += b
yield a
b = tr.next()

t = tetrahedral()
for i, n in enumerate(t):
if n == 169179692512835000:
print i
break
```
7. Mike said

Python 3 added accumulate() to the itertools library, which makes implementing
a linear search a snap.

```from itertools import accumulate, count

tgt = 169179692512835000

triangles = accumulate(count(1))
tetrahedrals = accumulate(triangles)
print(next(n for n,t in enumerate(tetrahedrals, 1) if t >= tgt))

```

Inspired by Graham’s binary search, I tried to find a way to use ‘bisect’
from the standard library. The Tetrahedrals class fakes being a list of
tetrahedral numbers by calculating the numbers on the fly based on the index
passed to __getitem__(). This is enough to get bisect() to work when using the
lo and hi parameters.

```
# binary search
from bisect import bisect

class Tetrahedrals:
def __getitem__(self, n):
return n * (n+1) * (n+2) // 6

tetras = Tetrahedrals()

tgt = 169179692512835000
n,nn = 0,1
while True:
t = bisect(tetras, tgt, n, nn)
if t < nn:
print(t - 1)
break

n, nn = nn, 2*nn

```
8. slabounty said

In Ruby …

```def linear(target, f)
n = 1
while ( f.call(n) != target)
n = n + 1
end
n
end

def binary(target, f)
low, high = 1, 2
while (f.call(high) < target) do high = high*2 end
mid = (high + low) / 2
while (fmid = f.call(mid)) != target do
fmid < target ?  (low, mid = mid, (mid + high) / 2) : (high, mid = mid, (low + mid) / 2)
end
mid
end

tetrahedral = lambda { |n| n * (n + 1) * (n + 2) / 6 }

1.upto(10) { |i| puts tetrahedral.call(i) }

puts linear(169179692512835000, tetrahedral)
puts binary(169179692512835000, tetrahedral)
```

We use a lambda to create a tetrahedral function we can pass around specifically to the linear (very slow in this case) and binary (much faster).

9. Axio said

;; With binary search. Quite naive.

```(define *number-of-balls* 169179692512835000)

(define (tetra n)
(/ (* n (1+ n) (+ 2 n)) 6))

;; bounds
(define (prax balls)
(let loop ((n 1) (last 1))
(let ((tn (tetra n)))
(cond
((= tn balls) n)
((< tn balls)
(loop (* n 2) n))
(else
(dicho balls last n))))))

;; dichotomic search
(define (dicho balls mi ma)
(let* ((mid (/ (+ ma mi) 2))
(mid-tn (tetra mid)))
(cond
((= balls mid-tn) mid)
((> balls mid-tn) (dicho balls mid ma))
(else (dicho balls mi mid)))) )

(define (main)
(pretty-print (prax *number-of-balls*)))

(main)
```
10. Axio said

Optimised/more clean version:

```(define *number-of-balls* 169179692512835000)

(define (tetra n)
(/ (* n (1+ n) (+ 2 n)) 6))

;; dichotomic search
(define (dicho balls mi ma)
(pp (list mi ma))
(let* ((mid (/ (+ ma mi) 2))
(mid-tn (tetra mid)))
(cond
((< (tetra ma) balls) (dicho balls ma (* 2 ma)))
((= balls mid-tn) mid)
((> balls mid-tn) (dicho balls mid ma))
(else (dicho balls mi mid)))) )

(define (main)
(pretty-print (dicho *number-of-balls* 1 1)))
```
11. Graham said

@Mike: very elegant. I was looking for a way to let bisect do the search for me, but didn’t think beyond having a very large list.

12. ardnew said

all brute force solutions?

here’s one using newton’s method:

```use strict;
use warnings;

my \$TOLERANCE = 0.1;

#
# the n'th tetrahedral number:
#
# T(n) = n(n + 1)(n + 2) / 6
#      = (n^3 + 3n^2 + 2n) / 6
#
sub t
{
my \$n = shift;
return (\$n * \$n * \$n + 3.0 * \$n * \$n + 2.0 * \$n) / 6.0;
}

#
# we are looking for the solution to:
#
# f(n) = T(n) = 169179692512835000
#      = T(n) - 169179692512835000 = 0
#
sub f
{
my \$n = shift;
return t(\$n) - 169179692512835000.0;

}

#
# the derivative of f(n) with respect to n:
#
# f'(n) = 3n^2 + 6n + 2
#
sub df
{
my \$n = shift;
return 3.0 * \$n * \$n + 6.0 * \$n + 2.0;
}

#
# perform newton's method to find the (unique!) root of f(n):
#
my \$xn   = 1000000.0; # initial guess
my \$xnp1 = 0.0;

do
{
\$xnp1 = \$xn;
\$xn = \$xn - f(\$xn) / df(\$xn);
#print "\$xn\n";
}
until (abs(\$xn - \$xnp1) <= \$TOLERANCE);

#
# naively round the solution
#
\$xn = int(\$xn + 0.5);

print "base number of balls = \$xn\n";
```
13. ardnew said

and since the theme seems to be keeping it as concise as possible, good luck reading this one

```do{\$b=\$a;\$a=\$a-(\$a**3+3*\$a**2+2*\$a-1015078155077010000)/(3*\$a**2+6*\$a+2)}until(abs(\$a-\$b)<=.1);

print "number of balls = ".int(\$a+.5)."\n";
```
14. THE G MAN said

(x^3+3*x^2+2*x)/6 = 169179692512835000; one real root: 1005000

Constant time. Thanks for playing.

15. GDean said

quick and dirty solution:

n<-500 #initial length of vector
#b<-20958500
b<-169179692512835000 #number of balls in tetrahedron of interest

for (y in 1:20) {

tri<-1
for (a in 2:n) {
tri<-c(tri,max(tri)+a)
}

dtri<-1

for (a in 2:n) {
dtri<-c(dtri,max(dtri)+tri[a])
}

c<-1
for (a in 1:n) {
if (dtri[a]==b) {
print(a)
c<-a
}
}

if (c==1) {
n<-n*10
}
else {
break
}
}