import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.symbolic.ring import SR
import sage.rings.polynomial.laurent_polynomial_ring
from sage.combinat.sf.sf import SymmetricFunctions
from sage.combinat.partition import Partitions
from math import factorial, prod
from collections import Counter


def orthogonal_weingarten(N, mu):
    k = sum(mu)
    sf = SymmetricFunctions(QQ)
    p, s, J = sf.power(), sf.schur(), sf.jack(t=QQ(2)).J()
    z = prod(r**m * factorial(m) for r, m in Counter(mu).items())
    total = QQ(0)
    for lam in Partitions(k):
        if len(lam) > N:
            continue
        f = s[tuple(2*r for r in lam)].scalar(p[[1]*(2*k)])
        c = p(J[lam]).coefficient(mu)
        d = prod(N + 2*j - i for i, row in enumerate(lam) for j in range(row))
        total += QQ(f*c) / QQ(d)
    return QQ(2**len(mu)*z) * total / QQ(factorial(2*k))


print(orthogonal_weingarten(13, (3, 1)))