def _prec_for_solve_diff_eqn_families(M, p):
    #UPDATE THIS with valuation of K[0]-1 and K[1]
    r"""
        A helper function for determining the (relative) precision of the input
        to solve_diff_eqn required in order obtain an answer with (relative)
        precision ``M``. The parameter ``p`` is the prime and ``k`` is the weight.

        Given input precision `M_\text{in}`, the output has precision

        .. MATH::

            M = M_\text{in} - \lceil\log_p(M_\text{in}) - 3.

    """
    # Do we need the weight?
    # A good guess to begin:
    if M < 1:
        raise ValueError("Desired precision M(=%s) must be at least 1."%(M))
    cp = (p - 2) / (p - 1)
    Min = ZZ(3 + M + ceil(ZZ(M).log(p)))
    # It looks like usually there are no iterations
    # For low M, there can be 1 or 2
    while M > Min*cp - ceil(log((Min * cp),p)) - 3: #THINK ABOUT THIS MORE
        Min += 1
        #print("An iteration in _prec_solve_diff_eqn")
    return Min

    def __init__(self, n, instance='key', m=None):
        """
        Construct LWE instance parameterised by security parameter ``n`` where
        all other parameters are chosen as in [CGW2013]_.

        INPUT:

        - ``n`` - security parameter (integer >= 89)
        - ``instance`` - one of

          - "key" - the LWE-instance that hides the secret key is generated
          - "encrypt" - the LWE-instance that hides the message is generated
            (default: ``key``)

        - ``m`` - number of allowed samples or ``None`` in which case ``m`` is
          chosen as in [CGW2013]_.  (default: ``None``)

        EXAMPLES::

            sage: from sage.crypto.lwe import UniformNoiseLWE
            sage: UniformNoiseLWE(89)
            LWE(89, 154262477, UniformSampler(0, 351), 'noise', 131)

            sage: UniformNoiseLWE(89, instance='encrypt')
            LWE(131, 154262477, UniformSampler(0, 497), 'noise', 181)
        """

        if n<89:
            raise TypeError("Parameter too small")

        n2 = n
        C  = 4/sqrt(2*pi)
        kk = floor((n2-2*log(n2, 2)**2)/5)
        n1 = floor((3*n2-5*kk)/2)
        ke = floor((n1-2*log(n1, 2)**2)/5)
        l  = floor((3*n1-5*ke)/2)-n2
        sk = ceil((C*(n1+n2))**(3/2))
        se = ceil((C*(n1+n2+l))**(3/2))
        q = next_prime(max(ceil((4*sk)**((n1+n2)/n1)), ceil((4*se)**((n1+n2+l)/(n2+l))), ceil(4*(n1+n2)*se*sk+4*se+1)))

        if kk<=0:
            raise TypeError("Parameter too small")

        if instance == 'key':
            D  = UniformSampler(0, sk-1)
            if m is None:
                m = n1
            LWE.__init__(self, n=n2, q=q, D=D, secret_dist='noise', m=m)
        elif instance == 'encrypt':
            D   = UniformSampler(0, se-1)
            if m is None:
                m = n2+l
            LWE.__init__(self, n=n1, q=q, D=D, secret_dist='noise', m=m)
        else:
            raise TypeError("Parameter instance=%s not understood."%(instance))

    def get_contents_bound_for_semi_definite_forms(self):
        r"""
        Return the bound for the contents of a semi positive definite
        form falling within this precision.

        EXAMPLES::

            sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision
            sage: prec = SiegelModularFormPrecision(7)
            sage: prec.get_contents_bound_for_semi_definite_forms()
            2
            sage: prec = SiegelModularFormPrecision(101)
            sage: prec.get_contents_bound_for_semi_definite_forms()
            26


        NOTE
            If (a,b,c) is semi positive definite, then it is GL(2,Z)-equivalent
            to a the form (0,0,g) where g is the contents (gcd) of a,b,c.
   
            I don't know what this does.  It seems to only do it for singular forms-- NR
        """
        if self.__type == 'infinity':
            return infinity
        elif self.__type == 'disc':
            return ceil((self.__prec+1)/4)
        elif self.__type == 'box':
            return self.__prec[2]
        else:
            raise RuntimeError, "Unexpected value of self.__type"        

    def _getitem_by_num(self, i):
        from sage.functions.other import ceil

        val = self._polygon(i)
        if val is Infinity:
            return self._base_exactprec
        else:
            return self._baseprec(ceil(val))

def griesmer_upper_bound(n,q,d,algorithm=None):
    r"""
    Returns the Griesmer upper bound.

    Returns the Griesmer upper bound for the number of elements in a
    largest linear code of minimum distance `d` in `\GF{q}^n`, cf. [HP2003]_.
    If the method is "gap", it wraps GAP's ``UpperBoundGriesmer``. 

    The bound states:

    .. MATH::

        `n\geq \sum_{i=0}^{k-1} \lceil d/q^i \rceil.`


    EXAMPLES:

    The bound is reached for the ternary Golay codes::

        sage: codes.bounds.griesmer_upper_bound(12,3,6)
        729
        sage: codes.bounds.griesmer_upper_bound(11,3,5)
        729

    ::

        sage: codes.bounds.griesmer_upper_bound(10,2,3)
        128
        sage: codes.bounds.griesmer_upper_bound(10,2,3,algorithm="gap")  # optional - gap_packages (Guava package)
        128

    TESTS::

        sage: codes.bounds.griesmer_upper_bound(11,3,6)
        243
        sage: codes.bounds.griesmer_upper_bound(11,3,6)
        243
    """
    _check_n_q_d(n, q, d)
    if algorithm=="gap":
        gap.load_package("guava")
        ans=gap.eval("UpperBoundGriesmer(%s,%s,%s)"%(n,d,q))
        return QQ(ans)
    else:
        #To compute the bound, we keep summing up the terms on the RHS
        #until we start violating the inequality.
        from sage.functions.other import ceil
        den = 1
        s = 0
        k = 0
        while s <= n:
            s += ceil(d/den)
            den *= q
            k = k + 1
        return q**(k-1)

    def _an_element_3d(self, x=0, y=0):
        r"""
        Returns an element in self.

        EXAMPLES::

            sage: from slabbe import DiscreteHyperplane
            sage: p = DiscreteHyperplane([1,pi,7], 1+pi+7, mu=10)
            sage: p._an_element_3d()
            (0, 0, 0)
        """
        a,b,c = self._v
        x_sqrt3 = ceil(x / sqrt(3))
        left  = ((a+b) * y + (a-b) * x_sqrt3 - self._mu) / (a+b+c)
        right = ((a+b) * y + (a-b) * x_sqrt3 - self._mu + self._omega) / (a+b+c)
        #print("left, right = ", left, right)
        #print("left, right = ", ceil(left), ceil(right)-1)
        # left <= z <= right
        znew = ceil(left)
        xnew = znew - y - x_sqrt3
        ynew = znew - y + x_sqrt3
        znew = ceil(right)-1
        #print(xnew, ynew, znew)
        #print(vector((xnew, ynew, znew)) in self)
        #print(vector((x,y,ceil(right)-1)) in self)
        v = vector((xnew, ynew, znew))
        if v in self:
            v.set_immutable()
            return v
        else:
            print("%s not in the plane" % v)
            print("trying similar points")
            v = vector((xnew, ynew, znew-1))
            if v in self:
                v.set_immutable()
                return v
            v = vector((xnew, ynew, znew+1))
            if v in self:
                v.set_immutable()
                return v
        raise ValueError("%s not in the plane" % v)

def _sympysage_ceiling(self):
    """
    EXAMPLES::

        sage: from sympy import Symbol, ceiling
        sage: assert ceil(x)._sympy_() == ceiling(Symbol('x'))
        sage: assert ceil(x) == ceiling(Symbol('x'))._sage_()
        sage: integrate(ceil(x), x, 0, infinity, algorithm='sympy')
        integrate(ceil(x), x, 0, +Infinity)
    """
    from sage.functions.other import ceil
    return ceil(self.args[0]._sage_())

def bench(p, start=2, stop=200, routine=test_irred):
    l = []
    i = ZZ(start)
    while i < stop:
	print i,
        try:
            l.append((i, routine(p, i, i+1)))
        except:
            print "failed"
	else:
	    print sum(l[-1][1])
        i += ceil(i/10)
    return l

def automorphy_factor_vector(p, a, c, k, chi, p_prec, var_prec, R):
    """
    EXAMPLES::
        
        sage: from sage.modular.pollack_stevens.families_util import automorphy_factor_vector
        sage: automorphy_factor_vector(3, 1, 3, 0, None, 4, 3, PowerSeriesRing(ZpCA(3), 'w'))
        [1 + O(3^20), O(3^21) + (3 + 3^2 + 2*3^3 + O(3^21))*w + (3^2 + 2*3^3 + O(3^22))*w^2, O(3^22) + (3^2 + 2*3^3 + O(3^22))*w + (2*3^2 + O(3^22))*w^2, O(3^22) + (3^2 + 3^3 + O(3^22))*w + (2*3^3 + O(3^23))*w^2]
        sage: automorphy_factor_vector(3, 1, 3, 2, None, 4, 3, PowerSeriesRing(ZpCA(3), 'w'))
        [1 + O(3^20),
         O(3^21) + (3 + 3^2 + 2*3^3 + O(3^21))*w + (3^2 + 2*3^3 + O(3^22))*w^2,
         O(3^22) + (3^2 + 2*3^3 + O(3^22))*w + (2*3^2 + O(3^22))*w^2,
         O(3^22) + (3^2 + 3^3 + O(3^22))*w + (2*3^3 + O(3^23))*w^2]
        sage: p, a, c, k, chi, p_prec, var_prec, R = 11, -3, 11, 0, None, 6, 4, PowerSeriesRing(ZpCA(11, 6), 'w')
        sage: automorphy_factor_vector(p, a, c, k, chi, p_prec, var_prec, R)
        [1 + O(11^6) + (7*11^2 + 4*11^3 + 11^4 + 9*11^5 + 3*11^6 + 10*11^7 + O(11^8))*w + (2*11^3 + 2*11^5 + 2*11^6 + 6*11^7 + 8*11^8 + O(11^9))*w^2 + (6*11^4 + 4*11^5 + 5*11^6 + 6*11^7 + 4*11^9 + O(11^10))*w^3,
         O(11^7) + (7*11 + 11^2 + 2*11^3 + 3*11^4 + 4*11^5 + 3*11^6 + O(11^7))*w + (2*11^2 + 4*11^3 + 5*11^4 + 10*11^5 + 10*11^6 + 2*11^7 + O(11^8))*w^2 + (6*11^3 + 6*11^4 + 2*11^5 + 9*11^6 + 9*11^7 + 2*11^8 + O(11^9))*w^3,
         O(11^8) + (3*11^2 + 4*11^4 + 2*11^5 + 6*11^6 + 10*11^7 + O(11^8))*w + (8*11^2 + 7*11^3 + 8*11^4 + 8*11^5 + 11^6 + O(11^8))*w^2 + (3*11^3 + 9*11^4 + 10*11^5 + 5*11^6 + 10*11^7 + 11^8 + O(11^9))*w^3,
         O(11^9) + (8*11^3 + 3*11^4 + 3*11^5 + 7*11^6 + 2*11^7 + 6*11^8 + O(11^9))*w + (10*11^3 + 6*11^5 + 7*11^6 + O(11^9))*w^2 + (4*11^3 + 11^4 + 8*11^8 + O(11^9))*w^3,
         O(11^10) + (2*11^4 + 9*11^5 + 8*11^6 + 5*11^7 + 4*11^8 + 6*11^9 + O(11^10))*w + (6*11^5 + 3*11^6 + 5*11^7 + 7*11^8 + 4*11^9 + O(11^10))*w^2 + (2*11^4 + 8*11^6 + 2*11^7 + 9*11^8 + 7*11^9 + O(11^10))*w^3,
         O(11^11) + (2*11^5 + 10*11^6 + 10*11^7 + 10*11^8 + 10*11^9 + 10*11^10 + O(11^11))*w + (5*11^5 + 10*11^6 + 10*11^7 + 10*11^8 + 10*11^9 + 10*11^10 + O(11^11))*w^2 + (2*11^5 + 10*11^6 + 10*11^7 + 10*11^8 + 10*11^9 + 10*11^10 + O(11^11))*w^3]
        sage: k = 2
        sage: automorphy_factor_vector(p, a, c, k, chi, p_prec, var_prec, R)
        [9 + 6*11^2 + 11^3 + 9*11^4 + 8*11^5 + O(11^6) + (8*11^2 + 8*11^3 + 10*11^4 + 6*11^5 + 5*11^7 + O(11^8))*w + (7*11^3 + 11^4 + 8*11^5 + 9*11^7 + 10*11^8 + O(11^9))*w^2 + (10*11^4 + 7*11^5 + 7*11^6 + 3*11^7 + 8*11^8 + 4*11^9 + O(11^10))*w^3,
         O(11^7) + (8*11 + 3*11^2 + 6*11^3 + 11^4 + 6*11^5 + 11^6 + O(11^7))*w + (7*11^2 + 4*11^3 + 5*11^4 + 10*11^6 + 5*11^7 + O(11^8))*w^2 + (10*11^3 + 3*11^4 + 4*11^5 + 7*11^6 + 10*11^7 + 3*11^8 + O(11^9))*w^3,
         O(11^8) + (5*11^2 + 2*11^3 + 10*11^4 + 3*11^5 + 8*11^6 + 7*11^7 + O(11^8))*w + (6*11^2 + 3*11^3 + 5*11^4 + 11^5 + 5*11^6 + 9*11^7 + O(11^8))*w^2 + (5*11^3 + 6*11^4 + 5*11^5 + 2*11^6 + 9*11^7 + 6*11^8 + O(11^9))*w^3,
         O(11^9) + (6*11^3 + 11^5 + 8*11^6 + 9*11^7 + 2*11^8 + O(11^9))*w + (2*11^3 + 8*11^4 + 4*11^5 + 6*11^6 + 11^7 + 8*11^8 + O(11^9))*w^2 + (3*11^3 + 11^4 + 3*11^5 + 11^6 + 5*11^7 + 6*11^8 + O(11^9))*w^3,
         O(11^10) + (7*11^4 + 5*11^5 + 3*11^6 + 10*11^7 + 10*11^8 + 11^9 + O(11^10))*w + (10*11^5 + 9*11^6 + 6*11^7 + 6*11^8 + 10*11^9 + O(11^10))*w^2 + (7*11^4 + 11^5 + 7*11^6 + 5*11^7 + 6*11^8 + 2*11^9 + O(11^10))*w^3,
         O(11^11) + (7*11^5 + 3*11^6 + 8*11^8 + 5*11^9 + 8*11^10 + O(11^11))*w + (11^5 + 6*11^6 + 7*11^7 + 11^8 + 2*11^10 + O(11^11))*w^2 + (7*11^5 + 3*11^6 + 8*11^8 + 5*11^9 + 8*11^10 + O(11^11))*w^3]
    """
    S = PolynomialRing(R, 'z')
    z = S.gens()[0]
    w = R.gen()
    aut = S(1)
    for n in range(1, var_prec):
        ## RP: I doubled the precision in "z" here to account for the loss of precision from plugging in arg in below
        ## This should be done better.
        LB = logpp_binom(n, p, ceil(p_prec * (p-1)/(p-2)))
        ta = ZZ(Qp(p, 2 * max(p_prec, var_prec)).teichmuller(a))
        arg = (a / ta - 1) / p + c / (p * ta) * z
        aut += LB(arg).truncate(p_prec) * (w ** n)
    aut *= (ta ** k)
    #if not (chi is None):
    #    aut *= chi(a)
    aut = aut.list()
    len_aut = len(aut)
    if len_aut == p_prec:
        return aut
    elif len_aut > p_prec:
        return aut[:p_prec]
    return aut + [R.zero_element()] * (p_prec - len_aut)

 def exp(self,workprec=Infinity):
     from sage.functions.other import ceil
     if workprec is Infinity:
         raise ApproximationError("unable to compute exp to infinite precision")
     parent = self.parent()
     pow = parent(1)
     res = parent(1)
     val = self.valuation()
     iter = ceil(workprec / (val - 1/(parent._p-1))) + 1
     for i in range(1,iter):
         pow *= self / parent(i)
         res += pow
         res = res.truncate(workprec)
     return res

    def _normalisation_factor_zz(self, tau=3):
        r"""
        This function returns an approximation of `∑_{x ∈ \ZZ^n}
        \exp(-|x|_2^2/(2σ²))`, i.e. the normalisation factor such that the sum
        over all probabilities is 1 for `\ZZⁿ`.

        If this ``self.B`` is not an identity matrix over `\ZZ` a
        ``NotImplementedError`` is raised.

        INPUT:

        - ``tau`` -- all vectors `v` with `|v|_∞ ≤ τ·σ` are enumerated
                     (default: ``3``).

        EXAMPLES::

            sage: from sage.stats.distributions.discrete_gaussian_lattice import DiscreteGaussianDistributionLatticeSampler
            sage: n = 3; sigma = 1.0; m = 1000
            sage: D = DiscreteGaussianDistributionLatticeSampler(ZZ^n, sigma)
            sage: f = D.f
            sage: c = D._normalisation_factor_zz(); c
            15.528...

            sage: l = [D() for _ in range(m)]
            sage: v = vector(ZZ, n, (0, 0, 0))
            sage: l.count(v), ZZ(round(m*f(v)/c))
            (57, 64)

        """
        if self.B != identity_matrix(ZZ, self.B.nrows()):
            raise NotImplementedError("This function is only implemented when B is an identity matrix.")

        f = self.f
        n = self.B.ncols()
        sigma = self._sigma
        return sum(f(x) for x in _iter_vectors(n, -ceil(tau * sigma),
                                               ceil(tau * sigma)))

    def next_height(self,N):
        """
        Return the next permissable height greater than or equal to N for
         curves in self's family.
        
        WARNING: This function my return a height for which only singular
                  curves exist. For example, in the short Weierstrass case
                  height 0 is permissable, as the curve Y^2 = X^3 (uniquely)
                  has height zero.

        INPUT:

            - ``N`` -- A non-negative integer

        OUTPUT:

            - A tuple consisting of three elements of the form
              (H, C, I) such that 
              H: The smallest height >= N
              C: A list of coefficients for curves of this height
              I: A list of indices indicating which of the above coefficients
              achieve this height. The remaining values in C  indicate the 
              max absolute value those coefficients are allowed to obtain
              without altering the height.

              For example, the tuple (4, [1, 2], [1]) for the short Weierstrass
              case denotes set of curves with height 4; these are all of the
              form Y^2 = X^3 + A*X + B, where B=2 and A ranges between -1 and 1.

        EXAMPLES::

            sage: from sage.schemes.elliptic_curves.curve_enumerator import *
            sage: C = CurveEnumerator(family="short_weierstrass")
            sage: C.next_height(4) 
            (4, [1, 2], [1])
            sage: C.next_height(60)
            (64, [4, 8], [0, 1])

            sage: C.next_height(-100)
            Traceback (most recent call last):
            ...
            AssertionError: Input must be non-negative integer.
        """

        assert N>=0, "Input must be non-negative integer."

        coeffs = [ceil(N**(1/n))-1 for n in self._pows]
        height = max([coeffs[i]**(self._pows[i]) for i in range(self._num_coeffs)])
        return self._height_increment(coeffs)

    def e_string_to_ground_state(self):
        r"""
        Returns a string of integers in the index set `(i_1,\ldots,i_k)` such that `e_{i_k} \cdots e_{i_1} self` is
        the ground state.

        INPUT:

        - ``self`` -- an element of a tensor product of Kirillov-Reshetikhin crystals.

        OUTPUT: a tuple of integers `(i_1,\ldots,i_k)`

        This method is only defined when ``self`` is an element of a tensor product of affine Kirillov-Reshetikhin crystals.
        It calculates a path from ``self`` to a ground state path using Demazure arrows as defined in
        Lemma 7.3 in [SchillingTingley2011]_.

        EXAMPLES::

            sage: K = KirillovReshetikhinCrystal(['A',2,1],1,1)
            sage: T = TensorProductOfCrystals(K,K)
            sage: t = T.module_generators[0]
            sage: t.e_string_to_ground_state()
            (0, 2)

            sage: K = KirillovReshetikhinCrystal(['C',2,1],1,1)
            sage: T = TensorProductOfCrystals(K,K)
            sage: t = T.module_generators[0]; t
            [[[1]], [[1]]]
            sage: t.e_string_to_ground_state()
            (0,)
            sage: x=t.e(0)
            sage: x.e_string_to_ground_state()
            ()
            sage: y=t.f_string([1,2,1,1,0]); y
            [[[2]], [[1]]]
            sage: y.e_string_to_ground_state()
            ()
        """
        assert self.parent().crystals[0].__module__ == 'sage.combinat.crystals.kirillov_reshetikhin', \
            "All crystals in the tensor product need to be Kirillov-Reshetikhin crystals"
        I = self.index_set()
        I.remove(0)
        ell = max(ceil(K.s()/K.cartan_type().c()[K.r()]) for K in self.parent().crystals)
        for i in I:
            if self.epsilon(i) > 0:
                return (i,) + (self.e(i)).e_string_to_ground_state()
        if self.epsilon(0) > ell:
            return (0,) + (self.e(0)).e_string_to_ground_state()
        return ()

    def calc_prec(self):
        if self.prec != None:
            return self.prec

        mp0 = MatrixPowerSexp(self.bsym,self.N-1,iprec=self.iprec,x0=self.x0sym)

        sexp_precision=RR(1)*log(abs(self.sexp_1(0.5)-mp0.sexp_1(0.5)),2.0)
        self.prec = (-sexp_precision).floor()
        print "sexp precision: " , self.prec

        cprec = self.prec+ceil(log(self.N)/log(2.0))

        #self.eigenvalues = [ ev.n(cprec) for ev in self.eigenvalues ]
        #self.IM = self.IM.n(cprec)
        #self.b = self.bsym.n(cprec)
        return self

def trunc(i):
    """
    Truncates to the integer closer to zero

    EXAMPLES::

        sage: from sage.combinat.crystals.tensor_product import trunc
        sage: trunc(-3/2), trunc(-1), trunc(-1/2), trunc(0), trunc(1/2), trunc(1), trunc(3/2)
        (-1, -1, 0, 0, 0, 1, 1)
        sage: isinstance(trunc(3/2), Integer)
        True
    """
    if i>= 0:
        return floor(i)
    else:
        return ceil(i)

    def __iter__(self):
        r"""
        Iterate over all GL(2,Z)-reduced  semi positive forms
        which are within the bounds of this precision.

        EXAMPLES::

            sage: from sage.modular.siegel.siegel_modular_form_prec import SiegelModularFormPrecision
            sage: prec = SiegelModularFormPrecision(11)
            sage: for k in prec.__iter__(): print k
            (0, 0, 0)
            (0, 0, 1)
            (0, 0, 2)
            (1, 0, 1)
            (1, 0, 2)
            (1, 1, 1)
            (1, 1, 2)


        NOTE
            The forms are enumerated in lexicographic order.
        """
        if self.__type == 'disc':
            bound = self.get_contents_bound_for_semi_definite_forms()
            for c in xrange(0, bound):
                yield (0,0,c)

            atop = isqrt(self.__prec // 3)
            if 3*atop*atop == self.__prec: atop -= 1
            for a in xrange(1,atop + 1):
                for b in xrange(a+1):
                    for c in xrange(a, ceil((b**2 + self.__prec)/(4*a))):
                        yield (a,b,c)

        elif 'box' == self.__type:
            (am, bm, cm) = self.__prec
            for a in xrange(am):
                for b in xrange( min(bm,a+1)):
                    for c in xrange(a, cm):
                        yield (a,b,c)

        else:
            raise RuntimeError, "Unexpected value of self.__type"
            
        raise StopIteration

def newton_iteration(G, n):
    r"""Returns a truncated series `y = y(x)` satisfying

    .. math::

        G(x,y(x)) \equiv 0 \bmod{x^r}

    where $r = \ceil{\log_2{n}}$. Based on the algorithm in [XXX].

    Parameters
    ----------
    G, x, y : polynomial
        A polynomial in `x` and `y`.
    n : int
        Requested degree of the series expansion.

    Notes
    -----
    This algorithm returns the series up to order :math:`2^r > n`. Any
    choice of order below :math:`2^r` will return the same series.

    """
    R = G.parent()
    x,y = R.gens()
    if n < 0:
        raise ValueError('Number of terms must be positive. (n=%d'%n)
    elif n == 0:
        return R(0)

    phi = G
    phiprime = phi.derivative(y)
    try:
        pi = R(x).polynomial(x)
        gi = R(0)
        si = R(phiprime(x,gi)).polynomial(x).inverse_mod(pi)
    except NotImplementedError:
        raise ValueError('Newton iteration for computing regular part of '
                         'Puiseux expansion failed. Curve is most likely '
                         'not regular at center.')

    r = ceil(log(n,2))
    for i in range(r):
        gi,si,pi = newton_iteration_step(phi,phiprime,gi,si,pi)
    return R(gi)

    def split_xor(self, monomial_list, equal_zero):
        """
        Split XOR chains into subchains.

        INPUT:

        - ``monomial_list`` - a list of monomials
        - ``equal_zero`` - is the constant coefficient zero?

        EXAMPLE::

            sage: from sage.sat.converters.polybori import CNFEncoder
            sage: from sage.sat.solvers.dimacs import DIMACS
            sage: B.<a,b,c,d,e,f> = BooleanPolynomialRing()
            sage: ce = CNFEncoder(DIMACS(), B, cutting_number=3)
            sage: ce.split_xor([1,2,3,4,5,6], False)
            [[[1, 7], False], [[7, 2, 8], True], [[8, 3, 9], True], [[9, 4, 10], True], [[10, 5, 11], True], [[11, 6], True]]

            sage: ce = CNFEncoder(DIMACS(), B, cutting_number=4)
            sage: ce.split_xor([1,2,3,4,5,6], False)
            [[[1, 2, 7], False], [[7, 3, 4, 8], True], [[8, 5, 6], True]]

            sage: ce = CNFEncoder(DIMACS(), B, cutting_number=5)
            sage: ce.split_xor([1,2,3,4,5,6], False)
            [[[1, 2, 3, 7], False], [[7, 4, 5, 6], True]]
        """
        c = self.cutting_number

        nm = len(monomial_list)
        step = ceil((c-2)/ZZ(nm) * nm)
        M = []

        new_variables = []
        for j in range(0, nm, step):
            m =  new_variables + monomial_list[j:j+step]
            if (j + step) < nm:
                new_variables = [self.var(None)]
                m += new_variables
            M.append([m, equal_zero])
            equal_zero = True
        return M

def compute_series_truncations(f, alpha):
    r"""Computes Puiseux series at :math:`x=\alpha` with necessary terms.

    The Puiseux series expansions of :math:`f = f(x,y)` centered at
    :math:`\alpha` are computed up to the number of terms needed for the
    integral basis algorithm to be successful. The expansion degree bounds are
    determined by :func:`compute_expansion_bounds`.

    Parameters
    ----------
    f : polynomial
    alpha : complex

    Returns
    -------
    list : PuiseuxXSeries
        A list of Puiseux series expansions centered at :math:`x = \alpha` with
        enough terms to compute integral bases as SymPy expressions.

    """
    # compute the parametric Puiseix series with the minimal number of terms
    # needed to distinguish them.
    pt = puiseux(f,alpha)
    px = [p for P in pt for p in P.xseries()]

    # compute the orders necessary for the integral basis algorithm. the orders
    # are on the Puiseux x-series (non-parametric) so scale by the ramification
    # index of each series
    N = compute_expansion_bounds(px)
    for i in range(len(N)):
        e = px[i].ramification_index
        N[i] = ceil(N[i]*e)

    order = max(N) + 1
    for pti in pt:
        pti.extend(order=order)

    # recompute the corresponding x-series with the extened terms
    px = [p for P in pt for p in P.xseries()]
    return px

def numeric_converter(value, cur=None):
    """
    Used for converting numeric values from Postgres to Python.

    INPUT:

    - ``value`` -- a string representing a decimal number.
    - ``cur`` -- a cursor, unused

    OUTPUT:

    - either a sage integer (if there is no decimal point) or a real number whose precision depends on the number of digits in value.
    """
    if value is None:
        return None
    if '.' in value:
        # The following is a good guess for the bit-precision,
        # but we use LmfdbRealLiterals to ensure that our number
        # prints the same as we got it.
        prec = ceil(len(value)*3.322)
        return LmfdbRealLiteral(RealField(prec), value)
    else:
        return Integer(value)