def _make_monic(self, f):
        r"""
        Let alpha be a root of f.  This function returns a monic
        integral polynomial g and an element d of the base field such
        that g(alpha*d) = 0.
        
        EXAMPLES::
        
            sage: K.<x> = FunctionField(QQ); R.<y> = K[];
            sage: L.<alpha> = K.extension(x^2*y^5 - 1/x)
            sage: g, d = L._make_monic(L.polynomial()); g,d
            (y^5 - x^12, x^3)
            sage: g(alpha*d)
            0
        """
        R = f.base_ring()
        if not isinstance(R, RationalFunctionField):
            raise NotImplementedError
        
        # make f monic
        n = f.degree()
        c = f.leading_coefficient()
        if c != 1:
            f = f / c

        # find lcm of denominators
        from sage.rings.arith import lcm
        # would be good to replace this by minimal...
        d = lcm([b.denominator() for b in f.list() if b])
        if d != 1:
            x = f.parent().gen()
            g = (d**n) * f(x/d)
        else:
            g = f
        return g, d

    def generalised_level(self):
        r"""
        Return the generalised level of self, i.e. the least common multiple of
        the widths of all cusps. 
        
        If self is *even*, Wohlfart's theorem tells us that this is equal to
        the (conventional) level of self when self is a congruence subgroup.
        This can fail if self is odd, but the actual level is at most twice the
        generalised level. See the paper by Kiming, Schuett and Verrill for
        more examples.
       
        EXAMPLE::

            sage: Gamma0(18).generalised_level()
            18
            sage: sage.modular.arithgroup.congroup_generic.CongruenceSubgroup(5).index()
            Traceback (most recent call last):
            ...
            NotImplementedError

        In the following example, the actual level is twice the generalised
        level. This is the group `G_2` from Example 17 of K-S-V.
        
        ::

            sage: G = CongruenceSubgroup(8, [ [1,1,0,1], [3,-1,4,-1] ])
            sage: G.level()
            8
            sage: G.generalised_level()
            4
        """
        return arith.lcm([self.cusp_width(c) for c in self.cusps()])

    def clear_denominators(self):
        r"""
        scales by the least common multiple of the denominators.

        OUTPUT: None.

        EXAMPLES::

            sage: R.<t>=PolynomialRing(QQ)
            sage: P.<x,y,z>=ProjectiveSpace(FractionField(R),2)
            sage: Q=P([t,3/t^2,1])
            sage: Q.clear_denominators(); Q
            (t^3 : 3 : t^2)

        ::

            sage: R.<x>=PolynomialRing(QQ)
            sage: K.<w>=NumberField(x^2-3)
            sage: P.<x,y,z>=ProjectiveSpace(K,2)
            sage: Q=P([1/w,3,0])
            sage: Q.clear_denominators(); Q
            (w : 9 : 0)

        ::

            sage: P.<x,y,z>=ProjectiveSpace(QQ,2)
            sage: X=P.subscheme(x^2-y^2);
            sage: Q=X([1/2,1/2,1]);
            sage: Q.clear_denominators(); Q
            (1 : 1 : 2)
        """
        self.scale_by(lcm([self[i].denominator() for i in range(self.codomain().ambient_space().dimension_relative())]))

    def _roots_univariate_polynomial(self, p, ring=None, multiplicities=None, algorithm=None):
        r"""
        Return a list of pairs ``(root,multiplicity)`` of roots of the polynomial ``p``.

        If the argument ``multiplicities`` is set to ``False`` then return the
        list of roots.

        .. SEEALSO::

            :meth:`_factor_univariate_polynomial`

        EXAMPLES::

            sage: R.<x> = PolynomialRing(GF(5),'x')
            sage: K = GF(5).algebraic_closure('t')

            sage: sorted((x^6 - 1).roots(K,multiplicities=False))
            [1, 4, 2*t2 + 1, 2*t2 + 2, 3*t2 + 3, 3*t2 + 4]
            sage: ((K.gen(2)*x - K.gen(3))**2).roots(K)
            [(3*t6^5 + 2*t6^4 + 2*t6^2 + 3, 2)]

            sage: for _ in xrange(10):
            ....:     p = R.random_element(degree=randint(2,8))
            ....:     for r in p.roots(K, multiplicities=False):
            ....:         assert p(r).is_zero(), "r={} is not a root of p={}".format(r,p)

        """
        from sage.rings.arith import lcm
        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

        # first build a polynomial over some finite field
        coeffs = [v.as_finite_field_element(minimal=True) for v in p.list()]
        l = lcm([c[0].degree() for c in coeffs])
        F, phi = self.subfield(l)
        P = p.parent().change_ring(F)

        new_coeffs = [self.inclusion(c[0].degree(), l)(c[1]) for c in coeffs]

        polys = [(g,m,l,phi) for g,m in P(new_coeffs).factor()]
        roots = []    # a list of pair (root,multiplicity)
        while polys:
            g,m,l,phi = polys.pop()
        
            if g.degree() == 1: # found a root
                r = phi(-g.constant_coefficient())
                roots.append((r,m))
            else: # look at the extension of degree g.degree() which contains at
                  # least one root of g
                ll = l * g.degree()
                psi = self.inclusion(l, ll)
                FF, pphi = self.subfield(ll)
                # note: there is no coercion from the l-th subfield to the ll-th
                # subfield. The line below does the conversion manually.
                g = PolynomialRing(FF, 'x')(map(psi, g))
                polys.extend((gg,m,ll,pphi) for gg,_ in g.factor())

        if multiplicities:
            return roots
        else:
            return [r[0] for r in roots]

    def ncusps(self):
        r"""
        Return the number of orbits of cusps (regular or otherwise) for this subgroup.

        EXAMPLE::

            sage: GammaH(33,[2]).ncusps()
            8
            sage: GammaH(32079, [21676]).ncusps()
            28800

        AUTHORS:

        - Jordi Quer

        """
        N = self.level()
        H = self._list_of_elements_in_H()
        c = ZZ(0)
        for d in [d for d in N.divisors() if d**2 <= N]:
            Nd = lcm(d,N//d)
            Hd = set([x % Nd for x in H])
            lenHd = len(Hd)
            if Nd-1 not in Hd: lenHd *= 2
            summand = euler_phi(d)*euler_phi(N//d)//lenHd
            if d**2 == N:
                c = c + summand
            else:
                c = c + 2*summand
        return c

    def _iterator_discriminant(self, epsilon, offset) :
        """
        Return iterator over all possible discriminants of a Fourier index given a content ``eps``
        and a discriminant offset modulo `D`.
        
        INPUT:
            - `\epsilon` -- Integer; Content of an index.
            - ``offset`` -- Integer; Offset modulo D of the discriminant.

        OUTPUT:
            - Iterator over Integers.

        TESTS::
            sage: fac = HermitianModularFormD2Factory(4,-3)
            sage: list(fac._iterator_discriminant(2, 3))
            [0, 12, 24]
            sage: list(fac._iterator_discriminant(2, 2))
            [8, 20]
            sage: list(fac._iterator_discriminant(1, 2))
            [2, 5, 8, 11, 14, 17, 20, 23, 26]
        """
        mD = -self.__D
        periode = lcm(mD, epsilon**2)
        offset = crt(offset, 0, mD, epsilon**2) % periode 
        
        return xrange( offset,
                       self.__precision._enveloping_discriminant_bound(),
                       periode )

    def _init_from_Vrepresentation(self, ambient_dim, vertices, rays, lines, minimize=True):
        """
        Construct polyhedron from V-representation data.

        INPUT:

        - ``ambient_dim`` -- integer. The dimension of the ambient space.
        
        - ``vertices`` -- list of point. Each point can be specified
           as any iterable container of
           :meth:`~sage.geometry.polyhedron.base.base_ring` elements.
        
        - ``rays`` -- list of rays. Each ray can be specified as any
          iterable container of
          :meth:`~sage.geometry.polyhedron.base.base_ring` elements.
        
        - ``lines`` -- list of lines. Each line can be specified as
          any iterable container of
          :meth:`~sage.geometry.polyhedron.base.base_ring` elements.

        EXAMPLES::

            sage: p = Polyhedron(backend='ppl')
            sage: from sage.geometry.polyhedron.backend_ppl import Polyhedron_QQ_ppl
            sage: Polyhedron_QQ_ppl._init_from_Vrepresentation(p, 2, [], [], [])
        """
        gs = Generator_System()
        if vertices is None: vertices = []
        for v in vertices:
            d = lcm([denominator(v_i) for v_i in v])
            dv = [ ZZ(d*v_i) for v_i in v ]
            gs.insert(point(Linear_Expression(dv, 0), d))
        if rays is None: rays = []
        for r in rays:
            d = lcm([denominator(r_i) for r_i in r])
            dr = [ ZZ(d*r_i) for r_i in r ]
            gs.insert(ray(Linear_Expression(dr, 0)))
        if lines is None: lines = []
        for l in lines:
            d = lcm([denominator(l_i) for l_i in l])
            dl = [ ZZ(d*l_i) for l_i in l ]
            gs.insert(line(Linear_Expression(dl, 0)))
        self._ppl_polyhedron = C_Polyhedron(gs)
        self._init_Vrepresentation_from_ppl(minimize)
        self._init_Hrepresentation_from_ppl(minimize)

def circle_drops(A,B):
    # Drops going around the unit circle for those A and B.
    # See http://user.math.uzh.ch/dehaye/thesis_students/Nicolas_Wider-Integrality_of_factorial_ratios.pdf
    # for longer description (not original, better references exist)    
    marks = lcm(lcm(A),lcm(B))    
    tmp = [0 for i in range(marks)]
    for a in A:
#        print tmp
        for i in range(a):
            if gcd(i, a) == 1:
                tmp[i*marks/a] -= 1
    for b in B:
#        print tmp
        for i in range(b):
            if gcd(i, b) == 1:
                tmp[i*marks/b] += 1
#    print tmp
    return [sum(tmp[:j]) for j in range(marks)]

 def __init__( self, q):
     """
     We initialize by a list of integers $[a_1,...,a_N]$. The
     lattice self is then $L = (L,\beta) = (G^{-1}\ZZ^n/\ZZ^n, G[x])$,
     where $G$ is the symmetric matrix
     $G=[a_1,...,a_n;*,a_{n+1},....;..;* ... * a_N]$,
     i.e.~the $a_j$ denote the elements above the diagonal of $G$.
     """
     self.__form = q
     # We compute the rank
     N = len(q)
     assert is_square( 1+8*N)
     n = Integer( (-1+isqrt(1+8*N))/2)
     self.__rank = n
     # We set up the Gram matrix
     self.__G = matrix( IntegerRing(), n, n)
     i = j = 0
     for a in q:
         self.__G[i,j] = self.__G[j,i] = Integer(a)
         if j < n-1:
             j += 1
         else:
             i += 1
             j = i
     # We compute the level
     Gi = self.__G**-1
     self.__Gi = Gi
     a = lcm( map( lambda x: x.denominator(), Gi.list()))
     I = Gi.diagonal()
     b =  lcm( map( lambda x: (x/2).denominator(), I))
     self.__level = lcm( a, b)
     # We define the undelying module and the ambient space
     self.__module = FreeModule( IntegerRing(), n)
     self.__space = self.__module.ambient_vector_space()
     # We compute a shadow vector
     self.__shadow_vector = self.__space([(a%2)/2 for a in self.__G.diagonal()])*Gi
     # We define a basis
     M = Matrix( IntegerRing(), n, n, 1)
     self.__basis = self.__module.basis()
     # We prepare a cache
     self.__dual_vectors = None
     self.__values = None
     self.__chi = {}

    def _roots_univariate_polynomial(self, p, ring=None, multiplicities=None, algorithm=None):
        r"""
        Return a list of pairs ``(root,multiplicity)`` of roots of the polynomial ``p``.

        If the argument ``multiplicities`` is set to ``False`` then return the
        list of roots.

        .. SEEALSO::

            :meth:`_factor_univariate_polynomial`

        EXAMPLES::

            sage: R.<x> = PolynomialRing(GF(5),'x')
            sage: K = GF(5).algebraic_closure('t')

            sage: sorted((x^6 - 1).roots(K,multiplicities=False))
            [1, 4, 2*t2 + 1, 2*t2 + 2, 3*t2 + 3, 3*t2 + 4]
            sage: ((K.gen(2)*x - K.gen(3))**2).roots(K)
            [(3*t6^5 + 2*t6^4 + 2*t6^2 + 3, 2)]

            sage: for _ in xrange(10):
            ....:     p = R.random_element(degree=randint(2,8))
            ....:     for r in p.roots(K, multiplicities=False):
            ....:         assert p(r).is_zero()

        """
        from sage.rings.arith import lcm
        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

        # first build a polynomial over some finite field
        coeffs = [v.as_finite_field_element(minimal=True) for v in p.list()]
        l = lcm([c[0].degree() for c in coeffs])
        F, phi = self.subfield(l)
        P = p.parent().change_ring(F)

        new_coeffs = [self.inclusion(c[0].degree(), l)(c[1]) for c in coeffs]

        roots = []    # a list of pair (root,multiplicity)
        for g, m in P(new_coeffs).factor():
            if g.degree() == 1:
                r = phi(-g.constant_coefficient())
                roots.append((r,m))
            else:
                ll = l * g.degree()
                psi = self.inclusion(l, ll)
                FF, pphi = self.subfield(ll)
                gg = PolynomialRing(FF, 'x')(map(psi, g))
                for r, _ in gg.roots():  # note: we know that multiplicity is 1
                    roots.append((pphi(r), m))

        if multiplicities:
            return roots
        else:
            return [r[0] for r in roots]

    def shadow_level( self):
        """
        Return the level of $L_ev$, where $L_ev$ is the kernel
        of the map $x\mapsto e(G[x]/2)$.

        REMARK
            Lemma: If $L$ is odd, if $s$ is a shadow vector, and if $N$ denotes the
            denominator of $G[s]/2$, then the level of $L_ev$
            equals the lcm of the order of $s$, $N$ and the level of $L$.

            (For the proof:  the requested level is the smallest integer $l$
            such that $lG[x]/2$ is integral for all $x$ in $L^*$ and $x$ in $s+L^*$
            since $L_ev^* = L^* \cup s+L^*$. This $l$ is the smallest integer
            such that $lG[s]/2$, $lG[x]/2$ and $ls^tGx$ are integral for all $x$ in $L^*$.)
        """
        if self.is_even():
            return self.level()
        s = self.a_shadow_vector()
        N = self.beta(s).denominator()
        h = lcm( [x.denominator() for x in s])
        return lcm( [h, N, self.level()])

    def _init_from_Hrepresentation(self, ambient_dim, ieqs, eqns, minimize=True):
        """
        Construct polyhedron from H-representation data.

        INPUT:

        - ``ambient_dim`` -- integer. The dimension of the ambient space.
        
        - ``ieqs`` -- list of inequalities. Each line can be specified
          as any iterable container of
          :meth:`~sage.geometry.polyhedron.base.base_ring` elements.

        - ``eqns`` -- list of equalities. Each line can be specified
          as any iterable container of
          :meth:`~sage.geometry.polyhedron.base.base_ring` elements.

        EXAMPLES::

            sage: p = Polyhedron(backend='ppl')
            sage: from sage.geometry.polyhedron.backend_ppl import Polyhedron_QQ_ppl
            sage: Polyhedron_QQ_ppl._init_from_Hrepresentation(p, 2, [], [])
        """
        cs = Constraint_System()
        if ieqs is None: ieqs = []
        for ieq in ieqs:
            d = lcm([denominator(ieq_i) for ieq_i in ieq])
            dieq = [ ZZ(d*ieq_i) for ieq_i in ieq ]
            b = dieq[0]
            A = dieq[1:]
            cs.insert(Linear_Expression(A, b) >= 0)
        if eqns is None: eqns = []
        for eqn in eqns:
            d = lcm([denominator(eqn_i) for eqn_i in eqn])
            deqn = [ ZZ(d*eqn_i) for eqn_i in eqn ]
            b = deqn[0]
            A = deqn[1:]
            cs.insert(Linear_Expression(A, b) == 0)
        self._ppl_polyhedron = C_Polyhedron(cs)
        self._init_Vrepresentation_from_ppl(minimize)
        self._init_Hrepresentation_from_ppl(minimize)

def _heckebasis(M):
    r"""
    Gives a basis of the hecke algebra of M as a ZZ-module

    INPUT:

    - ``M`` - a hecke module

    OUTPUT:

    - a list of hecke algebra elements represented as matrices

    EXAMPLES::

        sage: M = ModularSymbols(11,2,1)
        sage: sage.modular.hecke.algebra._heckebasis(M)
        [Hecke operator on Modular Symbols space of dimension 2 for Gamma_0(11) of weight 2 with sign 1 over Rational Field defined by:
        [1 0]
        [0 1],
        Hecke operator on Modular Symbols space of dimension 2 for Gamma_0(11) of weight 2 with sign 1 over Rational Field defined by:
        [0 1]
        [0 5]]
    """
    d = M.rank()
    VV = QQ ** (d ** 2)
    WW = ZZ ** (d ** 2)
    MM = MatrixSpace(QQ, d)
    MMZ = MatrixSpace(ZZ, d)
    S = []
    Denom = []
    B = []
    B1 = []
    for i in xrange(1, M.hecke_bound() + 1):
        v = M.hecke_operator(i).matrix()
        den = v.denominator()
        Denom.append(den)
        S.append(v)
    den = lcm(Denom)
    for m in S:
        B.append(WW((den * m).list()))
    UU = WW.submodule(B)
    B = UU.basis()
    for u in B:
        u1 = u.list()
        m1 = M.hecke_algebra()(MM(u1), check=False)
        # m1 = MM(u1)
        B1.append((1 / den) * m1)
    return B1

def summand(part, n):
    """
    Create the summand used in the Harrison count for a given partition.

    Args:
        part (tuple): A partition of `n` represented as a tuple.
        n (int): The integer for which `part` is a partition.

    Returns:
        int: The summand corresponding to the partition `part` of `n`.
    """

    t = 1
    count = list(cycle_count(part, n)) + (factorial(n)-n)*[0]
    for i in range(1,n+1):
        for j in range(1,n+1):
            s = sum([d*(count[d-1]) for d in divisors(lcm(i,j))])
            t = t*(s**(count[i-1]*count[j-1]*gcd(i,j)))
    t = t*factorial(n)/(prod(factorial(count[d-1])*(d**(count[d-1])) for d in range(1,n+1)))
    return t

 def exponent(self):
     """
     Return the exponent of this abelian group.
     
     EXAMPLES::
     
         sage: G = AbelianGroup([2,3,7]); G
         Multiplicative Abelian Group isomorphic to C2 x C3 x C7
         sage: G.exponent()
         42
         sage: G = AbelianGroup([2,4,6]); G
         Multiplicative Abelian Group isomorphic to C2 x C4 x C6
         sage: G.exponent()
         12
     """
     try:
         return self.__exponent
     except AttributeError:
         self.__exponent = e = lcm(self.invariants())
         return e

    def clear_denominators(self):
        r"""
        scales by the least common multiple of the denominators.

        OUTPUT: None.

        EXAMPLES::

            sage: R.<t> = PolynomialRing(QQ)
            sage: P.<x,y,z> = ProjectiveSpace(FractionField(R), 2)
            sage: Q = P([t, 3/t^2, 1])
            sage: Q.clear_denominators(); Q
            (t^3 : 3 : t^2)

        ::

            sage: R.<x> = PolynomialRing(QQ)
            sage: K.<w> = NumberField(x^2 - 3)
            sage: P.<x,y,z> = ProjectiveSpace(K, 2)
            sage: Q = P([1/w, 3, 0])
            sage: Q.clear_denominators(); Q
            (w : 9 : 0)

        ::

            sage: P.<x,y,z> = ProjectiveSpace(QQ, 2)
            sage: X = P.subscheme(x^2 - y^2);
            sage: Q = X([1/2, 1/2, 1]);
            sage: Q.clear_denominators(); Q
            (1 : 1 : 2)

        ::

            sage: PS.<x,y> = ProjectiveSpace(QQ, 1)
            sage: Q = PS.point([1, 2/3], False); Q
            (1 : 2/3)
            sage: Q.clear_denominators(); Q
            (3 : 2)
        """
        self.scale_by(lcm([t.denominator() for t in self]))

    def nregcusps(self):
        r"""
        Return the number of orbits of regular cusps for this subgroup. A cusp is regular
        if we may find a parabolic element generating the stabiliser of that
        cusp whose eigenvalues are both +1 rather than -1. If G contains -1,
        all cusps are regular.

        EXAMPLES::

            sage: GammaH(20, [17]).nregcusps()
            4
            sage: GammaH(20, [17]).nirregcusps()
            2
            sage: GammaH(3212, [2045, 2773]).nregcusps()
            1440
            sage: GammaH(3212, [2045, 2773]).nirregcusps()
            720

        AUTHOR:

        - Jordi Quer
        """
        if self.is_even():
            return self.ncusps()

        N = self.level()
        H = self._list_of_elements_in_H()

        c = ZZ(0)
        for d in [d for d in divisors(N) if d**2 <= N]:
            Nd = lcm(d,N//d)
            Hd = set([x%Nd for x in H])
            if Nd - 1 not in Hd:
                summand = euler_phi(d)*euler_phi(N//d)//(2*len(Hd))
                if d**2==N:
                    c = c + summand
                else:
                    c = c + 2*summand
        return c

#.........这里部分代码省略.........
            w = matrix(R, [[self.u0],[self.u1]])
            Fac = list(m.factor())
            Periods = {}

            #To compute the period mod m, we compute the least integer n such that A^n*w == w.  This necessarily
            #divides the order of A as a matrix in GL_2(Z/mZ).

            #We compute the period modulo all distinct prime powers dividing m, and combine via the lcm.
            #To compute the period mod p^e, we first compute the order mod p.  Then the period mod p^e
            #must divide p^{4e-4}*period(p), as the subgroup of matrices mod p^e, which reduce to
            #the identity mod p is of order (p^{e-1})^4.  So we compute the period mod p^e by successively
            #multiplying the period mod p by powers of p.

            for i in Fac:
                p = i[0]; e = i[1]
                #first compute the period mod p
                if p in self._period_dict:
                    perp = self._period_dict[p]
                else:
                    F = A.change_ring(GF(p))
                    v = w.change_ring(GF(p))
                    FF = F**(p-1)
                    p1fac = list((p-1).factor())

                    #The order of any matrix in GL_2(F_p) either divides p(p-1) or (p-1)(p+1).
                    #The order divides p-1 if it is diagaonalizable.  In any case, det(F^(p-1))=1,
                    #so if tr(F^(p-1)) = 2, then it must be triangular of the form [[1,a],[0,1]].
                    #The order of the subgroup of matrices of this form is p, so the order must divide
                    #p(p-1) -- in fact it must be a multiple of p.  If this is not the case, then the
                    #order divides (p-1)(p+1).  As the period divides the order of the matrix in GL_2(F_p),
                    #these conditions hold for the period as well.

                    #check if the order divides (p-1)
                    if FF*v == v:
                        M = p-1
                        Mfac = p1fac

                    #check if the trace is 2, then the order is a multiple of p dividing p*(p-1)
                    elif (FF).trace() == 2:
                        M = p-1
                        Mfac = p1fac
                        F = F**p        #replace F by F^p as now we only need to determine the factor dividing (p-1)

                    #otherwise it will divide (p+1)(p-1)
                    else :
                        M = (p+1)*(p-1)
                        p2fac = list((p+1).factor())        #factor the (p+1) and (p-1) terms seperately and then combine for speed
                        Mfac_dic = {}
                        for i in list(p1fac + p2fac):
                            if i[0] not in Mfac_dic:
                                Mfac_dic[i[0]] = i[1]
                            else :
                                Mfac_dic[i[0]] = Mfac_dic[i[0]] + i[1]
                        Mfac = [(i,Mfac_dic[i]) for i in Mfac_dic]

                    #Now use a fast order algorithm to compute the period.  We know that the period divides
                    #M = i_1*i_2*...*i_l where the i_j denote not necessarily distinct prime factors.  As
                    #F^M*v == v, for each i_j, if F^(M/i_j)*v == v, then the period divides (M/i_j).  After
                    #all factors have been iterated over, the result is the period mod p.

                    Mfac = list(Mfac)
                    C=[]

                    #expand the list of prime factors so every factor is with multiplicity 1

                    for i in xrange(len(Mfac)):
                        for j in xrange(Mfac[i][1]):
                            C.append(Mfac[i][0])

                    Mfac = C
                    n = M
                    for i in Mfac:
                        b = Integer(n/i)
                        if F**b*v == v:
                            n = b
                    perp = n

                #Now compute the period mod p^e by steping up by multiples of p
                F = A.change_ring(Integers(p**e))
                v = w.change_ring(Integers(p**e))
                FF = F**perp
                if FF*v == v:
                    perpe = perp
                else :
                    tries = 0
                    while True:
                        tries += 1
                        FF = FF**p
                        if FF*v == v:
                            perpe = perp*p**tries
                            break
                Periods[p] = perpe

            #take the lcm of the periods mod all distinct primes dividing m
            period = 1
            for p in Periods:
                period = lcm(Periods[p],period)

            self._period_dict[m] = period        #cache the period mod m
            return period

#.........这里部分代码省略.........
        specifies the algorithm to be used:

         - ``'qfsolve'`` -- Use Denis Simon's GP script ``qfsolve``
           (see ``sage.quadratic_forms.qfsolve.qfsolve``)

         - ``'rnfisnorm'`` -- Use PARI's function rnfisnorm
           (cannot be combined with ``obstruction = True``)

         - ``'local'`` -- Check if a local solution exists for all primes
           and infinite places of `\QQ` and apply the Hasse principle
           (cannot be combined with ``point = True``)

         - ``'default'`` -- Use ``'qfsolve'``

         - ``'magma'`` (requires Magma to be installed) --
           delegates the task to the Magma computer algebra
           system.

        EXAMPLES::

            sage: C = Conic(QQ, [1, 2, -3])
            sage: C.has_rational_point(point = True)
            (True, (1 : 1 : 1))
            sage: D = Conic(QQ, [1, 3, -5])
            sage: D.has_rational_point(point = True)
            (False, 3)
            sage: P.<X,Y,Z> = QQ[]
            sage: E = Curve(X^2 + Y^2 + Z^2); E
            Projective Conic Curve over Rational Field defined by X^2 + Y^2 + Z^2
            sage: E.has_rational_point(obstruction = True)
            (False, -1)

        The following would not terminate quickly with
        ``algorithm = 'rnfisnorm'`` ::

            sage: C = Conic(QQ, [1, 113922743, -310146482690273725409])
            sage: C.has_rational_point(point = True)
            (True, (-76842858034579/5424 : -5316144401/5424 : 1))
            sage: C.has_rational_point(algorithm = 'local', read_cache = False)
            True
            sage: C.has_rational_point(point=True, algorithm='magma', read_cache=False) # optional - magma
            (True, (30106379962113/7913 : 12747947692/7913 : 1))

        TESTS:

        Create a bunch of conics over `\QQ`, check if ``has_rational_point`` runs without errors
        and returns consistent answers for all algorithms. Check if all points returned are valid. ::

            sage: l = Sequence(cartesian_product_iterator([[-1, 0, 1] for i in range(6)]))
            sage: c = [Conic(QQ, a) for a in l if a != [0,0,0] and a != (0,0,0,0,0,0)]
            sage: d = []
            sage: d = [[C]+[C.has_rational_point(algorithm = algorithm, read_cache = False, obstruction = (algorithm != 'rnfisnorm'), point = (algorithm != 'local')) for algorithm in ['local', 'qfsolve', 'rnfisnorm']] for C in c[::10]] # long time: 7 seconds
            sage: assert all([e[1][0] == e[2][0] and e[1][0] == e[3][0] for e in d])
            sage: assert all([e[0].defining_polynomial()(Sequence(e[i][1])) == 0 for e in d for i in [2,3] if e[1][0]])
        """
        if read_cache:
            if self._rational_point is not None:
                if point or obstruction:
                    return True, self._rational_point
                else:
                    return True
            if self._local_obstruction is not None:
                if point or obstruction:
                    return False, self._local_obstruction
                else:
                    return False
            if (not point) and self._finite_obstructions == [] and \
               self._infinite_obstructions == []:
                if obstruction:
                    return True, None
                return True
        if self.has_singular_point():
            if point:
                return self.has_singular_point(point = True)
            if obstruction:
                return True, None
            return True
        if algorithm == 'default' or algorithm == 'qfsolve':
            M = self.symmetric_matrix()
            M *= lcm([ t.denominator() for t in M.list() ])
            pt = qfsolve(M)
            if pt in ZZ:
                if self._local_obstruction == None:
                    self._local_obstruction = pt
                if point or obstruction:
                    return False, pt
                return False
            pt = self.point([pt[0], pt[1], pt[2]])
            if point or obstruction:
                return True, pt
            return True
        ret = ProjectiveConic_number_field.has_rational_point( \
                                           self, point = point, \
                                           obstruction = obstruction, \
                                           algorithm = algorithm, \
                                           read_cache = read_cache)
        if point or obstruction:
            if is_RingHomomorphism(ret[1]):
                ret[1] = -1
        return ret

    def parametrization(self, point=None, morphism=True):
        r"""
        Return a parametrization `f` of ``self`` together with the
        inverse of `f`.

        If ``point`` is specified, then that point is used
        for the parametrization. Otherwise, use ``self.rational_point()``
        to find a point.

        If ``morphism`` is True, then `f` is returned in the form
        of a Scheme morphism. Otherwise, it is a tuple of polynomials
        that gives the parametrization.

        ALGORITHM:

        Uses Denis Simon's GP script ``qfparam``.
        See ``sage.quadratic_forms.qfsolve.qfparam``.

        EXAMPLES ::

            sage: c = Conic([1,1,-1])
            sage: c.parametrization()
            (Scheme morphism:
              From: Projective Space of dimension 1 over Rational Field
              To:   Projective Conic Curve over Rational Field defined by x^2 + y^2 - z^2
              Defn: Defined on coordinates by sending (x : y) to
                    (2*x*y : -x^2 + y^2 : x^2 + y^2),
             Scheme morphism:
              From: Projective Conic Curve over Rational Field defined by x^2 + y^2 - z^2
              To:   Projective Space of dimension 1 over Rational Field
              Defn: Defined on coordinates by sending (x : y : z) to
                    (1/2*x : 1/2*y + 1/2*z))

        An example with ``morphism = False`` ::

            sage: R.<x,y,z> = QQ[]
            sage: C = Curve(7*x^2 + 2*y*z + z^2)
            sage: (p, i) = C.parametrization(morphism = False); (p, i)
            ([-2*x*y, 7*x^2 + y^2, -2*y^2], [-1/2*x, -1/2*z])
            sage: C.defining_polynomial()(p)
            0
            sage: i[0](p) / i[1](p)
            x/y

        A ``ValueError`` is raised if ``self`` has no rational point ::

            sage: C = Conic(x^2 + 2*y^2 + z^2)
            sage: C.parametrization()
            Traceback (most recent call last):
            ...
            ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + 2*y^2 + z^2 has no rational points over Rational Field!

        A ``ValueError`` is raised if ``self`` is not smooth ::

            sage: C = Conic(x^2 + y^2)
            sage: C.parametrization()
            Traceback (most recent call last):
            ...
            ValueError: The conic self (=Projective Conic Curve over Rational Field defined by x^2 + y^2) is not smooth, hence does not have a parametrization.
        """
        if (not self._parametrization is None) and not point:
            par = self._parametrization
        else:
            if not self.is_smooth():
                raise ValueError("The conic self (=%s) is not smooth, hence does not have a parametrization." % self)
            if point == None:
                point = self.rational_point()
            point = Sequence(point)
            Q = PolynomialRing(QQ, 'x,y')
            [x, y] = Q.gens()
            gens = self.ambient_space().gens()
            M = self.symmetric_matrix()
            M *= lcm([ t.denominator() for t in M.list() ])
            par1 = qfparam(M, point)
            B = Matrix([[par1[i][j] for j in range(3)] for i in range(3)])
            # self is in the image of B and does not lie on a line,
            # hence B is invertible
            A = B.inverse()
            par2 = [sum([A[i,j]*gens[j] for j in range(3)]) for i in [1,0]]
            par = ([Q(pol(x/y)*y**2) for pol in par1], par2)
            if self._parametrization is None:
                self._parametrization = par
        if not morphism:
            return par
        P1 = ProjectiveSpace(self.base_ring(), 1, 'x,y')
        return P1.hom(par[0],self), self.Hom(P1)(par[1], check = False)