def _eval_(self, n, x):
        """
        EXAMPLES::

            sage: y=var('y')
            sage: bessel_I(y,x)
            bessel_I(y, x)
            sage: bessel_I(0.0, 1.0)
            1.26606587775201
            sage: bessel_I(1/2, 1)
            sqrt(2)*sinh(1)/sqrt(pi)
            sage: bessel_I(-1/2, pi)
            sqrt(2)*cosh(pi)/pi
        """
        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
                (is_inexact(n) or is_inexact(x))):
            coercion_model = get_coercion_model()
            n, x = coercion_model.canonical_coercion(n, x)
            return self._evalf_(n, x, parent(n))

        # special identities
        if n == Integer(1) / Integer(2):
            return sqrt(2 / (pi * x)) * sinh(x)
        elif n == -Integer(1) / Integer(2):
            return sqrt(2 / (pi * x)) * cosh(x)

        return None  # leaves the expression unevaluated

 def get_s(self):
     """Returns a single value of s as required for the Elligator map while there may even be two of them, in which case the returned values is the smaller of the two
     """ 
     try:
         return self._s
     except AttributeError:
         pass
     s = []
     d = self.__E.get_d()
     if is_square(-d):
         c = ( 2*(d-1) + 4*sqrt(-d) ) / ( 2*(d+1) )
         s_2 = 2/c
         if is_square(s_2):
             s.append(sqrt(s_2))
         
         c = ( 2*(d-1) - 4*sqrt(-d) ) / ( 2*(d+1) )
         s_2 = 2/c
         if is_square(s_2):
             s.append(sqrt(s_2))
     if len(s) == 0:
         return None
     elif len(s) == 1:
         return s[0]
     else:
         if s[0] < s[1]:
             return s[0]
         else: 
             return s[1]    

def platonic_dodecahedron():
    r"""Produce a triple consisting of a polyhedral version of the platonic dodecahedron,
    the associated cone surface, and a ConeSurfaceToPolyhedronMap from the surface
    to the polyhedron.

    EXAMPLES::

    sage: from flatsurf.geometry.polyhedra import platonic_dodecahedron
    sage: polyhedron,surface,surface_to_polyhedron = platonic_dodecahedron()
    sage: TestSuite(surface).run()
    r"""
    vertices=[]
    phi=AA(1+sqrt(5))/2
    F=NumberField(phi.minpoly(),"phi",embedding=phi)
    phi=F.gen()
    for x in xrange(-1,3,2):
        for y in xrange(-1,3,2):
            for z in xrange(-1,3,2):
                vertices.append(vector(F,(x,y,z)))
    for x in xrange(-1,3,2):
        for y in xrange(-1,3,2):
            vertices.append(vector(F,(0,x*phi,y/phi)))
            vertices.append(vector(F,(y/phi,0,x*phi)))
            vertices.append(vector(F,(x*phi,y/phi,0)))
    scale=AA(2/sqrt(1+(phi-1)**2+(1/phi-1)**2))
    p=Polyhedron(vertices=vertices)
    s,m = polyhedron_to_cone_surface(p,scaling_factor=scale)
    return p,s,m

    def _eval_(self, a, z):
        """
        EXAMPLES::

            sage: struve_H(0,0)
            0
            sage: struve_H(pi,0)
            0
            sage: struve_H(-1/2,x)
            sqrt(2)*sqrt(1/(pi*x))*sin(x)
            sage: struve_H(1/2,-1)
            -sqrt(2)*sqrt(-1/pi)*(cos(1) - 1)
            sage: struve_H(1/2,pi)
            2*sqrt(2)/pi
            sage: struve_H(2,x)
            struve_H(2, x)
            sage: struve_H(-3/2,x)
            -bessel_J(3/2, x)
        """
        from sage.symbolic.ring import SR
        if z.is_zero() \
                and (SR(a).is_numeric() or SR(a).is_constant()) \
                and a.real() >= -1:
                return ZZ(0)
        if a == -Integer(1)/2:
            from sage.functions.trig import sin
            return sqrt(2/(pi*z)) * sin(z)
        if a == Integer(1)/2:
            from sage.functions.trig import cos
            return sqrt(2/(pi*z)) * (1-cos(z))
        if a < 0 and not SR(a).is_integer() and SR(2*a).is_integer():
            from sage.rings.rational_field import QQ
            n = (a*(-2) - 1)/2
            return Integer(-1)**n * bessel_J(n+QQ(1)/2, z)

    def __eq__(self, other):
        r"""
        Return ``True`` if the isometries are the same and ``False`` otherwise.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: A = UHP.get_isometry(identity_matrix(2))
            sage: B = UHP.get_isometry(-identity_matrix(2))
            sage: A == B
            True

            sage: HM = HyperbolicPlane().HM()
            sage: A = HM.random_isometry()
            sage: A == A
            True
        """
        if not isinstance(other, HyperbolicIsometry):
            return False
        test_matrix = bool((self.matrix() - other.matrix()).norm() < EPSILON)
        if self.domain().is_isometry_group_projective():
            A,B = self.matrix(), other.matrix() # Rename for simplicity
            m = self.matrix().ncols()
            A = A / sqrt(A.det(), m) # Normalized to have determinant 1
            B = B / sqrt(B.det(), m)
            test_matrix = ((A - B).norm() < EPSILON
                           or (A + B).norm() < EPSILON)
        return self.domain() is other.domain() and test_matrix

def SO21_to_SL2R(M):
    r"""
    A homomorphism from `SO(2, 1)` to `SL(2, \RR)`.

    Note that this is not the only homomorphism, but it is the only one
    that works in the context of the implemented 2D hyperbolic geometry
    models.

    EXAMPLES::

        sage: from sage.geometry.hyperbolic_space.hyperbolic_coercion import SO21_to_SL2R
        sage: (SO21_to_SL2R(identity_matrix(3)) - identity_matrix(2)).norm() < 10**-4
        True
    """
    ####################################################################
    # SL(2,R) is the double cover of SO (2,1)^+, so we need to choose  #
    # a lift.  I have formulas for the absolute values of each entry   #
    # a,b ,c,d of the lift matrix(2,[a,b,c,d]), but we need to choose  #
    # one entry to be positive.  I choose d for no particular reason,  #
    # unless d = 0, then we choose c > 0.  The basic strategy for this #
    # function is to find the linear map induced by the SO(2,1)        #
    # element on the Lie algebra sl(2, R).  This corresponds to the    #
    # Adjoint action by a matrix A or -A in SL(2,R).  To find which    #
    # matrix let X,Y,Z be a basis for sl(2,R) and look at the images   #
    # of X,Y,Z as well as the second and third standard basis vectors  #
    # for 2x2 matrices (these are traceless, so are in the Lie         #
    # algebra).  These corresponds to AXA^-1 etc and give formulas     #
    # for the entries of A.                                            #
    ####################################################################
    (m_1,m_2,m_3,m_4,m_5,m_6,m_7,m_8,m_9) = M.list()
    d = sqrt(Integer(1)/Integer(2)*m_5 - Integer(1)/Integer(2)*m_6 -
             Integer(1)/Integer(2)*m_8 + Integer(1)/Integer(2)*m_9)
    if M.det() > 0:  # EPSILON?
        det_sign = 1
    elif M.det() < 0:  # EPSILON?
        det_sign = -1
    if d > 0:  # EPSILON?
        c = (-Integer(1)/Integer(2)*m_4 + Integer(1)/Integer(2)*m_7)/d
        b = (-Integer(1)/Integer(2)*m_2 + Integer(1)/Integer(2)*m_3)/d
        ad = det_sign*1 + b*c  # ad - bc = pm 1
        a = ad/d
    else:  # d is 0, so we make c > 0
        c = sqrt(-Integer(1)/Integer(2)*m_5 - Integer(1)/Integer(2)*m_6 +
                 Integer(1)/Integer(2)*m_8 + Integer(1)/Integer(2)*m_9)
        d = (-Integer(1)/Integer(2)*m_4 + Integer(1)/Integer(2)*m_7)/c
            # d = 0, so ad - bc = -bc = pm 1.
        b = - (det_sign*1)/c
        a = (Integer(1)/Integer(2)*m_4 + Integer(1)/Integer(2)*m_7)/b
    A = matrix(2, [a, b, c, d])
    return A

 def _0f1(b, z):
     F12 = cosh(2 * sqrt(z))
     F32 = sinh(2 * sqrt(z)) / (2 * sqrt(z))
     if 2 * b == 1:
         return F12
     if 2 * b == 3:
         return F32
     if 2 * b > 3:
         return ((b - 2) * (b - 1) / z * (_0f1(b - 2, z) -
                 _0f1(b - 1, z)))
     if 2 * b < 1:
         return (_0f1(b + 1, z) + z / (b * (b + 1)) *
                 _0f1(b + 2, z))
     raise ValueError

def RankingScale_USAU_2013():
    r"""
    EXAMPLES::

        sage: from slabbe.ranking_scale import RankingScale_USAU_2013
        sage: R = RankingScale_USAU_2013()
        sage: R.table()
          Position   Serie 1500   Serie 1000   Serie 500   Serie 250
          1          1500         1000         500         250
          2          1366         911          455         228
          3          1252         835          417         209
          4          1152         768          384         192
          5          1061         708          354         177
          6          979          653          326         163
          7          903          602          301         151
          8          832          555          278         139
          9          767          511          256         128
          10         706          471          236         118
          11         648          432          217         109
          12         595          397          199         100
          13         544          363          182         91
          14         497          332          166         83
          15         452          302          151         76
          16         410          274          137         69
          17         371          248          124         62
          18         334          223          112         56
          19         299          199          100         50
          20         266          177          89          45
          21         235          157          79          40
          22         206          137          69          35
          23         178          119          60          30
          24         152          102          51          26
          25         128          86           43          22
          26         106          71           36          18
          27         85           57           29          15
          28         66           44           22          12
          29         47           32           16          9
          30         31           21           11          6
          31         15           10           6           3
          32         1            1            1           1
    """
    L1500 = [0] + discrete_curve(32, 1500, K=1, R=sqrt(2), base=e)
    L1000 = [0] + discrete_curve(32, 1000, K=1, R=sqrt(2), base=e)
    L500 = [0] + discrete_curve(32, 500, K=1, R=sqrt(2), base=e)
    L250 = [0] + discrete_curve(32, 250, K=1, R=sqrt(2), base=e)

    scales = L1500, L1000, L500, L250
    names = ['Serie 1500', 'Serie 1000', 'Serie 500', 'Serie 250']
    return RankingScale(names, scales)

def platonic_icosahedron():
    r"""Produce a triple consisting of a polyhedral version of the platonic icosahedron,
    the associated cone surface, and a ConeSurfaceToPolyhedronMap from the surface
    to the polyhedron.

    EXAMPLES::

    sage: from flatsurf.geometry.polyhedra import platonic_icosahedron
    sage: polyhedron,surface,surface_to_polyhedron = platonic_icosahedron()
    sage: TestSuite(surface).run()
    r"""
    vertices=[]
    phi=AA(1+sqrt(5))/2
    F=NumberField(phi.minpoly(),"phi",embedding=phi)
    phi=F.gen()
    for i in xrange(3):
        for s1 in xrange(-1,3,2):
            for s2 in xrange(-1,3,2):
                p=3*[None]
                p[i]=s1*phi
                p[(i+1)%3]=s2
                p[(i+2)%3]=0
                vertices.append(vector(F,p))
    p=Polyhedron(vertices=vertices)
    
    s,m = polyhedron_to_cone_surface(p)
    return p,s,m

    def _test_representation(self, **options):
        """
        Check (on some elements) that ``self`` is a representation of the
        given semigroup.

        EXAMPLES::

            sage: G = groups.permutation.Dihedral(4)
            sage: R = G.regular_representation()
            sage: R._test_representation()

            sage: G = CoxeterGroup(['A',4,1], base_ring=ZZ)
            sage: M = CombinatorialFreeModule(QQ, ['v'])
            sage: from sage.modules.with_basis.representation import Representation
            sage: on_basis = lambda g,m: M.term(m, (-1)**g.length())
            sage: R = Representation(G, M, on_basis, side="right")
            sage: R._test_representation(max_runs=500)
        """
        from sage.functions.other import sqrt
        tester = self._tester(**options)
        S = tester.some_elements()
        L = []
        max_len = int(sqrt(tester._max_runs)) + 1
        for i,x in enumerate(self._semigroup):
            L.append(x)
            if i >= max_len:
                break
        for x in L:
            for y in L:
                for elt in S:
                    if self._left_repr:
                        tester.assertEqual(x*(y*elt), (x*y)*elt)
                    else:
                        tester.assertEqual((elt*y)*x, elt*(y*x))

 def random_point(self):
     while(1):
         x = self.base_ring().random_element()
         test = (x**2 - 1)/(self._d*(x**2) - 1)
         if is_square(test):
             y = sqrt(test)
             return self([x,y,1])    

def solve_degree2_to_integer_range(a,b,c):
    r"""
    Returns the greatest integer range `[i_1, i_2]` such that
    `i_1 > x_1` and `i_2 < x_2` where `x_1, x_2` are the two zeroes of the equation in `x`:
    `ax^2+bx+c=0`.

    If there is no real solution to the equation, it returns an empty range with negative coefficients.

    INPUT:

    - ``a``, ``b`` and ``c`` -- coefficients of a second degree equation, ``a`` being the coefficient of
      the higher degree term.

    EXAMPLES::

        sage: from sage.coding.guruswami_sudan.utils import solve_degree2_to_integer_range
        sage: solve_degree2_to_integer_range(1, -5, 1)
        (1, 4)

    If there is no real solution::

        sage: solve_degree2_to_integer_range(50, 5, 42)
        (-2, -1)
    """
    D = b**2 - 4*a*c
    if D < 0:
        return (-2,-1)
    sD = float(sqrt(D))
    minx, maxx = (-b-sD)/2.0/a , (-b+sD)/2.0/a
    mini, maxi = (ligt(minx), gilt(maxx))
    if mini > maxi:
        return (-2,-1)
    else:
        return (mini,maxi)

 def cmp_ir(self,z):
     """ 
     returns -1 for left, 0 for in, and 1 for right from initial region
     cut line is on the north ray from L.
     """
     L = self.L
     x0 = self.x0
     if x0 > 0.5:
         if real(z) > real(L) and abs(z) < abs(L):
             return 0
         if real(z) < real(L):
             return -1
         if real(z) > real(L):
             return 1
     else:
         if imag(z) > imag(L):
             if real(z) > real(L):
                 return 1
             if real(z) < real(L):
                 return -1
         if real(z) < real(L) and real(z) > log(real(L)) + log(sqrt(1+tan(imag(z))**2)):
             return 0
         if real(z) > real(L):
             return 1
         if real(z) < real(L):
             return -1

    def cmp_ir(self,z):
        """ 
        returns -1 for left, 0 for in, and 1 for right from initial region
        cut line is on the north ray from pfp.

        Works only for real x0.
        """
        pfp = self.pfp
        x0 = self.x0
        if x0 > 0.5:
            print z,abs(z)
            if real(z) >= real(pfp) and abs(z) < abs(pfp):
                return 0
            if real(z) < real(pfp):
                return -1
            if real(z) > real(pfp):
                return 1
        else:
            if imag(z) > imag(pfp):
                if real(z) > real(pfp):
                    return 1
                if real(z) < real(pfp):
                    return -1
            if real(z) < real(pfp) and real(z) > log(real(pfp)) + log(sqrt(1+tan(imag(z))**2)):
                return 0
            if real(z) > real(pfp):
                return 1
            if real(z) < real(pfp):
                return -1

    def iter_positive_forms_with_content(self) :
        if self.__disc is infinity :
            raise ValueError, "infinity is not a true filter index"
        
        
        if self.__reduced :        
            for a in xrange(1,isqrt(self.__disc // 3) + 1) :
                for b in xrange(a+1) :
                    g = gcd(a, b)
                    for c in xrange(a, (b**2 + (self.__disc - 1))//(4*a) + 1) :
                        yield (a,b,c), gcd(g,c)
        else :
            maxtrace = floor(5*self.__disc / 15 + sqrt(self.__disc)/2)
            for a in xrange(1, maxtrace + 1) :
                for c in xrange(1, maxtrace - a + 1) :
                    g = gcd(a,c)
                    
                    Bu = isqrt(4*a*c - 1)

                    di = 4*a*c - self.__disc
                    if di >= 0 :
                        Bl = isqrt(di) + 1 
                    else :
                        Bl = 0
                    
                    for b in xrange(-Bu, -Bl + 1) :
                        yield (a,b,c), gcd(g,b)
                    for b in xrange(Bl, Bu + 1) :
                        yield (a,b,c), gcd(g,b)
        #! if self.__reduced

        raise StopIteration

    def show(self, show_hyperboloid=True, **graphics_options):
        r"""
        Plot ``self``.

        EXAMPLES::

            sage: from sage.geometry.hyperbolic_space.hyperbolic_geodesic import *
            sage: g = HyperbolicPlane().HM().random_geodesic()
            sage: g.show()
            Graphics3d Object
        """
        x = SR.var('x')
        opts = self.graphics_options()
        opts.update(graphics_options)
        v1, u2 = [vector(k.coordinates()) for k in self.endpoints()]
        # Lorentzian Gram Shmidt.  The original vectors will be
        # u1, u2 and the orthogonal ones will be v1, v2.  Except
        # v1 = u1, and I don't want to declare another variable,
        # hence the odd naming convention above.
        # We need the Lorentz dot product of v1 and u2.
        v1_ldot_u2 = u2[0]*v1[0] + u2[1]*v1[1] - u2[2]*v1[2]
        v2 = u2 + v1_ldot_u2 * v1
        v2_norm = sqrt(v2[0]**2 + v2[1]**2 - v2[2]**2)
        v2 = v2 / v2_norm
        v2_ldot_u2 = u2[0]*v2[0] + u2[1]*v2[1] - u2[2]*v2[2]
        # Now v1 and v2 are Lorentz orthogonal, and |v1| = -1, |v2|=1
        # That is, v1 is unit timelike and v2 is unit spacelike.
        # This means that cosh(x)*v1 + sinh(x)*v2 is unit timelike.
        hyperbola = cosh(x)*v1 + sinh(x)*v2
        endtime = arcsinh(v2_ldot_u2)
        from sage.plot.plot3d.all import parametric_plot3d
        pic = parametric_plot3d(hyperbola, (x, 0, endtime), **graphics_options)
        if show_hyperboloid:
            pic += self._model.get_background_graphic()
        return pic

    def perpendicular_bisector(self): #UHP
        r"""
        Return the perpendicular bisector of the hyperbolic geodesic ``self``
        if that geodesic has finite length.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: g = UHP.random_geodesic()
            sage: h = g.perpendicular_bisector()
            sage: c = lambda x: x.coordinates()
            sage: bool(c(g.intersection(h)[0]) - c(g.midpoint()) < 10**-9)
            True

        Infinite geodesics cannot be bisected::

            sage: UHP.get_geodesic(0, 1).perpendicular_bisector()
            Traceback (most recent call last):
            ...
            ValueError: the length must be finite
        """
        if self.length() == infinity:
            raise ValueError("the length must be finite")
        start = self._start.coordinates()
        d = self._model._dist_points(start, self._end.coordinates()) / 2
        S = self.complete()._to_std_geod(start)
        T1 = matrix([[exp(d/2), 0], [0, exp(-d/2)]])
        s2 = sqrt(2) * 0.5
        T2 = matrix([[s2, -s2], [s2, s2]])
        isom_mtrx = S.inverse() * (T1 * T2) * S # We need to clean this matrix up.
        if (isom_mtrx - isom_mtrx.conjugate()).norm() < 5*EPSILON: # Imaginary part is small.
            isom_mtrx = (isom_mtrx + isom_mtrx.conjugate()) / 2 # Set it to its real part.
        H = self._model.get_isometry(isom_mtrx)
        return self._model.get_geodesic(H(self._start), H(self._end))

def RankingScale_CQU4_2011():
    r"""
    EXEMPLES::
        
        sage: from slabbe.ranking_scale import RankingScale_CQU4_2011
        sage: R = RankingScale_CQU4_2011()
        sage: R.table()
          Position   Grand Chelem   Mars Attaque   La Flotte   Petit Chelem
          1          1000           1000           800         400
          2          938            938            728         355
          3          884            884            668         317
          4          835            835            614         285
          5          791            791            566         256
          6          750            750            522         230
          7          711            711            482         206
          8          675            675            444         184
          9          641            641            409         164
          10         609            609            377         146
        ...
          95         0              11             0           0
          96         0              10             0           0
          97         0              9              0           0
          98         0              8              0           0
          99         0              7              0           0
          100        0              6              0           0
          101        0              4              0           0
          102        0              3              0           0
          103        0              2              0           0
          104        0              1              0           0
    """
    from sage.rings.integer_ring import ZZ
    serieA = [0] + discrete_curve(50, 1000, K=1, R=sqrt(2), base=e) #Movember, Bye Bye, Cdf, MA
    la_flotte = [0] + discrete_curve(32, 800, K=1, R=sqrt(2), base=e)  # la flotte
    serieB = [0] + discrete_curve(24, 400, K=1, R=sqrt(2), base=e) # october fest, funenuf, la viree

    nb_ma = 104
    pivot_x = 40
    pivot_y = serieA[pivot_x]
    slope = (1-pivot_y) / (nb_ma-pivot_x)
    L = [pivot_y + slope * (p-pivot_x) for p in range(pivot_x, nb_ma+1)]
    L = [ZZ(round(_)) for _ in L]
    mars_attaque = serieA[:pivot_x] + L

    scales = serieA, mars_attaque, la_flotte, serieB
    names = ["Grand Chelem", "Mars Attaque", "La Flotte", "Petit Chelem"]
    return RankingScale(names, scales)

            def _2f1(a, b, c, z):
                """
                Evaluation of 2F1(a, b; c; z), assuming a, b, c positive
                integers or half-integers
                """
                if b == c:
                    return (1 - z) ** (-a)
                if a == c:
                    return (1 - z) ** (-b)
                if a == 0 or b == 0:
                    return Integer(1)
                if a > b:
                    a, b = b, a
                if b >= 2:
                    F1 = _2f1(a, b - 1, c, z)
                    F2 = _2f1(a, b - 2, c, z)
                    q = (b - 1) * (z - 1)
                    return (((c - 2 * b + 2 + (b - a - 1) * z) * F1 +
                            (b - c - 1) * F2) / q)
                if c > 2:
                    # how to handle this case?
                    if a - c + 1 == 0 or b - c + 1 == 0:
                        raise NotImplementedError
                    F1 = _2f1(a, b, c - 1, z)
                    F2 = _2f1(a, b, c - 2, z)
                    r1 = (c - 1) * (2 - c - (a + b - 2 * c + 3) * z)
                    r2 = (c - 1) * (c - 2) * (1 - z)
                    q = (a - c + 1) * (b - c + 1) * z
                    return (r1 * F1 + r2 * F2) / q

                if (a, b, c) == (R12, 1, 2):
                    return (2 - 2 * sqrt(1 - z)) / z
                if (a, b, c) == (1, 1, 2):
                    return -log(1 - z) / z
                if (a, b, c) == (1, R32, R12):
                    return (1 + z) / (1 - z) ** 2
                if (a, b, c) == (1, R32, 2):
                    return 2 * (1 / sqrt(1 - z) - 1) / z
                if (a, b, c) == (R32, 2, R12):
                    return (1 + 3 * z) / (1 - z) ** 3
                if (a, b, c) == (R32, 2, 1):
                    return (2 + z) / (2 * (sqrt(1 - z) * (1 - z) ** 2))
                if (a, b, c) == (2, 2, 1):
                    return (1 + z) / (1 - z) ** 3
                raise NotImplementedError

 def n_random_points(self,n):
     v = []
     while(len(v) < n):
         x = self.base_ring().random_element()
         test = (x**2 - 1)/(self._d*(x**2) - 1)
         if is_square(test):
             y = sqrt(test)
             v.append(self([x,y,1]))
     v.sort()
     return v