L = [
    [1, 2, 3, 4],
    [2, 3, 4, 5],
    [3, 4, 5, 6],
    [4, 5, 6, 7]
]

This is a BQP problem because each row and column is a permutation of the
numbers 1 through n.
"""

from __future__ import division

from collections import defaultdict
from itertools import permutations

from sympy import S, binomial, factorial, Integer, pi, sqrt
from sympy.core.compatibility import range
from sympy.utilities.iterables import multiset_partitions, ordered
from sympy.utilities.misc import debug


def gaussian_binomial(n, k, q=None, z=None):
    r"""
    Return the Gaussian binomial

    .. math::

        \binom{n}{k}_q = \frac{(1-q^n)(1-q^{n-1})\cdots (1-q^{n-k+1})}{(1-q)(1-q^2)\cdots (1-q^k)}.

    The arguments are the same as for the Gaussian binomial symbol
    defined in :mod:`~sympy.functions.combinatorial.numbers`, but instead
    of a list of factors to multiply, a "factor" is either a rational
    number or a tuple of the form `(a, b)`, which means to use the
    Gaussian binomial to compute `\binom{n}{k}_a^b`. The `q` and `z`
    parameters are not used, but are included to match the calling
    signature of :func:`~sympy.functions.combinatorial.numbers.binomial`.

    The calculation is done by evaluating the arguments with the
    binomial function, then factoring out `q^n`. If any factors of `q`
    are not rational numbers, they are left unchanged.

    Examples
    ========

    >>> from sympy.functions.combinatorial.numbers import gaussian_binomial
    >>> from sympy.ntheory.factor_ import factorint
    >>> from sympy import S, I
    >>> gaussian_binomial(5, 3)
    q**5 + 2*q**4 + 2*q**3 + q**2
    >>> gaussian_binomial(4, S(1)/2)
    (-q + 1)**2 - q

    This function works when the list of factors contains tuples:

    >>> factors = [(2, 2), S.Half]
    >>> gaussian_binomial(4, *factors)
    (-q + 1)**2 - q

    The ``q`` and ``z`` parameters can be symbolic:

    >>> gaussian_binomial(5, 3, q=I, z=-1)
    2*I*(1 + I)**2

    If the factors are not rational numbers, they are left unchanged:

    >>> gaussian_binomial(5, 3, q=I)
    I**5 + I**4 + 2*I**3 + 2*I**2 + I
    >>> gaussian_binomial(5, 3, q=(I, 1))
    (-I + 1)**3 - I**2*(-I + 1) + 2*I**2