def main():
    size = 20
    # size 'rights' + size 'downs'
    # combinations formula
    steps = size * 2
    result = f(steps) / (f(size) * f(steps - size))
    print('Result: ', int(result))

def fixation_prob(n, m, a):
    matrix = [[0 for i in range(len(a))] for j in range(m)]
    
    # Population is diploid so number of alleles is 2N    
    n = 2*n
    
    for i in range(len(a)):
        # Frequency of each number of recessive alleles in the current pop.
        curr_gen = [0 for x in range(n+1)]
        curr_gen[a[i]] = 1
        
        # For each generatjon...
        for gen in range(m):
            next_gen = [0 for x in range(n+1)]
            
            for j in range(n+1): # frequency in next generation...
                for k in range(n+1): # frequency in current geneneration...
                    p1 = f(n) / f(j) / f(n-j) 
                    p2 = (k/n)**j * (1-(k/n))**(n-j) 
                    p3 = curr_gen[k]
                    next_gen[j] += p1 * p2 * p3
                    
            curr_gen = next_gen
            
            # Store the common logarithm of the answer.
            matrix[gen][i] = log10(curr_gen[0])

    return matrix

    def handle(self, *args, **options):
        school_dict = dict()
        for m in Project.objects.filter(
            category="micro_sumo",is_confirmed=True):
            if school_dict.has_key(m.manager.school):
                school_dict[m.manager.school] += 1
            else:
                school_dict[m.manager.school] = 1
        school_pairs = school_dict.items()
        school_pairs.sort(key=lambda tup: tup[1])
        school_list = map(lambda x: x[0], school_pairs)
        robot_count = sum(school_dict.values())

        # generate groups
        if robot_count % 4 == 1:
            group_count = robot_count / 4
        elif robot_count % 4 == 2:
            group_count = robot_count / 4 + 1
        elif robot_count % 4 == 3:
            group_count = robot_count / 4 + 1
        elif robot_count % 4 == 0:
            group_count = robot_count / 4

        for i in range(1, group_count+1):
            SumoGroup.objects.create(order=i, is_final=False)

        group_list = SumoGroup.objects.filter(is_final=False)
        active_groups = list(group_list)
        passive_groups = list()
        shuffle(active_groups)

        for school in school_list:
            for m in Project.objects.filter(
                    category="micro_sumo",
                    manager__school=school, is_confirmed=True):
                if len(active_groups) == 0:
                    active_groups = passive_groups[:]
                    shuffle(active_groups)
                    passive_groups = list()
                current_group = active_groups.pop()
                passive_groups.append(current_group)
                SumoGroupTeam.objects.create(
                    group=current_group, robot=m)
        self.stdout.write('Sumo groups generated.')

        for group in SumoGroup.objects.all():
            order = 1
            teams = SumoGroupTeam.objects.filter(group=group)
            team_list = list(teams)
            count = len(team_list)
            for j in range(f(count) / f(2) / f(count-2) / 2):
                for i in range(0, len(team_list) / 2):
                    SumoGroupMatch.objects.create(
                        home=team_list[i].robot,
                        away=team_list[len(team_list) - i - 1].robot,
                        group=group, order=order)
                    order += 1
                hold = team_list.pop()
                team_list.insert(1, hold)
        self.stdout.write("Fixtures generated.")

def perms(alist):
    acounter = Counter(alist)
    n = len(alist)
    ans = f(n)
    for key in acounter:
        ans /= f(acounter[key])
    return ans

def pure_py(xyz, Snlm, Tnlm, nmax, lmax):
    from scipy.special import lpmv, gegenbauer, eval_gegenbauer, gamma
    from math import factorial as f

    Plm = lambda l,m,costh: lpmv(m, l, costh)
    Ylmth = lambda l,m,costh: np.sqrt((2*l+1)/(4 * np.pi) * f(l-m)/f(l+m)) * Plm(l,m,costh)

    twopi = 2*np.pi
    sqrtpi = np.sqrt(np.pi)
    sqrt4pi = np.sqrt(4*np.pi)

    r = np.sqrt(np.sum(xyz**2, axis=0))
    X = xyz[2]/r # cos(theta)
    sinth = np.sqrt(1 - X**2)
    phi = np.arctan2(xyz[1], xyz[0])
    xsi = (r - 1) / (r + 1)

    density = 0
    potenti = 0
    gradien = np.zeros_like(xyz)
    sph_gradien = np.zeros_like(xyz)
    for l in range(lmax+1):
        r_term1 = r**l / (r*(1+r)**(2*l+3))
        r_term2 = r**l / (1+r)**(2*l+1)
        for m in range(l+1):
            for n in range(nmax+1):
                Cn = gegenbauer(n, 2*l+3/2)
                Knl = 0.5 * n * (n+4*l+3) + (l+1)*(2*l+1)
                rho_nl = Knl / twopi * sqrt4pi * r_term1 * Cn(xsi)
                phi_nl = -sqrt4pi * r_term2 * Cn(xsi)

                density += rho_nl * Ylmth(l,m,X) * (Snlm[n,l,m]*np.cos(m*phi) +
                                                    Tnlm[n,l,m]*np.sin(m*phi))
                potenti += phi_nl * Ylmth(l,m,X) * (Snlm[n,l,m]*np.cos(m*phi) +
                                                    Tnlm[n,l,m]*np.sin(m*phi))

                # derivatives
                dphinl_dr = (2*sqrtpi*np.power(r,-1 + l)*np.power(1 + r,-3 - 2*l)*
                              (-2*(3 + 4*l)*r*eval_gegenbauer(-1 + n,2.5 + 2*l,(-1 + r)/(1 + r)) +
                              (1 + r)*(l*(-1 + r) + r)*eval_gegenbauer(n,1.5 + 2*l,(-1 + r)/(1 + r))))
                sph_gradien[0] += dphinl_dr * Ylmth(l,m,X) * (Snlm[n,l,m]*np.cos(m*phi) +
                                                              Tnlm[n,l,m]*np.sin(m*phi))

                A = np.sqrt((2*l+1) / (4*np.pi)) * np.sqrt(gamma(l-m+1) / gamma(l+m+1))
                dYlm_dth = A / sinth * (l*X*Plm(l,m,X) - (l+m)*Plm(l-1,m,X))
                sph_gradien[1] += (1/r) * dYlm_dth * phi_nl * (Snlm[n,l,m]*np.cos(m*phi) +
                                                               Tnlm[n,l,m]*np.sin(m*phi))

                sph_gradien[2] += (m/(r*sinth)) * phi_nl * Ylmth(l,m,X) * (-Snlm[n,l,m]*np.sin(m*phi) +
                                                                           Tnlm[n,l,m]*np.cos(m*phi))

    cosphi = np.cos(phi)
    sinphi = np.sin(phi)
    gradien[0] = sinth*cosphi*sph_gradien[0] + X*cosphi*sph_gradien[1] - sinphi*sph_gradien[2]
    gradien[1] = sinth*sinphi*sph_gradien[0] + X*sinphi*sph_gradien[1] + cosphi*sph_gradien[2]
    gradien[2] = X*sph_gradien[0] - sinth*sph_gradien[1]

    return density, potenti, gradien

def solve(fname):
    v,c = read_connections(fname)
    c = sorted(c,key=lambda x:x[2])
    from math import factorial as f
    n = len(c)
    print(f(n)/f(n-v+1)/f(v-1))
    for subset in itertools.combinations(c,v-1):
        if fully_connected(v,subset):
            return weight(c)-weight(subset)

def surprise_value(ss):
    from dist import simple_dist
    from math import factorial as f
    if len(ss) < 2: return 0

    l = len(ss)
    nC2 = f(l) / 2 / f(l-2)
    total_dist = sum(simple_dist(a,b) for (a,b) in itertools.combinations(ss,2))
    return total_dist / nC2

 def cntMatrix(self, A):
     # determine how far columns are sorted
     N, count = len(A), 0
     for i in xrange(1, N + 1):
         if N % i == 0:
             # check if all i chunks of array are sorted
             if self.chunkSortedChecker(i, A):
                 print 'For i: ', i,  pow(f(i), N / i, mod)
                 count += pow(f(i), N / i, mod)
     return count % mod

def bp_prob(s, a, b, k):
    n = len(s)
    N = count_kmers(s, 1)
    nn = [(i[1], i[0] in [a,b]) for i in N.items()]
    
    if k > min(N[a],N[b]) or k < 0: 
        return 0 
    
    ns = sum([( f(n-j) / ( reduce(lambda x,y: x*( f(y[0] - (j if y[1] else 0) ) ) , nn, 1) * f(j-k) * f(k) ) ) * (-1)**(j-k)\
              for j in range(k,min(N[a],N[b])+1)])
    
    return (ns / ( f(n) / ( reduce(lambda x,y: x*f(y[0]) , nn, 1) ) / 10**100   )   )*(10.0**(-100))

 def getPermutation(self, n, k):
     from math import factorial as f
     nums = range(1, n + 1)
     result = ''
     for i in range(n):
         q, r = k / f(n - 1 - i), k % f(n - 1 - i)
         if r > 0:
             q += 1
         result += str(nums[q - 1])
         nums.remove(nums[q - 1])
         k = r
     return result

def solve_equation(a, b, c):
    d = b * b - 4 * a * c
    if d > 0:
        x1 = (-b + f(d)) / (2 * a)
        x2 = (-b - f(d)) / (2 * a)
        return x1, x2
    elif d == 0:
        x = -b / (2 * a)
        return x, x
    else:
        x1 = complex(-b / (2 * a), f(-d) / (2 * a))
        x2 = complex(-b / (2 * a), -f(-d) / (2 * a))
        return x1, x2

def factorial_digits(x):
    x = str(x)
    total = 0
    for c in x:
        total += f(int(c))

    return total

def foo(num):
	a=[num]
	while 1:
		tmp=sum([f(int(i)) for i in list(str(a[-1]))])
		if tmp in a: break
		a.append(tmp)
	return len(a)

def ncr(n,r):
    a = max(n-r,r)
    b = min(n-r,r)
    fac = 1
    while n>a:
        fac *= n
        n -= 1
    return fac/f(b)

def main():
    # the max that we can reach is 9! * the number of digits
    i = 1
    lim = f(9)
    while lim * i > pow(10, i):
        i += 1
    lim *= i

    return sum([n for n in xrange(lim + 1) if is_factorian(n)])

def foo2(num):
	a=[num]
	while 1:
		if a[-1] in b: tmp=b[a[-1]]
		else:
			tmp=sum([f(int(i)) for i in list(str(a[-1]))])
			b[a[-1]]=tmp
		if tmp in a: break
		a.append(tmp)
	return len(a)

def probability(n, m, g, k):  
    # Population is diploid so number of alleles is 2N    
    n = 2*n
    
    # Hold the frequency of each number of dominant alleles in the current pop.
    curr_gen = [0 for i in range(n+1)]
    curr_gen[m] = 1
    
    # For each generation...
    for gen in range(g):
        next_gen = [0 for i in range(n+1)]
        
        for i in range(n+1): # next generation...
            for j in range(n+1): # current geneneration...
                a = f(n) / f(i) / f(n-i) 
                b = (j/n)**i * (1-(j/n))**(n-i) 
                c = curr_gen[j]
                next_gen[i] += a * b * c
                
        curr_gen = next_gen

    return sum(curr_gen[:-k])

def main(): # how many lattice-routes (right, down) are there through a 20x20 grid?
	
	n = (20,20)
	path_length = n[0]

	'''
	# n! ouch!
	one_route = [0 for i in xrange(n[0])] + [1 for i in xrange(n[0])] # one-route right:0, down:1
	all_routes = set(list(permutations(one_route)))
	return len(all_routes)
	
	# this method assumes symmetry for grid
	# calculating combinations is much faster than permutations... like n^2 ?
	return len(list(combinations(range(path_length*2),path_length)))
	'''

	# I think we can do better...
	# oh duh... we don't need a list of all the paths,
	# just get the number -> use math:
	# n! / ( r! (n - r)! )
	n = path_length*2
	r = path_length
	ret = f(n)/( f(r) * f(n-r) )
	return ret

def solve():
    m = f(9)
    for i in range(1,10):
        if i*m < 10**(i-1):
            lim = 10**(i-1)
            break
    n = 10
    s = 0
    while n < lim: 
        p = f_sum(n)
        if p == n:
            s += n
        if p > n:
            n = (n//10 + 1)*10
        else:
            n += 1
    return s

def main():
    # Seek to find the number of non-repeating chains below one million
    # that have 60 terms

    # Upper limit
    lim = 1000000

    # First populate a cache of the factorials so we don't repeat them
    factorials = dict((str(n), f(n)) for n in range(10))

    # Cache chain lengths
    chain_lengths = {}

    for n in xrange(lim):
        current_chain = []
        current_chain.append(n)

        new = n
        while True:
            new = sum(factorials[i] for i in str(new))

            if new in chain_lengths:
                # We know how long this chain is
                break

            if new in current_chain:
                # repeated digit, chain is over
                break
            current_chain.append(new)

        # Find the additional contribution from the last sum
        additional_size = chain_lengths.get(new, 0)

        for i, num in enumerate(current_chain):
            if num not in chain_lengths:
                chain_lengths[num] = additional_size + len(current_chain) - i

    # Find all chains that are 60 sections long before they repeat
    # Note: in python, True is treated as 1 and False as 0
    solution = sum(value == 60 for value in chain_lengths.values())

    return solution