def _eval_line(self, line, allow_use_file=True, wait_for_prompt=True, restart_if_needed=False):
        """
        EXAMPLES::

            sage: gp._eval_line('2+2')
            '4'

        TESTS:

        We verify that trac 11617 is fixed::

            sage: gp._eval_line('a='+str(range(2*10^5)))[:70]
            '[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,'
        """
        line = line.strip()
        if len(line) == 0:
            return ''
        a = Expect._eval_line(self, line,
                              allow_use_file=allow_use_file,
                              wait_for_prompt=wait_for_prompt)
        if a.find("the PARI stack overflows") != -1:
            verbose("automatically doubling the PARI stack and re-executing current input line")
            b = self.eval("allocatemem()")
            if b.find("Warning: not enough memory") != -1:
                raise RuntimeError(a)
            return self._eval_line(line, allow_use_file=allow_use_file,
                                   wait_for_prompt=wait_for_prompt)
        else:
            return a

    def __call__(self, x):
        r"""
        Evaluate this character at an element of `\ZZ_p^\times`.

        EXAMPLES::

            sage: kappa = pAdicWeightSpace(23)(1 + 23^2 + O(23^20), 4, False)
            sage: kappa(2)
            16 + 7*23 + 7*23^2 + 16*23^3 + 23^4 + 20*23^5 + 15*23^7 + 11*23^8 + 12*23^9 + 8*23^10 + 22*23^11 + 16*23^12 + 13*23^13 + 4*23^14 + 19*23^15 + 6*23^16 + 7*23^17 + 11*23^19 + O(23^20)
            sage: kappa(-1)
            1 + O(23^20)
            sage: kappa(23)
            0
            sage: kappa(2 + 2*23 + 11*23^2 + O(23^3))
            16 + 7*23 + O(23^3)
        """

        if not isinstance(x, pAdicGenericElement):
            x = Qp(self._p)(x)
        if x.valuation() != 0:
            return 0

        teich = x.parent().teichmuller(x)
        xx = x / teich
        if (xx - 1).valuation() <= 0:
            raise ArithmeticError
        verbose("Normalised element is %s" % xx)

        e = xx.log() / self.parent()._param.log()
        verbose("Exponent is %s" % e)

        return teich**(self.t) * (self.w.log() * e).exp()

    def hecke_bound(self):
        r"""
        Return an integer B such that the Hecke operators `T_n`, for `n\leq B`,
        generate the full Hecke algebra as a module over the base ring. Note
        that we include the `n` with `n` not coprime to the level.

        At present this returns an unproven guess for non-cuspidal spaces which
        appears to be valid for `M_k(\Gamma_0(N))`, where k and N are the
        weight and level of self. (It is clearly valid for *cuspidal* spaces
        of any fixed character, as a consequence of the Sturm bound theorem.)
        It returns a hopelessly wrong answer for spaces of full level
        `\Gamma_1`.

        TODO: Get rid of this dreadful bit of code.

        EXAMPLE::

            sage: ModularSymbols(17, 4).hecke_bound()
            15
            sage: ModularSymbols(Gamma1(17), 4).hecke_bound() # wrong!
            15
        """
        try:
            if self.is_cuspidal():
                return Gamma0(self.level()).sturm_bound(self.weight())
        except AttributeError:
            pass
        misc.verbose("WARNING: ambient.py -- hecke_bound; returning unproven guess.")
        return Gamma0(self.level()).sturm_bound(self.weight()) + 2*Gamma0(self.level()).dimension_eis(self.weight()) + 5

    def _eval_line(self, line, reformat=True, allow_use_file=False,
                   wait_for_prompt=True, restart_if_needed=False):
        """
        EXAMPLES::

            sage: print octave._eval_line('2+2')  #optional - octave
              ans =  4
        """
        if not wait_for_prompt:
            return Expect._eval_line(self, line)
        if line == '':
            return ''
        if self._expect is None:
            self._start()
        if allow_use_file and len(line)>3000:
            return self._eval_line_using_file(line)
        try:
            E = self._expect
            verbose("in = '%s'"%line,level=3)
            E.sendline(line)
            E.expect(self._prompt)
            out = E.before
            verbose("out = '%s'"%out,level=3)
        except EOF:
            if self._quit_string() in line:
                return ''
        except KeyboardInterrupt:
            self._keyboard_interrupt()
        if reformat:
            if '>>> ' in out and 'syntax error' in out:
                raise SyntaxError(out)
        out = "\n".join(out.splitlines()[1:])
        return out

    def _next_var_name(self):
        """
        Return the name of the next unused interface variable name.

        EXAMPLES::
        
            sage: g = Gp()
            sage: g._next_var_name()
            'sage[1]'
            sage: g(2)^2
            4
            sage: g._next_var_name()
            'sage[5]'
        """
        self.__seq += 1
        self.__seq %= 1000
        #print 'wtf: %s' % self.__seq
        if self.__seq >= self.__var_store_len:
            if self.__var_store_len == 0:
                self.eval('sage=vector(%s,k,0);'%self.__init_list_length)
                self.__var_store_len = self.__init_list_length
            else:
                self.eval('sage=concat(sage, vector(%s,k,0));'%self.__var_store_len)
                self.__var_store_len *= 2
                verbose("doubling PARI/sage object vector: %s"%self.__var_store_len)
        return 'sage[%s]'%self.__seq

def add_row(A, b, pivots, include_zero_rows):
    """
    The add row procedure.

    INPUT:
        A -- a matrix in Hermite normal form with n column
        b -- an n x 1 row matrix
        pivots -- sorted list of integers; the pivot positions of A.

    OUTPUT:
        H -- the Hermite normal form of A.stack(b).
        new_pivots -- the pivot columns of H.

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: A = matrix(ZZ, 2, 3, [-21, -7, 5, 1,20,-7])
        sage: b = matrix(ZZ, 1,3, [-1,1,-1])
        sage: hnf.add_row(A, b, A.pivots(), True)
        (
        [ 1  6 29]
        [ 0  7 28]
        [ 0  0 46], [0, 1, 2]
        )
        sage: A.stack(b).echelon_form()
        [ 1  6 29]
        [ 0  7 28]
        [ 0  0 46]
    """
    t = verbose('add hnf row')
    v = b.row(0)
    H, pivs = A._add_row_and_maintain_echelon_form(b.row(0), pivots)
    if include_zero_rows and H.nrows() != A.nrows()+1:
        H = H.matrix_from_rows(range(A.nrows()+1))
    verbose('finished add hnf row', t)
    return H, pivs

    def cuspidal_submodule_q_expansion_basis(self, weight, prec=None):
        r"""
        Calculate a basis of `q`-expansions for the space of cusp forms of
        weight ``weight`` for this group.

        INPUT:

        - ``weight`` (integer) -- the weight
        - ``prec`` (integer or None) -- precision of `q`-expansions to return

        ALGORITHM: Uses the method :meth:`cuspidal_ideal_generators` to
        calculate generators of the ideal of cusp forms inside this ring. Then
        multiply these up to weight ``weight`` using the generators of the
        whole modular form space returned by :meth:`q_expansion_basis`.

        EXAMPLES::

            sage: R = ModularFormsRing(Gamma0(3))
            sage: R.cuspidal_submodule_q_expansion_basis(20)
            [q - 8532*q^6 - 88442*q^7 + O(q^8), q^2 + 207*q^6 + 24516*q^7 + O(q^8), q^3 + 456*q^6 + O(q^8), q^4 - 135*q^6 - 926*q^7 + O(q^8), q^5 + 18*q^6 + 135*q^7 + O(q^8)]

        We compute a basis of a space of very large weight, quickly (using this
        module) and slowly (using modular symbols), and verify that the answers
        are the same. ::

            sage: A = R.cuspidal_submodule_q_expansion_basis(80, prec=30)  # long time (1s on sage.math, 2013)
            sage: B = R.modular_forms_of_weight(80).cuspidal_submodule().q_expansion_basis(prec=30)  # long time (19s on sage.math, 2013)
            sage: A == B # long time
            True
        """
        d = self.modular_forms_of_weight(weight).cuspidal_submodule().dimension()
        if d == 0: return []

        minprec = self.modular_forms_of_weight(weight).sturm_bound()
        if prec is None:
            prec = working_prec = minprec
        else:
            working_prec = max(prec, minprec)

        gen_weight = min(6, weight)

        while 1:
            verbose("Trying to generate the %s-dimensional cuspidal submodule at weight %s using generators of weight up to %s" % (d, weight, gen_weight))
            G = self.cuspidal_ideal_generators(maxweight=gen_weight, prec=working_prec)

            flist = []
            for (j, f, F) in G:
                for g in self.q_expansion_basis(weight - j, prec=working_prec):
                    flist.append(g*f)

            A = self.base_ring() ** working_prec
            W = A.span([A(f.padded_list(working_prec)) for f in flist])
            if W.rank() == d and (self.base_ring().is_field() or W.index_in_saturation() == 1):
                break
            else:
                gen_weight += 1
                verbose("Need more generators: trying again with generators of weight up to %s" % gen_weight)

        R = G[0][1].parent()
        return [R(list(x), prec=prec) for x in W.gens()]

    def Tq_eigenvalue(self, q, p=None, M=None, check=True):
        r"""
        Eigenvalue of `T_q` modulo `p^M`

        INPUT:

        - ``q`` -- prime of the Hecke operator
        - ``p`` -- prime we are working modulo (default: None)
        - ``M`` -- degree of accuracy of approximation (default: None)
        - ``check`` --

        OUTPUT:

        - Constant `c` such that `self|T_q - c * self` has valuation greater than
          or equal to `M` (if it exists), otherwise raises ValueError

        EXAMPLES::

            sage: E = EllipticCurve('11a')
            sage: from sage.modular.pollack_stevens.space import ps_modsym_from_elliptic_curve
            sage: phi = ps_modsym_from_elliptic_curve(E)
            sage: phi.values()
            [-1/5, 3/2, -1/2]
            sage: phi_ord = phi.p_stabilize(p = 3, ap = E.ap(3), M = 10, ordinary = True)
            sage: phi_ord.Tq_eigenvalue(2,3,10) + 2
            O(3^10)

            sage: phi_ord.Tq_eigenvalue(3,3,10)
            2 + 3^2 + 2*3^3 + 2*3^4 + 2*3^6 + 3^8 + 2*3^9 + O(3^10)
            sage: phi_ord.Tq_eigenvalue(3,3,100)
            Traceback (most recent call last):
            ...
            ValueError: not a scalar multiple
        """
        qhecke = self.hecke(q)
        gens = self.parent().source().gens()
        if p is None:
            p = self.parent().prime()
        i = 0
        g = gens[i]
        verbose("Computing eigenvalue")
        while self._map[g].is_zero(p, M):
            if not qhecke._map[g].is_zero(p, M):
                raise ValueError("not a scalar multiple")
            i += 1
            try:
                g = gens[i]
            except IndexError:
                raise ValueError("self is zero")
        aq = self._map[g].find_scalar(qhecke._map[g], p, M, check)
        if check:
            verbose("Checking that this is actually an eigensymbol")
            if p is None or M is None:
                for g in gens[1:]:
                    if qhecke._map[g] != aq * self._map[g]:
                        raise ValueError("not a scalar multiple")
            elif (qhecke - aq * self).valuation(p) < M:
                raise ValueError("not a scalar multiple")
        return aq

def extract_ones_data(H, pivots):
    """
    Compute ones data and corresponding submatrices of H.  This is
    used to optimized the add_row function.

    INPUT:

    - H -- a matrix in HNF
    - pivots -- list of all pivot column positions of H

    OUTPUT:

    C, D, E, onecol, onerow, non_onecol, non_onerow
    where onecol, onerow, non_onecol, non_onerow are as for
    the ones function, and C, D, E are matrices:

    - C -- submatrix of all non-onecol columns and onecol rows
    - D -- all non-onecol columns and other rows
    - E -- inverse of D

    If D isn't invertible or there are 0 or more than 2 non onecols,
    then C, D, and E are set to None.

    EXAMPLES::

        sage: H = matrix(ZZ, 3, 4, [1, 0, 0, 7, 0, 1, 5, 2, 0, 0, 6, 6])
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf
        sage: matrix_integer_dense_hnf.extract_ones_data(H, [0,1,2])
        (
        [0]
        [5], [6], [1/6], [0, 1], [0, 1], [2], [2]
        )

    Here we get None's since the (2,2) position submatrix is not invertible.
        sage: H = matrix(ZZ, 3, 5, [1, 0, 0, 45, -36, 0, 1, 0, 131, -107, 0, 0, 0, 178, -145]); H
        [   1    0    0   45  -36]
        [   0    1    0  131 -107]
        [   0    0    0  178 -145]
        sage: import sage.matrix.matrix_integer_dense_hnf as matrix_integer_dense_hnf
        sage: matrix_integer_dense_hnf.extract_ones_data(H, [0,1,3])
        (None, None, None, [0, 1], [0, 1], [2], [2])
    """
    onecol, onerow, non_onecol, non_onerow = ones(H, pivots)
    verbose('extract_ones -- got submatrix of size %s'%len(non_onecol))
    if len(non_onecol) in [1,2]:
        # Extract submatrix of all non-onecol columns and onecol rows
        C = H.matrix_from_rows_and_columns(onerow, non_onecol)
        # Extract submatrix of all non-onecol columns and other rows
        D = H.matrix_from_rows_and_columns(non_onerow, non_onecol).transpose()
        tt = verbose("extract ones -- INVERT %s x %s"%(len(non_onerow), len(non_onecol)), level=1)
        try:
            E = D**(-1)
        except ZeroDivisionError:
            C = D = E = None
        verbose("done inverting", tt, level=1)
        return C, D, E, onecol, onerow, non_onecol, non_onerow
    else:
        return None, None, None, onecol, onerow, non_onecol, non_onerow

    def __call__(self, z, prec=None):
        r"""
        Evaluate ``self`` at a point `z \in X_0(N)` where `z` is given by a
        representative in the upper half plane.

        All computations are done with ``prec``
        bits of precision. If ``prec`` is not given, use the precision of `z`.

        EXAMPLES::

            sage: E = EllipticCurve('37a')
            sage: phi = E.modular_parametrization()
            sage: phi((sqrt(7)*I - 17)/74, 53)
            (...e-16 - ...e-16*I : ...e-16 + ...e-16*I : 1.00000000000000)

        Verify that the mapping is invariant under the action of `\Gamma_0(N)`
        on the upper half plane::

            sage: E = EllipticCurve('11a')
            sage: phi = E.modular_parametrization()
            sage: tau = CC((1+1j)/5)
            sage: phi(tau)
            (-3.92181329652811 - 12.2578555525366*I : 44.9649874434872 + 14.3257120944681*I : 1.00000000000000)
            sage: phi(tau+1)
            (-3.92181329652810 - 12.2578555525366*I : 44.9649874434872 + 14.3257120944681*I : 1.00000000000000)
            sage: phi((6*tau+1) / (11*tau+2))
            (-3.9218132965285... - 12.2578555525369*I : 44.964987443489... + 14.325712094467...*I : 1.00000000000000)

        We can also apply the modular parametrization to a Heegner point on `X_0(N)`::

            sage: H = heegner_points(389,-7,5); H
            All Heegner points of conductor 5 on X_0(389) associated to QQ[sqrt(-7)]
            sage: x = H[0]; x
            Heegner point 5/778*sqrt(-7) - 147/778 of discriminant -7 and conductor 5 on X_0(389)
            sage: E = EllipticCurve('389a'); phi = E.modular_parametrization()
            sage: phi(x)
            Heegner point of discriminant -7 and conductor 5 on elliptic curve of conductor 389
            sage: phi(x).quadratic_form()
            389*x^2 + 147*x*y + 14*y^2


        ALGORITHM:

            Integrate the modular form attached to this elliptic curve from
            `z` to `\infty` to get a point on the lattice representation of
            `E`, then use the Weierstrass `\wp` function to map it to the
            curve itself.
        """
        if isinstance(z, heegner.HeegnerPointOnX0N):
            return z.map_to_curve(self.curve())
        # Map to the CC of CC/PeriodLattice.
        tm = verbose("Evaluating modular parameterization to precision %s bits" % prec)
        w = self.map_to_complex_numbers(z, prec=prec)
        # Map to E via Weierstrass P
        z = self._E.elliptic_exponential(w)
        verbose("Finished evaluating modular parameterization", tm)
        return z

    def upper_bound_on_elliptic_factors(self, p=None, ellmax=2):
        r"""
        Return an upper bound (provably correct) on the number of
        elliptic curves of conductor equal to the level of this
        supersingular module.

        INPUT:

        - ``p`` - (default: 997) prime to work modulo

        ALGORITHM: Currently we only use `T_2`.  Function will be
        extended to use more Hecke operators later.

        The prime p is replaced by the smallest prime that doesn't
        divide the level.

        EXAMPLE::

            sage: SupersingularModule(37).upper_bound_on_elliptic_factors()
            2

        (There are 4 elliptic curves of conductor 37, but only 2 isogeny
        classes.)
        """
        # NOTE: The heuristic runtime is *very* roughly `p^2/(2\cdot 10^6)`.
        # ellmax -- (default: 2) use Hecke operators T_ell with ell <= ellmax
        if p is None:
            p = 997

        while self.level() % p == 0:
            p = rings.next_prime(p)

        ell = 2
        t = self.hecke_matrix(ell).change_ring(rings.GF(p))

        # TODO: temporarily try using sparse=False
        # turn this off when sparse rank is optimized.
        t = t.dense_matrix()

        B = 2 * math.sqrt(ell)
        bnd = 0
        lower = -int(math.floor(B))
        upper = int(math.floor(B)) + 1
        for a in range(lower, upper):
            tm = verbose("computing T_%s - %s" % (ell, a))
            if a == lower:
                c = a
            else:
                c = 1
            for i in range(t.nrows()):
                t[i, i] += c
            tm = verbose("computing kernel", tm)
            # dim = t.kernel().dimension()
            dim = t.nrows() - t.rank()
            bnd += dim
            verbose("got dimension = %s; new bound = %s" % (dim, bnd), tm)
        return bnd

    def _find_alpha(self, p, k, M=None, ap=None, new_base_ring=None, ordinary=True, check=True, find_extraprec=True):
        """
        Finds `alpha`, a `U_p` eigenvalue, which is found as a root of
        the polynomial `x^2 - ap * x + p^(k+1)`.
        
        INPUT:

        - ``p`` -- prime
        - ``k`` -- Pollack-Stevens weight
        - ``M`` -- precision (default = None) of `Q_p`
        - ``ap`` -- Hecke eigenvalue at p (default = None)
        - ``new_base_ring`` -- field of definition of `alpha` (default = None)
        - ``ordinary`` -- True if the prime is ordinary (default = True)
        - ``check`` -- check to see if the prime is ordinary (default = True)
        - ``find_extraprec`` -- setting this to True finds extra precision (default = True)

        OUTPUT:

        - ``alpha`` --  `U_p` eigenvalue
        - ``new_base_ring`` -- field of definition of `alpha` with precision at least `newM`
        - ``newM`` -- new precision
        - ``eisenloss`` -- loss of precision
        - ``q`` -- a prime not equal to p which was used to find extra precision
        - ``aq`` -- the Hecke eigenvalue `aq` corresponding to `q`

        EXAMPLES::

            sage: from sage.modular.pollack_stevens.space import ps_modsym_from_elliptic_curve
            sage: E = EllipticCurve('11a')
            sage: p = 5
            sage: M = 10
            sage: k = 0
            sage: phi = ps_modsym_from_elliptic_curve(E)
            sage: phi._find_alpha(p,k,M)
            (1 + 4*5 + 3*5^2 + 2*5^3 + 4*5^4 + 4*5^5 + 4*5^6 + 3*5^7 + 2*5^8 + 3*5^9 + 3*5^10 + 3*5^12 + O(5^13), 5-adic Field with capped relative precision 13, 12, 1, None, None)

        """
        if ap is None:
            ap = self.Tq_eigenvalue(p, check=check)
        if check and ap.valuation(p) > 0:
            raise ValueError("p is not ordinary")
        poly = PolynomialRing(ap.parent(), 'x')([p**(k+1), -ap, 1])
        if new_base_ring is None:
            # These should actually be completions of disc.parent()
            if p == 2:
                # is this the right precision adjustment for p=2?
                new_base_ring = Qp(2, M+1)
            else:
                new_base_ring = Qp(p, M)
            set_padicbase = True
        else:
            set_padicbase = False
        try:
            verbose("finding alpha: rooting %s in %s"%(poly, new_base_ring))
            (v0,e0),(v1,e1) = poly.roots(new_base_ring)
        except TypeError, ValueError:
            raise ValueError("new base ring must contain a root of x^2 - ap * x + p^(k+1)")

def modS_relations(syms):
    """
    Compute quotient of Manin symbols by the S relations.

    Here S is the 2x2 matrix [0, -1; 1, 0].

    INPUT:

    -  ``syms`` - manin_symbols.ManinSymbols

    OUTPUT:

    -  ``rels`` - set of pairs of pairs (j, s), where if
       mod[i] = (j,s), then x_i = s\*x_j (mod S relations)

    EXAMPLES::

        sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_gamma0
        sage: from sage.modular.modsym.relation_matrix import modS_relations

    ::

        sage: syms = ManinSymbolList_gamma0(2, 4); syms
        Manin Symbol List of weight 4 for Gamma0(2)
        sage: modS_relations(syms)
        set([((3, -1), (4, 1)), ((5, -1), (5, 1)), ((1, 1), (6, 1)), ((0, 1), (7, 1)), ((3, 1), (4, -1)), ((2, 1), (8, 1))])

    ::

        sage: syms = ManinSymbolList_gamma0(7, 2); syms
        Manin Symbol List of weight 2 for Gamma0(7)
        sage: modS_relations(syms)
        set([((3, 1), (4, 1)), ((2, 1), (7, 1)), ((5, 1), (6, 1)), ((0, 1), (1, 1))])

    Next we do an example with Gamma1::

        sage: from sage.modular.modsym.manin_symbols import ManinSymbolList_gamma1
        sage: syms = ManinSymbolList_gamma1(3,2); syms
        Manin Symbol List of weight 2 for Gamma1(3)
        sage: modS_relations(syms)
        set([((3, 1), (6, 1)), ((0, 1), (5, 1)), ((0, 1), (2, 1)), ((3, 1), (4, 1)), ((6, 1), (7, 1)), ((1, 1), (2, 1)), ((1, 1), (5, 1)), ((4, 1), (7, 1))])
    """
    if not isinstance(syms, manin_symbols.ManinSymbolList):
        raise TypeError, "syms must be a ManinSymbolList"
    tm = misc.verbose()
    # We will fill in this set with the relations x_i + s*x_j = 0,
    # where the notation is as in _sparse_2term_quotient.
    rels = set()
    for i in xrange(len(syms)):
        j, s = syms.apply_S(i)
        assert j != -1
        if i < j:
            rels.add( ((i,1),(j,s)) )
        else:
            rels.add( ((j,s),(i,1)) )
    misc.verbose("finished creating S relations",tm)
    return rels

def double_det (A, b, c, proof):
    """
    Compute the determinants of the stacked integer matrices
    A.stack(b) and A.stack(c).

    INPUT:

    - A -- an (n-1) x n matrix
    - b -- an 1 x n matrix
    - c -- an 1 x n matrix
    - proof -- whether or not to compute the det modulo enough times to
      provably compute the determinant.

    OUTPUT:

    - a pair of two integers.

    EXAMPLES::

        sage: from sage.matrix.matrix_integer_dense_hnf import double_det
        sage: A = matrix(ZZ, 2, 3, [1,2,3, 4,-2,5])
        sage: b = matrix(ZZ, 1, 3, [1,-2,5])
        sage: c = matrix(ZZ, 1, 3, [8,2,10])
        sage: A.stack(b).det()
        -48
        sage: A.stack(c).det()
        42
        sage: double_det(A, b, c, False)
        (-48, 42)
    """
    # We use the "two for the price of one" algorithm, which I made up. (William Stein)

    # This is a clever trick!  First we transpose everything.  Then
    # we use that if [A|b]*v = c then [A|c]*w = b with w easy to write down!
    # In fact w is got from v by dividing all entries by -v[n], where n is the
    # number of rows of v, and *also* dividing the last entry of w by v[n] again.
    # See this as an algebra exercise where you have to think of matrix vector
    # multiply as "linear combination of columns".
    A = A.transpose()
    b = b.transpose()
    c = c.transpose()
    t = verbose('starting double det')
    B = A.augment(b)
    v = B.solve_right(-c)

    db = det_given_divisor(B, v.denominator(), proof=proof)

    n = v.nrows()
    vn = v[n-1,0]
    w = (-1/vn)*v
    w[n-1] = w[n-1]/vn
    dc = det_given_divisor(A.augment(c), w.denominator(), proof=proof)

    verbose('finished double det', t)

    return (db, dc)

    def _lift_to_OMS(self, p, M, new_base_ring, check):
        """
        Returns a (`p`-adic) overconvergent modular symbol with `M` moments which lifts self up to an Eisenstein error

        Here the Eisenstein error is a symbol whose system of Hecke eigenvalues equals `ell+1` for `T_ell` when `ell`
        does not divide `Np` and 1 for `U_q` when `q` divides `Np`.

        INPUT:

        - ``p`` -- prime
        - ``M`` -- integer equal to the number of moments
        - ``new_base_ring`` -- new base ring

        OUTPUT:

        - An overconvergent modular symbol whose specialization equals self up to some Eisenstein error.

        EXAMPLES::


        """
        D = {}
        manin = self.parent().source()
        MSS = self.parent()._lift_parent_space(p, M, new_base_ring)
        verbose("Naive lifting: newM=%s, new_base_ring=%s"%(M, MSS.base_ring()))
        half = ZZ(1) / ZZ(2)
        for g in manin.gens()[1:]:
            twotor = g in manin.reps_with_two_torsion
            threetor = g in manin.reps_with_three_torsion
            if twotor:
                # See [PS] section 4.1
                gam = manin.two_torsion[g]
                mu = self._map[g].lift(p, M, new_base_ring)
                D[g] = (mu * gam - mu) * half
            elif threetor:
                # See [PS] section 4.1
                gam = manin.three_torsion[g]
                mu = self._map[g].lift(p, M, new_base_ring)
                D[g] = (2 * mu - mu * gam - mu * (gam**2)) * half
            else:
                # no two or three torsion
                D[g] = self._map[g].lift(p, M, new_base_ring)

        t = self.parent().coefficient_module().lift(p, M, new_base_ring).zero_element()
        for h in manin[2:]:
            R = manin.relations(h)
            if len(R) == 1:
                c, A, g = R[0]
                if c == 1:
                    t += self._map[h].lift(p, M, new_base_ring)
                elif A is not Id:
                    # rules out extra three torsion terms
                    t += c * self._map[g].lift(p, M, new_base_ring) * A
        D[manin.gen(0)] = t.solve_diff_eqn()  ###### Check this!
        return MSS(D)

    def _find_extraprec(self, p, M, alpha, check):
        q, aq, eisenloss = self._find_aq(p, M, check)
        newM = M + eisenloss

        # We also need to add precision to account for denominators appearing while solving the difference equation.
        eplog = (newM -1).exact_log(p)
        while eplog < (newM + eplog).exact_log(p):
            eplog = (newM + eplog).exact_log(p)
            verbose("M = %s, newM = %s, eplog=%s"%(M, newM, eplog), level=2)
        newM += eplog
        return newM, eisenloss, q, aq

def p_saturation(A, p, proof=True):
    """
    INPUT:

    - A -- a matrix over ZZ
    - p -- a prime
    - proof -- bool (default: True)

    OUTPUT:

    The p-saturation of the matrix A, i.e., a new matrix in Hermite form
    whose row span a ZZ-module that is p-saturated.

    EXAMPLES::

        sage: from sage.matrix.matrix_integer_dense_saturation import p_saturation
        sage: A = matrix(ZZ, 2, 2, [3,2,3,4]); B = matrix(ZZ, 2,3,[1,2,3,4,5,6])
        sage: A.det()
        6
        sage: C = A*B; C
        [11 16 21]
        [19 26 33]
        sage: C2 = p_saturation(C, 2); C2
        [ 1  8 15]
        [ 0  9 18]
        sage: C2.index_in_saturation()
        9
        sage: C3 = p_saturation(C, 3); C3
        [ 1  0 -1]
        [ 0  2  4]
        sage: C3.index_in_saturation()
        2
    """
    tm = verbose("%s-saturating a %sx%s matrix"%(p, A.nrows(), A.ncols()))
    H = A.hermite_form(include_zero_rows=False, proof=proof)
    while True:
        if p == 2:
            A = H.change_ring(GF(p))
        else:
            try:
                # Faster than change_ring
                A = H._reduce(p)
            except OverflowError:
                # fall back to generic GF(p) matrices
                A = H.change_ring(GF(p))
        assert A.nrows() <= A.ncols()
        K = A.kernel()
        if K.dimension() == 0:
            verbose("done saturating", tm)
            return H
        B = K.basis_matrix().lift()
        C = ((B * H) / p).change_ring(ZZ)
        H = H.stack(C).hermite_form(include_zero_rows=False, proof=proof)
    verbose("done saturating", tm)

def multiply_forms_to_weight(forms, weight, stop_dim=None):
    r"""
    Given a list of pairs ``(k,f)``, where `k` is an integer and `f` is a power
    series, and a weight l, return all weight l forms obtained by multiplying
    together the given forms. 

    INPUT:

    - forms -- list of pairs (k, f) with k an integer and f a power series
    - weight -- an integer
    - stop_dim -- integer (optional): if set to an integer and we find that the
      series so far span a space of at least this dimension, then stop
      multiplying more forms together.

    EXAMPLES::

        sage: import sage.modular.modform.find_generators as f
        sage: forms = [(4, 240*eisenstein_series_qexp(4,5)), (6,504*eisenstein_series_qexp(6,5))]
        sage: f.multiply_forms_to_weight(forms, 12)
        [(12, 1 - 1008*q + 220752*q^2 + 16519104*q^3 + 399517776*q^4 + O(q^5)), (12, 1 + 720*q + 179280*q^2 + 16954560*q^3 + 396974160*q^4 + O(q^5))]
        sage: f.multiply_forms_to_weight(forms, 24)
        [(24, 1 - 2016*q + 1457568*q^2 - 411997824*q^3 + 16227967392*q^4 + O(q^5)), (24, 1 - 288*q - 325728*q^2 + 11700864*q^3 + 35176468896*q^4 + O(q^5)), (24, 1 + 1440*q + 876960*q^2 + 292072320*q^3 + 57349833120*q^4 + O(q^5))]
        sage: dimension_modular_forms(SL2Z,24)
        3
    """
    verbose('multiplying forms up to weight %s'%weight)
    # Algorithm: run through the subsets of forms and for each check
    # whether or not the sum of the weights (with coefficients -- i.e.,
    # account for multiplicities) of the forms equals weight.
    # If so, multiply those together and append them to the output
    # list v

    # The answer list
    v = []
    n = len(forms)

    # List of weights
    from sage.combinat.integer_vector_weighted import WeightedIntegerVectors
    wts = WeightedIntegerVectors(weight, [f[0] for f in forms])
    
    for c in wts:
        if sum(c[i]*forms[i][0] for i in xrange(n) if c[i]) != weight:
            raise ArithmeticError, "Can't get here!"
        g = prod(forms[i][1]**c[i] for i in xrange(n))
        v.append((weight, g))
        if stop_dim and len(v) >= stop_dim:
            z = span_of_series([f for _, f in v]).dimension()                
            if z >= stop_dim:
                return v
    return v

 def _lmul_(self, right):
     """
     Scalar multiplication self*right.
     """
     #RH: mostly "copied" from dist.pyx but then changed
     ans = CoeffMod_OMS_element(None, self.parent(), None, False)
     #p = self.parent().prime()
     if right.is_zero():
         verbose("right is zero: %s"%(right), level=2)
         ans._moments = self.parent().approx_module(0)([])
         ans.ordp = min(self.parent().precision_cap(), right.valuation()+self.ordp)
     else:
         v, u = right.val_unit()
         ans._moments = self._moments * u
         ans.ordp = self.ordp + v
     return ans

def hnf_square(A, proof):
    """
    INPUT:
        a nonsingular n x n matrix A over the integers.
    OUTPUT:
        the Hermite normal form of A.

    EXAMPLES:
        sage: import sage.matrix.matrix_integer_dense_hnf as hnf
        sage: A = matrix(ZZ, 3, [-21, -7, 5, 1,20,-7, -1,1,-1])
        sage: hnf.hnf_square(A, False)
        [ 1  6 29]
        [ 0  7 28]
        [ 0  0 46]
        sage: A.echelon_form()
        [ 1  6 29]
        [ 0  7 28]
        [ 0  0 46]
    """
    n = A.nrows()
    m = A.ncols()
    if n != m:
        raise ValueError("A must be square.")

    # Small cases -- don't use this algorithm
    if n <= 3:
        return A.echelon_form(algorithm="pari")

    if A.rank() < A.nrows():
        raise ValueError("matrix must have full rank")
                           
    

    t = verbose("starting slicings")
    B = A.matrix_from_rows(range(m-2)).matrix_from_columns(range(n-1))
    c = A.matrix_from_rows([m-2]).matrix_from_columns (range(n-1))
    d = A.matrix_from_rows([m-1]).matrix_from_columns (range(n-1))
    b = A.matrix_from_columns([n-1]).matrix_from_rows(range(m-2))
    verbose("done slicing", t)

    try:
        (d1,d2) = double_det (B,c,d, proof=proof)
    except (ValueError, ZeroDivisionError), msg:
        d1 = B.stack(c).det(proof=proof)
        d2 = B.stack(d).det(proof=proof)