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)))))))
And the answer is:
> (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.
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
}
}