```# By Kirby Urner, 4D Solutions

import math, primes
from functions import Fraction

# useful constants (used in fibo3)

root5   = 5.0**0.5
gold    = (1+root5)/2.0 # same as phi, named to not conflict with phi()

tau     = 1/gold

def tri(n):
# triangular num n = sum of n consecutive counting nos

if n<=1: return n
else: return n + tri(n-1)

def tri2(n):
"""triangular num n = sum of n consecutive counting nos"""
return n*(n+1)/2

def sqr(n):
"""square num = sum of 2 consecutive triangular nos"""
return tri(n) + tri(n-1)

def tetra(n):
"""tetrahedral num = sum of n consecutive triangular nos"""
if n <= 1: return n
else: return tetra(n-1) + tri(n)

def tetra2(n):
"""tetrahedral num = sum of n consecutive triangular nos"""
sum = 0
for i in range(1,n+1):
sum = sum + tri(i)
return sum

def tetra3(n):
"""tetrahedral num = sum of n consecutive triangular nos"""
return (1.0/6.0)*n*(n+1)*(n+2)

def tetra4(n):
"""tetrahedral num = sum of n consecutive triangular nos"""
f=n-1
return (1.0/6.0)*f**3 + f**2 + (11.0/6.0)*f + 1

def hoct(n):
"""half octahedral num = sum n of consecutive squares"""
if n <= 1: return n
else: return hoct(n-1) + sqr(n)

def icoshell(n):
"""
number of spheres in an icosahedral shell
if n spheres along an edge
"""
if n==1: return 1
elif n<1: return 0
f= n-1
return 10*f**2 + 2

def icoshell2(n):
"""alternative function for nth icosa shell number"""
return 20*tri(n)-30*(n-2)-48

def cubocta(n):
"""
cumulative number of spheres in cuboctahedral packing
out to n-spheres per edge
"""
if n==1: return 1
elif n<1: return 0
return icoshell(n) + cubocta(n-1)

def cubocta2(n):
"""
cumulative number of spheres in cuboctahedral packing
out to n-spheres per edge
"""
return (10.0/3.0)*n**3 - 5*n**2 + (11.0/3.0)*n - 1

def cubocta3(n):
"""
cumulative number of spheres in cuboctahedral packing
out to n-spheres per edge
"""
sum = 0
for i in range(1,n+1):
sum = sum + icoshell(i)
return sum

def fact(n):
"""
factorial of n (recursive)
"""
if n <= 1: return long(1)
else: return long(n) * fact(n-1)

def fact2(n):
"""non-recursive factorial"""
return reduce(mul,[1L]+range(1,n+1))

def kfact(n,k):
"""n!/(n-k)!"""
return reduce(mul,[1L]+range(1,n+1)[n-k:])

def pascal(n,k):
"""row n, column k"""
return kfact(n,k)/fact2(k)

def prow(n):
"""return a row of Pascal's Triangle"""
row = []
for c in range(0,n+1):
row.append(pascal(n,c))
return row

def sumprow(n):
"""return the sum of the terms in a row of Pascal's Triangle"""

def showpascal(n):
"""print rows 0-n of Pascal's Triangle"""
for r in range(0,n+1):
row = prow(r)
print map(int,row)

def pasctet(n,r,c):  # thanks to David Chako

"""
lookup entry (r = row, c = column)
in level n of Pascal's Tetrahedron
"""
i,j,k = getijk(n,r,c)
return fact(n)/(fact(i)*fact(j)*fact(k))

def getijk(n,r,c):
i = n-r; j = r-c; k = n-i-j
return [i,j,k]

def pasctetrow(n,r):
"""
return a row of Pascal's Tetrahedron at level n
"""
row = []
for c in range(0,r+1):
row.append(pasctet(n,r,c))
return row

def pasctetlev(n):
"""print level n of Pascal's Tetrahedron"""
for r in range(n+1):
row = pasctetrow(n,r)
print row

def ptri(n):
"""lookup nth triangular number in pascal's triangle"""
return pascal(n+1,n-1)

def ptetra(n):
"""lookup nth tetrahedral number in pascal's triangle"""
return pascal(n+2,n-1)

fibcache = {}          # thanks to Tim Peters

def fibo(n):
"""return nth fibonacci number (recursive, cached)"""
if n == 0: return 0L
elif n == 1: return 1L
elif _fibcache.has_key(n): return _fibcache[n]
else:
result = fibo(n-1) + fibo(n-2)
_fibcache[n] = result
return result

def fibo1(n):
"""
slightly modififed from Python tutorial -- non-recursive
so you can get fibo1(10000) without hitting upper limits
on recursive depth
"""
i, a, b = 0, 0, 1L
while i < n:
a, b = b, a+b
i = i + 1  # or i += 1 in Python 2.0

return a

def fibo2(n):
"""
uses cut through Pascal's triangle to return
nth fibonacci number, fibonacci(5) returns 5L
"""
rval=0
if n%2==0: k = n/2
else: k = (n+1)/2
n=n-1
for i in range(k):
rval = rval + pascal(n-i,i)
return rval

def fibo3(n): # thanks David Zachmann

"""return nth fibonacci number as a floating point"""
return (gold**n - (-tau)**n)/root5

_bern = [Fraction(1,1),Fraction(-1,2)]

def bernoulli(n):
"""
Using pascal's triangle in series.py:

Bn = Bernoulli #       prow(n)
===============================
1                       row 0
1 1                     row 1
1+2B1               =0  row 2
1+3B1+ 3B2          =0  row 3
1+4B1+ 6B2+4B3      =0  row 4
1+5B1+10B2+10B3+5B4 =0  row 5
...
"""
if n<len(_bern):
return _bern[n]
for k in range(len(_bern)+1,n+2):
row = prow(k)
_bern.append(nxtbern)
return _bern[n]

def coeffs(x):
"""
See: Chapter 20, The Bernoulli Era, in
'The History of Mathematics' by Carl Boyer
(2nd edition, John Wiley and Sons, 1991), pg. 419.
"""
A = [(1,x+1),(1,2)]
if x<2:  return A
top,bot,k = long(x),2L,2L
while x>1:
b = bernoulli(k)
num,den = top*b.p,bot*b.q
f = primes.gcd(num,den)
if f>1: num,den = num/f,den/f
A.append((num,den))
top = top*(x-1)*(x-2)
bot = bot*(k+1)*(k+2)
k = k+2
x = x-2
return A

def poly(n):
"""
Generate the terms of a polynomial giving the sum of the
first m consecutive integers each raised to the nth power
"""
A = coeffs(n)
k = 0
terms = []
for c in A:
terms.append( "("+str(c[0])+".0/"+str(c[1])+".0)*m**"+str(n+1-k))
if k>1: k = k+2
else: k=k+1
return terms

def evalpoly(t,m):
"""
evaluate a polynomial given a list of terms from poly,
and a value for m
"""
m = long(m)

def sigma(m,n):
"""sum consecutive integers 1...m each to the nth power"""
sum = 0L
for i in range(m+1):
sum = sum + pow(long(i),n)
return sum

def phi(depth):
"""use continued fraction of user-supplied depth to refine phi"""
if depth>0:
return 1.0 + 1.0/phi(depth - 1)
else:
return 1.0

def pifib(n):
"""
return approximation for pi using
arctan(1) = 4 * sum of arctan(1.0/Fib(2i+1)) for i = 1...n
"""
sum = 0
for i in range(1,n+1):
sum = sum + math.atan(1.0/fibo(2*i+1))
return 4*sum

def euler(depth):
"""use continued fraction of user-supplied depth to refine e"""
n = 1.0
return 2.0 + 1.0/efract(n,depth)

def efract(n,depth):
if n<depth:
return n + n/efract(n+1,depth)
else: return n

def pi(depth):
"""use continued fraction of user-supplied depth to refine pi"""
n = 1.0
return 4.0 * 1.0/(1 + 1.0/pifract(n,depth))

def pifract(n,depth):
if n<depth:
return 2 + ((2*n+1)**2)/pifract(n+1,depth)
else: return 2.0

```
`# code highlighted using py2html.py version 0.8`