Tetrahedral Numbers
September 13, 2011
Before I started writing my own programming exercises, I enjoyed solving other programming exercises available on the internet, many of them mathematical in nature. Today’s exercise comes from http://open-cs.net/problems.php?id=18 and concerns triangular and tetrahedral numbers.
A triangular number tells the number of ways that balls can be stacked in a triangle. The first triangular number is 1, the second is 3 (a row of 1 plus a row of 2), the third triangular number is 6 (the first two rows plus a row of 3), the fourth triangular number is 10 (adding a row of 4), the fifth triangular number is 15 (think of the 15 balls on a pool table), and so on.
A tetrahedral number is the three-dimensional equivalent of a triangular number; think of cannonballs stacked in a three-sided pyramid. The top layer has one ball, the second layer has 3 balls (the second triangular number) so the second tetrahedral number is 1 + 3 = 4, the third layer has 6 balls giving a total of 10 balls in the tetrahedron, and so on.
Your task is to find the base of the tetrahedron that contains 169179692512835000 balls. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
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
My Haskell solution (see http://bonsaicode.wordpress.com/2011/09/13/programming-praxis-tetrahedral-numbers/ for a version with comments):
import Data.List triangular :: [Integer] triangular = scanl1 (+) [1..] tetrahedral :: [Integer] tetrahedral = scanl1 (+) triangular main :: IO () main = print . maybe 0 succ $ findIndex (== 169179692512835000) tetrahedral[…] today’s Programming Praxis, our goal is to find the base of the three-sided pyramid that has […]
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)
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)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 breakPython 3 added accumulate() to the itertools library, which makes implementing
a linear search a snap.
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*nnIn 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).
;; 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)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)))@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.
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";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";(x^3+3*x^2+2*x)/6 = 169179692512835000; one real root: 1005000
Constant time. Thanks for playing.
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
}
}