def METHOD_NAME(a, n):
    f = lambda j, i: reduce(mul, ((a[i] - a[k] if k != i else 1) for k in range(0, n - 1)), 1)
    Ff = np.fromfunction(np.vectorize(f), (1, n - 1), dtype=int)
    f = (
        lambda i, nth: (
            reduce(mul, [(np.poly1d([1, -a[k]]) if k != i else 1) for k in range(0, n - 1)], 1)
        ).coef[n - 1 - nth - 1]
        / Ff[0, i]
    )
    F = np.fromfunction(np.vectorize(f), (n - 1, n - 1), dtype=int)
    f = lambda i, j: -a[i] ** (n - 1)
    t = np.fromfunction(np.vectorize(f), (n - 1, 1), dtype=int)
    T = np.append(np.eye(n - 1), t, axis=1)
    return np.append(F.T.dot(T), np.array([np.eye(n)[n - 1]]), axis=0)