# By Kirby Urner, 4D Solutions
# Last modified: Oct 14, 2000

import math, primes
from functions import Fraction
from operator import add, mul

# 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"""
    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):  # 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)
       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