def HalfPeriodsInTermsOfLambdas(L1, L2, L3, lvec_and_Mlist = None, HP0 = None, prec = None, max_iters = 20):
    K = L1.parent()
    L0 = Matrix(K, 3, 1, [L1, L2, L3])
    if lvec_and_Mlist is None:
        assert prec is not None
        lvec, Mlist = load_lvec_and_Mlist(prec)
    else:
        lvec, Mlist = lvec_and_Mlist
    # Evaluates a matrix M with entries in Z[[x,y,z]] at points x0,y0,z0
    def ev(x0,y0,z0):
        return [lvec[0](x0,y0,z0), lvec[1](x0,y0,z0), ((1-z0)/(1+z0))**2 * lvec[2](x0,y0,z0)]
    if HP0 is None:
        HP0 = [K(0),K(0),K(0)]
    Pn = Matrix(K,3,1,HP0) # 0th approximation
    n_iters = 0
    while n_iters < 20:
        n_iters += 1
        Jinv = evaluate_twisted_jacobian_matrix(*Pn.list(),Mlist = Mlist).inverse()
        FPn = matrix(3,1, ev(*Pn.list()))
        Pnn = Pn - Jinv * (FPn - L0)
        print('(%s)'%n_iters, [(u-v).valuation() for u,v in zip(Pn.list(), Pnn.list())])
        if all([u == v for u,v in zip(Pn.list(), Pnn.list())]):
            return Pn
        Pn = Pnn
    raise RuntimeError,"Does not converge"

def p_adic_l_invariant_additive(logA,logB, alpha, beta, Tmatrix):
    K = logA.parent()
    logA, logB = K(logA), K(logB)
    x,y,z,t = Tmatrix.change_ring(K).list()
    M = Matrix(K,2,2,[alpha,x*alpha+z*beta,beta, y*alpha+t*beta])
    n  = Matrix(K,2,1,[logA, logB])
    return M.solve_right(n).list()

    def discriminant(self):
        """
        Return the discriminant of this ring, which is the discriminant of
        the trace pairing.

        .. note::

           One knows that for modular abelian varieties, the
           endomorphism ring should be isomorphic to an order in a
           number field. However, the discriminant returned by this
           function will be `2^n` ( `n =`
           self.dimension()) times the discriminant of that order,
           since the elements are represented as 2d x 2d
           matrices. Notice, for example, that the case of a one
           dimensional abelian variety, whose endomorphism ring must
           be ZZ, has discriminant 2, as in the example below.

        EXAMPLES::

            sage: J0(33).endomorphism_ring().discriminant()
            -64800
            sage: J0(46).endomorphism_ring().discriminant()  # long time (6s on sage.math, 2011)
            24200000000
            sage: J0(11).endomorphism_ring().discriminant()
            2
        """
        g = self.gens()
        M = Matrix(ZZ,len(g), [ (g[i]*g[j]).trace()
                                for i in range(len(g)) for j in range(len(g)) ])
        return M.determinant()

    def space(self):
        r'''
        Calculates the space of cocyles modulo coboundaries, as a Z-module.

        TESTS:

        sage: from darmonpoints.sarithgroup import *
        sage: from darmonpoints.cohomology_abstract import *
        sage: from darmonpoints.ocmodule import *
        sage: GS = BigArithGroup(5, 6,1,use_shapiro=False,outfile='/tmp/darmonpoints.tmp') #  optional - magma
        sage: G = GS.large_group() #  optional - magma
        sage: V = OCVn(5,1)     #  optional - magma
        sage: Coh = CohomologyGroup(G,V,trivial_action = False) #  optional - magma
        '''
        verb = get_verbose()
        set_verbose(0)
        V = self.coefficient_module()
        R = V.base_ring()
        Vdim = V.dimension()
        G = self.group()
        gens = G.gens()
        ambient = R**(Vdim * len(gens))
        # Now find the subspace of cocycles
        A = Matrix(R, Vdim * len(gens), 0)
        for r in G.get_relation_words():
            Alist = self.fox_gradient(r)
            newA = block_matrix(Alist, nrows = 1)
            A = A.augment(newA.transpose())
        A = A.transpose()
        cocycles = ambient.submodule([ambient(o) for o in A.right_kernel_matrix().rows()])
        gmat = block_matrix([self._gen_pows[i][1] - 1 for i in range(len(G.gens()))], nrows = len(G.gens()))
        coboundaries = cocycles.submodule([ambient(o) for o in gmat.columns()])
        ans = cocycles.quotient(coboundaries)
        set_verbose(verb)
        return ans

    def __init__(self, A, gens=None, given_by_matrix=False):
        """
        EXAMPLES::

            sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1, 0], [0, 1]]), Matrix([[0, 1], [0, 0]])])
            sage: I = A.ideal(A([0,1]))
            sage: TestSuite(I).run(skip="_test_category") # Currently ideals are not using the category framework
        """
        k = A.base_ring()
        n = A.degree()
        if given_by_matrix:
            self._basis_matrix = gens
            gens = gens.rows()
        elif gens is None:
            self._basis_matrix = Matrix(k, 0, n)
        elif isinstance(gens, (list, tuple)):
            B = [FiniteDimensionalAlgebraIdeal(A, x).basis_matrix() for x in gens]
            B = reduce(lambda x, y: x.stack(y), B, Matrix(k, 0, n))
            self._basis_matrix = B.echelon_form().image().basis_matrix()
        elif is_Matrix(gens):
            gens = FiniteDimensionalAlgebraElement(A, gens)
        elif isinstance(gens, FiniteDimensionalAlgebraElement):
            gens = gens.vector()
            B = Matrix([gens * b for b in A.table()])
            self._basis_matrix = B.echelon_form().image().basis_matrix()
        Ideal_generic.__init__(self, A, gens)

    def diagonal_matrix(self):
        r"""
        Returns a diagonal matrix `D` and a matrix `T` such that `T^t A T = D`
        holds, where `(x, y, z) A (x, y, z)^t` is the defining polynomial
        of the conic ``self``.

        EXAMPLES:

        ::

            sage: c = Conic(QQ, [1,2,3,4,5,6])
            sage: d, t = c.diagonal_matrix(); d, t
            (
            [    1     0     0]  [   1   -1 -7/6]
            [    0     3     0]  [   0    1 -1/3]
            [    0     0 41/12], [   0    0    1]
            )
            sage: t.transpose()*c.symmetric_matrix()*t
            [    1     0     0]
            [    0     3     0]
            [    0     0 41/12]

        Diagonal matrices are only defined in characteristic different
        from `2`:

        ::

            sage: c = Conic(GF(4, 'a'), [0, 1, 1, 1, 1, 1])
            sage: c.is_smooth()
            True
            sage: c.diagonal_matrix()
            Traceback (most recent call last):
            ...
            ValueError: The conic self (= Projective Conic Curve over Finite Field in a of size 2^2 defined by x*y + y^2 + x*z + y*z + z^2) has no symmetric matrix because the base field has characteristic 2
        """
        A = self.symmetric_matrix()
        B = self.base_ring()
        basis = [vector(B,{2:0,i:1}) for i in range(3)]
        for i in range(3):
            zerovalue = (basis[i]*A*basis[i].column()== 0)
            if zerovalue:
                for j in range(i+1,3):
                    if basis[j]*A*basis[j].column() != 0:
                        b = basis[i]
                        basis[i] = basis[j]
                        basis[j] = b
                        zerovalue = False
            if zerovalue:
                for j in range(i+1,3):
                    if basis[i]*A*basis[j].column() != 0:
                        basis[i] = basis[i]+basis[j]
                        zerovalue = False
            if not zerovalue:
                l = (basis[i]*A*basis[i].column())
                for j in range(i+1,3):
                    basis[j] = basis[j] - \
                               (basis[i]*A*basis[j].column())/l * basis[i]
        T = Matrix(basis).transpose()
        return T.transpose()*A*T, T

 def LocalizedMomentMatrix(self, idx):
     """
     Returns localized moment matrix corresponding to $g_{idx}$
     """
     
     tmp_vec = self.MonomialsVec(self.Relaxation-self.HalfDegs[idx])
     m = Matrix(1, len(tmp_vec), tmp_vec)
     return self.Constraints[idx]* (m.transpose() * m)

    def linear_approximation_matrix(self):
        """
        Return linear approximation matrix ``A`` for this S-box.

        Let ``i_b`` be the ``b``-th bit of ``i`` and ``o_b`` the
        ``b``-th bit of ``o``. Then ``v = A[i,o]`` encodes the bias of
        the equation ``sum( i_b * x_i ) = sum( o_b * y_i )`` if
        ``x_i`` and ``y_i`` represent the input and output variables
        of the S-box.

        See [He2002]_ for an introduction to linear cryptanalysis.

        EXAMPLES::

            sage: from sage.crypto.sbox import SBox
            sage: S = SBox(7,6,0,4,2,5,1,3)
            sage: S.linear_approximation_matrix()
            [ 4  0  0  0  0  0  0  0]
            [ 0  0  0  0  2  2  2 -2]
            [ 0  0 -2 -2 -2  2  0  0]
            [ 0  0 -2  2  0  0 -2 -2]
            [ 0  2  0  2 -2  0  2  0]
            [ 0 -2  0  2  0  2  0  2]
            [ 0 -2 -2  0  0 -2  2  0]
            [ 0 -2  2  0 -2  0  0 -2]

        According to this matrix the first bit of the input is equal
        to the third bit of the output 6 out of 8 times::

            sage: for i in srange(8): print(S.to_bits(i)[0] == S.to_bits(S(i))[2])
            False
            True
            True
            True
            False
            True
            True
            True
        """
        m = self.m
        n = self.n

        nrows = 1<<m
        ncols = 1<<n

        B = BooleanFunction(self.m)
        L = []
        for j in range(ncols):
            for i in range(nrows):
                B[i] = ZZ(self(i)&j).popcount()
            L.append(B.walsh_hadamard_transform())

        A = Matrix(ZZ, ncols, nrows, L)
        A = -A.transpose()/2
        A.set_immutable()

        return A

    def __init__(self, A, elt=None, check=True):
        """
        TESTS::

            sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
            sage: A(QQ(4))
            Traceback (most recent call last):
            ...
            TypeError: elt should be a vector, a matrix, or an element of the base field

            sage: B = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [-1,0]])])
            sage: elt = B(Matrix([[1,1], [-1,1]])); elt
            e0 + e1
            sage: TestSuite(elt).run()
            sage: B(Matrix([[0,1], [1,0]]))
            Traceback (most recent call last):
            ...
            ValueError: matrix does not define an element of the algebra
        """
        AlgebraElement.__init__(self, A)
        k = A.base_ring()
        n = A.degree()
        if elt is None:
            self._vector = vector(k, n)
            self._matrix = Matrix(k, n)
        else:
            if isinstance(elt, int):
                elt = Integer(elt)
            elif isinstance(elt, list):
                elt = vector(elt)
            if A == elt.parent():
                self._vector = elt._vector.base_extend(k)
                self._matrix = elt._matrix.base_extend(k)
            elif k.has_coerce_map_from(elt.parent()):
                e = k(elt)
                if e == 0:
                    self._vector = vector(k, n)
                    self._matrix = Matrix(k, n)
                elif A.is_unitary():
                    self._vector = A._one * e
                    self._matrix = Matrix.identity(k, n) * e
                else:
                    raise TypeError("algebra is not unitary")
            elif is_Vector(elt):
                self._vector = elt.base_extend(k)
                self._matrix = Matrix(k, sum([elt[i] * A.table()[i] for i in xrange(n)]))
            elif is_Matrix(elt):
                if not A.is_unitary():
                    raise TypeError("algebra is not unitary")
                self._vector = A._one * elt
                if not check or sum([self._vector[i]*A.table()[i] for i in xrange(n)]) == elt:
                    self._matrix = elt
                else:
                    raise ValueError("matrix does not define an element of the algebra")
            else:
                raise TypeError("elt should be a vector, a matrix, " +
                                "or an element of the base field")

 def apply_matrix(self,m,in_place=True, mapping=False):
     r"""
     Carry out the GL(2,R) action of m on this surface and return the result.
     
     If in_place=True, then this is done in place and changes the surface. 
     This can only be carried out if the surface is finite and mutable.
     
     If mapping=True, then we return a GL2RMapping between this surface and its image. 
     In this case in_place must be False.
     
     If in_place=False, then a copy is made before the deformation.
     """
     if mapping==True:
         assert in_place==False, "Can not modify in place and return a mapping."
         return GL2RMapping(self, m)
     if not in_place:
         if self.is_finite():
             from sage.structure.element import get_coercion_model
             cm = get_coercion_model()
             field = cm.common_parent(self.base_ring(), m.base_ring())
             s=self.copy(mutable=True, new_field=field)
             return s.apply_matrix(m)
         else:
             return m*self
     else:
         # Make sure m is in the right state
         from sage.matrix.constructor import Matrix
         m=Matrix(self.base_ring(), 2, 2, m)
         assert m.det()!=self.base_ring().zero(), "Can not deform by degenerate matrix."
         assert self.is_finite(), "In place GL(2,R) action only works for finite surfaces."
         us=self.underlying_surface()
         assert us.is_mutable(), "In place changes only work for mutable surfaces."
         for label in self.label_iterator():
             us.change_polygon(label,m*self.polygon(label))
         if m.det()<self.base_ring().zero():
             # Polygons were all reversed orientation. Need to redo gluings.
             
             # First pass record new gluings in a dictionary.
             new_glue={}
             seen_labels=set()
             for p1 in self.label_iterator():
                 n1=self.polygon(p1).num_edges()
                 for e1 in xrange(n1):
                     p2,e2=self.opposite_edge(p1,e1)
                     n2=self.polygon(p2).num_edges()
                     if p2 in seen_labels:
                         pass
                     elif p1==p2 and e1>e2:
                         pass
                     else:
                         new_glue[(p1, n1-1-e1)]=(p2, n2-1-e2)
                 seen_labels.add(p1)
             # Second pass: reassign gluings
             for (p1,e1),(p2,e2) in new_glue.iteritems():
                 us.change_edge_gluing(p1,e1,p2,e2)
         return self

 def matrix_rep(self,B=None):
     r"""
     Returns a matrix representation of ``self``.
     """
     #Express the element in terms of the basis B
     if B is None:
         B = self._parent.basis()
     A=Matrix(self._parent._R,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)])
     tmp=A.solve_right(self._val)
     return tmp

 def decompose(self):
     """
         Gives an SOS decomposition of f if exists as a list of polynomials
         of at most half degree of f.
         This method also fills the 'Info' as 'minimize' does. In addition,
         returns Info['is sos'] which is Boolean depending on the status of 
         sdp solver.
     """
     n = self.NumVars
     N0 = self.NumMonomials(n, self.MainPolyHlfDeg)
     N1 = self.NumMonomials(n, self.MainPolyTotDeg)
     self.MatSize = [N0, N1]
     
     vec = self.MonomialsVec(self.MainPolyHlfDeg)
     m = Matrix(1, N0, vec)
     Mmnt = m.transpose() * m
     Blck = [[] for i in range(N1)]
     C = []
     
     h = Matrix(self.Field, N0, N0, 0)
     C.append(h)
     decomp = []
     
     for i in range(N1):
         p = self.Monomials[i]
         A = self.Calpha(p, Mmnt)
         Blck[i].append(A)
     
     from SDP import SDP
     sos_sdp = SDP.sdp(solver = self.solver, Settings = {'detail':self.detail})
     sos_sdp.solve(C, self.PolyCoefFullVec(), Blck)
     
     if sos_sdp.Info['Status'] == 'Optimal':
         self.Info['status'] = 'Feasible'
         GramMtx = Matrix(sos_sdp.Info['X'][0])
         try:
             self.Info['Message'] = "A SOS decomposition of the polynomial were found."
             self.Info['is sos'] = True
             H1 = GramMtx.cholesky();
             tmpM = Matrix(1, N0, vec)
             decomp = list(tmpM*H1)[0]
             self.Info['Wall'] = sos_sdp.Info['Wall']
             self.Info['CPU'] = sos_sdp.Info['CPU']
         except:
             self.Info['Message'] = "The given polynomial seems to be a sum of squares, but no SOS decomposition were extracted."
             self.Info['is sos'] = False
             self.Info['Wall'] = sos_sdp.Info['Wall']
             self.Info['CPU'] = sos_sdp.Info['CPU']
     else:
         self.Info['Message'] = "The given polynomial is not a sum of squares."
         self.Info['status'] = 'Infeasible'
         self.Info['is sos']= False
     
     self.Info["Size"] = self.MatSize
     return decomp

    def _get_powers_and_mult(self,a,b,c,d,lambd,vect):
        r"""
        Compute the action of a matrix on the basis elements.

        EXAMPLES:

        ::

        """
        try:
            xnew = self._powers[(a,b,c,d)]
        except KeyError:
            R=self._PowerSeries
            r=R([b,a])
            s=R([d,c])
            n=self._n

            if self._depth == n+1:
                rpows=[R(1)]
                spows=[R(1)]
                for ii in range(n):
                    rpows.append(r*rpows[ii])
                    spows.append(s*spows[ii])
                x=Matrix(self._Rmod,n+1,n+1,0)
                for ii in range(n+1):
                    y=rpows[ii]*spows[n-ii]
                    for jj in range(self._depth):
                        x[ii,jj]=y[jj]
            else:
                ratio = r*(s**(-1))
                y = s**n
                x = Matrix(self._Rmod,self._depth,self._depth,0)
                for jj in range(self._depth):
                    x[0,jj] = y[jj]
                for ii in range(1,self._depth):
                    y *= ratio
                    for jj in range(self._depth):
                        x[ii,jj] = y[jj]

            xnew = x.change_ring(self._R) #ZZ)
            # if self._Rmod is self._R:
            #     xnew = x
            # else:
            #     #xnew = x.change_ring(self._R)
            #     xnew = x.change_ring(ZZ)

            self._powers[(a,b,c,d)] = xnew

        tmp = xnew*vect
        return self._R(lambd)*tmp #.change_ring(self._R)

    class _FiniteBasisConverter:
        def __init__(self, P, comb_mod, basis):
            r"""
            Basis should be a finite set of polynomials
            """
            self._poly_ring = P
            self._module = comb_mod
            self._basis = basis

            max_deg = max([self._poly_ring(b).degree() for b in self._basis])
            monoms = []
            for b in self._basis:
                poly = self._poly_ring(b)
                monoms += poly.monomials()
            monoms_list = tuple(Set(monoms))

            # check if the basis represented in terms of Monomials is efficient
            degs = [self._poly_ring(m).degree() for m in monoms]
            min_deg, max_deg = min(degs), max(degs)
            monoms_obj = Monomials(self._poly_ring, (min_deg, max_deg + 1))

            if monoms_obj.cardinality() < 2 * len(monoms_list):
                computational_basis = monoms_obj
            else:
                computational_basis = monoms_list
            self._monomial_module = PolynomialFreeModule(
                P=self._poly_ring, basis=computational_basis)
            cols = [self._monomial_module(b).to_vector() for b in self._basis]
            self._basis_mat = Matrix(cols).transpose()
            if self._basis_mat.ncols() > self._basis_mat.rank():
                raise ValueError(
                    "Basis polynomials are not linearly independent")

        def convert(self, p):
            r"""
            Algorithm is to convert all polynomials into monomials and use
            linear algebra to solve for the appropriate coefficients in this
            common basis.
            """
            try:
                p_vect = self._monomial_module(p).to_vector()
                decomp = self._basis_mat.solve_right(p_vect)
            except ValueError:
                raise ValueError(
                    "Value %s is not spanned by the basis polynomials" % p)
            polys = [v[1] * self._module.monomial(v[0])
                     for v in zip(self._basis, decomp)]
            module_p = sum(polys, self._module.zero())
            return module_p

def poly_dual_basis(P, poly_basis):
    r"""
    Return a collection of polynomials which are dual under the differential
    bilinear form to a given homogeneous collection

    INPUT:

    - ``P`` -- a polynomial ring
    - ``poly_basis`` -- a collection of polynomials in ``P`` which are
      homogeneous and linearly independent

    OUTPUT:

    - the dual basis of the polynomials in ``poly_basis`` in their span

    EXAMPLES:

        sage: P.<x, y> = PolynomialRing(QQ)
        sage: poly_basis = (1, x, x+y)
        sage: poly_dual_basis(P, poly_basis)
        [1, x - y, y]
        sage: poly_basis = (1, 2*x - y, x^2, x^2 + x*y)
        sage: poly_dual_basis(P, poly_basis)
        [1, 2/5*x - 1/5*y, 1/2*x^2 - x*y, x*y]
    """
    # recast poly_basis to ensure elements are all from P
    poly_basis = [P(p) for p in poly_basis]
    # compute max degree of basis polynomials for linear algebra computations
    deg = max([p.degree() for p in poly_basis])
    # construct polynomial free module for linear algebra computations
    monoms = Monomials(P, (0, deg))
    poly_module = PolynomialFreeModule(P, basis=monoms)
    # compute the values of the bilinear form <m|m> for basis monomials m
    bilinear_form_coeffs = []
    for b in poly_module.basis().keys():
        # each b is a monomial in P of degree at most deg
        b = P(b)
        bilinear_form_coeffs.append(prod(map(factorial, b.degrees())))
    # compute dual basis
    A = Matrix([poly_module(p).to_vector() for p in poly_basis])
    D = Matrix.diagonal(bilinear_form_coeffs, sparse=False)
    B = (A * D * A.transpose()).inverse()
    # reconstruct dual basis polynomials from corresponding vectors
    dual_basis = []
    for col in B.columns():
        q = sum([coeff * p for coeff, p in zip(col, poly_basis)])
        dual_basis.append(q)
    return dual_basis

    def __init__(self,segment,rho,rhoexp):
        """
        Initializes self.

        See ``AssociatedFactor`` for full documentation.

        """
        self.segment = segment
        self.rho = rho
        self.rhoexp = rhoexp
        self.Fplus = self.rho.degree()

        if self.segment.frame.is_first():
            # In the first frame, so FFbase is the residue class field of O
            self.FFbase = self.segment.frame.R
        else:
            # Not the first frame
            self.FFbase = self.segment.frame.prev.FF

        if self.Fplus == 1:
            self.FF = self.FFbase
            self.FFz = PolynomialRing(self.FF,'z'+str(self.segment.frame.depth))
            # rho is linear delta is the root of rho
            self.delta = self.rho.roots()[0][0]
        else:
            self.FF = GF(self.FFbase.order()**self.Fplus,'a'+str(self.segment.frame.depth))
            self.FFz = PolynomialRing(self.FF,'z'+str(self.segment.frame.depth))
            self.FFbase_gamma = (self.FFz(self.FFbase.modulus())).roots()[0][0]
            FFrho = self.FFz([self.FFbase_elt_to_FF(a) for a in list(rho)])
            self.gamma = FFrho.roots()[0][0]
            basis = [(self.gamma**j*self.FFbase_gamma**i).polynomial() for j in range(0,self.Fplus) for i in range(0,self.FFbase.degree())]
            self.basis_trans_mat = Matrix([self.FF(b)._vector_() for b in basis])

    def matrix_rep(self,B=None):
        r"""

        EXAMPLES:

        This example illustrates ...

        ::

        """
        #Express the element in terms of the basis B
        if(B is None):
            B=self._parent.basis()
        A=Matrix(self._parent._R,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)])
        tmp=A.solve_right(self._val)
        return tmp

 def acting_matrix(self, g, prec = None):
     dim = len(self.basis())
     ans = Matrix(self._V.base_ring(),dim, 0)
     for v in self.basis():
         gv = g * v
         gvlist = []
         for w in gv._val:
             try:
                 wlist = list(w)
             except TypeError:
                 wlist = list(w._moments)
             if prec is None:
                 gvlist.extend(wlist)
             else:
                 gvlist.extend(wlist[:prec])
         ans = ans.augment(Matrix(self._V.base_ring(),dim,1,gvlist))
     return ans

 def linear_relation(self, List, Psi, verbose=True, prec=None):
     for Phi in List:
         assert Phi.valuation() >= 0, "Symbols must be integral"
     assert Psi.valuation() >= 0
     R = self.base()
     Rbase = R.base()
     w = R.gen()
     d = len(List)
     if d == 0:
         if Psi.is_zero():
             return [None, R(1)]
         else:
             return [None, 0]
     if prec is None:
         M, var_prec = Psi.precision_absolute()
     else:
         M, var_prec = prec
     p = self.prime()
     V = R**d
     RSR = LaurentSeriesRing(Rbase, R.variable_name())
     VSR = RSR**self.source().ngens()
     List_TMs = [VSR(Phi.list_of_total_measures()) for Phi in List]
     Psi_TMs = VSR(Psi.list_of_total_measures())
     A = Matrix(RSR, List_TMs).transpose()
     try:
         sol = V([vv.power_series() for vv in A.solve_right(Psi_TMs)])
     except ValueError:
         #try "least squares"
         if verbose:
             print "Trying least squares."
         sol = (A.transpose() * A).solve_right(A.transpose() * Psi_TMs)
         #check precision (could make this better, checking each power of w)
         p_prec = M
         diff = Psi_TMs - sum([sol[i] * List_TMs[i] for i in range(len(List_TMs))])
         for i in diff:
             for j in i.list():
                 if p_prec > j.valuation():
                     p_prec = j.valuation()
         if verbose:
             print "p-adic precision is now", p_prec
         #Is this right?
         sol = V([R([j.add_bigoh(p_prec) for j in i.list()]) for i in sol])
     return [sol, R(-1)]

def solve_coefficient_system(Q, equations, vars):
    # NOTE: to make things easier (and uniform) in the univariate case a dummy
    # variable is added to the polynomial ring. See compute_bd()
    a = Q.gens()[:-1]
    B = Q.base_ring()

    # construct the coefficient system and right-hand side
    system = [[e.coefficient({ai:1}) for ai in a] for e in equations]
    rhs = [-e.constant_coefficient() for e in equations]
    system = Matrix(B, system)
    rhs = Matrix(B, rhs).transpose()

    # we only allow unique solutions. return None if there are infinitely many
    # solutions or if no solution exists. Sage will raise a ValueError in both
    # circumstances
    try:
        sol = system.solve_right(rhs)
    except ValueError:
        return None
    return sol