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.

About these ads

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

  2. 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
    
  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
    }
    }

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 600 other followers

%d bloggers like this: