import itertools
import math
import functools
import collections
import fractions

[docs]def nCk(n,k):
r'''

Compute the binomial coefficients :math:\binom nk = \frac{n!}{k!(n-k)!}

>>> nCk(0,0)
1
>>> nCk(1,0)
1
>>> nCk(4, 2)
6
>>> [[nCk(n,k) for k in range(n+1)] for n in range(7)] == [
...          [1],
...        [1, 1],
...       [1, 2, 1],
...      [1, 3, 3, 1],
...     [1, 4, 6, 4, 1],
...    [1, 5, 10, 10, 5, 1],
...  [1, 6, 15, 20, 15, 6, 1]                              ]
True
'''

f = math.factorial
return f(n) // f(k) // f(n-k)

def tostr(expr):
'''

>>> q =  ([((('x', 0),), 1), ((('x', 1),), 1)], fractions.Fraction(1, 2))
>>> tostr(q)
'1/2 (x0 + x1)'
'''
expression, multiple = expr
terms = []
for t,c in expression:
s = ''
if c > 1:
s += str(c)
s += ' '
s += ''.join(map(str, sum(t,())))
terms.append(s)
s = ' + '.join(terms)
if not multiple == 1:
s = '{}/{} ({})'.format(multiple.numerator, multiple.denominator, s)
return s

'''

Compute the quadrature of a homogenous polynomial given by term='xxy'
over the vertices verts=(0,1,2) of a simplex.

([((), 1)], Fraction(1, 1))
([((('x', 0),), 1), ((('x', 1),), 1)], Fraction(1, 2))
([((('x', 0),), 1), ((('x', 1),), 1), ((('x', 2),), 1)], Fraction(1, 3))
...    [((('x', 0), ('x', 0)), 1),
...     ((('x', 0), ('x', 1)), 1),
...     ((('x', 1), ('x', 1)), 1)], fractions.Fraction(1, 3))
True
...    [((('x', 0), ('y', 0)), 2),
...      ((('x', 0), ('y', 1)), 1),
...      ((('x', 0), ('y', 2)), 1),
...      ((('x', 1), ('y', 0)), 1),
...      ((('x', 1), ('y', 1)), 2),
...      ((('x', 1), ('y', 2)), 1),
...      ((('x', 2), ('y', 0)), 1),
...      ((('x', 2), ('y', 1)), 1),
...      ((('x', 2), ('y', 2)), 2)], fractions.Fraction(1, 12))
True
'''

k = len(verts)-1    # dimension of simplex
q = len(term)       # order of term
P = itertools.permutations(term)
V = itertools.combinations_with_replacement(verts, q)
Xp = itertools.product(P, V)
X = (tuple(sorted(zip(*x))) for x in Xp)
counter = collections.Counter(X)

gcd = functools.reduce(fractions.gcd, counter.values())

return (sorted([(c, counter[c] // gcd) for c in counter ]),
fractions.Fraction(gcd, math.factorial(q)*nCk(k+q, q)))

def q(term, verts):

[docs]def test_line_segment():
'''

>>> q('x')
'1/2 (x0 + x1)'
>>> q('xx')
'1/3 (x0x0 + x0x1 + x1x1)'
>>> q('xxx')
'1/4 (x0x0x0 + x0x0x1 + x0x1x1 + x1x1x1)'
>>> q('xy')
'1/6 (2 x0y0 + x0y1 + x1y0 + 2 x1y1)'
>>> q('xxy')
'1/12 (3 x0x0y0 + x0x0y1 + 2 x0x1y0 + 2 x0x1y1 + x1x1y0 + 3 x1x1y1)'
'''
pass

[docs]def test_triangle():
'''

>>> q('x')
'1/3 (x0 + x1 + x2)'
>>> q('xx')
'1/6 (x0x0 + x0x1 + x0x2 + x1x1 + x1x2 + x2x2)'
>>> q('xxx')
'1/10 (x0x0x0 + x0x0x1 + x0x0x2 + x0x1x1 + x0x1x2 + x0x2x2 + x1x1x1 + x1x1x2 + x1x2x2 + x2x2x2)'
>>> q('xy')
'1/12 (2 x0y0 + x0y1 + x0y2 + x1y0 + 2 x1y1 + x1y2 + x2y0 + x2y1 + 2 x2y2)'
'''
pass

[docs]def test_tetrahedron():
'''