import math

import numpy
import sympy

from ..helpers import article
from ._helpers import TetrahedronScheme

citation = article(
    authors=["P. Silvester"],
    title="Symmetric quadrature formulae for simplexes",
    journal="Math. Comp.",
    volume="24",
    pages="95-100",
    year="1970",
    url="https://doi.org/10.1090/S0025-5718-1970-0258283-6",
)


def _newton_cotes(n, point_fun):
    degree = n

    # points
    idx = numpy.array(
        [
            [i, j, k, n - i - j - k]
            for i in range(n + 1)
            for j in range(n + 1 - i)
            for k in range(n + 1 - i - j)
        ]
    )
    points = point_fun(idx, n)

    # weights
    if n == 0:
        weights = numpy.ones(1)
        return weights, points, weights, degree

    def get_poly(t, m, n):
        return sympy.prod(
            [
                sympy.poly((t - point_fun(k, n)) / (point_fun(m, n) - point_fun(k, n)))
                for k in range(m)
            ]
        )

    weights = numpy.empty(len(points))
    idx = 0
    for i in range(n + 1):
        for j in range(n + 1 - i):
            for k in range(n + 1 - i - j):
                L = n - i - j - k
                # Compute weight.
                # Define the polynomial which to integrate over the
                # tetrahedron.
                t = sympy.DeferredVector("t")
                g = (
                    get_poly(t[0], i, n)
                    * get_poly(t[1], j, n)
                    * get_poly(t[2], k, n)
                    * get_poly(t[3], L, n)
                )
                # The integral of monomials over a tetrahedron are well-known,
                # see Silvester.
                weights[idx] = numpy.sum(
                    [
                        c
                        * numpy.prod([math.factorial(k) for k in m])
                        * 6.0
                        / math.factorial(numpy.sum(m) + 3)
                        for m, c in zip(g.monoms(), g.coeffs())
                    ]
                )
                idx += 1
    return weights, points, degree, citation


def newton_cotes_closed(n):
    return TetrahedronScheme(
        f"Newton-Cotes (closed, {n})", *_newton_cotes(n, lambda k, n: k / float(n))
    )


def newton_cotes_open(n):
    scheme = TetrahedronScheme(
        f"Newton-Cotes (open, {n})",
        *_newton_cotes(n, lambda k, n: (k + 1) / float(n + 4)),
    )
    if n == 0:
        scheme.degree = 1
    return scheme
