def anf_coeffs(truthvalues):
    """
    Convert a list of truth values of some boolean expression
    to the list of coefficients of the polynomial mod 2 (exclusive
    disjunction) representing the boolean expression in ANF
    (i.e., the "Zhegalkin polynomial").

    There are 2^n possible Zhegalkin monomials in n variables, since
    each monomial is fully specified by the presence or absence of
    each variable.

    We can enumerate all the monomials. For example, boolean
    function with four variables (a, b, c, d) can contain
    up to 2^4 = 16 monomials. The 13-th monomial is the
    product a & b & d, because 13 in binary is 1, 1, 0, 1.

    A given monomial's presence or absence in a polynomial corresponds
    to that monomial's coefficient being 1 or 0 respectively.

    Examples
    ========
    >>> from sympy.logic.boolalg import anf_coeffs, bool_monomial, Xor
    >>> from sympy.abc import a, b, c
    >>> truthvalues = [0, 1, 1, 0, 0, 1, 0, 1]
    >>> coeffs = anf_coeffs(truthvalues)
    >>> coeffs
    [0, 1, 1, 0, 0, 0, 1, 0]
    >>> polynomial = Xor(*[
    ...     bool_monomial(k, [a, b, c])
    ...     for k, coeff in enumerate(coeffs) if coeff == 1
    ... ])
    >>> polynomial
    b ^ c ^ (a & b)

    """

    s = '{0:b}'.format(len(truthvalues))
    n = len(s) - 1

    if len(truthvalues) != 2**n:
        raise ValueError("The number of truth values must be a power of two, "
                         "got %d" % len(truthvalues))

    coeffs = [[v] for v in truthvalues]

    for i in range(n):
        tmp = []
        for j in range(2 ** (n-i-1)):
            tmp.append(coeffs[2*j] +
                list(map(lambda x, y: x^y, coeffs[2*j], coeffs[2*j+1])))
        coeffs = tmp

    return coeffs[0]