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 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.arithgroup_perm.HsuExample18().generalised_level()
            24

        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 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 __common_minimal_basering(chi, psi):
    """
    Find the smallest basering over which chi and psi are valued, and
    return new chi and psi valued in that ring.

    EXAMPLES::

        sage: sage.modular.modform.eis_series.__common_minimal_basering(DirichletGroup(1)[0], DirichletGroup(1)[0])
        (Dirichlet character modulo 1 of conductor 1, Dirichlet character modulo 1 of conductor 1)

        sage: sage.modular.modform.eis_series.__common_minimal_basering(DirichletGroup(3).0, DirichletGroup(5).0)
        (Dirichlet character modulo 3 of conductor 3 mapping 2 |--> -1, Dirichlet character modulo 5 of conductor 5 mapping 2 |--> zeta4)

        sage: sage.modular.modform.eis_series.__common_minimal_basering(DirichletGroup(12).0, DirichletGroup(36).0)
        (Dirichlet character modulo 12 of conductor 4 mapping 7 |--> -1, 5 |--> 1, Dirichlet character modulo 36 of conductor 4 mapping 19 |--> -1, 29 |--> 1)
    """
    chi = chi.minimize_base_ring()
    psi = psi.minimize_base_ring()
    n = lcm(chi.base_ring().zeta().multiplicative_order(),
                  psi.base_ring().zeta().multiplicative_order())
    if n <= 2:
        K = QQ
    else:
        K = CyclotomicField(n)
    chi = chi.change_ring(K)
    psi = psi.change_ring(K)
    return chi, psi

    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
        # 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 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 _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 range(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.arith.all 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')([psi(_) for _ in 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 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))

def _heckebasis(M):
    r"""
    Return 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 range(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 _upper_bound_for_longest_cycle(self, s):
        """
        EXAMPLES::

            sage: cis = species.PartitionSpecies().cycle_index_series()
            sage: cis._upper_bound_for_longest_cycle([4])
            4
            sage: cis._upper_bound_for_longest_cycle([3,1])
            3
            sage: cis._upper_bound_for_longest_cycle([2,2])
            2
            sage: cis._upper_bound_for_longest_cycle([2,1,1])
            2
            sage: cis._upper_bound_for_longest_cycle([1,1,1,1])
            1
        """
        if s == []:
            return 1
        return min(self._card(sum(s)), lcm(list(s)))

    def exponent(self):
        """
        Return the exponent of this finite abelian group.

        OUTPUT: Integer

        EXAMPLES::

            sage: t = J0(33).hecke_operator(7)
            sage: G = t.kernel()[0]; G
            Finite subgroup with invariants [2, 2, 2, 2, 4, 4] over QQ of Abelian variety J0(33) of dimension 3
            sage: G.exponent()
            4
        """
        try:
            return self.__exponent
        except AttributeError:
            e = lcm(self.invariants())
            self.__exponent = e
            return e

    def additive_order(self):
        """
        Return the additive order of this element.

        EXAMPLES::

            sage: V = span([[1/2,1,1],[3/2,2,1],[0,0,1]],ZZ); W = V.span([2*V.0+4*V.1, 9*V.0+12*V.1, 4*V.2])
            sage: Q = V/W; Q
            Finitely generated module V/W over Integer Ring with invariants (4, 12)
            sage: Q.0.additive_order()
            4
            sage: Q.1.additive_order()
            12
            sage: (Q.0+Q.1).additive_order()
            12
            sage: V = span([[1/2,1,1],[3/2,2,1],[0,0,1]],ZZ); W = V.span([2*V.0+4*V.1, 9*V.0+12*V.1])
            sage: Q = V/W; Q
            Finitely generated module V/W over Integer Ring with invariants (12, 0)
            sage: Q.0.additive_order()
            12
            sage: type(Q.0.additive_order())
            <type 'sage.rings.integer.Integer'>
            sage: Q.1.additive_order()
            +Infinity
        """
        Q = self.parent()
        I = Q.invariants()
        v = self.vector()

        from sage.rings.all import infinity, Mod, Integer
        from sage.arith.all import lcm
        n = Integer(1)
        for i, a in enumerate(I):
            if a == 0:
                if v[i] != 0:
                    return infinity
            else:
                n = lcm(n, Mod(v[i],a).additive_order())
        return n

    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

def find_gens_list(klist, r = 0, omax = Infinity, use_lucas = True, verbose = False):
    '''
    Use Rain's method to find a generator of a subfield of degree r
    within the ambient fields in klist.

    The algorithm is done in two steps:

    1. Find a small integer m, such that the m-th roots of unity
       generate an extension of the subfields, together with some
       additional constraints detailed below. This is done by
       `find_root_order`.

    2. Find a uniquely defined Galois orbit of m-th Gaussian periods.
       Return arbitrary elements of these orbits.
       This is done by `find_unique_orbit`.

    Read the docstrings of `find_root_order` and `find_unique_orbit`
    to find out more on the algorithm.
    '''
    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)

    # Find a suitable m, see doc below
    o, G = find_root_order(p, [ngcd, nlcm], r, verbose=verbose)
    if o > omax:
        raise RuntimeError

    # Construct extension of the fields, if needed
    #
    # Note: `find_unique_orbit` is horribly slow if k is not a
    # field. Some composita tricks would be welcome.
    return tuple(find_gen_with_data(k, r, o, G, use_lucas=use_lucas, verbose=verbose) for k in klist)

def carmichael_lambda(n):
    r"""
    Return the Carmichael function of a positive integer ``n``.

    The Carmichael function of `n`, denoted `\lambda(n)`, is the smallest
    positive integer `k` such that `a^k \equiv 1 \pmod{n}` for all
    `a \in \ZZ/n\ZZ` satisfying `\gcd(a, n) = 1`. Thus, `\lambda(n) = k`
    is the exponent of the multiplicative group `(\ZZ/n\ZZ)^{\ast}`.

    INPUT:

    - ``n`` -- a positive integer.

    OUTPUT:

    - The Carmichael function of ``n``.

    ALGORITHM:

    If `n = 2, 4` then `\lambda(n) = \varphi(n)`. Let `p \geq 3` be an odd
    prime and let `k` be a positive integer. Then
    `\lambda(p^k) = p^{k - 1}(p - 1) = \varphi(p^k)`. If `k \geq 3`, then
    `\lambda(2^k) = 2^{k - 2}`. Now consider the case where `n > 3` is
    composite and let `n = p_1^{k_1} p_2^{k_2} \cdots p_t^{k_t}` be the
    prime factorization of `n`. Then

    .. MATH::

        \lambda(n)
        = \lambda(p_1^{k_1} p_2^{k_2} \cdots p_t^{k_t})
        = \text{lcm}(\lambda(p_1^{k_1}), \lambda(p_2^{k_2}), \dots, \lambda(p_t^{k_t}))

    EXAMPLES:

    The Carmichael function of all positive integers up to and including 10::

        sage: from sage.crypto.util import carmichael_lambda
        sage: list(map(carmichael_lambda, [1..10]))
        [1, 1, 2, 2, 4, 2, 6, 2, 6, 4]

    The Carmichael function of the first ten primes::

        sage: list(map(carmichael_lambda, primes_first_n(10)))
        [1, 2, 4, 6, 10, 12, 16, 18, 22, 28]

    Cases where the Carmichael function is equivalent to the Euler phi
    function::

        sage: carmichael_lambda(2) == euler_phi(2)
        True
        sage: carmichael_lambda(4) == euler_phi(4)
        True
        sage: p = random_prime(1000, lbound=3, proof=True)
        sage: k = randint(1, 1000)
        sage: carmichael_lambda(p^k) == euler_phi(p^k)
        True

    A case where `\lambda(n) \neq \varphi(n)`::

        sage: k = randint(1, 1000)
        sage: carmichael_lambda(2^k) == 2^(k - 2)
        True
        sage: carmichael_lambda(2^k) == 2^(k - 2) == euler_phi(2^k)
        False

    Verifying the current implementation of the Carmichael function using
    another implementation. The other implementation that we use for
    verification is an exhaustive search for the exponent of the
    multiplicative group `(\ZZ/n\ZZ)^{\ast}`. ::

        sage: from sage.crypto.util import carmichael_lambda
        sage: n = randint(1, 500)
        sage: c = carmichael_lambda(n)
        sage: def coprime(n):
        ....:     return [i for i in range(n) if gcd(i, n) == 1]
        sage: def znpower(n, k):
        ....:     L = coprime(n)
        ....:     return list(map(power_mod, L, [k]*len(L), [n]*len(L)))
        sage: def my_carmichael(n):
        ....:     for k in range(1, n):
        ....:         L = znpower(n, k)
        ....:         ones = [1] * len(L)
        ....:         T = [L[i] == ones[i] for i in xrange(len(L))]
        ....:         if all(T):
        ....:             return k
        sage: c == my_carmichael(n)
        True

    Carmichael's theorem states that `a^{\lambda(n)} \equiv 1 \pmod{n}`
    for all elements `a` of the multiplicative group `(\ZZ/n\ZZ)^{\ast}`.
    Here, we verify Carmichael's theorem. ::

        sage: from sage.crypto.util import carmichael_lambda
        sage: n = randint(1, 1000)
        sage: c = carmichael_lambda(n)
        sage: ZnZ = IntegerModRing(n)
        sage: M = ZnZ.list_of_elements_of_multiplicative_group()
        sage: ones = [1] * len(M)
        sage: P = [power_mod(a, c, n) for a in M]
        sage: P == ones
#.........这里部分代码省略.........

#.........这里部分代码省略.........

                F = GF(ell)
                a0 = F(self.u0); a1 = F(self.u1) #a0 and a1 are variables for terms in sequence
                bf, cf = F(self.b), F(self.c)

                for n in range(Bound): # n is the index of the a0

                    #Check whether a0 is a perfect power mod ell
                    if _is_p_power_mod(a0, p, ell) :
                        #if a0 is a perfect power mod ell, check if nth term is ppower
                        if _is_p_power(self(n), p):
                            powers.append(n)

                    a0, a1 = a1, bf*a1 + cf*a0        #step up the variables

            else :

                powers = []        #documents the indices of the sequence that provably correspond to pth powers
                cong = [0]        #list of necessary congruences on the index for it to correspond to pth powers
                Possible_count = {}    #keeps track of the number of rounds a congruence lasts in cong

                #These parameters are involved in how we choose primes to increase the modulus
                qqold = 1        #we believe that we know complete information coming from primes good by qqold
                M1 = 1            #we have congruences modulo M1, this may not be the tightest list
                M2 = p            #we want to move to have congruences mod M2
                qq = 1            #the largest prime power divisor of M1 is qq

                #This loop ups the modulus.
                while True:

                    #Try to get good data mod M2

                    #patience of how long we should search for a "good prime"
                    patience = 0.01 * _estimated_time(lcm(M2,p*next_prime_power(qq)), M1, len(cong), p)
                    tries = 0

                    #This loop uses primes to get a small set of congruences mod M2.
                    while True:

                        #only proceed if took less than patience time to find the next good prime
                        ell = _next_good_prime(p, self, qq, patience, qqold)
                        if ell:

                            #gather congruence data for the sequence mod ell, which will be mod period(ell) = modu
                            cong1, modu = _find_cong1(p, self, ell)

                            CongNew = []        #makes a new list from cong that is now mod M = lcm(M1, modu) instead of M1
                            M = lcm(M1, modu)
                            for k in range(M // M1):
                                for i in cong:
                                    CongNew.append(k * M1 + i)
                            cong = set(CongNew)

                            M1 = M

                            killed_something = False        #keeps track of when cong1 can rule out a congruence in cong

                            #CRT by hand to gain speed
                            for i in list(cong):
                                if not (i % modu in cong1):        #congruence in cong is inconsistent with any in cong1
                                    cong.remove(i)            #remove that congruence
                                    killed_something = True

                            if M1 == M2:
                                if not killed_something:
                                    tries += 1

#.........这里部分代码省略.........
            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 diagonalizable.  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 separately 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 range(len(Mfac)):
                        for j in range(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 stepping 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

def row_reduced_form(M,transformation=False):
    """
    This function computes a row reduced form of a matrix over a rational
    function field `k(x)`, for `k` a field.

    INPUT:

     - `M` - a matrix over `k(x)` or `k[x]` for `k` a field.
     - `transformation` - A boolean (default: `False`). If this boolean is set to `True` a second matrix is output (see OUTPUT).
     
    OUTPUT:

    If `transformation` is `False`, the output is `W`, a row reduced form of `M`.
    
    If `transformation` is `True`, this function will output a pair `(W,N)` consisting of two matrices over `k(x)`:

    1. `W` - a row reduced form of `M`.
    2. `N` - an invertible matrix over `k(x)` satisfying `NW = M`.

    EXAMPLES:

    The fuction expects matrices over the rational function field, but
    other examples below show how one can provide matrices over the ring
    of polynomials (whose quotient field is the rational function field).

    ::

        sage: R.<t> = GF(3)['t']
        sage: K = FractionField(R)
        sage: import sage.matrix.matrix_misc
        sage: sage.matrix.matrix_misc.row_reduced_form(matrix([[(t-1)^2/t],[(t-1)]]))
        [(2*t + 1)/t]
        [          0]

    The last example shows the usage of the transformation parameter.
        
    ::
        sage: Fq.<a> = GF(2^3)
        sage: Fx.<x> = Fq[]
        sage: A = matrix(Fx,[[x^2+a,x^4+a],[x^3,a*x^4]])
        sage: from sage.matrix.matrix_misc import row_reduced_form
        sage: row_reduced_form(A,transformation=True)
        (
        [(a^2 + 1)*x^3 + x^2 + a                       a]  [      1 a^2 + 1]
        [                    x^3                   a*x^4], [      0                 1]
        )
            
    NOTES:

    See docstring for row_reduced_form method of matrices for
    more information.
    """

    # determine whether M has polynomial or rational function coefficients
    R0 = M.base_ring()

    #Compute the base polynomial ring
    if R0 in _Fields:
        R = R0.base()
    else:
        R = R0
    from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
    if not is_PolynomialRing(R) or not R.base_ring().is_field():
        raise TypeError("the coefficients of M must lie in a univariate polynomial ring over a field")

    t = R.gen()

    # calculate least-common denominator of matrix entries and clear
    # denominators. The result lies in R
    from sage.arith.all import lcm
    from sage.matrix.constructor import matrix
    from sage.misc.functional import numerator
    if R0 in _Fields:
        den = lcm([a.denominator() for a in M.list()])
        num = matrix([[numerator(_) for _ in v] for v in (M*den).rows()])
    else:
        # No need to clear denominators
        num = M

    r = [list(v) for v in num.rows()]

    if transformation:
        N = matrix(num.nrows(), num.nrows(), R(1)).rows()


    rank = 0
    num_zero = 0
    if M.is_zero():
        num_zero = len(r)
    while rank != len(r) - num_zero:
        # construct matrix of leading coefficients
        v = []
        for w in map(list, r):
            # calculate degree of row (= max of degree of entries)
            d = max([e.numerator().degree() for e in w])

            # extract leading coefficients from current row
            x = []
            for y in w:
                if y.degree() >= d and d >= 0:   x.append(y.coefficients(sparse=False)[d])
#.........这里部分代码省略.........

    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 the PARI/GP function ``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, x^2 + 7*y^2, -2*x^2], [-1/2*x, 1/7*y + 1/14*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 is None:
                point = self.rational_point()
            point = Sequence(point)
            Q = PolynomialRing(QQ, '