##########################################################################
#  Copyright (C) 2006 Jaap Spies, jaapspies@gmail.com
#  Copyright (C) 2006 William Stein, wstein@gmail.com
#
#  Distributed under the terms of the GNU General Public License (GPL):
#
#                  http://www.gnu.org/licenses/
##########################################################################

def A079909(n):
    r"""
    function returns solutions to the Dancing School problem with 4 girls

    The value is $per(B)$, the permanent of the (0,1)-matrix $B$
    of size $4 \times 4+n$ with $b(i,j)=1$ if and only if $i \le j \le i+n$.

    We use precalculated values for $n = 0, 1$ For $n > 1$ we use the
    polynomial n^4 - 2*n^3 + 9*n^2 - 8*n + 6. Calculated with the function
    dance.

    See NAW 5/7 nr. 4 december 2006 p. 285

    INPUT:
        n -- non negative number

    OUTPUT
        integer

    EXAMPLES:

        sage: A079909(20)
        147446

        sage: A079909(10)
        8826



    AUTHOR:
        - Jaap Spies (2006-12-10)
    """
    if n == 0:
        return 1
    if n == 1:
        return 5
    else:
        return n^4 - 2*n^3 + 9*n^2 - 8*n + 6

#    return A079909_list(100)[n]
#    we will use precalculated values.


def A079909_list(n):
    r"""
    function returns list of solutions to the Dancing School problem with 4 girls

    INPUT:
        n -- non negative number

    OUTPUT
        list of integers

    EXAMPLES:
        sage: print A079909_list(10)
        [1, 5, 26, 90, 246, 566, 1146, 2106, 3590, 5766, 8826]

        sage: print A079909_list(1)
        [1, 5]

        sage: print A079909_list(0)
        [1]


    AUTHOR:
        - Jaap Spies (2006-12-10)
    """
    if n == 0:
        return [1]
    elif n == 1:
        return [1, 5]
    else:
        pol = dance(4)
        a = [pol(i) for i in range(2, n+1)]
        a.insert(0,1)
        a.insert(1,5)
        return a




# use variable 'h' in the polynomial ring over the rationals

h = QQ['h'].gen()

def dance(m):
    """
        Generates the polynomial solutions of the Dancing School Problem
        Based on a modification of theorem 7.2.1 from Brualdi and Ryser,
        Combinatorial Matrix Theory
        See NAW 5/7 nr. 4 december 2006 p. 285

        INPUT:
            m -- integer 

        OUTPUT:
            polynomial in 'h'

        EXAMPLES:
            sage: dance(4)
            h^4 - 2*h^3 + 9*h^2 - 8*h + 6

            sage: dance(5)
            h^5 - 5*h^4 + 25*h^3 - 55*h^2 + 80*h - 46


        AUTHOR:
            - Jaap Spies (2006)
    """    
    n = 2*m-2
    M = MatrixSpace(QQ, m, n)
    A = M([0 for i in range(m*n)])
    for i in range(m):
        for j in range(n):
            if i > j or j > i + n - m:
                A[i,j] = 1
    rv = A.rook_vector()
    s = sum([(-1)^k*rv[k]*falling_factorial(m+h-k, m-k) for k in range(m+1)])
    return s