'Compute the coefficients of a polynomial in the Bernstein basis
given the values and derivatives at the edges.
Return the coefficients of a polynomial in the Bernstein basis
defined on `[xa, xb]` and having the values and derivatives at the
endpoints ``xa`` and ``xb`` as specified by ``ya`` and ``yb``.
The polynomial constructed is of the minimal possible degree, i.e.,
if the lengths of ``ya`` and ``yb`` are ``na`` and ``nb``, the degree
of the polynomial is ``na + nb - 1``.
Parameters
xa : float
Left-hand end point of the interval
xb : float
Right-hand end point of the interval
ya : array_like
Derivatives at ``xa``. ``ya[0]`` is the value of the function, and
``ya[i]`` for ``i > 0`` is the value of the ``i``-th derivative.
yb : array_like
Derivatives at ``xb``.
Returns
array
coefficient array of a polynomial having specified derivatives
Notes
This uses several facts from life of Bernstein basis functions.
First of all,
.. math:: b\'_{a, n} = n (b_{a-1, n-1} - b_{a, n-1})
If B(x) is a linear combination of the form
.. math:: B(x) = \sum_{a=0}^{n} c_a b_{a, n},
then :math: B\'(x) = n \sum_{a=0}^{n-1} (c_{a+1} - c_{a}) b_{a, n-1}.
Iterating the latter one, one finds for the q-th derivative
.. math:: B^{q}(x) = n!/(n-q)! \sum_{a=0}^{n-q} Q_a b_{a, n-q},
with
.. math:: Q_a = \sum_{j=0}^{q} (-)^{j+q} comb(q, j) c_{j+a}
This way, only `a=0` contributes to :math: `B^{q}(x = xa)`, and
`c_q` are found one by one by iterating `q = 0, ..., na`.
At `x = xb` it\'s the same with `a = n - q`.'

@staticmethod
def _construct_from_derivatives(xa, xb, ya, yb):

    (ya, yb) = (np.asarray(ya), np.asarray(yb))
    if (ya.shape[1:] != yb.shape[1:]):
        raise ValueError('ya    and    yb    have    incompatible    dimensions.')
    (dta, dtb) = (ya.dtype, yb.dtype)
    if (np.issubdtype(dta, np.complexfloating) or np.issubdtype(dtb, np.complexfloating)):
        dt = np.complex_
    else:
        dt = np.float_
    (na, nb) = (len(ya), len(yb))
    n = (na + nb)
    c = np.empty((((na + nb),) + ya.shape[1:]), dtype=dt)
    for q in range(0, na):
        c[q] = ((ya[q] / spec.poch((n - q), q)) * ((xb - xa) ** q))
        for j in range(0, q):
            c[q] -= ((((-1) ** (j + q)) * comb(q, j)) * c[j])
    for q in range(0, nb):
        c[((- q) - 1)] = (((yb[q] / spec.poch((n - q), q)) * ((-1) ** q)) * ((xb - xa) ** q))
        for j in range(0, q):
            c[((- q) - 1)] -= ((((-1) ** (j + 1)) * comb(q, (j + 1))) * c[((- q) + j)])
    return c
