def iter_positive_forms_with_content(self) :
        if self.__disc is infinity :
            raise ValueError, "infinity is not a true filter index"
        
        
        if self.__reduced :        
            for a in xrange(1,isqrt(self.__disc // 3) + 1) :
                for b in xrange(a+1) :
                    g = gcd(a, b)
                    for c in xrange(a, (b**2 + (self.__disc - 1))//(4*a) + 1) :
                        yield (a,b,c), gcd(g,c)
        else :
            maxtrace = floor(5*self.__disc / 15 + sqrt(self.__disc)/2)
            for a in xrange(1, maxtrace + 1) :
                for c in xrange(1, maxtrace - a + 1) :
                    g = gcd(a,c)
                    
                    Bu = isqrt(4*a*c - 1)

                    di = 4*a*c - self.__disc
                    if di >= 0 :
                        Bl = isqrt(di) + 1 
                    else :
                        Bl = 0
                    
                    for b in xrange(-Bu, -Bl + 1) :
                        yield (a,b,c), gcd(g,b)
                    for b in xrange(Bl, Bu + 1) :
                        yield (a,b,c), gcd(g,b)
        #! if self.__reduced

        raise StopIteration

    def _find_cusps(self):
        r"""
        Calculate the reduced representatives of the equivalence classes of
        cusps for this group. Adapted from code by Ron Evans.

        EXAMPLES::

            sage: Gamma(8).cusps() # indirect doctest
            [0, 1/4, 1/3, 3/8, 1/2, 2/3, 3/4, 1, 4/3, 3/2, 5/3, 2, 7/3, 5/2, 8/3, 3, 7/2, 11/3, 4, 14/3, 5, 6, 7, Infinity]
        """
        n = self.level()
        C = [QQ(x) for x in range(n)]

        n0=n//2
        n1=(n+1)//2

        for r in range(1, n1):
            if r > 1 and gcd(r,n)==1:
                C.append(ZZ(r)/ZZ(n))
            if n0==n/2 and gcd(r,n0)==1:
                C.append(ZZ(r)/ZZ(n0))

        for s in range(2,n1):
            for r in range(1, 1+n):
                if GCD_list([s,r,n])==1:
                    # GCD_list is ~40x faster than gcd, since gcd wastes loads
                    # of time initialising a Sequence type.
                    u,v = _lift_pair(r,s,n)
                    C.append(ZZ(u)/ZZ(v))

        return [Cusp(x) for x in sorted(C)] + [Cusp(1,0)]

    def _iter_positive_forms_with_content_and_discriminant(self) :
        if self.__disc is infinity :
            raise ValueError, "infinity is not a true filter index"
        
        if self.__reduced :
            for (l, (u,x)) in enumerate(self.__p1list) :
                if u == 0 :
                    for a in xrange(self.__level, isqrt(self.__disc // 4) + 1, self.__level) :
                        frpa = 4 * a
                        for b in xrange(a+1) :
                            g = gcd(a // self.__level,b)
                            bsq = b**2

                            for c in xrange(a, (b**2 + (self.__disc - 1))//(4*a) + 1) :
                                yield (((a,b,c), l), gcd(g,c), frpa*c - bsq)
                else :
                    for a in xrange(1, isqrt(self.__disc // 3) + 1) :
                        frpa = 4 * a
                        for b in xrange(a+1) :
                            g = gcd(a, b)
                            bsq = b**2
                                                        
                            ## We need x**2 * a + x * b + c % N == 0
                            h = (-((x**2 + 1) * a + x * b)) % self.__level
                            for c in xrange( a + h,
                                             (b**2 + (self.__disc - 1))//(4*a) + 1, self.__level ) :
                                yield (((a,b,c), l), gcd(g,(x**2 * a + x * b + c) // self.__level), frpa*c - bsq)
        #! if self.__reduced
        else :
            raise NotImplementedError
        
        raise StopIteration

def normalize_coefficients(self, c):
    r"""
    If our coefficient ring is the field of fractions over a univariate
    polynomial ring over the rationals, then we should clear both the
    numerator and denominator of the denominators of their
    coefficients.

    INPUT:

    - ``self`` -- a Jack basis of the symmetric functions
    - ``c`` -- a coefficient in the base ring of ``self``

    OUTPUT:

    - divide numerator and denominator by the greatest common divisor

    EXAMPLES::

        sage: JP = SymmetricFunctions(FractionField(QQ['t'])).jack().P()
        sage: t = JP.base_ring().gen()
        sage: a = 2/(1/2*t+1/2)
        sage: JP._normalize_coefficients(a)
        4/(t + 1)
        sage: a = 1/(1/3+1/6*t)
        sage: JP._normalize_coefficients(a)
        6/(t + 2)
        sage: a = 24/(4*t^2 + 12*t + 8)
        sage: JP._normalize_coefficients(a)
        6/(t^2 + 3*t + 2)
    """
    BR = self.base_ring()
    if is_FractionField(BR) and BR.base_ring() == QQ:
        denom = c.denominator()
        numer = c.numerator()

        # Clear the denominators
        a = lcm([i.denominator() for i in denom.coefficients(sparse=False)])
        b = lcm([i.denominator() for i in numer.coefficients(sparse=False)])
        l = Integer(a).lcm(Integer(b))
        denom *= l
        numer *= l

        # Divide through by the gcd of the numerators
        a = gcd([i.numerator() for i in denom.coefficients(sparse=False)])
        b = gcd([i.numerator() for i in numer.coefficients(sparse=False)])
        l = Integer(a).gcd(Integer(b))

        denom = denom // l
        numer = numer // l

        return c.parent()(numer, denom)
    else:
        return c

    def _find_cusps(self):
        r"""
        Return an ordered list of inequivalent cusps for self, i.e. a
        set of representatives for the orbits of self on
        `\mathbf{P}^1(\QQ)`.  These are returned in a reduced
        form; see self.reduce_cusp for the definition of reduced.

        ALGORITHM:
            Lemma 3.2 in Cremona's 1997 book shows that for the action
            of Gamma1(N) on "signed projective space"
            `\Q^2 / (\Q_{\geq 0}^+)`, we have `u_1/v_1 \sim u_2 / v_2`
            if and only if `v_1 = v_2 \bmod N` and `u_1 = u_2 \bmod
            gcd(v_1, N)`. It follows that every orbit has a
            representative `u/v` with `v \le N` and `0 \le u \le
            gcd(v, N)`.  We iterate through all pairs `(u,v)`
            satisfying this.

            Having found a set containing at least one of every
            equivalence class modulo Gamma1(N), we can be sure of
            picking up every class modulo GammaH(N) since this
            contains Gamma1(N); and the reduce_cusp call does the
            checking to make sure we don't get any duplicates.

        EXAMPLES::

            sage: Gamma1(5)._find_cusps()
            [0, 2/5, 1/2, Infinity]
            sage: Gamma1(35)._find_cusps()
            [0, 2/35, 1/17, 1/16, 1/15, 1/14, 1/13, 1/12, 3/35, 1/11, 1/10, 1/9, 4/35, 1/8, 2/15, 1/7, 1/6, 6/35, 1/5, 3/14, 8/35, 1/4, 9/35, 4/15, 2/7, 3/10, 11/35, 1/3, 12/35, 5/14, 13/35, 2/5, 3/7, 16/35, 17/35, 1/2, 8/15, 4/7, 3/5, 9/14, 7/10, 5/7, 11/14, 4/5, 6/7, 9/10, 13/14, Infinity]
            sage: Gamma1(24)._find_cusps() == Gamma1(24).cusps(algorithm='modsym')
            True
            sage: GammaH(24, [13,17])._find_cusps() == GammaH(24,[13,17]).cusps(algorithm='modsym')
            True
        """

        s = []
        hashes = []
        N = self.level()

        for d in range(1, 1+N):
            w = N.gcd(d)
            M = int(w) if w > 1 else 2
            for a in range(1,M):
                if gcd(a, w) != 1:
                    continue
                while gcd(a, d) != 1:
                    a += w
                c = self.reduce_cusp(Cusp(a,d))
                h = hash(c)
                if not h in hashes:
                    hashes.append(h)
                    s.append(c)
        return sorted(s)

def DuadicCodeOddPair(F,S1,S2):
    """
    Constructs the "odd pair" of duadic codes associated to the
    "splitting" S1, S2 of n.

    .. warning::

       Maybe the splitting should be associated to a sum of
       q-cyclotomic cosets mod n, where q is a *prime*.

    EXAMPLES::

        sage: from sage.coding.code_constructions import _is_a_splitting
        sage: n = 11; q = 3
        sage: C = Zmod(n).cyclotomic_cosets(q); C
        [[0], [1, 3, 4, 5, 9], [2, 6, 7, 8, 10]]
        sage: S1 = C[1]
        sage: S2 = C[2]
        sage: _is_a_splitting(S1,S2,11)
        True
        sage: codes.DuadicCodeOddPair(GF(q),S1,S2)
        ([11, 6] Cyclic Code over GF(3),
         [11, 6] Cyclic Code over GF(3))

    This is consistent with Theorem 6.1.3 in [HP2003]_.
    """
    from .cyclic_code import CyclicCode
    n = len(S1) + len(S2) + 1
    if not _is_a_splitting(S1,S2,n):
        raise TypeError("%s, %s must be a splitting of %s."%(S1,S2,n))
    q = F.order()
    k = Mod(q,n).multiplicative_order()
    FF = GF(q**k,"z")
    z = FF.gen()
    zeta = z**((q**k-1)/n)
    P1 = PolynomialRing(FF,"x")
    x = P1.gen()
    g1 = prod([x-zeta**i for i in S1+[0]])
    g2 = prod([x-zeta**i for i in S2+[0]])
    j = sum([x**i/n for i in range(n)])
    P2 = PolynomialRing(F,"x")
    x = P2.gen()
    coeffs1 = [_lift2smallest_field(c)[0] for c in (g1+j).coefficients(sparse=False)]
    coeffs2 = [_lift2smallest_field(c)[0] for c in (g2+j).coefficients(sparse=False)]
    gg1 = P2(coeffs1)
    gg2 = P2(coeffs2)
    gg1 = gcd(gg1, x**n - 1)
    gg2 = gcd(gg2, x**n - 1)
    C1 = CyclicCode(length = n, generator_pol = gg1)
    C2 = CyclicCode(length = n, generator_pol = gg2)
    return C1,C2

    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 eval(self, expansion, weight = None) :
     if weight is None :
         try :
             weight = expansion.weight()
         except AttributeError :
             raise ValueError, "weight must be defined for the Hecke action"
     
     precision = expansion.precision()
     if precision.is_infinite() :
         precision = expansion._bounding_precision()
     else :
         precision = precision._hecke_operator(self.__l)
     characters = expansion.non_zero_components()
     expansion_level = precision.level()
     
     if gcd(expansion_level, self.__l) != 1 :
         raise ValueError, "Level of expansion and of Hecke operator must be coprime."
     
     hecke_expansion = dict()
     for ch in characters :
         hecke_expansion[ch] = dict( (k, self.hecke_coeff(expansion, ch, k, weight, expansion_level))
                                     for k in precision )
     
     result = expansion.parent()._element_constructor_(hecke_expansion)
     result._set_precision(expansion.precision()._hecke_operator(self.__l))
     
     return result

    def get_representatives( self, t, N) :
        r"""
        A helper function used in hecke_coeff that computes the right
        coset representatives of $\Gamma^0(t) \cap \Gamma_0(N) \ Gamma_0(N)$ where
        $\Gamma^0(t)$ is the subgroup of $SL(2,Z)$ where the upper right hand
        corner is divisible by $t$.

        NOTE
            We use the bijection $\Gamma^0(t)\SL(2,Z) \rightarrow P^1(\Z/t\Z)$
            given by $A \mapsto [1:0]A$.
        """
        if t == 1 : return [(1,0,0,1)]
        
        rep_list = []
        
        for (x,y) in P1List(t):
            ## we know that (N, x, y) = 1
            N1 = gcd(N, x)
            if N1 != 1 :
                x = x + t * N // N1

            ## we calculate a pair c,d satisfying a minimality condition
            ## to make later multiplications cheaper
            (_, d, c) = Integer(x)._xgcd(Integer(N * y), minimal=True)
            
            #print (x, y, -N * c, d)
            rep_list.append((x, y, -N * c, d))
                
        return rep_list

def find_gens_list(klist, r = 0, bound = None, verbose = True):
    """
    Use an elliptic curve variation of Rain's method to find a generator
    of a subfield of degree r within the ambient fields in klist.
    """
    p = klist[0].characteristic()
    ngcd = gcd([k.degree() for k in klist])
    nlcm = lcm([k.degree() for k in klist])
    if r == 0:
        r = ngcd
    assert all(k.degree() % r == 0 for k in klist)

    # This might be useful if elliptic curves defined over an 
    # extension of F_p are used.
    kb = k.base_ring()

    # List of candidates for l.
    lT = find_l(kb, r, bound)
    if lT is None:
        raise RuntimeError, "no suitable l found"

    # Find an elliptic curve with the given trace. 
    E = find_elliptic_curve(kb, lT)
    if E is None:
        raise RuntimeError, "no suitable elliptic curve found"

    return tuple(find_unique_orbit(E.change_ring(k), E.cardinality(extension_degree=k.degree()), lT[0], r) for k in klist)

    def order_at_cusp(self, cusp):
        r"""
        Return the order of vanishing of self at the given cusp.

        INPUT:


        -  ``cusp`` -  a CuspFamily object


        OUTPUT:


        - an integer


        EXAMPLES::

            sage: e = EtaProduct(2, {2:24, 1:-24})
            sage: e.order_at_cusp(CuspFamily(2, 1)) # cusp at infinity
            1
            sage: e.order_at_cusp(CuspFamily(2, 2)) # cusp 0
            -1
        """
        if not isinstance(cusp, CuspFamily):
            raise TypeError("Argument (=%s) should be a CuspFamily" % cusp)
        if cusp.level() != self.level():
            raise ValueError("Cusp not on right curve!")
        return 1/ZZ(24)/gcd(cusp.width(), self.level()//cusp.width()) * sum( [ell*self.r(ell)/cusp.width() * (gcd(cusp.width(), self.level()//ell))**2  for ell in self._keys] )

def _normalize_H(H, level):
    """
    Normalize representatives for a given subgroup H of the units
    modulo level.

    .. NOTE::

        This function does *not* make any attempt to find a minimal
        set of generators for H. It simply normalizes the inputs for use
        in hashing.

    EXAMPLES::

        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([23], 10)
        [3]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([1,5], 7)
        [5]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([4,18], 14)
        Traceback (most recent call last):
        ...
        ArithmeticError: The generators [4, 18] must be units modulo 14
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([3,17], 14)
        [3]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([-1,7,9], 10)
        [7, 9]
    """
    H = [ZZ(h) for h in H]
    for h in H:
        if gcd(h, level) > 1:
            raise ArithmeticError('The generators %s must be units modulo %s'%(H, level))
    H = sorted(set([h%level for h in H]))
    if 1 in H:
        H.remove(1)
    return H

    def __init__(self, base_field, length, designed_distance,
                 primitive_root=None, offset=1, jump_size=1, b=0):
        """
        TESTS:

        ``designed_distance`` must be between 1 and ``length`` (inclusive),
        otherwise an exception is raised::

            sage: C = codes.BCHCode(GF(2), 15, 16)
            Traceback (most recent call last):
            ...
            ValueError: designed_distance must belong to [1, n]
        """
        if not (0 < designed_distance <= length):
            raise ValueError("designed_distance must belong to [1, n]")

        if base_field not in Fields or not base_field.is_finite():
            raise ValueError("base_field has to be a finite field")

        q = base_field.cardinality()
        s = Zmod(length)(q).multiplicative_order()
        if gcd(jump_size, q ** s - 1) != 1:
            raise ValueError("jump_size must be coprime with the order of "
                             "the multiplicative group of the splitting field")

        D = [(offset + jump_size * i) % length
             for i in range(designed_distance - 1)]

        super(BCHCode, self).__init__(field=base_field, length=length,
                                      D=D, primitive_root=primitive_root)
        self._default_decoder_name = "UnderlyingGRS"
        self._jump_size = jump_size
        self._offset = offset
        self._designed_distance = designed_distance

def _lift_pair(U,V,N):
    r"""
    Utility function. Given integers ``U, V, N``, with `N \ge 1` and `{\rm
    gcd}(U, V, N) = 1`, return a pair `(u, v)` congruent to `(U, V) \bmod N`,
    such that `{\rm gcd}(u,v) = 1`, `u, v \ge 0`, `v` is as small as possible,
    and `u` is as small as possible for that `v`.

    *Warning*: As this function is for internal use, it does not do a
    preliminary sanity check on its input, for efficiency. It will recover
    reasonably gracefully if ``(U, V, N)`` are not coprime, but only after
    wasting quite a lot of cycles!

    EXAMPLES::

        sage: from sage.modular.arithgroup.congroup_gamma import _lift_pair
        sage: _lift_pair(2,4,7)
        (9, 4)
        sage: _lift_pair(2,4,8) # don't do this
        Traceback (most recent call last):
        ...
        ValueError: (U, V, N) must be coprime
    """
    u = U % N
    v = V % N
    if v == 0:
        if u == 1:
            return (1,0)
        else:
            v = N
    while gcd(u, v) > 1:
        u = u+N
        if u > N*v: raise ValueError("(U, V, N) must be coprime")
    return (u, v)

def get_overconvergent_class_matrices(p,E,prec, sign_at_infinity,use_ps_dists = False,use_sage_db = False,parallelize = False,progress_bar = False):
    # If the moments are pre-calculated, will load them. Otherwise, calculate and
    # save them to disk.
    if use_ps_dists == False:
        raise NotImplementedError, 'Must use distributions from Pollack-Stevens code in the split case'

    sgninfty = 'plus' if sign_at_infinity == 1 else 'minus'
    dist_type = 'ps' if use_ps_dists == True else 'fm'
    fname = 'moments_%s_%s_%s_%s_%s.sobj'%(p,E.cremona_label(),sgninfty,prec,dist_type)
    if use_sage_db:
        try:
            Phi = db(fname)
            return Phi
        except IOError: pass
    phi0 = E.pollack_stevens_modular_symbol()
    if sign_at_infinity == 1:
        phi0 = phi0.plus_part()
    elif sign_at_infinity == -1:
        phi0 = phi0.minus_part()
    else:
        assert sign_at_infinity == 0
        phi0 = phi0.plus_part() + phi0.minus_part()
    phi0 = 1 / gcd([val.moment(0) for val in phi0.values()]) * phi0
    Phi = phi0.lift(p,M = prec - 1,algorithm = 'stevens',eigensymbol = True)
    Phi._liftee = phi0
    return Phi

def num_cusps_of_width(N, d):
    r"""
    Return the number of cusps on `X_0(N)` of width d.

    INPUT:


    -  ``N`` - (integer): the level

    -  ``d`` - (integer): an integer dividing N, the cusp
       width


    EXAMPLES::

        sage: [num_cusps_of_width(18,d) for d in divisors(18)]
        [1, 1, 2, 2, 1, 1]
        sage: num_cusps_of_width(4,8)
        Traceback (most recent call last):
        ...
        ValueError: N and d must be positive integers with d|N
    """
    N = ZZ(N)
    d = ZZ(d)
    if N <= 0 or d <= 0 or (N % d) != 0:
        raise ValueError("N and d must be positive integers with d|N")

    return euler_phi(gcd(d, N//d))

                def rational_catalan_number(self, p, polynomial=False):
                    r"""
                    Return the ``p``-th rational Catalan number
                    associated to ``self``.

                    It is defined by

                    .. MATH::

                        \prod_{i = 1}^n \frac{p + (p(d_i-1)) \mod h)}{d_i},

                    where `d_1, \ldots, d_n` are the degrees and
                    `h` is the Coxeter number. See [STW2016]_
                    for this formula.

                    INPUT:

                    - ``polynomial`` -- optional boolean (default ``False``)
                      if ``True``, return instead the `q`-analogue as a
                      polynomial in `q`

                    REFERENCES:

                    .. [STW2016] C. Stump, H. Thomas, N. Williams.
                       *Cataland II*, in preparation, 2016.

                    EXAMPLES::

                        sage: W = ColoredPermutations(1,3)
                        sage: [W.rational_catalan_number(p) for p in [5,7,8]]
                        [7, 12, 15]

                        sage: W = ColoredPermutations(2,2)
                        sage: [W.rational_catalan_number(p) for p in [7,9,11]]
                        [10, 15, 21]

                    TESTS::

                        sage: W = ColoredPermutations(1,4)
                        sage: W.rational_catalan_number(3, polynomial=True)
                        q^6 + q^4 + q^3 + q^2 + 1
                    """
                    from sage.arith.all import gcd
                    from sage.combinat.q_analogues import q_int

                    h = self.coxeter_number()
                    if not gcd(h,p) == 1:
                        raise ValueError("parameter p = %s is not coprime to the Coxeter number %s" % (p, h))

                    if polynomial:
                        f = q_int
                    else:
                        f = lambda n: n

                    num = prod(f(p + (p * (deg - 1)) % h)
                               for deg in self.degrees())
                    den = prod(f(deg) for deg in self.degrees())
                    return num // den

def blift(LF, Li, p, S=None):
    r"""
    Search for a solution to the given list of inequalities.

    If found, lift the solution to
    an appropriate valuation. See Lemma 3.3.6 in [Molnar]_

    INPUT:

    - ``LF`` -- a list of integer polynomials in one variable (the normalized coefficients).

    - ``Li`` -- an integer, the bound on coefficients.

    - ``p`` -- a prime.

    OUTPUT:

    - Boolean -- whether or not the lift is successful.

    - integer -- the lift.

    EXAMPLES::

        sage: R.<b> = PolynomialRing(QQ)
        sage: from sage.schemes.projective.endPN_minimal_model import blift
        sage: blift([8*b^3 + 12*b^2 + 6*b + 1, 48*b^2 + 483*b + 117, 72*b + 1341, -24*b^2 + 411*b + 99, -144*b + 1233, -216*b], 2, 3)
        (True, 4)
    """

    P = LF[0].parent()
    #Determine which inequalities are trivial, and scale the rest, so that we only lift
    #as many times as needed.
    keepScaledIneqs = [scale(P(coeff),Li,p) for coeff in LF if coeff != 0]
    keptVals = [i[2] for i in keepScaledIneqs if i[0]]
    if keptVals != []:
        #Determine the valuation to lift until.
        liftval = max(keptVals)
    else:
        #All inequalities are satisfied.
        return True,1
    if S is None:
        S = PolynomialRing(Zmod(p),'b')
    keptScaledIneqs = [S(i[1]) for i in keepScaledIneqs if i[0]]
    #We need a solution for each polynomial on the left hand side of the inequalities,
    #so we need only find a solution for their gcd.
    g = gcd(keptScaledIneqs)
    rts = g.roots(multiplicities = False)
    for r in rts:
        #Recursively try to lift each root
        r_initial = QQ(r)
        newInput = P([r_initial, p])
        LG = [F(newInput) for F in LF]
        lift,lifted = blift(LG,Li,p,S=S)
        if lift:
            #Lift successful.
            return True,r_initial + p*lifted
    #Lift non successful.
    return False,0

    def primitive(self, signed=True):
        """
        Return hyperplane defined by primitive equation.

        INPUT:

        - ``signed`` -- boolean (optional, default: ``True``); whether
          to preserve the overall sign

        OUTPUT:

        Hyperplane whose linear expression has common factors and
        denominators cleared. That is, the same hyperplane (with the
        same sign) but defined by a rescaled equation. Note that
        different linear expressions must define different hyperplanes
        as comparison is used in caching.

        If ``signed``, the overall rescaling is by a positive constant
        only.

        EXAMPLES::

            sage: H.<x,y> = HyperplaneArrangements(QQ)
            sage: h = -1/3*x + 1/2*y - 1;  h
            Hyperplane -1/3*x + 1/2*y - 1
            sage: h.primitive()
            Hyperplane -2*x + 3*y - 6
            sage: h == h.primitive()
            False
            sage: (4*x + 8).primitive()
            Hyperplane x + 0*y + 2

            sage: (4*x - y - 8).primitive(signed=True)   # default
            Hyperplane 4*x - y - 8
            sage: (4*x - y - 8).primitive(signed=False)
            Hyperplane -4*x + y + 8
        """
        from sage.arith.all import lcm, gcd

        coeffs = self.coefficients()
        try:
            d = lcm([x.denom() for x in coeffs])
            n = gcd([x.numer() for x in coeffs])
        except AttributeError:
            return self
        if not signed:
            for x in coeffs:
                if x > 0:
                    break
                if x < 0:
                    d = -d
                    break
        parent = self.parent()
        d = parent.base_ring()(d)
        n = parent.base_ring()(n)
        if n == 0:
            n = parent.base_ring().one()
        return parent(self * d / n)

 def arith_prod_of_partitions(l1, l2):
     # Given two partitions `l_1` and `l_2`, we construct a new partition `l_1 \\boxtimes l_2` by
     # the following procedure: each pair of parts `a \\in l_1` and `b \\in l_2` contributes
     # `\\gcd (a, b)`` parts of size `\\lcm (a, b)` to `l_1 \\boxtimes l_2`. If `l_1` and `l_2`
     # are partitions of integers `n` and `m`, respectively, then `l_1 \\boxtimes l_2` is a
     # partition of `nm`.
     term_iterable = chain.from_iterable(repeat(lcm(pair), gcd(pair))
                                         for pair in product(l1, l2))
     return Partition(sorted(term_iterable, reverse=True))