import math, primes
from functions import Fraction
from operator import add, mul
root5 = 5.0**0.5
gold = (1+root5)/2.0
tau = 1/gold
def tri(n):
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"""
return reduce(add,prow(n))
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):
"""
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 = {}
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
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):
"""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)
nxtbern = -reduce(add,map(mul,_bern[:k-1],row[:k-1]))/row[-2]
_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)
return round((reduce(add,map(eval,t))),0)
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