def module_composition_factors(self, algorithm=None):
        r"""
        Return a list of triples consisting of [base field, dimension,
        irreducibility], for each of the Meataxe composition factors
        modules. The ``algorithm="verbose"`` option returns more information,
        but in Meataxe notation.

        EXAMPLES::

            sage: F=GF(3);MS=MatrixSpace(F,4,4)
            sage: M=MS(0)
            sage: M[0,1]=1;M[1,2]=1;M[2,3]=1;M[3,0]=1
            sage: G = MatrixGroup([M])
            sage: G.module_composition_factors()
            [(Finite Field of size 3, 1, True),
             (Finite Field of size 3, 1, True),
             (Finite Field of size 3, 2, True)]
            sage: F = GF(7); MS = MatrixSpace(F,2,2)
            sage: gens = [MS([[0,1],[-1,0]]),MS([[1,1],[2,3]])]
            sage: G = MatrixGroup(gens)
            sage: G.module_composition_factors()
            [(Finite Field of size 7, 2, True)]

        Type ``G.module_composition_factors(algorithm='verbose')`` to get a
        more verbose version.

        For more on MeatAxe notation, see
        http://www.gap-system.org/Manuals/doc/ref/chap69.html
        """
        from sage.misc.sage_eval import sage_eval
        F = self.base_ring()
        if not(F.is_finite()):
            raise NotImplementedError("Base ring must be finite.")
        q = F.cardinality()
        gens = self.gens()
        n = self.degree()
        MS = MatrixSpace(F,n,n)
        mats = [] # initializing list of mats by which the gens act on self
        W = self.matrix_space().row_space()
        for g in gens:
            p = MS(g.matrix())
            m = p.rows()
            mats.append(m)
        mats_str = str(gap([[list(r) for r in m] for m in mats]))
        gap.eval("M:=GModuleByMats("+mats_str+", GF("+str(q)+"))")
        gap.eval("MCFs := MTX.CompositionFactors( M )")
        N = eval(gap.eval("Length(MCFs)"))
        if algorithm == "verbose":
            print(gap.eval('MCFs') + "\n")
        L = []
        for i in range(1,N+1):
            gap.eval("MCF := MCFs[%s]"%i)
            L.append(tuple([sage_eval(gap.eval("MCF.field")),
                            eval(gap.eval("MCF.dimension")),
                            sage_eval(gap.eval("MCF.IsIrreducible")) ]))
        return sorted(L)

    def _element_constructor_(self, x):
        """
        Construct a complex number.

        EXAMPLES::

            sage: CC((1,2)) # indirect doctest
            1.00000000000000 + 2.00000000000000*I
        """
        if not isinstance(x, (real_mpfr.RealNumber, tuple)):
            if isinstance(x, complex_double.ComplexDoubleElement):
                return complex_number.ComplexNumber(self, x.real(), x.imag())
            elif isinstance(x, str):
                # TODO: this is probably not the best and most
                # efficient way to do this.  -- Martin Albrecht
                return complex_number.ComplexNumber(self,
                            sage_eval(x.replace(' ',''), locals={"I":self.gen(),"i":self.gen()}))

            late_import()
            if isinstance(x, NumberFieldElement_quadratic) and list(x.parent().polynomial()) == [1, 0, 1]:
                (re, im) = list(x)
                return complex_number.ComplexNumber(self, re, im)

            try:
                return self(x.sage())
            except (AttributeError, TypeError):
                pass
            try:
                return x._complex_mpfr_field_( self )
            except AttributeError:
                pass
        return complex_number.ComplexNumber(self, x)

def gen_laguerre(n, a, x):
    """
    Returns the generalized Laguerre polynomial for integers `n > -1`.
    Typically, `a = 1/2` or `a = -1/2`.

    REFERENCES:

    - Table on page 789 in [ASHandbook]_.

    EXAMPLES::

        sage: x = PolynomialRing(QQ, 'x').gen()
        sage: gen_laguerre(2,1,x)
        1/2*x^2 - 3*x + 3
        sage: gen_laguerre(2,1/2,x)
        1/2*x^2 - 5/2*x + 15/8
        sage: gen_laguerre(2,-1/2,x)
        1/2*x^2 - 3/2*x + 3/8
        sage: gen_laguerre(2,0,x)
        1/2*x^2 - 2*x + 1
        sage: gen_laguerre(3,0,x)
        -1/6*x^3 + 3/2*x^2 - 3*x + 1
    """
    _init()
    return sage_eval(maxima.eval("gen_laguerre(%s,%s,x)" % (ZZ(n), a)), locals={"x": x})

    def _element_constructor_(self, x):
        r"""
        Convert ``x`` into ``self``.

        EXAMPLES::

            sage: R = ShuffleAlgebra(QQ,'xy')
            sage: x, y = R.gens()
            sage: R(3) # indirect doctest
            3*B[word: ]
            sage: R(x)
            B[word: x]
        """
        P = x.parent()
        if isinstance(P, ShuffleAlgebra):
            if P is self:
                return x
            if not (P is self.base_ring()):
                return self.element_class(self, x.monomial_coefficients())
        # ok, not a shuffle algebra element (or should not be viewed as one).
        if isinstance(x, basestring):
            from sage.misc.sage_eval import sage_eval
            return sage_eval(x,locals=self.gens_dict())
        R = self.base_ring()
        # coercion via base ring
        x = R(x)
        if x == 0:
            return self.element_class(self,{})
        else:
            return self.from_base_ring_from_one_basis(x)

def ultraspherical(n, a, x):
    """
    Returns the ultraspherical (or Gegenbauer) polynomial for integers
    `n > -1`.

    Computed using Maxima.

    REFERENCE:

    - AS 22.5.27

    EXAMPLES::

        sage: x = PolynomialRing(QQ, 'x').gen()
        sage: ultraspherical(2,3/2,x)
        15/2*x^2 - 3/2
        sage: ultraspherical(2,1/2,x)
        3/2*x^2 - 1/2
        sage: ultraspherical(1,1,x)
        2*x
        sage: t = PolynomialRing(RationalField(),"t").gen()
        sage: gegenbauer(3,2,t)
        32*t^3 - 12*t
    """
    _init()
    return sage_eval(maxima.eval("ultraspherical(%s,%s,x)" % (ZZ(n), a)), locals={"x": x})

def spherical_bessel_Y(n,var, algorithm="maxima"):
    r"""
    Returns the spherical Bessel function of the second kind for
    integers n -1.
    
    Reference: AS 10.1.9 page 437 and AS 10.1.15 page 439.
    
    EXAMPLES::
    
        sage: x = PolynomialRing(QQ, 'x').gen()
        sage: spherical_bessel_Y(2,x)
        -((3/x^2 - 1)*cos(x) + 3*sin(x)/x)/x
    """
    if algorithm=="scipy":
        import scipy.special
        ans = str(scipy.special.sph_yn(int(n),float(var)))
        ans = ans.replace("(","")
        ans = ans.replace(")","")
        ans = ans.replace("j","*I")
        return sage_eval(ans)
    elif algorithm == 'maxima':
        _init()
        return meval("spherical_bessel_y(%s,%s)"%(ZZ(n),var))
    else:
        raise ValueError, "unknown algorithm '%s'"%algorithm

    def eval(self, s, globals=None, locals=None):
        r"""
        EXAMPLES::

            sage: html.eval('<hr>')
            <html><font color='black'><hr></font></html>
            ''
        """
        if globals is None:
            globals = {}
        if locals is None:
            locals = {}
        s = str(s)
        s = math_parse(s)
        t = ''
        while len(s) > 0:
            i = s.find('<sage>')
            if i == -1:
                 t += s
                 break
            j = s.find('</sage>')
            if j == -1:
                 t += s
                 break
            t += s[:i] + '<script type="math/tex">%s</script>'%\
                     latex(sage_eval(s[6+i:j], locals=locals))
            s = s[j+7:]
        print "<html><font color='black'>%s</font></html>"%t
        return ''

 def eigenvalues( self, index_list):
     evs = DB.find( { 'owner_id': self.__id,
                      'data_type': 'ev',
                      'index': { '$in': [ str(l) for l in index_list]}
                      })
     loc_f = self.__field.gens_dict() 
     return dict( (eval(ev['index']),sage_eval( ev['data'], locals = loc_f)) for  ev in evs)

def legendre_P(n, x):
    """
    Returns the Legendre polynomial of the first kind for integers
    `n > -1`.

    REFERENCE:

    - AS 22.5.35 page 779.

    EXAMPLES::

        sage: P.<t> = QQ[]
        sage: legendre_P(2,t)
        3/2*t^2 - 1/2
        sage: legendre_P(3, 1.1)
        1.67750000000000
        sage: legendre_P(3, 1 + I)
        7/2*I - 13/2
        sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
        [-179  242]
        [-484  547]
        sage: legendre_P(3, GF(11)(5))
        8
    """
    _init()
    return sage_eval(maxima.eval("legendre_p(%s,x)" % ZZ(n)), locals={"x": x})

        def right_coset_representatives(self):
            r"""
            Return the right coset representatives of ``self``.

            EXAMPLES::

                sage: W = ReflectionGroup(['A',2])                      # optional - gap3
                sage: for w in W:                                       # optional - gap3
                ....:     rcr = w.right_coset_representatives()         # optional - gap3
                ....:     print("%s %s"%(w.reduced_word(),              # optional - gap3
                ....:                    [v.reduced_word() for v in rcr]))  # optional - gap3
                [] [[], [2], [1], [2, 1], [1, 2], [1, 2, 1]]
                [2] [[], [2], [1]]
                [1] [[], [1], [1, 2]]
                [1, 2] [[]]
                [2, 1] [[]]
                [1, 2, 1] [[], [2], [2, 1]]
            """
            from sage.combinat.root_system.reflection_group_complex import _gap_return

            W = self.parent()
            T = W.reflections()
            T_fix = [i + 1 for i in T.keys() if self.fix_space().is_subspace(T[i].fix_space())]
            S = str(
                gap3(
                    "ReducedRightCosetRepresentatives(%s,ReflectionSubgroup(%s,%s))"
                    % (W._gap_group._name, W._gap_group._name, T_fix)
                )
            )
            return sage_eval(_gap_return(S, coerce_obj="W"), locals={"self": self, "W": W})

def hermite(n, x):
    """
    Returns the Hermite polynomial for integers `n > -1`.

    REFERENCE:

    - AS 22.5.40 and 22.5.41, page 779.

    EXAMPLES::

        sage: x = PolynomialRing(QQ, 'x').gen()
        sage: hermite(2,x)
        4*x^2 - 2
        sage: hermite(3,x)
        8*x^3 - 12*x
        sage: hermite(3,2)
        40
        sage: S.<y> = PolynomialRing(RR)
        sage: hermite(3,y)
        8.00000000000000*y^3 - 12.0000000000000*y
        sage: R.<x,y> = QQ[]
        sage: hermite(3,y^2)
        8*y^6 - 12*y^2
        sage: w = var('w')
        sage: hermite(3,2*w)
        8*(8*w^2 - 3)*w
    """
    _init()
    return sage_eval(maxima.eval("hermite(%s,x)" % ZZ(n)), locals={"x": x})

def gen_laguerre(n, a, x):
    """
    Returns the generalized Laguerre polynomial for integers `n > -1`.
    Typically, a = 1/2 or a = -1/2.

    REFERENCE:

    - table on page 789 in AS.

    EXAMPLES::

        sage: x = PolynomialRing(QQ, 'x').gen()
        sage: gen_laguerre(2,1,x)
        1/2*x^2 - 3*x + 3
        sage: gen_laguerre(2,1/2,x)
        1/2*x^2 - 5/2*x + 15/8
        sage: gen_laguerre(2,-1/2,x)
        1/2*x^2 - 3/2*x + 3/8
        sage: gen_laguerre(2,0,x)
        1/2*x^2 - 2*x + 1
        sage: gen_laguerre(3,0,x)
        -1/6*x^3 + 3/2*x^2 - 3*x + 1
    """
    from sage.functions.all import sqrt

    _init()
    return sage_eval(maxima.eval("gen_laguerre(%s,%s,x)" % (ZZ(n), a)), locals={"x": x})

def best_known_covering_design_www(v, k, t, verbose=False):
    r"""
    Gives the best known `(v,k,t)` covering design, using the database
    available at `<http://www.ccrwest.org/>`_

    INPUeT:

    - ``v`` -- integer, the size of the point set for the design
    - ``k`` -- integer, the number of points per block
    - ``t`` -- integer, the size of sets covered by the blocks
    - ``verbose`` -- bool (default=``False``), print verbose message

    OUTPUT:

    A :class:`CoveringDesign` object representing the ``(v,k,t)``-covering
    design with smallest number of blocks available in the database.

    EXAMPLES::

        sage: from sage.combinat.designs.covering_design import best_known_covering_design_www
        sage: C = best_known_covering_design_www(7, 3, 2)   # optional - internet
        sage: print(C)                                     # optional - internet
        C(7,3,2) = 7
        Method: lex covering
        Submitted on: 1996-12-01 00:00:00
        0  1  2
        0  3  4
        0  5  6
        1  3  5
        1  4  6
        2  3  6
        2  4  5

    This function raises a ValueError if the ``(v,k,t)`` parameters are not
    found in the database.
    """
    # import compatible with py2 and py3
    from six.moves.urllib.request import urlopen

    from sage.misc.sage_eval import sage_eval

    v = int(v)
    k = int(k)
    t = int(t)

    param = ("?v=%s&k=%s&t=%s"%(v,k,t))

    url = "http://www.ccrwest.org/cover/get_cover.php"+param
    if verbose:
        print("Looking up the bounds at %s" % url)
    f = urlopen(url)
    s = f.read()
    f.close()

    if 'covering not in database' in s:   #not found
        str = "no (%d,%d,%d) covering design in database\n"%(v,k,t)
        raise ValueError(str)

    return sage_eval(s)

    def __init__( self, doc):

        self.__collection = doc.get( 'collection')
        self.__name = doc.get( 'name')

        weight = doc.get( 'weight')
        self.__weight = sage_eval( weight) if weight else weight
       
        field = doc.get( 'field')
        R = PolynomialRing( IntegerRing(), name = 'x')
        self.__field = sage_eval( field, locals = R.gens_dict()) if field else field

        self.__explicit_formula = doc.get( 'explicit_formula')
        self.__type = doc.get( 'type')
        self.__is_eigenform = doc.get( 'is_eigenform')
        self.__is_integral = doc.get( 'is_integral')
        self.__representation = doc.get( 'representation')
        self.__id = doc.get( '_id')

def bessel_Y(nu,z,algorithm="maxima", prec=53):
    r"""
    Implements the "Y-Bessel function", or "Bessel function of the 2nd
    kind", with index (or "order") nu and argument z.
    
    .. note::

       Currently only prec=53 is supported.
    
    Defn::
    
                    cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
                   -------------------------------------------------
                                     sin(nu*pi)
    
    if nu is not an integer and by taking a limit otherwise.
    
    Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.
    
    This is computed using Maxima by default.
    
    EXAMPLES::
    
        sage: bessel_Y(2,1.1,"scipy")
        -1.4314714939...
        sage: bessel_Y(2,1.1)   
        -1.4314714939590...
        sage: bessel_Y(3.001,2.1) 
        -1.0299574976424...

    TESTS::

        sage: bessel_Y(2,1.1, algorithm="pari")
        Traceback (most recent call last):
        ...
        NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms
    """
    if algorithm=="scipy":
        if prec != 53:
            raise ValueError, "for the scipy algorithm the precision must be 53"
        import scipy.special
        ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))
        ans = ans.replace("(","")
        ans = ans.replace(")","")
        ans = ans.replace("j","*I")
        ans = sage_eval(ans)
        return real(ans) if z in RR else ans
    elif algorithm == "maxima":
        if prec != 53:
            raise ValueError, "for the maxima algorithm the precision must be 53"
        return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
    elif algorithm == "pari":
        raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"
    else:
        raise ValueError, "unknown algorithm '%s'"%algorithm

 def Fourier_coefficients( self, det_list):
     fcs = DB.find( { 'owner_id': self.__id,
                      'data_type': 'fc',
                      'det': { '$in': [ str(d) for d in det_list]}
                      })
     P = PolynomialRing( self.__field, names = 'x,y')
     loc = P.gens_dict()
     loc.update ( self.__field.gens_dict())
     return dict( (Integer(fcd['det']),
                   dict( (tuple( eval(f)), sage_eval( fcd['data'][f], locals = loc))
                         for f in fcd['data'] ))
                  for fcd in fcs)

def number_of_automorphisms__souvigner(self):
    """
    Uses the Souvigner code to compute the number of automorphisms.

    EXAMPLES::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1,1])
        sage: Q.number_of_automorphisms__souvigner()                           # optional -- souvigner
        3840
        sage: 2^5 * factorial(5)
        3840

    """
    ## Write an input text file
    F_filename = '/tmp/tmp_isom_input' + str(random()) + ".txt"
    F = open(F_filename, 'w')
    #F = tempfile.NamedTemporaryFile(prefix='tmp_auto_input', suffix=".txt")   ## This fails because the Souvigner code doesn't like random filenames (hyphens are bad...)!
    F.write("#1 \n")
    n = self.dim()
    F.write(str(n) + "x0 \n")      ## Use the lower-triangular form
    for i in range(n):
        for j in range(i+1):
            if i == j:
                F.write(str(2 * self[i,j]) + " ")
            else:
                F.write(str(self[i,j]) + " ")
        F.write("\n")
    F.flush()
    #print "Input filename = ", F.name
    #os.system("less " + F.name)

    ## Call the Souvigner automorphism code
    souvigner_auto_path = os.path.join(SAGE_LOCAL,'bin','Souvigner_AUTO')                 ## FIX THIS LATER!!!
    G1 = tempfile.NamedTemporaryFile(prefix='tmp_auto_ouput', suffix=".txt")
    #print "Output filename = ", G1.name
    os.system(souvigner_auto_path + " '" +  F.name + "' > '" + G1.name +"'")


    ## Read the output
    G2 = open(G1.name, 'r')
    for line in G2:
        if line.startswith("|Aut| = "):
            num_of_autos = sage_eval(line.lstrip("|Aut| = "))
            F.close()
            G1.close()
            G2.close()
            os.system("rm -f " + F_filename)
            #os.system("rm -f " + G1.name)
            return num_of_autos

    ## Raise and error if we're here:
    raise RuntimeError("Oops! There is a problem...")

    def _element_constructor_(self, x, y=1, coerce=True):
        """
        Construct an element of this fraction field.

        EXAMPLES::

            sage: F = QQ['x,y'].fraction_field()
            sage: F._element_constructor_(1)
            1
            sage: F._element_constructor_(F.gen(0)/F.gen(1))
            x/y
            sage: F._element_constructor_('1 + x/y')
            (x + y)/y

        ::

            sage: K = ZZ['x,y'].fraction_field()
            sage: x,y = K.gens()

        ::

            sage: F._element_constructor_(x/y)
            x/y
        """
        Element = self._element_class
        if isinstance(x, Element):
            if x.parent() is self:
                return x
            else:
                return Element(self, x.numerator(), x.denominator())
        elif isinstance(x, basestring):
            try:
                from sage.misc.sage_eval import sage_eval
                x = sage_eval(x, self.gens_dict_recursive())
                y = sage_eval(str(y), self.gens_dict_recursive())
                return Element(self, x, y, coerce=coerce, reduce=True)
            except NameError, e:
                raise TypeError, "unable to convert string"

def repr_short_to_parent(s):
    r"""
    Helper method for the growth group factory, which converts a short
    representation string to a parent.

    INPUT:

    - ``s`` -- a string, short representation of a parent.

    OUTPUT:

    A parent.

    The possible short representations are shown in the examples below.

    EXAMPLES::

        sage: from sage.rings.asymptotic.misc import repr_short_to_parent
        sage: repr_short_to_parent('ZZ')
        Integer Ring
        sage: repr_short_to_parent('QQ')
        Rational Field
        sage: repr_short_to_parent('SR')
        Symbolic Ring
        sage: repr_short_to_parent('NN')
        Non negative integer semiring

    TESTS::

        sage: repr_short_to_parent('abcdef')
        Traceback (most recent call last):
        ...
        ValueError: Cannot create a parent out of 'abcdef'.
        > *previous* NameError: name 'abcdef' is not defined
    """
    from sage.misc.sage_eval import sage_eval
    try:
        P = sage_eval(s)
    except Exception as e:
        raise combine_exceptions(
            ValueError("Cannot create a parent out of '%s'." % (s,)), e)

    from sage.misc.lazy_import import LazyImport
    if type(P) is LazyImport:
        P = P._get_object()

    from sage.structure.parent import is_Parent
    if not is_Parent(P):
        raise ValueError("'%s' does not describe a parent." % (s,))
    return P

def hypergeometric_U(alpha, beta, x, algorithm="pari", prec=53):
    r"""
    Default is a wrap of PARI's hyperu(alpha,beta,x) function.
    Optionally, algorithm = "scipy" can be used.

    The confluent hypergeometric function `y = U(a,b,x)` is
    defined to be the solution to Kummer's differential equation

    .. math::

             xy'' + (b-x)y' - ay = 0.

    This satisfies `U(a,b,x) \sim x^{-a}`, as
    `x\rightarrow \infty`, and is sometimes denoted
    ``x^{-a}2_F_0(a,1+a-b,-1/x)``. This is not the same as Kummer's
    `M`-hypergeometric function, denoted sometimes as
    ``_1F_1(alpha,beta,x)``, though it satisfies the same DE that
    `U` does.

    .. warning::

       In the literature, both are called "Kummer confluent
       hypergeometric" functions.

    EXAMPLES::

        sage: hypergeometric_U(1,1,1,"scipy")
        0.596347362323...
        sage: hypergeometric_U(1,1,1)
        0.59634736232319...
        sage: hypergeometric_U(1,1,1,"pari",70)
        0.59634736232319407434...
    """
    if algorithm == "scipy":
        if prec != 53:
            raise ValueError("for the scipy algorithm the precision must be 53")
        import scipy.special

        ans = str(scipy.special.hyperu(float(alpha), float(beta), float(x)))
        ans = ans.replace("(", "")
        ans = ans.replace(")", "")
        ans = ans.replace("j", "*I")
        return sage_eval(ans)
    elif algorithm == "pari":
        from sage.libs.pari.all import pari

        R = RealField(prec)
        return R(pari(R(alpha)).hyperu(R(beta), R(x), precision=prec))
    else:
        raise ValueError("unknown algorithm '%s'" % algorithm)