def _derivative_(self, u, m, diff_param):
        """
        EXAMPLES::

            sage: x,m = var('x,m')
            sage: elliptic_eu(x,m).diff(x)
            sqrt(-m*jacobi_sn(x, m)^2 + 1)*jacobi_dn(x, m)
            sage: elliptic_eu(x,m).diff(m)
            1/2*(elliptic_eu(x, m)
             - elliptic_f(jacobi_am(x, m), m))/m
             - 1/2*(m*jacobi_cn(x, m)*jacobi_sn(x, m)
             - (m - 1)*x
             - elliptic_eu(x, m)*jacobi_dn(x, m))*sqrt(-m*jacobi_sn(x, m)^2 + 1)/((m - 1)*m)
        """
        from sage.functions.jacobi import jacobi, jacobi_am
        if diff_param == 0:
            return (sqrt(-m * jacobi('sn', u, m) ** Integer(2) +
                         Integer(1)) * jacobi('dn', u, m))
        elif diff_param == 1:
            return (Integer(1) / Integer(2) *
                    (elliptic_eu(u, m) - elliptic_f(jacobi_am(u, m), m)) / m -
                    Integer(1) / Integer(2) * sqrt(-m * jacobi('sn', u, m) **
                    Integer(2) + Integer(1)) * (m * jacobi('sn', u, m) *
                    jacobi('cn', u, m) - (m - Integer(1)) * u -
                    elliptic_eu(u, m) * jacobi('dn', u, m)) /
                    ((m - Integer(1)) * m))

    def Weyl_law_N(self,T,T1=None):
        r"""
        The counting function for this space. N(T)=#{disc. ev.<=T}
        
        INPUT:
        
        -  ``T`` -- double


        EXAMPLES::

            sage: M=MaassWaveForms(MySubgroup(Gamma0(1))
            sage: M.Weyl_law_N(10)
            0.572841337202191
            
        """
        (c1,c2,c3,c4,c5)=self._Weyl_law_const
        cc1=RR(c1); cc2=RR(c2); cc3=RR(c3); cc4=RR(c4); cc5=RR(c5)
        #print "c1,c2,c3,c4,c5=",cc1,cc2,cc3,cc4,cc5
        t=sqrt(T*T+0.25)
        try: 
            lnt=ln(t)
        except TypeError:
            lnt=mpmath.ln(t)
        #print "t,ln(t)=",t,lnt
        NT=cc1*t*t-cc2*t*lnt+cc3*t+cc4*t+cc5
        if(T1<>None):
            t=sqrt(T1*T1+0.25)
            NT1=cc1*(T1*T1+0.25)-cc2*t*ln(t)+cc3*t+cc4*t+cc5
            return RR(abs(NT1-NT))
        else:
            return RR(NT)

    def _eval_(self, n, m, theta, phi, **kwargs):
        r"""
        TESTS::

            sage: x, y = var('x y')
            sage: spherical_harmonic(1, 2, x, y)
            0
            sage: spherical_harmonic(1, -2, x, y)
            0
            sage: spherical_harmonic(1/2, 2, x, y)
            spherical_harmonic(1/2, 2, x, y)
            sage: spherical_harmonic(3, 2, x, y)
            1/8*sqrt(30)*sqrt(7)*cos(x)*e^(2*I*y)*sin(x)^2/sqrt(pi)
            sage: spherical_harmonic(3, 2, 1, 2)
            1/8*sqrt(30)*sqrt(7)*cos(1)*e^(4*I)*sin(1)^2/sqrt(pi)
            sage: spherical_harmonic(3 + I, 2., 1, 2)
            -0.351154337307488 - 0.415562233975369*I

        Check that :trac:`20939` is fixed::

            sage: ex = spherical_harmonic(3,2,1,2*pi/3)
            sage: QQbar(ex * sqrt(pi)/cos(1)/sin(1)^2).minpoly()
            x^4 + 105/32*x^2 + 11025/1024
        """
        if n in ZZ and m in ZZ and n > -1:
            if abs(m) > n:
                return ZZ(0)
            if m == 0 and theta.is_zero():
                return sqrt((2*n+1)/4/pi)
            from sage.arith.misc import factorial
            from sage.functions.trig import cos
            from sage.functions.orthogonal_polys import gen_legendre_P
            return (sqrt(factorial(n-m) * (2*n+1) / (4*pi * factorial(n+m))) *
                    exp(I*m*phi) * gen_legendre_P(n, m, cos(theta)) *
                    (-1)**m).simplify_trig()

def plot_slope_field(f, xrange, yrange, **kwds):
    r"""
    ``plot_slope_field`` takes a function of two variables xvar and yvar
    (for instance, if the variables are `x` and `y`, take `f(x,y)`), and at
    representative points `(x_i,y_i)` between xmin, xmax, and ymin, ymax
    respectively, plots a line with slope `f(x_i,y_i)` (see below).

    ``plot_slope_field(f, (xvar, xmin, xmax), (yvar, ymin, ymax))``

    EXAMPLES:

    A logistic function modeling population growth::

        sage: x,y = var('x y')
        sage: capacity = 3 # thousand
        sage: growth_rate = 0.7 # population increases by 70% per unit of time
        sage: plot_slope_field(growth_rate*(1-y/capacity)*y, (x,0,5), (y,0,capacity*2))
        Graphics object consisting of 1 graphics primitive

    Plot a slope field involving sin and cos::

        sage: x,y = var('x y')
        sage: plot_slope_field(sin(x+y)+cos(x+y), (x,-3,3), (y,-3,3))
        Graphics object consisting of 1 graphics primitive

    Plot a slope field using a lambda function::

        sage: plot_slope_field(lambda x,y: x+y, (-2,2), (-2,2))
        Graphics object consisting of 1 graphics primitive

    TESTS:

    Verify that we're not getting warnings due to use of headless quivers
    (trac #11208)::

        sage: x,y = var('x y')
        sage: import numpy # bump warnings up to errors for testing purposes
        sage: old_err = numpy.seterr('raise')
        sage: plot_slope_field(sin(x+y)+cos(x+y), (x,-3,3), (y,-3,3))
        Graphics object consisting of 1 graphics primitive
        sage: dummy_err = numpy.seterr(**old_err)
    """
    slope_options = {'headaxislength': 0, 'headlength': 1e-9, 'pivot': 'middle'}
    slope_options.update(kwds)

    from sage.functions.all import sqrt
    from inspect import isfunction
    if isfunction(f):
        norm_inverse=lambda x,y: 1/sqrt(f(x,y)**2+1)
        f_normalized=lambda x,y: f(x,y)*norm_inverse(x,y)
    else:
        norm_inverse = 1/sqrt((f**2+1))
        f_normalized=f*norm_inverse
    return plot_vector_field((norm_inverse, f_normalized), xrange, yrange, **slope_options)

def is_triangular_number(n):
    """
    Determines if the integer n is a triangular number.
    (I.e. determine if n = a*(a+1)/2 for some natural number a.)
    If so, return the number a, otherwise return False.

    Note: As a convention, n=0 is considered triangular for the
    number a=0 only (and not for a=-1).

    WARNING: Any non-zero value will return True, so this will test as
    True iff n is triangular and not zero.  If n is zero, then this
    will return the integer zero, which tests as False, so one must test

        if is_triangular_number(n) != False:

    instead of

        if is_triangular_number(n):

    to get zero to appear triangular.

    INPUT:
        an integer

    OUTPUT:
        either False or a non-negative integer

    EXAMPLES:
        sage: is_triangular_number(3)
        2
        sage: is_triangular_number(1)
        1
        sage: is_triangular_number(2)
        False
        sage: is_triangular_number(0)
        0
        sage: is_triangular_number(-1)
        False
        sage: is_triangular_number(-11)
        False
        sage: is_triangular_number(-1000)
        False
        sage: is_triangular_number(-0)
        0
        sage: is_triangular_number(10^6 * (10^6 +1)/2)
        1000000
    """
    if n < 0:
        return False
    elif n == 0:
        return ZZ(0)
    else:
        from sage.functions.all import sqrt
        ## Try to solve for the integer a
        try:
            disc_sqrt = ZZ(sqrt(1+8*n))
            a = ZZ( (ZZ(-1) + disc_sqrt) / ZZ(2) )
            return a
        except StandardError:
            return False

    def _2x2_matrix_entries(self, beta):
        r"""
        Young's representations are constructed by combining
        `2\times2`-matrices that depend on ``beta`` For the orthogonal
        representation, this is the following matrix::

            ``[     -beta       sqrt(1-beta^2) ]``
            ``[ sqrt(1-beta^2)       beta      ]``

        EXAMPLES::

            sage: from sage.combinat.symmetric_group_representations import YoungRepresentation_Orthogonal
            sage: orth = YoungRepresentation_Orthogonal([2,1])
            sage: orth._2x2_matrix_entries(1/2)
            (-1/2, 1/2*sqrt(3), 1/2*sqrt(3), 1/2)
        """
        return (-beta, sqrt(1-beta**2), sqrt(1-beta**2), beta)

    def set_params(lam, k):
        n = pow(2, ceil(log(lam**2 * k)/log(2))) # dim of poly ring, closest power of 2 to k(lam^2)
        q = next_prime(ZZ(2)**(8*k*lam) * n**k, proof=False) # prime modulus

        sigma = int(sqrt(lam * n))
        sigma_prime = lam * int(n**(1.5))

        return (n, q, sigma, sigma_prime, k)

    def _derivative_(self, x, diff_param=None):
        """
        Derivative of inverse erf function.

        EXAMPLES::

            sage: erfinv(x).diff(x)
            1/2*sqrt(pi)*e^(erfinv(x)^2)
        """
        return sqrt(pi)*exp(erfinv(x)**2)/2

    def __lalg__(self,D):
        r"""
        For positive `D`, this function evaluates the quotient
        `L(E_D,1)\cdot \sqrt(D)/\Omega_E` where `E_D` is the twist of
        `E` by `D`, `\Omega_E` is the least positive period of `E`.
        For negative `E`, it is the quotient
        `L(E_D,1)\cdot \sqrt(-D)/\Omega^{-}_E`
        where `\Omega^{-}_E` is the least positive imaginary part of a
        non-real period of `E`.

        EXAMPLES::

            sage: E = EllipticCurve('11a1')
            sage: m = E.modular_symbol(sign=+1, implementation='sage')
            sage: m.__lalg__(1)
            1/5
            sage: m.__lalg__(3)
            5/2

        """
        from sage.functions.all import sqrt
        # the computation of the L-value could take a lot of time,
        # but then the conductor is so large
        # that the computation of modular symbols for E took even longer

        E = self._E
        ED = E.quadratic_twist(D)
        lv = ED.lseries().L_ratio() # this is L(ED,1) divided by the Néron period omD of ED
        lv *= ED.real_components() # now it is by the least positive period
        omD = ED.period_lattice().basis()[0]
        if D > 0 :
            om = E.period_lattice().basis()[0]
            q = sqrt(D)*omD/om * 8
        else :
            om = E.period_lattice().basis()[1].imag()
            if E.real_components() == 1:
                om *= 2
            q = sqrt(-D)*omD/om*8

        # see padic_lseries.pAdicLeries._quotient_of_periods_to_twist
        # for the explanation of the second factor
        verbose('real approximation is %s'%q)
        return lv/8 * QQ(int(round(q)))

def mrrw1_bound_asymp(delta,q):
    """
    Computes the first asymptotic McEliese-Rumsey-Rodemich-Welsh bound
    for the information rate, provided `0 < \delta < 1-1/q`.
    
    EXAMPLES::
    
        sage: mrrw1_bound_asymp(1/4,2)
        0.354578902665
    """
    return RDF(entropy((q-1-delta*(q-2)-2*sqrt((q-1)*delta*(1-delta)))/q,q))

def elias_bound_asymp(delta,q):
    """
    Computes the asymptotic Elias bound for the information rate,
    provided `0 < \delta 1-1/q`.
    
    EXAMPLES::
    
        sage: elias_bound_asymp(1/4,2)
        0.39912396330...
    """
    r = 1-1/q
    return RDF((1-entropy(r-sqrt(r*(r-delta)), q)))

def mrrw1_bound_asymp(delta,q):
    """
    The first asymptotic McEliese-Rumsey-Rodemich-Welsh bound.

    This only makes sense when `0 < \delta < 1-1/q`.

    EXAMPLES::

        sage: codes.bounds.mrrw1_bound_asymp(1/4,2)   # abs tol 4e-16
        0.3545789026652697
    """
    return RDF(entropy((q-1-delta*(q-2)-2*sqrt((q-1)*delta*(1-delta)))/q,q))

    def __init__(self, params, asym=False):

        (self.n, self.q, sigma, self.sigma_prime, self.k) = params

        S, x = PolynomialRing(ZZ, 'x').objgen()
        self.R = S.quotient_ring(S.ideal(x**self.n + 1))

        Sq = PolynomialRing(Zmod(self.q), 'x')
        self.Rq = Sq.quotient_ring(Sq.ideal(x**self.n + 1))

        # draw z_is uniformly from Rq and compute its inverse in Rq
        if asym:
            z = [self.Rq.random_element() for i in range(self.k)]
            self.zinv = [z_i**(-1) for z_i in z]
        else: # or do symmetric version
            z = self.Rq.random_element()
            zinv = z**(-1)
            z, self.zinv = zip(*[(z,zinv) for i in range(self.k)])

        # set up some discrete Gaussians
        DGSL_sigma = DGSL(ZZ**self.n, sigma)
        self.D_sigma = lambda: self.Rq(list(DGSL_sigma()))

        # discrete Gaussian in ZZ^n with stddev sigma_prime, yields random level-0 encodings
        DGSL_sigmap_ZZ = DGSL(ZZ**self.n, self.sigma_prime)
        self.D_sigmap_ZZ = lambda: self.Rq(list(DGSL_sigmap_ZZ()))

        # draw g repeatedly from a Gaussian distribution of Z^n (with param sigma)
        # until g^(-1) in QQ[x]/<x^n + 1> is small (< n^2)
        Sk = PolynomialRing(QQ, 'x')
        K = Sk.quotient_ring(Sk.ideal(x**self.n + 1)) 
        while True:
            l = self.D_sigma()
            ginv_K = K(mod_near_poly(l, self.q))**(-1)
            ginv_size = vector(ginv_K).norm()

            if ginv_size < self.n**2:
                g = self.Rq(l)
                self.ginv = g**(-1)
                break

        # discrete Gaussian in I = <g>, yields random encodings of 0
        short_g = vector(ZZ, mod_near_poly(g,self.q))
        DGSL_sigmap_I = DGSL(short_g, self.sigma_prime)
        self.D_sigmap_I = lambda: self.Rq(list(DGSL_sigmap_I()))

        # compute zero-testing parameter p_zt
        # randomly draw h (in Rq) from a discrete Gaussian with param q^(1/2)
        self.h = self.Rq(list(DGSL(ZZ**self.n, round(sqrt(self.q)))()))

        # create p_zt
        self.p_zt = self.ginv * self.h * prod(z)

def plot_slope_field(f, xrange, yrange, **kwds):
    r"""
    ``plot_slope_field`` takes a function of two variables xvar and yvar
    (for instance, if the variables are `x` and `y`, take `f(x,y)`), and at
    representative points `(x_i,y_i)` between xmin, xmax, and ymin, ymax
    respectively, plots a line with slope `f(x_i,y_i)` (see below).

    ``plot_slope_field(f, (xvar, xmin, xmax), (yvar, ymin, ymax))``

    EXAMPLES: 

    A logistic function modeling population growth::

        sage: x,y = var('x y')
        sage: capacity = 3 # thousand
        sage: growth_rate = 0.7 # population increases by 70% per unit of time
        sage: plot_slope_field(growth_rate*(1-y/capacity)*y, (x,0,5), (y,0,capacity*2))

    Plot a slope field involving sin and cos::

        sage: x,y = var('x y')
        sage: plot_slope_field(sin(x+y)+cos(x+y), (x,-3,3), (y,-3,3))

    Plot a slope field using a lambda function::

        sage: plot_slope_field(lambda x,y: x+y, (-2,2), (-2,2))
    """
    slope_options = {'headaxislength': 0, 'headlength': 0, 'pivot': 'middle'}
    slope_options.update(kwds)

    from sage.functions.all import sqrt
    from inspect import isfunction
    if isfunction(f):
        norm_inverse=lambda x,y: 1/sqrt(f(x,y)**2+1)
        f_normalized=lambda x,y: f(x,y)*norm_inverse(x,y)
    else:
        norm_inverse = 1/sqrt((f**2+1))
        f_normalized=f*norm_inverse
    return plot_vector_field((norm_inverse, f_normalized), xrange, yrange, **slope_options)

def elias_bound_asymp(delta,q):
    """
    The asymptotic Elias bound for the information rate.

    This only makes sense when `0 < \delta < 1-1/q`.

    EXAMPLES::

        sage: codes.bounds.elias_bound_asymp(1/4,2)
        0.39912396330...
    """
    r = 1-1/q
    return RDF((1-entropy(r-sqrt(r*(r-delta)), q)))

def gen_legendre_Q(n,m,x):
    """
    Returns the generalized (or associated) Legendre function of the
    second kind for integers `n>-1`, `m>-1`.

    Maxima restricts m = n. Hence the cases m n are computed using the
    same recursion used for gen_legendre_P(n,m,x) when m is odd and
    1.

    EXAMPLES::

        sage: P.<t> = QQ[]
        sage: gen_legendre_Q(2,0,t)
        3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1))
        sage: gen_legendre_Q(2,0,t) - legendre_Q(2, t)
        0
        sage: gen_legendre_Q(3,1,0.5)
        2.49185259170895
        sage: gen_legendre_Q(0, 1, x)
        -1/sqrt(-x^2 + 1)
        sage: gen_legendre_Q(2, 4, x).factor()
        48*x/((x - 1)^2*(x + 1)^2)
    """
    from sage.functions.all import sqrt
    if m <= n:
        _init()
        return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x})
    if m == n + 1 or n == 0:
        if m.mod(2).is_zero():
            denom = (1 - x**2)**(m/2)
        else:
            denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2)
        if m == n + 1:
            return (-1)**m*(m-1).factorial()*2**n/denom
        else:
            return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom)
    else:
        return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-(n+m-1)*gen_legendre_Q(n-1,m-1,x))/sqrt(1-x**2)

    def standard_deviation(self):
        r"""
        The standard deviation of the discrete random variable.

        Let `S` be the probability space of `X` = self,
        with probability function `p`, and `E(X)` be the
        expectation of `X`. Then the standard deviation of
        `X` is defined to be

        .. MATH::

                     \sigma(X) = \sqrt{ \sum_{x \in S} p(x) (X(x) - E(x))^2}
        """
        return sqrt(self.variance())

 def _quotient_of_periods_to_twist(self,D):
     r"""
     For a fundamental discriminant `D` of a quadratic number field
     this computes the constant `\eta` such that
     `\sqrt{D}\cdot\Omega_{E_D}^{+} =\eta\cdot
     \Omega_E^{sign(D)}`. As in [MTT]_ page 40.  This is either 1
     or 2 unless the condition on the twist is not satisfied,
     e.g. if we are 'twisting back' to a semi-stable curve.
     
     REFERENCES:
     
     - [MTT] B. Mazur, J. Tate, and J. Teitelbaum,
       On `p`-adic analogues of the conjectures of Birch and 
       Swinnerton-Dyer, Invertiones mathematicae 84, (1986), 1-48.        
     
     .. note: No check on precision is made, so this may fail for huge `D`.
     
     EXAMPLES::
     
     """
     from sage.functions.all import sqrt
     # This funciton does not depend on p and could be moved out of this file but it is needed only here
     
     # Note that the number of real components does not change by twisting.
     if D == 1:
         return 1
     if D > 1:
         Et = self._E.quadratic_twist(D)
         qt = Et.period_lattice().basis()[0]/self._E.period_lattice().basis()[0]
         qt *= sqrt(qt.parent()(D))
     else:
         Et = self._E.quadratic_twist(D)
         qt = Et.period_lattice().basis()[0]/self._E.period_lattice().basis()[1].imag()
         qt *= sqrt(qt.parent()(-D))
     verbose('the real approximation is %s'%qt)
     # we know from MTT that the result has a denominator 1
     return QQ(int(round(8*qt)))/8

    def translation_standard_deviation(self, map):
        r"""
        The standard deviation of the translated discrete random variable
        `X \circ e`, where `X` = self and `e` =
        map.

        Let `S` be the probability space of `X` = self,
        with probability function `p`, and `E(X)` be the
        expectation of `X`. Then the standard deviation of
        `X` is defined to be

        .. MATH::

                     \sigma(X) = \sqrt{ \sum_{x \in S} p(x) (X(x) - E(x))^2}
        """
        return sqrt(self.translation_variance(map))

    def _derivative_(self, n, m, theta, phi, diff_param):
        r"""
        TESTS::

            sage: n, m, theta, phi = var('n m theta phi')
            sage: spherical_harmonic(n, m, theta, phi).diff(theta)
            m*cot(theta)*spherical_harmonic(n, m, theta, phi)
             + sqrt(-(m + n + 1)*(m - n))*e^(-I*phi)*spherical_harmonic(n, m + 1, theta, phi)
            sage: spherical_harmonic(n, m, theta, phi).diff(phi)
            I*m*spherical_harmonic(n, m, theta, phi)
        """
        if diff_param == 2:
            return m * cot(theta) * spherical_harmonic(n, m, theta, phi) + sqrt((n - m) * (n + m + 1)) * exp(
                -I * phi
            ) * spherical_harmonic(n, m + 1, theta, phi)
        if diff_param == 3:
            return I * m * spherical_harmonic(n, m, theta, phi)

        raise ValueError("only derivative with respect to theta or phi" " supported")