def Torus(r=2, R=3, name="Torus"):
        r"""
        Returns a torus obtained by revolving a circle of radius ``r`` around
        a coplanar axis ``R`` units away from the center of the circle. The
        parametrization used is

        .. MATH::

            \begin{aligned}
              x(u, v) & = (R + r \cos(v)) \cos(u); \\
              y(u, v) & = (R + r \cos(v)) \sin(u); \\
              z(u, v) & = r \sin(v).
            \end{aligned}

        INPUT:

        - ``r``, ``R`` -- Minor and major radius of the torus.

        - ``name`` -- string. Name of the surface.

        EXAMPLES::

            sage: torus = surfaces.Torus(); torus
            Parametrized surface ('Torus') with equation ((2*cos(v) + 3)*cos(u), (2*cos(v) + 3)*sin(u), 2*sin(v))
            sage: torus.plot()
            Graphics3d Object

        """

        u, v = var("u, v")
        torus_eq = [(R + r * cos(v)) * cos(u), (R + r * cos(v)) * sin(u), r * sin(v)]
        coords = ((u, 0, 2 * pi), (v, 0, 2 * pi))

        return ParametrizedSurface3D(torus_eq, coords, name)

    def Klein(r=1, name="Klein bottle"):
        r"""
        Returns the Klein bottle, in the figure-8 parametrization given by

        .. MATH::

            \begin{aligned}
              x(u, v) & = (r + \cos(u/2)\cos(v) - \sin(u/2)\sin(2v)) \cos(u); \\
              y(u, v) & = (r + \cos(u/2)\cos(v) - \sin(u/2)\sin(2v)) \sin(u); \\
              z(u, v) & = \sin(u/2)\cos(v) + \cos(u/2)\sin(2v).
            \end{aligned}

        INPUT:

        - ``r`` -- radius of the "figure-8" circle.

        - ``name`` -- string. Name of the surface.

        EXAMPLES::

            sage: klein = surfaces.Klein(); klein
            Parametrized surface ('Klein bottle') with equation (-(sin(1/2*u)*sin(2*v) - cos(1/2*u)*sin(v) - 1)*cos(u), -(sin(1/2*u)*sin(2*v) - cos(1/2*u)*sin(v) - 1)*sin(u), cos(1/2*u)*sin(2*v) + sin(1/2*u)*sin(v))
            sage: klein.plot()
            Graphics3d Object

        """

        u, v = var("u, v")
        x = (r + cos(u / 2) * sin(v) - sin(u / 2) * sin(2 * v)) * cos(u)
        y = (r + cos(u / 2) * sin(v) - sin(u / 2) * sin(2 * v)) * sin(u)
        z = sin(u / 2) * sin(v) + cos(u / 2) * sin(2 * v)
        klein_eq = [x, y, z]
        coords = ((u, 0, 2 * pi), (v, 0, 2 * pi))

        return ParametrizedSurface3D(klein_eq, coords, name)

    def Dini(a=1, b=1, name="Dini's surface"):
        r"""
        Returns Dini's surface, with parametrization

        .. MATH::

            \begin{aligned}
              x(u, v) & = a \cos(u)\sin(v); \\
              y(u, v) & = a \sin(u)\sin(v); \\
              z(u, v) & = u + \log(\tan(v/2)) + \cos(v).
            \end{aligned}

        INPUT:

        - ``a, b`` -- surface parameters.

        - ``name`` -- string. Name of the surface.

        EXAMPLES::

            sage: dini = surfaces.Dini(a=3, b=4); dini
            Parametrized surface ('Dini's surface') with equation (3*cos(u)*sin(v), 3*sin(u)*sin(v), 4*u + 3*cos(v) + 3*log(tan(1/2*v)))
            sage: dini.plot()  # not tested -- known bug (see #10132)

        """

        u, v = var("u, v")
        dini_eq = [a * cos(u) * sin(v), a * sin(u) * sin(v), a * (cos(v) + log(tan(v / 2))) + b * u]
        coords = ((u, 0, 2 * pi), (v, 0, 2 * pi))

        return ParametrizedSurface3D(dini_eq, coords, name)

    def Crosscap(r=1, name="Crosscap"):
        r"""
        Returns a crosscap surface, with parametrization

        .. MATH::

            \begin{aligned}
              x(u, v) & = r(1 + \cos(v)) \cos(u); \\
              y(u, v) & = r(1 + \cos(v)) \sin(u); \\
              z(u, v) & = - r\tanh(u - \pi) \sin(v).
            \end{aligned}

        INPUT:

        - ``r`` -- surface parameter.

        - ``name`` -- string. Name of the surface.

        EXAMPLES::

            sage: crosscap = surfaces.Crosscap(); crosscap
            Parametrized surface ('Crosscap') with equation ((cos(v) + 1)*cos(u), (cos(v) + 1)*sin(u), -sin(v)*tanh(-pi + u))
            sage: crosscap.plot()
            Graphics3d Object

        """

        u, v = var("u, v")
        crosscap_eq = [r * (1 + cos(v)) * cos(u), r * (1 + cos(v)) * sin(u), -tanh(u - pi) * r * sin(v)]
        coords = ((u, 0, 2 * pi), (v, 0, 2 * pi))

        return ParametrizedSurface3D(crosscap_eq, coords, name)

    def show(self, boundary=True, **options):
        r"""
        Plot ``self``.

        EXAMPLES::

            sage: HyperbolicPlane().PD().get_geodesic(0, 1).show()
            Graphics object consisting of 2 graphics primitives
        """
        opts = dict([('axes', False), ('aspect_ratio', 1)])
        opts.update(self.graphics_options())
        opts.update(options)
        end_1, end_2 = [CC(k.coordinates()) for k in self.endpoints()]
        bd_1, bd_2 = [CC(k.coordinates()) for k in self.ideal_endpoints()]
        # Check to see if it's a line
        if bool(real(bd_1)*imag(bd_2) - real(bd_2)*imag(bd_1))**2 < EPSILON:
            from sage.plot.line import line
            pic = line([(real(bd_1),imag(bd_1)),(real(bd_2),imag(bd_2))],
                       **opts)
        else:
            # If we are here, we know it's not a line
            # So we compute the center and radius of the circle
            center = (1/(real(bd_1)*imag(bd_2) - real(bd_2)*imag(bd_1)) *
                ((imag(bd_2)-imag(bd_1)) + (real(bd_1)-real(bd_2))*I))
            radius = RR(abs(bd_1 - center)) # abs is Euclidean distance
            # Now we calculate the angles for the parametric plot
            theta1 = CC(end_1 - center).arg()
            theta2 = CC(end_2 - center).arg()
            if theta2 < theta1:
                theta1, theta2 = theta2, theta1
            from sage.calculus.var import var
            from sage.plot.plot import parametric_plot
            x = var('x')
            mid = (theta1 + theta2)/2.0
            if (radius*cos(mid) + real(center))**2 + \
               (radius*sin(mid) + imag(center))**2 > 1.0:
                # Swap theta1 and theta2
                tmp = theta1 + 2*pi
                theta1 = theta2
                theta2 = tmp
                pic = parametric_plot((radius*cos(x) + real(center),
                                       radius*sin(x) + imag(center)),
                                      (x, theta1, theta2), **opts)

            else:
                pic = parametric_plot((radius*cos(x) + real(center),
                                   radius*sin(x) + imag(center)),
                                  (x, theta1, theta2), **opts)
        if boundary:
            bd_pic = self._model.get_background_graphic()
            pic = bd_pic + pic
        return pic

    def transform(self, radius=None, azimuth=None, elevation=None):
        """
        A spherical elevation coordinates transform.

        EXAMPLE::

            sage: T = SphericalElevation('radius', ['azimuth', 'elevation'])
            sage: T.transform(radius=var('r'), azimuth=var('theta'), elevation=var('phi'))
            (r*cos(phi)*cos(theta), r*cos(phi)*sin(theta), r*sin(phi))
        """
        return (radius * cos(elevation) * cos(azimuth),
                radius * cos(elevation) * sin(azimuth),
                radius * sin(elevation))

        def plot_arc(radius, p, q, **opts):
            # TODO: THIS SHOULD USE THE EXISTING PLOT OF ARCS!
            # plot the arc from p to q differently depending on the type of self
            p = ZZ(p)
            q = ZZ(q)
            t = var('t')
            if p - q in [1, -1]:
                def f(t):
                    return (radius * cos(t), radius * sin(t))
                (p, q) = sorted([p, q])
                angle_p = vertex_to_angle(p)
                angle_q = vertex_to_angle(q)
                return parametric_plot(f(t), (t, angle_q, angle_p), **opts)
            if self.type() == 'A':
                angle_p = vertex_to_angle(p)
                angle_q = vertex_to_angle(q)
                if angle_p < angle_q:
                    angle_p += 2 * pi
                internal_angle = angle_p - angle_q
                if internal_angle > pi:
                    (angle_p, angle_q) = (angle_q + 2 * pi, angle_p)
                    internal_angle = angle_p - angle_q
                angle_center = (angle_p+angle_q) / 2
                hypotenuse = radius / cos(internal_angle / 2)
                radius_arc = hypotenuse * sin(internal_angle / 2)
                center = (hypotenuse * cos(angle_center),
                          hypotenuse * sin(angle_center))
                center_angle_p = angle_p + pi / 2
                center_angle_q = angle_q + 3 * pi / 2

                def f(t):
                    return (radius_arc * cos(t) + center[0],
                            radius_arc * sin(t) + center[1])
                return parametric_plot(f(t), (t, center_angle_p,
                                              center_angle_q), **opts)
            elif self.type() == 'D':
                if p >= q:
                    q += self.r()
                px = -2 * pi * p / self.r() + pi / 2
                qx = -2 * pi * q / self.r() + pi / 2
                arc_radius = (px - qx) / 2
                arc_center = qx + arc_radius

                def f(t):
                    return exp(I * ((cos(t) + I * sin(t)) *
                                    arc_radius + arc_center)) * radius
                return parametric_plot((real_part(f(t)), imag_part(f(t))),
                                       (t, 0, pi), **opts)

    def Helicoid(h=1, name="Helicoid"):
        r"""
        Returns a helicoid surface, with parametrization

        .. MATH::

            \begin{aligned}
              x(\rho, \theta) & = \rho \cos(\theta); \\
              y(\rho, \theta) & = \rho \sin(\theta); \\
              z(\rho, \theta) & = h\theta/(2\pi).
            \end{aligned}

        INPUT:

        - ``h`` -- distance along the z-axis between two
          successive turns of the helicoid.

        - ``name`` -- string. Name of the surface.

        EXAMPLES::

            sage: helicoid = surfaces.Helicoid(h=2); helicoid
            Parametrized surface ('Helicoid') with equation (rho*cos(theta), rho*sin(theta), theta/pi)
            sage: helicoid.plot()
            Graphics3d Object

        """

        rho, theta = var("rho, theta")
        helicoid_eq = [rho * cos(theta), rho * sin(theta), h * theta / (2 * pi)]
        coords = ((rho, -2, 2), (theta, 0, 2 * pi))

        return ParametrizedSurface3D(helicoid_eq, [rho, theta], name)

    def Catenoid(c=1, name="Catenoid"):
        r"""
        Returns a catenoid surface, with parametric representation

        .. MATH::

            \begin{aligned}
              x(u, v) & = c \cosh(v/c) \cos(u); \\
              y(u, v) & = c \cosh(v/c) \sin(u); \\
              z(u, v) & = v.
            \end{aligned}

        INPUT:

        - ``c`` -- surface parameter.

        - ``name`` -- string. Name of the surface.


        EXAMPLES::

            sage: cat = surfaces.Catenoid(); cat
            Parametrized surface ('Catenoid') with equation (cos(u)*cosh(v), cosh(v)*sin(u), v)
            sage: cat.plot()
            Graphics3d Object

        """
        u, v = var("u, v")
        catenoid_eq = [c * cosh(v / c) * cos(u), c * cosh(v / c) * sin(u), v]
        coords = ((u, 0, 2 * pi), (v, -1, 1))

        return ParametrizedSurface3D(catenoid_eq, coords, name)

    def _i_rotation(self, z, alpha):
        r"""
        Return the resulting point after applying a hyperbolic
        rotation centered at `0 + i` and angle ``alpha`` to ``z``.

        INPUT:

        - ``z``-- point in the upper complex halfplane to which
          apply the isometry

        - ``alpha``-- angle of rotation (radians,counterwise)

        OUTPUT:

        - rotated point in the upper complex halfplane

        TESTS::

            sage: from sage.plot.hyperbolic_regular_polygon import HyperbolicRegularPolygon
            sage: P = HyperbolicRegularPolygon(4, pi/4, 1+I, {})
            sage: P._i_rotation(2+I, pi/2)
            I - 2
        """
        _a = alpha / 2
        _c = cos(_a)
        _s = sin(_a)
        G = matrix([[_c, _s], [-_s, _c]])
        return (G[0][0] * z + G[0][1]) / (G[1][0] * z + G[1][1])

    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 __init__(self, coxeter_matrix, base_ring, index_set):
        """
        Initialize ``self``.

        EXAMPLES::

            sage: W = CoxeterGroup([[1,3,2],[3,1,3],[2,3,1]])
            sage: TestSuite(W).run() # long time
            sage: W = CoxeterGroup([[1,3,2],[3,1,4],[2,4,1]], base_ring=QQbar)
            sage: TestSuite(W).run() # long time
            sage: W = CoxeterGroup([[1,3,2],[3,1,6],[2,6,1]])
            sage: TestSuite(W).run(max_runs=30) # long time
            sage: W = CoxeterGroup([[1,3,2],[3,1,-1],[2,-1,1]])
            sage: TestSuite(W).run(max_runs=30) # long time
        """
        self._matrix = coxeter_matrix
        self._index_set = index_set
        n = ZZ(coxeter_matrix.nrows())
        MS = MatrixSpace(base_ring, n, sparse=True)
        # FIXME: Hack because there is no ZZ \cup \{ \infty \}: -1 represents \infty
        if base_ring is UniversalCyclotomicField():
            val = lambda x: base_ring.gen(2*x) + ~base_ring.gen(2*x) if x != -1 else base_ring(2)
        else:
            from sage.functions.trig import cos
            from sage.symbolic.constants import pi
            val = lambda x: base_ring(2*cos(pi / x)) if x != -1 else base_ring(2)
        gens = [MS.one() + MS({(i, j): val(coxeter_matrix[i, j])
                               for j in range(n)})
                for i in range(n)]
        FinitelyGeneratedMatrixGroup_generic.__init__(self, n, base_ring,
                                                      gens,
                                                      category=CoxeterGroups())

    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 _circ_arc(t0, t1, c, r, num_pts=500):
    r""" Circular arc
    INPUTS:
    - ''t0'' -- starting parameter
    - ''t1'' -- ending parameter
    - ''c''  -- center point of the circle
    - ''r''  -- radius of circle
    - ''num_pts''  -- (default 100) number of points on polygon
    OUTPUT:
    - ''ca'' -- a polygonal approximation of a circular arc centered
    at c and radius r, starting at t0 and ending at t1


    EXAMPLES::

        sage: ca=_circ_arc(0.1,0.2,0.0,1.0,100)

    """
    from sage.plot.plot import parametric_plot
    from sage.functions.trig import cos, sin
    from sage.all import var

    t00 = t0
    t11 = t1
    ## To make sure the line is correct we reduce all arguments to the same branch,
    ## e.g. [0,2pi]
    pi = RR.pi()
    while t00 < 0.0:
        t00 = t00 + RR(2.0 * pi)
    while t11 < 0:
        t11 = t11 + RR(2.0 * pi)
    while t00 > 2 * pi:
        t00 = t00 - RR(2.0 * pi)
    while t11 > 2 * pi:
        t11 = t11 - RR(2.0 * pi)

    xc = CC(c).real()
    yc = CC(c).imag()
    # L0 =
    # [[RR(r*cos(t00+i*(t11-t00)/num_pts))+xc,RR(r*sin(t00+i*(t11-t00)/num_pts))+yc]
    # for i in range(0 ,num_pts)]
    t = var("t")
    if t11 > t00:
        ca = parametric_plot((r * cos(t) + xc, r * sin(t) + yc), (t, t00, t11))
    else:
        ca = parametric_plot((r * cos(t) + xc, r * sin(t) + yc), (t, t11, t00))
    return ca

def _draw_funddom(coset_reps,format="S"):
    r""" Draw a fundamental domain for G.
    
    INPUT:
    
    - ``format``  -- (default 'Disp') How to present the f.d.
    -   ``S`` -- Display directly on the screen
    
    EXAMPLES::        


    sage: G=MySubgroup(Gamma0(3))
    sage: G._draw_funddom()
        
    """
    pi=RR.pi()
    pi_3 = pi / RR(3.0)
    from sage.plot.plot import (Graphics,line)
    from sage.functions.trig import (cos,sin)
    g=Graphics()
    x1=RR(-0.5) ; y1=RR(sqrt(3 )/2 )
    x2=RR(0.5) ; y2=RR(sqrt(3 )/2 )
    xmax=RR(20.0) 
    l1 = line([[x1,y1],[x1,xmax]])
    l2 = line([[x2,y2],[x2,xmax]])
    l3 = line([[x2,xmax],[x1,xmax]]) # This is added to make a closed contour
    c0=_circ_arc(RR(pi/3.0) ,RR(2.0*pi)/RR(3.0) ,0 ,1 ,100 )
    tri=c0+l1+l3+l2
    g=g+tri
    for A in coset_reps:
        [a,b,c,d]=A
        if(a==1  and b==0  and c==0  and d==1 ):
            continue
        if(a<0 ):
            a=RR(-a); b=RR(-b); c=RR(-c); d=RR(-d) 
        else:
            a=RR(a); b=RR(b); c=RR(c); d=RR(d) 
        if(c==0 ): # then this is easier
            L0 = [[cos(pi_3*RR(i/100.0))+b,sin(pi_3*RR(i/100.0))] for i in range(100 ,201 )]
            L1 = [[x1+b,y1],[x1+b,xmax]]
            L2 = [[x2+b,y2],[x2+b,xmax]]
            L3 = [[x2+b,xmax],[x1+b,xmax]]
            c0=line(L0); l1=line(L1); l2=line(L2); l3=line(L3)
            tri=c0+l1+l3+l2
            g=g+tri
        else:
            den=(c*x1+d)**2 +c**2 *y1**2 
            x1_t=(a*c*(x1**2 +y1**2 )+(a*d+b*c)*x1+b*d)/den
            y1_t=y1/den
            den=(c*x2+d)**2 +c**2 *y2**2 
            x2_t=(a*c*(x2**2 +y2**2 )+(a*d+b*c)*x2+b*d)/den
            y2_t=y2/den
            inf_t=a/c
            c0=_geodesic_between_two_points(x1_t,y1_t,x2_t,y2_t)
            c1=_geodesic_between_two_points(x1_t,y1_t,inf_t,0. )
            c2=_geodesic_between_two_points(x2_t,y2_t,inf_t,0.0)
            tri=c0+c1+c2
            g=g+tri
    return g

 def plot(self):
     from sage.functions.trig import sin, cos
     from sage.plot.circle import circle
     from sage.plot.point import point
     from sage.plot.text import text
     p = circle((0,0),1)
     for i in range(self._dimension):
         a = self._alpha[i]
         p += point([cos(2*pi*a),sin(2*pi*a)], color='blue', size=100)
         p += text(r"$\alpha_%i$"%(self._i_alpha[i]+1),
                   [1.2*cos(2*pi*a),1.2*sin(2*pi*a)],fontsize=40)
     for i in range(self._dimension):
         b = self._beta[i]
         p += point([cos(2*pi*b),sin(2*pi*b)], color='red', size=100)
         p += text(r"$\beta_%i$"%(self._i_beta[i]+1),
                   [1.2*cos(2*pi*b),1.2*sin(2*pi*b)],fontsize=40)
     p.show(axes=False, xmin=-1, xmax=1, ymin=-1, ymax=1)

    def plot(self, show_box=False, colors=["white","lightgray","darkgray"]):
        r"""
        Return a plot of ``self``.

        INPUT:

        - ``show_box`` -- boolean (default: ``False``); if ``True``,
          also shows the visible tiles on the `xy`-, `yz`-, `zx`-planes

        - ``colors`` -- (default: ``["white", "lightgray", "darkgray"]``)
          list ``[A, B, C]`` of 3 strings representing colors

        EXAMPLES::

            sage: PP = PlanePartition([[4,3,3,1],[2,1,1],[1,1]])
            sage: PP.plot()
            Graphics object consisting of 27 graphics primitives
        """
        x = self._max_x
        y = self._max_y
        z = self._max_z
        from sage.functions.trig import cos, sin
        from sage.plot.polygon import polygon
        from sage.symbolic.constants import pi
        from sage.plot.plot import plot
        Uside = [[0,0], [cos(-pi/6),sin(-pi/6)], [0,-1], [cos(7*pi/6),sin(7*pi/6)]]
        Lside = [[0,0], [cos(-pi/6),sin(-pi/6)], [cos(pi/6),sin(pi/6)], [0,1]]
        Rside = [[0,0], [0,1], [cos(5*pi/6),sin(5*pi/6)], [cos(7*pi/6),sin(7*pi/6)]]
        Xdir = [cos(7*pi/6), sin(7*pi/6)]
        Ydir = [cos(-pi/6), sin(-pi/6)]
        Zdir = [0, 1]
        def move(side, i, j, k):
            return [[P[0]+i*Xdir[0]+j*Ydir[0]+k*Zdir[0],
                     P[1]+i*Xdir[1]+j*Ydir[1]+k*Zdir[1]]
                    for P in side]
        def add_topside(i, j, k):
            return polygon(move(Uside,i,j,k), edgecolor="black", color=colors[0])
        def add_leftside(i, j, k):
            return polygon(move(Lside,i,j,k), edgecolor="black", color=colors[1])
        def add_rightside(i, j, k):
            return polygon(move(Rside,i,j,k), edgecolor="black", color=colors[2])
        TP = plot([])
        for r in range(len(self.z_tableau())):
            for c in range(len(self.z_tableau()[r])):
                if self.z_tableau()[r][c] > 0 or show_box:
                    TP += add_topside(r, c, self.z_tableau()[r][c])
        for r in range(len(self.y_tableau())):
            for c in range(len(self.y_tableau()[r])):
                if self.y_tableau()[r][c] > 0 or show_box:
                    TP += add_rightside(c, self.y_tableau()[r][c], r)
        for r in range(len(self.x_tableau())):
            for c in range(len(self.x_tableau()[r])):
                if self.x_tableau()[r][c] > 0 or show_box:
                    TP += add_leftside(self.x_tableau()[r][c], r, c)
        TP.axes(show=False)
        return TP

    def __init__(self, sides, i_angle, center, options):
        """
        Initialize HyperbolicRegularPolygon.

        EXAMPLES::

            sage: from sage.plot.hyperbolic_regular_polygon import HyperbolicRegularPolygon
            sage: print(HyperbolicRegularPolygon(5,pi/2,I, {}))
            Hyperbolic regular polygon (sides=5, i_angle=1/2*pi, center=1.00000000000000*I)
        """
        self.center = CC(center)
        if self.center.imag() <= 0 :
            raise ValueError("center: %s is not a valid point in the upper half plane model of the hyperbolic plane"%(self.center))
        if sides < 3 :
            raise ValueError("degenerated polygons (sides<=2) are not supported")
        if i_angle <=0 or i_angle >= pi:
            raise ValueError("interior angle %s must be in (0, pi) interval"%(i_angle))
        if pi*(sides-2) - sides*i_angle <= 0 :
            raise ValueError("there exists no hyperbolic regular compact polygon, for sides=%s the interior angle must be less than %s"%(sides, pi * (sides-2) / sides))
        self.sides = sides
        self.i_angle = i_angle
        beta = 2 * pi / self.sides # compute the rotation angle to be used ahead
        alpha = self.i_angle / Integer(2)
        I = CC(0, 1)
        # compute using cosine theorem the radius of the circumscribed circle
        # using the triangle formed by the radius and the three known angles
        r = arccosh(cot(alpha) * (1 + cos(beta)) / sin(beta))

        # The first point will be always on the imaginary axis limited
        # to 8 digits for efficiency in the subsequent calculations.
        z_0 = [I*(e**r).n(digits=8)]

        # Compute the dilation isometry used to move the center
        # from I to the imaginary part of the given center.
        scale = self.center.imag()

        # Compute the parabolic isometry to move the center to the
        # real part of the given center.
        h_disp = self.center.real()

        d_z_k = [z_0[0]*scale + h_disp]  #d_k has the points for the polygon in the given center
        z_k = z_0                      #z_k has the Re(z)>0 vertices for the I centered polygon 
        r_z_k = []                     #r_z_k has the Re(z)<0 vertices
        if is_odd(self.sides):
            vert = (self.sides - 1) / 2
        else:
            vert = self.sides / 2 - 1
        for k in range(0, vert):
            # Compute with 8 digits to accelerate calculations
            new_z_k = self._i_rotation(z_k[-1], beta).n(digits=8)
            z_k = z_k + [new_z_k]
            d_z_k = d_z_k + [new_z_k * scale + h_disp]
            r_z_k=[-(new_z_k).conjugate() * scale + h_disp] + r_z_k
        if is_odd(self.sides):
            HyperbolicPolygon.__init__(self, d_z_k + r_z_k, options)
        else:
            z_opo = [I * (e**(-r)).n(digits=8) * scale + h_disp]
            HyperbolicPolygon.__init__(self, d_z_k + z_opo + r_z_k, options)

    def bilinear_form(self, R=None):
        """
        Return the bilinear form over ``R`` associated to ``self``.

        INPUT:

        - ``R`` -- (default: universal cyclotomic field) a ring used to
          compute the bilinear form

        EXAMPLES::

            sage: CoxeterType(['A', 2, 1]).bilinear_form()
            [   1 -1/2 -1/2]
            [-1/2    1 -1/2]
            [-1/2 -1/2    1]
            sage: CoxeterType(['H', 3]).bilinear_form()
            [                      1                    -1/2                       0]
            [                   -1/2                       1 1/2*E(5)^2 + 1/2*E(5)^3]
            [                      0 1/2*E(5)^2 + 1/2*E(5)^3                       1]
            sage: C = CoxeterMatrix([[1,-1,-1],[-1,1,-1],[-1,-1,1]])
            sage: C.bilinear_form()
            [ 1 -1 -1]
            [-1  1 -1]
            [-1 -1  1]
        """

        n = self.rank()
        mat = self.coxeter_matrix()._matrix
        base_ring = mat.base_ring()

        from sage.rings.universal_cyclotomic_field import UniversalCyclotomicField

        UCF = UniversalCyclotomicField()
        if UCF.has_coerce_map_from(base_ring):
            R = UCF
        else:
            R = base_ring
        # Compute the matrix with entries `- \cos( \pi / m_{ij} )`.
        if R is UCF:
            val = lambda x: (R.gen(2 * x) + ~R.gen(2 * x)) / R(-2) if x > -1 else R.one() * x
        else:
            from sage.functions.trig import cos
            from sage.symbolic.constants import pi

            val = lambda x: -R(cos(pi / SR(x))) if x > -1 else x

        MS = MatrixSpace(R, n, sparse=True)
        MC = MS._get_matrix_class()

        bilinear = MC(
            MS,
            entries={(i, j): val(mat[i, j]) for i in range(n) for j in range(n) if mat[i, j] != 2},
            coerce=True,
            copy=True,
        )
        bilinear.set_immutable()
        return bilinear

    def plot(self, **kwargs):
        """
        Return a graphics object representing the Kontsevich graph.

        INPUT:

        - ``edge_labels`` (boolean, default True) -- if True, show edge labels.
        - ``indices`` (boolean, default False) -- if True, show indices as
          edge labels instead of L and R; see :meth:`._latex_`.
        - ``upright`` (boolean, default False) -- if True, try to plot the
          graph with the ground vertices on the bottom and the rest on top.
        """
        if not 'edge_labels' in kwargs:
            kwargs['edge_labels'] = True        # show edge labels by default
        if 'indices' in kwargs:
            del kwargs['indices']
            KG = DiGraph(self)
            for (k,e) in enumerate(self.edges()):
                KG.delete_edge(e)
                KG.add_edge((e[0], e[1], chr(97 + k)))
            return KG.plot(**kwargs)
        if len(self.ground_vertices()) == 2 and 'upright' in kwargs:
            del kwargs['upright']
            kwargs['save_pos'] = True
            DiGraph.plot(self, **kwargs)
            positions = self.get_pos()
            # translate F to origin:
            F_pos = vector(positions[self.ground_vertices()[0]])
            for p in positions:
                positions[p] = list(vector(positions[p]) - F_pos)
            # scale F - G distance to 1:
            G_len = abs(vector(positions[self.ground_vertices()[1]]))
            for p in positions:
                positions[p] = list(vector(positions[p])/G_len)
            # rotate the vector F - G to (1,0)
            from math import atan2
            theta = -atan2(positions[self.ground_vertices()[1]][1], positions[self.ground_vertices()[1]][0])
            for p in positions:
                positions[p] = list(matrix([[cos(theta),-(sin(theta))],[sin(theta),cos(theta)]]) * vector(positions[p]))
            # flip if most things are below the x-axis:
            if len([(x,y) for (x,y) in positions.values() if y < 0])/len(self.internal_vertices()) > 0.5:
                for p in positions:
                    positions[p] = [positions[p][0], -positions[p][1]]
        return DiGraph.plot(self, **kwargs)