def __init__(self, k, p=None, prec_cap=20, base=None, character=None, tuplegen=None, act_on_left=False):
     """
     - ``character`` --
       - None (default)
       - (chi, None)
       - (None, n) (n integral)
       - (chi, n)
       - lambda (for n half-integral use this form)
     """
     if p is not None:
         p = ZZ(p)
     if base is None:
         if p is None: raise ValueError("specify p or a base")
         base = ZpCA(p,prec_cap)
     elif isinstance(base, pAdicGeneric):
         if base.prime() != p: raise ValueError("p must be the same as the prime of base")
         if base.precision_cap() != prec_cap: raise ValueError("prec_cap must match the precision cap of base")
     elif prec_cap > k+1: # non-classical
         if p is None or not p.is_prime(): raise ValueError("p must be prime for non-classical weight")
     from sage.rings.padics.pow_computer import PowComputer_long
     # should eventually be the PowComputer on ZpCA once that uses longs.
     Dist, WeightKAction = get_dist_classes(p, prec_cap, base)
     self.Element = Dist
     if Dist is Dist_long:
         self.prime_pow = PowComputer_long(p, prec_cap, prec_cap, prec_cap, 0)
     Parent.__init__(self, base)
     self._k = k
     self._p = p
     self._prec_cap = prec_cap
     act = WeightKAction(self, character, tuplegen, act_on_left)
     self._act = act
     self._populate_coercion_lists_(action_list=[act])

    def random_element(self, bound=100, *args, **kwds):
        r"""
        Return a random element of `{\rm SL}_2(\ZZ)` with entries whose
        absolute value is strictly less than bound (default 100).
        Additional arguments and keywords are passed to the random_element
        method of ZZ.

        (Algorithm: Generate a random pair of integers at most bound. If they
        are not coprime, throw them away and start again. If they are, find an
        element of `{\rm SL}_2(\ZZ)` whose bottom row is that, and
        left-multiply it by `\begin{pmatrix} 1 & w \\ 0 & 1\end{pmatrix}` for
        an integer `w` randomly chosen from a small enough range that the
        answer still has entries at most bound.)

        It is, unfortunately, not true that all elements of SL2Z with entries <
        bound appear with equal probability; those with larger bottom rows are
        favoured, because there are fewer valid possibilities for w.

        EXAMPLES::

            sage: SL2Z.random_element()
            [60 13]
            [83 18]
            sage: SL2Z.random_element(5)
            [-1  3]
            [ 1 -4]

        Passes extra positional or keyword arguments through::

            sage: SL2Z.random_element(5, distribution='1/n')
            [ 1 -4]
            [ 0  1]
        """
        if bound <= 1: raise ValueError("bound must be greater than 1")
        c = ZZ.random_element(1-bound, bound, *args, **kwds)
        d = ZZ.random_element(1-bound, bound, *args, **kwds)
        if gcd(c,d) != 1: # try again
            return self.random_element(bound, *args, **kwds)
        else:
            a,b,c,d = lift_to_sl2z(c,d,0)
            whi = bound
            wlo = bound
            if c > 0:
                whi = min(whi, ((bound - a)/ZZ(c)).ceil())
                wlo = min(wlo, ((bound + a)/ZZ(c)).ceil())
            elif c < 0:
                whi = min(whi, ((bound + a)/ZZ(-c)).ceil())
                wlo = min(wlo, ((bound - a)/ZZ(-c)).ceil())

            if d > 0:
                whi = min(whi, ((bound - b)/ZZ(d)).ceil())
                wlo = min(wlo, ((bound + b)/ZZ(d)).ceil())
            elif d < 0:
                whi = min(whi, ((bound + b)/ZZ(-d)).ceil())
                wlo = min(wlo, ((bound - b)/ZZ(-d)).ceil())

            w = ZZ.random_element(1-wlo, whi, *args, **kwds)
            a += c*w
            b += d*w
            return self([a,b,c,d])

    def random_element(self, algorithm='default'):
        r"""
        Returns a random element of self, optionally using the
        algorithm argument to decide how it generates the
        element. Algorithms currently implemented:

        - default: Choose `a_i`, `i >= 0`, randomly between `0` and
          `p-1` until a nonzero choice is made. Then continue choosing
          `a_i` randomly between `0` and `p-1` until we reach
          precision_cap, and return `\sum a_i p^i`.

        EXAMPLES::

            sage: Zp(5,6).random_element()
            3 + 3*5 + 2*5^2 + 3*5^3 + 2*5^4 + 5^5 + O(5^6)
            sage: ZpCA(5,6).random_element()
            4*5^2 + 5^3 + O(5^6)
            sage: ZpFM(5,6).random_element()
            2 + 4*5^2 + 2*5^4 + 5^5 + O(5^6)
        """
        if (algorithm == 'default'):
            if self.is_capped_relative():
                i = 0
                a_i = ZZ.random_element(self.prime())
                while a_i.is_zero():
                    i += 1
                    a_i = ZZ.random_element(self.prime())
                return self((self.prime()**i)*(a_i + self.prime()*ZZ.random_element(self.prime_pow.pow_Integer_Integer(self.precision_cap()-1))))
            else:
                return self(ZZ.random_element(self.prime_pow.pow_Integer_Integer(self.precision_cap())))
        else:
            raise NotImplementedError("Don't know %s algorithm"%algorithm)

    def _validate(self, n):
        """
        Verify that n is positive and has at most 4095 digits.

        INPUT:

        - ``n`` -- integer.

        OUTPUT:

        The integer as a Sage integer.  This function raises a
        ValueError if the two conditions listed above are not both
        satisfied.  It is here because GMP-ECM silently ignores all
        digits of input after the 4095th!

        EXAMPLES::

            sage: ecm = ECM()
            sage: ecm._validate(3)
            3
            sage: ecm._validate(0)
            Traceback (most recent call last):
            ...
            ValueError: n must be positive
            sage: ecm._validate(10^5000)
            Traceback (most recent call last):
            ...
            ValueError: n must have at most 4095 digits
        """
        n = ZZ(n)
        if n <= 0:
            raise ValueError("n must be positive")
        if n.ndigits() > 4095:
            raise ValueError("n must have at most 4095 digits")
        return n

    def __init__(self, data, V=None):
        if isinstance(data, (tuple, list)):
            self._dict = {}
            for mon, mult in data:
                mon = V(mon)
                if mon.is_zero():
                    raise ValueError('zero in denominator')
                mult = ZZ(mult)
                if mult <= 0:
                    raise ValueError('non-positive multiplicity in denominator')
                mon.set_immutable()
                self._dict[mon] = mult

        elif isinstance(data, dict):
            self._dict = data

        elif isinstance(data, FactoredDenominator):
            self._dict = data._dict
            self._tuple = data._tuple
            return

        elif isinstance(data, Vector_integer_dense):
            data = V(data)
            data.set_immutable()
            self._dict = {data: ZZ.one()}
            self._tuple = ((data, ZZ.one()),)
            return

        else:
            raise TypeError('invalid data of type {} to initialized a FactoredDenominator'.format(type(data)))

        self._tuple = tuple(sorted(self._dict.items()))

    def lift(self, P, prec=20):
        r"""
        Given a point `P` in the formal group of the elliptic curve `E` with split multiplicative reduction,
        this produces an element `u` in `\QQ_p^{\times}` mapped to the point `P` by the Tate parametrisation.
        The algorithm return the unique such element in `1+p\ZZ_p`.

        INPUT:

        - ``P`` - a point on the elliptic curve.

        - ``prec`` - the `p`-adic precision, default is 20.

        EXAMPLES::

            sage: e = EllipticCurve('130a1')
            sage: eq = e.tate_curve(5)
            sage: P = e([-6,10])
            sage: l = eq.lift(12*P, prec=10); l
            1 + 4*5 + 5^3 + 5^4 + 4*5^5 + 5^6 + 5^7 + 4*5^8 + 5^9 + O(5^10)

        Now we map the lift l back and check that it is indeed right.::

            sage: eq.parametrisation_onto_original_curve(l)
            (4*5^-2 + 2*5^-1 + 4*5 + 3*5^3 + 5^4 + 2*5^5 + 4*5^6 + O(5^7) : 2*5^-3 + 5^-1 + 4 + 4*5 + 5^2 + 3*5^3 + 4*5^4 + O(5^6) : 1 + O(5^20))
            sage: e5 = e.change_ring(Qp(5,9))
            sage: e5(12*P)
            (4*5^-2 + 2*5^-1 + 4*5 + 3*5^3 + 5^4 + 2*5^5 + 4*5^6 + O(5^7) : 2*5^-3 + 5^-1 + 4 + 4*5 + 5^2 + 3*5^3 + 4*5^4 + O(5^6) : 1 + O(5^9))
        """
        p = self._p
        R = Qp(self._p, prec)
        if not self._E == P.curve():
            raise ValueError("The point must lie on the original curve.")
        if not self.is_split():
            raise ValueError("The curve must have split multiplicative reduction.")
        if P.is_zero():
            return R.one()
        if P[0].valuation(p) >= 0:
            raise ValueError("The point must lie in the formal group.")

        Eq = self.curve(prec=prec)
        C, r, s, t = self._isomorphism(prec=prec)
        xx = r + C ** 2 * P[0]
        yy = t + s * C ** 2 * P[0] + C ** 3 * P[1]
        try:
            Eq([xx, yy])
        except Exception:
            raise RuntimeError("Bug : Point %s does not lie on the curve " %
                               (xx, yy))

        tt = -xx / yy
        eqhat = Eq.formal()
        eqlog = eqhat.log(prec + 3)
        z = eqlog(tt)
        u = ZZ.one()
        fac = ZZ.one()
        for i in range(1, 2 * prec + 1):
            fac *= i
            u += z ** i / fac
        return u

def is_hyperbolic(self, p):
    r"""
    Check if the quadratic form is a sum of hyperbolic planes over
    the `p`-adic numbers `\QQ_p` or over the real numbers `\RR`.

    REFERENCES:

    This criteria follows from Cassels's "Rational Quadratic Forms":

    - local invariants for hyperbolic plane (Lemma 2.4, p58)
    - direct sum formulas (Lemma 2.3, p58)

    INPUT:

    - `p` -- a prime number > 0 or `-1` for the infinite place

    OUTPUT:

    boolean

    EXAMPLES::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1])
        sage: Q.is_hyperbolic(-1)
        False
        sage: Q.is_hyperbolic(2)
        False
        sage: Q.is_hyperbolic(3)
        False
        sage: Q.is_hyperbolic(5)     ## Here -1 is a square, so it's true.
        True
        sage: Q.is_hyperbolic(7)
        False
        sage: Q.is_hyperbolic(13)    ## Here -1 is a square, so it's true.
        True
    """
    ## False for odd-dim'l forms
    if self.dim() % 2:
        return False

    ## True for the zero form
    if not self.dim():
        return True

    ## Compare local invariants
    ## Note: since the dimension is even, the extra powers of 2 in
    ##        self.det() := Det(2*Q) don't affect the answer!
    m = ZZ(self.dim() // 2)
    if p == -1:
        return self.signature() == 0

    if p == 2:
        return (QQ(self.det() * (-1) ** m).is_padic_square(p) and
                self.hasse_invariant(p) ==
                (-1) ** m.binomial(2))  # here -1 is hilbert_symbol(-1,-1,2)

    return (QQ(self.det() * (-1) ** m).is_padic_square(p) and
            self.hasse_invariant(p) == 1)

def iterator_fast(n, l):
    """
    Iterate over all ``l`` weighted integer vectors with total weight ``n``.

    INPUT:

    - ``n`` -- an integer
    - ``l`` -- the weights in weakly decreasing order

    EXAMPLES::

        sage: from sage.combinat.integer_vector_weighted import iterator_fast
        sage: list(iterator_fast(3, [2,1,1]))
        [[1, 1, 0], [1, 0, 1], [0, 3, 0], [0, 2, 1], [0, 1, 2], [0, 0, 3]]
        sage: list(iterator_fast(2, [2]))
        [[1]]

    Test that :trac:`20491` is fixed::

        sage: type(list(iterator_fast(2, [2]))[0][0])
        <type 'sage.rings.integer.Integer'>
    """
    if n < 0:
        return

    zero = ZZ.zero()
    one = ZZ.one()

    if not l:
        if n == 0:
            yield []
        return
    if len(l) == 1:
        if n % l[0] == 0:
            yield [n // l[0]]
        return

    k = 0
    cur = [n // l[k] + one]
    rem = n - cur[-1] * l[k] # Amount remaining
    while cur:
        cur[-1] -= one
        rem += l[k]
        if rem == zero:
            yield cur + [zero] * (len(l) - len(cur))
        elif cur[-1] < zero or rem < zero:
            rem += cur.pop() * l[k]
            k -= 1
        elif len(l) == len(cur) + 1:
            if rem % l[-1] == zero:
                yield cur + [rem // l[-1]]
        else:
            k += 1
            cur.append(rem // l[k] + one)
            rem -= cur[-1] * l[k]

            def rank(self, x):
                r"""
                Return the rank of an element of this Cartesian product.

                The *rank* of ``x`` is its position in the
                enumeration. It is an integer between ``0`` and
                ``n-1`` where ``n`` is the cardinality of this set.

                .. SEEALSO::

                    - :meth:`EnumeratedSets.ParentMethods.rank`
                    - :meth:`unrank`

                EXAMPLES::

                    sage: C = cartesian_product([GF(2), GF(11), GF(7)])
                    sage: C.rank(C((1,2,5)))
                    96
                    sage: C.rank(C((0,0,0)))
                    0

                    sage: for c in C: print(C.rank(c))
                    0
                    1
                    2
                    3
                    4
                    5
                    ...
                    150
                    151
                    152
                    153

                    sage: F1 = FiniteEnumeratedSet('abcdefgh')
                    sage: F2 = IntegerRange(250)
                    sage: F3 = Partitions(20)
                    sage: C = cartesian_product([F1, F2, F3])
                    sage: c = C(('a', 86, [7,5,4,4]))
                    sage: C.rank(c)
                    54213
                    sage: C.unrank(54213)
                    ('a', 86, [7, 5, 4, 4])
                """
                from builtins import zip
                from sage.rings.integer_ring import ZZ
                x = self(x)
                b = ZZ.one()
                rank = ZZ.zero()
                for f, c in zip(reversed(x.cartesian_factors()),
                                reversed(self.cartesian_factors())):
                    rank += b * c.rank(f)
                    b *= c.cardinality()
                return rank

 def getcoeff(self, key, **kwds) :
     (ch, k) = key
     if ch != self.__ch :
         return ZZ.zero()
     
     if self.__series is None :
         self.__series = \
          self.__factory.by_taylor_expansion( self.__fs, self.__weight )
     
     try :
         return self.__series[k]
     except KeyError :
         return ZZ.zero()

    def Heisenberg(p):
        r"""
        Return the Heisenberg origami.

        EXAMPLES::

            sage: from surface_dynamics.all import *

            sage: h2 = origamis.Heisenberg(2)
            sage: h2.stratum_component()
            H_3(1^4)^c
            sage: h2.veech_group().index()
            3
            sage: h2.sum_of_lyapunov_exponents()
            2

            sage: h3 = origamis.Heisenberg(3)
            sage: h3.stratum_component()
            H_10(2^9)^even
            sage: h3.veech_group().index()
            1
            sage: h3.sum_of_lyapunov_exponents()
            5

            sage: h5 = origamis.Heisenberg(5)
            sage: h5.stratum_component()
            H_51(4^25)^even
            sage: h5.veech_group().index()
            1
            sage: h5.sum_of_lyapunov_exponents()
            15
        """
        p = ZZ(p)
        if not p.is_prime():
            raise ValueError("p (={}) must be prime".format(p))
        N = p**3  # map (a,b,d) -> a*p**2 + b*p + d
        pp = p*p
        r = [None]*N
        u = [None]*N
        for a in xrange(p):
            for b in xrange(p):
                for d in xrange(p):
                    n = a*pp + b*p + d
                    r[n] = ((a+1)%p)*pp + b*p + d
                    u[n] = a*pp + ((b+1)%p)*p + ((a+d)%p)
        return Origami(r,u,
                as_tuple=True,
                name="Heisenberg origami on GF(%d)"%p)

        def number_of_reflections(self):
            r"""
            Return the number of reflections of ``self``.

            For real groups, this coincides with the number of
            reflection hyperplanes.

            This implementation uses that it is given by the sum of
            the degrees of ``self`` minus its rank.

            .. SEEALSO:: :meth:`number_of_reflection_hyperplanes`

            EXAMPLES::

                sage: [SymmetricGroup(i).number_of_reflections() for i in range(int(8))]
                [0, 0, 1, 3, 6, 10, 15, 21]

                sage: W = ColoredPermutations(1,3)
                sage: W.number_of_reflections()
                3
                sage: W = ColoredPermutations(2,3)
                sage: W.number_of_reflections()
                9
                sage: W = ColoredPermutations(4,3)
                sage: W.number_of_reflections()
                21
                sage: W = ReflectionGroup((4,2,3))         # optional - gap3
                sage: W.number_of_reflections()            # optional - gap3
                15
            """
            from sage.rings.all import ZZ
            return ZZ.sum(deg - 1 for deg in self.degrees())

def build_Lambdalist_from_AB(A,B,T, scaling):
    # T is the matrix of Hecke acting on Homology
    try:
        x,y,z,t = T.list()
    except AttributeError:
        x,y,z,t = T
    alpha = ZZ(z)/ZZ(y)
    beta = ZZ(t-x)/ZZ(y)
    d = alpha.denominator().lcm(beta.denominator())
    alpha, beta = ZZ(d*alpha), ZZ(d*beta)
    ans = []
    K = A.parent()
    for A0, B0 in product(our_nroot(A,scaling,return_all=True),our_nroot(B,scaling,return_all=True)):
        for B1 in our_nroot(B0, d,return_all=True):
            ans.append(Matrix(K,2,2,[A0,B0,B1**alpha,A0*B1**beta]))
    return ans

    def _kf(self, i):
        r"""
        Return the value `k_f` with respect to ``i`` and ``self``.

        INPUT:

        - ``i`` -- an element of the index set

        EXAMPLES::

            sage: M = crystals.infinity.NakajimaMonomials(['F',4,1])
            sage: m = M.module_generators[0].f_string([0,1,4,3])
            sage: [m._kf(i) for i in M.index_set()]
            [0, 0, 2, 0, 0]
        """
        if all(i != x[0] for x in self._Y):
            return ZZ.zero()

        d = copy(self._Y)
        K = max(key[1] for key in d if key[0] == i)
        for a in range(K):
            if (i,a) in d:
                continue
            else:
                d[(i,a)] = 0
        S = sorted((x for x in six.iteritems(d) if x[0][0] == i), key=lambda x: x[0][1])
        sum = 0
        phi = self.phi(i)
        for var,exp in S:
            sum += exp
            if sum == phi:
                return var[1]

        def character_ring(self, modules_base_ring = None, base_ring = None, side = 'right', q = None):
            """
            Returns the character ring of ``self``

            EXAMPLES::

                sage: S = Semigroups().Finite().example()
                sage: S.rename("S")

                sage: A = S.character_ring(); A
                The right-character ring of S over Rational Field

                sage: A.category()
                Category of character rings of finite semigroups over Integer Ring with realizations

                sage: QQ.<q> = QQ[]
                sage: A = S.character_ring(q = q); A
                The right-character ring of S over Rational Field
            """
            from sage.rings.rational_field import QQ
            if modules_base_ring is None:
                modules_base_ring = QQ
            from sage.rings.integer_ring import ZZ
            if q is None:
                q = ZZ.one()
            if base_ring is None:
                base_ring = q.parent()
            from sage.categories.rings import Rings
            assert base_ring in Rings()
            assert q in base_ring
            from sage_semigroups.monoids.character_ring import AbstractCharacterRing
            return AbstractCharacterRing(self, modules_base_ring, base_ring, side = side, q = q)

    def phi(self, i):
        r"""
        Return the value of `\varphi_i` on ``self``.

        INPUT:

        - ``i`` -- an element of the index set

        EXAMPLES::

            sage: M = crystals.infinity.NakajimaMonomials(['D',4,3])
            sage: m = M.module_generators[0].f(1)
            sage: [m.phi(i) for i in M.index_set()]
            [1, -1, 1]

            sage: M = crystals.infinity.NakajimaMonomials(['C',4,1])
            sage: m = M.module_generators[0].f_string([4,2,3])
            sage: [m.phi(i) for i in M.index_set()]
            [0, 1, -1, 2, -1]
        """
        if i not in self.parent().index_set():
            raise ValueError("i must be an element of the index set")
        if not self._Y or all(x[0] != i for x in self._Y):
            return ZZ.zero()

        d = copy(self._Y)
        K = max(x[1] for x in d if x[0] == i)
        for a in range(K):
            if (i,a) in d:
                continue
            else:
                d[(i,a)] = 0
        S = sorted((x for x in six.iteritems(d) if x[0][0] == i), key=lambda x: x[0][1])
        return max(sum(S[k][1] for k in range(s)) for s in range(1,len(S)+1))

    def degree(self):
        r"""
        EXAMPLES::

            sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator

            sage: V = ZZ**3
            sage: FactoredDenominator([], None).degree()
            0
            sage: FactoredDenominator([((1,0,0), 2)], V).degree()
            2
            sage: FactoredDenominator([((0,1,2), 3), ((1,1,1), 1)], V).degree()
            4
            sage: FactoredDenominator([((0,-1,2), 1), ((1,0,0), 1), ((0,0,2), 1)], V).degree()
            3

        TESTS::

            sage: parent(FactoredDenominator([], None).degree())
            Integer Ring
            sage: parent(FactoredDenominator([((1,0,0), 2)], V).degree())
            Integer Ring
        """
        if not self._tuple:
            return ZZ.zero()
        return sum(m for _,m in self._tuple)

        def number_of_reflection_hyperplanes(self):
            r"""
            Return the number of reflection hyperplanes of ``self``.

            This is also the number of distinguished reflections.  For
            real groups, this coincides with the number of
            reflections.

            This implementation uses that it is given by the sum of
            the codegrees of ``self`` plus its rank.

            .. SEEALSO:: :meth:`number_of_reflections`

            EXAMPLES::

                sage: W = ColoredPermutations(1,3)
                sage: W.number_of_reflection_hyperplanes()
                3
                sage: W = ColoredPermutations(2,3)
                sage: W.number_of_reflection_hyperplanes()
                9
                sage: W = ColoredPermutations(4,3)
                sage: W.number_of_reflection_hyperplanes()
                15
                sage: W = ReflectionGroup((4,2,3))          # optional - gap3
                sage: W.number_of_reflection_hyperplanes()  # optional - gap3
                15
            """
            from sage.rings.all import ZZ
            return ZZ.sum(codeg + 1 for codeg in self.codegrees())

def number_of_rooted_trees(n):
    r"""
    Return the number of rooted trees with `n` nodes.

    Compute the number `a(n)` of rooted trees with `n` nodes using the
    recursive formula ([SL000081]_):

    .. MATH::

        a(n+1) = \frac{1}{n} \sum_{k=1}^{n} \left( \sum_{d|k} d a(d) \right) a(n-k+1)

    EXAMPLES::

        sage: from sage.combinat.rooted_tree import number_of_rooted_trees
        sage: [number_of_rooted_trees(i) for i in range(10)]
        [0, 1, 1, 2, 4, 9, 20, 48, 115, 286]

    REFERENCES:

    .. [SL000081] Sloane's :oeis:`A000081`
    """
    if n == 0:
        return Integer(0)
    if n == 1:
        return Integer(1)
    n = Integer(n)
    return sum(sum(d * number_of_rooted_trees(d) for d in k.divisors()) *
               number_of_rooted_trees(n - k)
               for k in ZZ.range(1, n)) // (n - 1)

    def __contains__(self, x):
        r"""
        Tests if ``x`` is an element of ``self``.

        INPUT:

        - ``x`` -- matrix

        EXAMPLES::

            sage: from sage.combinat.integer_matrices import IntegerMatrices
            sage: IM = IntegerMatrices([4], [1,2,1])
            sage: matrix([[1, 2, 1]]) in IM
            True
            sage: matrix(QQ, [[1, 2, 1]]) in IM
            True
            sage: matrix([[2, 1, 1]]) in IM
            False

        TESTS::

            sage: from sage.combinat.integer_matrices import IntegerMatrices
            sage: IM = IntegerMatrices([4], [1,2,1])
            sage: [1, 2, 1] in IM
            False
            sage: matrix([[-1, 3, 1]]) in IM
            False

        """
        from sage.matrix.matrix import Matrix

        if not isinstance(x, Matrix):
            return False
        row_sums = [ZZ.zero()] * x.nrows()
        col_sums = [ZZ.zero()] * x.ncols()
        for i in range(x.nrows()):
            for j in range(x.ncols()):
                x_ij = x[i, j]
                if x_ij not in ZZ or x_ij < 0:
                    return False
                row_sums[i] += x_ij
                col_sums[j] += x_ij
            if row_sums[i] != self._row_sums[i]:
                return False
        if col_sums != self._col_sums:
            return False
        return True