def __init__(self, projection_point):
        """
        Create a stereographic projection function.

        INPUT:

        - ``projection_point`` -- a list of coordinates in the
          appropriate dimension, which is the point projected from.

        EXAMPLES::

            sage: from sage.geometry.polyhedron.plot import ProjectionFuncStereographic
            sage: proj = ProjectionFuncStereographic([1.0,1.0])
            sage: proj.__init__([1.0,1.0])
            sage: proj.house
            [-0.7071067811...  0.7071067811...]
            [ 0.7071067811...  0.7071067811...]
            sage: TestSuite(proj).run(skip='_test_pickling')
        """
        self.projection_point = vector(projection_point)
        self.dim = self.projection_point.degree()

        pproj = vector(RDF,self.projection_point)
        self.psize = norm(pproj)
        if (self.psize).is_zero():
            raise ValueError, "projection direction must be a non-zero vector."
        v = vector(RDF, [0.0]*(self.dim-1) + [self.psize]) - pproj
        polediff = matrix(RDF,v).transpose()
        denom = RDF((polediff.transpose()*polediff)[0][0])
        if denom.is_zero():
            self.house = identity_matrix(RDF,self.dim)
        else:
            self.house = identity_matrix(RDF,self.dim) \
            - 2*polediff*polediff.transpose()/denom   # Householder reflector

    def __init__(self, v, omega, mu=0, prec=None):
        r"""
        EXAMPLES::

            sage: from slabbe import DiscreteHyperplane
            sage: p = DiscreteHyperplane([1,pi,7], 1+pi+7, mu=0)
            sage: p
            Set of points x in ZZ^3 satisfying: 0 <= (1, pi, 7) . x + 0 < pi + 8

        ::

            sage: p = DiscreteHyperplane([1,pi,7], 1+pi+7, mu=20)
            sage: vector((0,0,0)) in p
            False
            sage: p = DiscreteHyperplane([1,pi,7], 1+pi+7, mu=0)
            sage: vector((0,0,0)) in p
            True
        """
        if prec is None:
            self._v = vector(v)
            self._omega = omega
            self._mu = mu
        else:
            RF = RealField(prec=prec)
            self._v = vector(RF, v)
            self._omega = RF(omega)
            self._mu = RF(mu)
        def contain(p):
            #print("est-ce proche : ", self._v.dot_product(p) + self._mu)
            return  0 <= self._v.dot_product(p) + self._mu < self._omega
        DiscreteSubset.__init__(self, dimension=len(self._v), predicate=contain)

    def compatibility_degree(self, alpha, beta):
        if self.is_finite():
            tube_contribution = -1
        elif self.is_affine():
            gck = self.gamma().associated_coroot()
            if any([gck.scalar(alpha) != 0, gck.scalar(beta) != 0]):
                tube_contribution = -1
            else:
                sup_a = self._tube_support(alpha)
                sup_b = self._tube_support(beta)
                if all([x in sup_b for x in sup_a]) or all([x in sup_a for x in sup_b]):
                    tube_contribution = -1
                else:
                    nbh_a = self._tube_nbh(alpha)
                    tube_contribution = len([ x for x in nbh_a if x in sup_b ])
        else:
            raise ValueError("compatibility degree is implemented only for finite and affine types")
        
        initial = self.initial_cluster()
        if alpha in initial:
            return max(beta[initial.index(alpha)],0)

        alphacheck = alpha.associated_coroot()

        if beta in initial:
            return max(alphacheck[initial.index(beta)],0)

        Ap = -matrix(self.rk, map(lambda x: max(x,0), self.b_matrix().list() ) )
        Am =  matrix(self.rk, map(lambda x: min(x,0), self.b_matrix().list() ) )

        a = vector(alphacheck)
        b = vector(beta)

        return max( -a*b-a*Am*b, -a*b-a*Ap*b, tube_contribution )

def semi_norm_v(M, v,  p=2, verbose=False):
    r"""
    Return the semi norm on the hyperplane orthogonal to v.

    EXAMPLES::

        sage: from slabbe.matrix_cocycle import semi_norm_v
        sage: A1 = matrix(3, [1,-1,-1, 0,1,0, 0,0,1]).inverse()
        sage: semi_norm_v(A1, vector( (1,1,1)))      # tolerance 0.0001
        0.9999999999890247
        sage: semi_norm_v(A1, vector( (1,1,1)), p=1)   # tolerance 0.0001
        0.9999394820959548
        sage: semi_norm_v(A1, vector( (1,1,1)), p=oo)   # tolerance 0.0001
        1.0

    """
    from sage.modules.free_module_element import vector
    from sage.numerical.optimize import minimize_constrained
    def func(z):
        vz = vector(z)
        return - (M*vz).norm(p) / vz.norm(p)
    cons = [lambda z: v * vector(z),
            lambda z: - v * vector(z)]
    x0 = range(len(v))
    x0[0] = v[1]
    x0[1] = -v[0]
    rep = minimize_constrained(func, cons, x0)
    if verbose:
        print(rep, rep.norm(), rep*v)
    return -func(rep)

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 relative_field_representation(self, b):
        r"""
        Returns a vector representation of the field element ``b`` in the basis
        of the absolute field over the relative field.

        INPUT:

        - ``b`` -- an element of the absolute field

        EXAMPLES::

            sage: from sage.coding.relative_finite_field_extension import *
            sage: Fqm.<aa> = GF(16)
            sage: Fq.<a> = GF(4)
            sage: FE = RelativeFiniteFieldExtension(Fqm, Fq)
            sage: b = aa^3 + aa^2 + aa + 1
            sage: FE.relative_field_representation(b)
            (1, a + 1)
        """
        if not b in self.absolute_field():
            raise ValueError("The input has to be an element of the absolute field")
        s = self.relative_field_degree()
        if s == 1:
            return vector(b)
        else:
            Fq = self.relative_field()
            vect = self._flattened_relative_field_representation(b)
            sm = self.absolute_field_degree()
            list_elts = []
            for i in range(0, sm, s):
                list_elts.append(Fq(vect[i:i+s]))
            return vector(Fq, list_elts)

    def kernel_vector(self, way='LLL', verbose=False):
        r"""
        todo: clean this

        EXAMPLES::

            sage: from slabbe import ChristoffelGraph
            sage: C = ChristoffelGraph((2,5,7))
            sage: C.kernel_vector()
            [(-1, -1, 1), (3, -4, 0)]

        """
        from sage.arith.misc import gcd
        if way == 'vect_gcd':
            a,b,c = self._v
            gcd_ac = gcd(a,c)
            gcd_bc = gcd(b,c)
            U = ua,ub,uc = vector((c,0,-a)) / gcd(a,c)
            V = va,vb,vc = vector((0,c,-b)) / gcd(b,c)
            rows = U,V
        elif way == 'echelon':
            a,b,c = self._v
            m = matrix(ZZ, 4, [1,1,1,c,0,-a,0,c,-b,b,-a,0])
            me = m.echelon_form()
            if verbose:
                print(me)
            rows = me[1],me[2]
        elif way == 'LLL':
            dim = self.dimension()
            if dim == 3:
                a,b,c = self._v
                M = matrix(ZZ, 4, [1,1,1,c,0,-a,0,c,-b,b,-a,0])
            elif dim == 4:
                a,b,c,d = self._v
                M = matrix(ZZ, 7, (1,1,1,1,b,-a,0,0,c,0,-a,0,0,c,-b,0,d,0,0,-a,0,d,0,-b,0,0,d,-c))
            else:
                raise ValueError("dimension (=%s) must be 3 or 4" % dim)
            rows = M.LLL().rows()
            VS = rows[0].parent()
            zero = VS(0)
            un = VS((1,)*dim)
            assert zero in rows, "(0,0,0) not in LLL result"
            assert un in rows, "(1,1,1) not in LLL result"
            while zero in rows: rows.remove(zero)
            while un in rows: rows.remove(un)
        elif way == 'vect':
            a,b,c = self._v
            U = ua,ub,uc = vector((c,0,-a))
            V = va,vb,vc = vector((0,c,-b))
            rows = U,V
        else:
            raise ValueError("unknown way")
        R = matrix(rows)
        if sum(map(abs, R.minors(dim-1))) != sum(map(abs,self._v)):
            print(R)
            print(R.minors(dim-1))
            print(sum(map(abs, R.minors(dim))))
            print(sum(map(abs,self._v)))
            raise Exception("The result (=%s) is false " % rows)
        return rows

 def space(self):
     r'''
     Calculates the homology space as a Z-module.
     '''
     verb = get_verbose()
     set_verbose(0)
     V = self.coefficient_module()
     R = V.base_ring()
     Vdim = V.dimension()
     G = self.group()
     gens = G.gens()
     ambient = R**(Vdim * len(gens))
     if self.trivial_action():
         cycles = ambient
     else:
         # Now find the subspace of cycles
         A = Matrix(R, Vdim, 0)
         for g in gens:
             for v in V.gens():
                 A = A.augment(matrix(R,Vdim,1,list(vector(g**-1 * v - v))))
         K = A.right_kernel_matrix()
         cycles = ambient.submodule([ambient(list(o)) for o in K.rows()])
     boundaries = []
     for r in G.get_relation_words():
         grad = self.twisted_fox_gradient(G(r).word_rep)
         for v in V.gens():
             boundaries.append(cycles(ambient(sum([list(a * vector(v)) for a in grad],[]))))
     boundaries = cycles.submodule(boundaries)
     ans = cycles.quotient(boundaries)
     set_verbose(verb)
     return ans

    def __init__(self, v, start=(0,0,0)):
        r"""
        EXAMPLES::

            sage: from slabbe import BilliardCube
            sage: b = BilliardCube((1,pi,sqrt(2)))
            sage: b
            Cubic billiard of direction (1, pi, sqrt(2))

        TESTS::

            sage: vector((0,0,0)) in b
            True
            sage: vector((0,0,1)) in b
            False
            sage: vector((0,1,0)) in b
            True
            sage: vector((1,0,0)) in b
            False
            sage: vector((0,-1,0)) in b
            True

        """
        a,b,c = self._v = vector(v)
        sx,sy,sz = self._start = vector(start)
        px = DiscretePlane([0,c,-b], b+c, mu=(b+c)/2 - sy*c + sz*b)
        py = DiscretePlane([c,0,-a], a+c, mu=(a+c)/2 - sx*c + sz*a)
        pz = DiscretePlane([b,-a,0], a+b, mu=(a+b)/2 - sx*b + sy*a)
        Intersection.__init__(self, (px,py,pz))

def _arc(p,q,s,**kwds):
    #rewrite this to use polar_plot and get points to do filled triangles
    from sage.misc.functional import det
    from sage.plot.line import line
    from sage.misc.functional import norm
    from sage.symbolic.all import pi
    from sage.plot.arc import arc
     
    p,q,s = map( lambda x: vector(x), [p,q,s])
     
    # to avoid running into division by 0 we set to be colinear vectors that are
    # almost colinear
    if abs(det(matrix([p-s,q-s])))<0.01:
        return line((p,q),**kwds)
     
    (cx,cy)=var('cx','cy')
    equations=[
            2*cx*(s[0]-p[0])+2*cy*(s[1]-p[1]) == s[0]**2+s[1]**2-p[0]**2-p[1]**2,
            2*cx*(s[0]-q[0])+2*cy*(s[1]-q[1]) == s[0]**2+s[1]**2-q[0]**2-q[1]**2
            ]
    c = vector( [solve( equations, (cx,cy), solution_dict=True )[0][i] for i in [cx,cy]] )
     
    r = norm(p-c)
     
    a_p,a_q,a_s = map( _to_angle, [p-c,q-c,s-c])
    angles = [a_p,a_q,a_s]
    angles.sort()
     
    if a_s == angles[0]:
        return arc( c, r, angle=angles[2], sector=(0,2*pi-angles[2]+angles[1]), **kwds)
    if a_s == angles[1]:
        return arc( c, r, angle=angles[0], sector=(0,angles[2]-angles[0]), **kwds)
    if a_s == angles[2]:
        return arc( c, r, angle=angles[1], sector=(0,2*pi-angles[1]+angles[0]), **kwds)

    def decode_to_code(self, y):
        r"""
        Corrects the errors in ``word`` and returns a codeword.

        EXAMPLES::

            sage: C = codes.GeneralizedReedSolomonCode(GF(16, 'aa').list()[:13], 5)
            sage: Cs = codes.SubfieldSubcode(C, GF(4, 'a'))
            sage: D = codes.decoders.SubfieldSubcodeOriginalCodeDecoder(Cs)
            sage: Chan = channels.StaticErrorRateChannel(Cs.ambient_space(), D.decoding_radius())
            sage: c = Cs.random_element()
            sage: y = Chan(c)
            sage: c == D.decode_to_code(y)
            True
        """
        C = self.code()
        D = self.original_decoder()
        FE = C.embedding()
        phi = FE.embedding()
        y_or = vector([phi(i) for i in y])
        c_or = D.decode_to_code(y_or)
        if 'list-decoder' in self.decoder_type():
            result = []
            for c in c_or:
                if all(FE.is_in_relative_field(x) for x in c):
                    result.append(vector(map(FE.cast_into_relative_field, c)))
            return result
        else:
            if all(FE.is_in_relative_field(x) for x in c_or):
                return vector([FE.cast_into_relative_field(i, check=False)
                               for i in c_or])
            else:
                raise DecodingError("Original decoder does not output a "
                "subfield codeword. You may have exceeded the decoding radius.")

def plot_fan_stereographically(rays, walls, northsign=1, north=vector((-1,-1,-1)), right=vector((1,0,0)), colors=None, thickness=None):
    from sage.plot.graphics import Graphics
    from sage.plot.point import point
    from sage.misc.flatten import flatten
    from sage.plot.line import line
    from sage.misc.functional import norm
    
    if colors == None:
        colors = dict([('walls','black'),('rays','red')])

    if thickness == None:
        thickness = dict([('walls',0.5),('rays',20)])


    G = Graphics()
    
    for (u,v) in walls:
        G += _stereo_arc(vector(u),vector(v),vector(u+v),north=northsign*north,right=right,color=colors['walls'],thickness=thickness['walls'],zorder=len(G))
   
    for v in rays: 
        G += point(_stereo_coordinates(vector(v),north=northsign*north,right=right),color=colors['rays'],zorder=len(G),size=thickness['rays'])
    
    G.set_aspect_ratio(1)
    G._show_axes = False
    return G

    def first_coordinate_plane(self):
        """
        Restrict to the first coordinate plane.

        OUTPUT:

        A new double description pair with the constraint `x_0 = 0`
        added.

        EXAMPLES::

            sage: A = matrix([(1, 1), (-1, 1)])
            sage: from sage.geometry.polyhedron.double_description import StandardAlgorithm
            sage: DD, _ = StandardAlgorithm(A).initial_pair()
            sage: DD
            Double description pair (A, R) defined by
            A = [ 1  1],   R = [ 1/2 -1/2]
                [-1  1]        [ 1/2  1/2]
            sage: DD.first_coordinate_plane()
            Double description pair (A, R) defined by
                [ 1  1]
            A = [-1  1],   R = [  0]
                [-1  0]        [1/2]
                [ 1  0]
        """
        R = self.problem.base_ring()
        d = self.problem.dim()
        a_neg = vector(R, [-self.one] + [self.zero] * (d - 1))
        a_pos = vector(R, [+self.one] + [self.zero] * (d - 1))
        new = self._make_new(self.A, self.R)
        new.add_inequality(a_neg)
        new.add_inequality(a_pos)
        return new

    def parallelotope(self, generators):
        r"""
        Return the parallelotope spanned by the generators.

        INPUT:

        - ``generators`` -- an iterable of anything convertible to vector
          (for example, a list of vectors) such that the vectors all
          have the same dimension.

        OUTPUT:

        The parallelotope. This is the multi-dimensional
        generalization of a parallelogram (2 generators) and a
        parallelepiped (3 generators).

        EXAMPLES::

            sage: polytopes.parallelotope([ (1,0), (0,1) ])
            A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices
            sage: polytopes.parallelotope([[1,2,3,4],[0,1,0,7],[3,1,0,2],[0,0,1,0]])
            A 4-dimensional polyhedron in QQ^4 defined as the convex hull of 16 vertices
        """
        try:
            generators = [ vector(QQ,v) for v in generators ]
            base_ring = QQ
        except TypeError:
            generators = [ vector(RDF,v) for v in generators ]
            base_ring = RDF

        from sage.combinat.combination import Combinations
        par =  [ 0*generators[0] ]
        par += [ sum(c) for c in Combinations(generators) if c!=[] ]
        return Polyhedron(vertices=par, base_ring=base_ring)

    def plot_cluster_fan_stereographically(self, northsign=1, north=None, right=None, colors=None):
        from sage.plot.graphics import Graphics
        from sage.plot.point import point
        from sage.misc.flatten import flatten
        from sage.plot.line import line
        from sage.misc.functional import norm

        if self.rk !=3:
            raise ValueError("Can only stereographically project fans in 3d.")
        if not self.is_finite() and self._depth == infinity:
            raise ValueError("For infinite algebras you must specify the depth.")

        if north == None:
            if self.is_affine():
                north = vector(self.delta())
            else:
                north = vector( (-1,-1,-1) )
        if right == None:
            if self.is_affine():
                right = vector(self.gamma())
            else:
                right = vector( (1,0,0) )
        if colors == None:
            colors = dict([(0,'red'),(1,'green'),(2,'blue'),(3,'cyan'),(4,'yellow')])
        G = Graphics()

        roots = list(self.g_vectors())
        compatible = []
        while roots:
            x = roots.pop()
            for y in roots:
                if self.compatibility_degree(x,y) == 0:
                    compatible.append((x,y))
        for (u,v) in compatible:
            G += _stereo_arc(vector(u),vector(v),vector(u+v),north=northsign*north,right=right,thickness=0.5,color='black')

        for i in range(3):
            orbit = self.ith_orbit(i)
            for j in orbit:
                G += point(_stereo_coordinates(vector(orbit[j]),north=northsign*north,right=right),color=colors[i],zorder=len(G))

        if self.is_affine():
            tube_vectors = map(vector,flatten(self.affine_tubes()))
            for v in tube_vectors:
                G += point(_stereo_coordinates(v,north=northsign*north,right=right),color=colors[3],zorder=len(G))
            if north != vector(self.delta()):
                G += _stereo_arc(tube_vectors[0],tube_vectors[1],vector(self.delta()),north=northsign*north,right=right,thickness=2,color=colors[4],zorder=0)
            else:
                # FIXME: refactor this before publishing
                tube_projections = [
                        _stereo_coordinates(v,north=northsign*north,right=right)
                        for v in tube_vectors ]
                t=min((G.get_minmax_data()['xmax'],G.get_minmax_data()['ymax']))
                G += line([tube_projections[0],tube_projections[0]+t*(_normalize(tube_projections[0]-tube_projections[1]))],thickness=2,color=colors[4],zorder=0)
                G += line([tube_projections[1],tube_projections[1]+t*(_normalize(tube_projections[1]-tube_projections[0]))],thickness=2,color=colors[4],zorder=0)
        G.set_aspect_ratio(1)
        G._show_axes = False
        return G

    def closest_vector(self, t):
        """
        Compute the closest vector in the embedded lattice to a given vector.

        INPUT:

        - ``t`` -- the target vector to compute the closest vector to

        OUTPUT:

        The vector in the lattice closest to ``t``.

        EXAMPLES::

            sage: from sage.modules.free_module_integer import IntegerLattice
            sage: L = IntegerLattice([[1, 0], [0, 1]])
            sage: L.closest_vector((-6, 5/3))
            (-6, 2)

        ALGORITHM:

        Uses the algorithm from [Mic2010]_.

        REFERENCES:

        .. [Mic2010] D. Micciancio, P. Voulgaris. *A Deterministic Single
           Exponential Time Algorithm for Most Lattice Problems based on
           Voronoi Cell Computations*.
           Proceedings of the 42nd ACM Symposium Theory of Computation, 2010.
        """
        voronoi_cell = self.voronoi_cell()

        def projection(M, v):
            Mt = M.transpose()
            P = Mt * (M * Mt) ** (-1) * M
            return P * v

        t = projection(matrix(self.reduced_basis), vector(t))

        def CVPP_2V(t, V, voronoi_cell):
            t_new = t
            while not voronoi_cell.contains(t_new.list()):
                v = max(V, key=lambda v: t_new * v / v.norm() ** 2)
                t_new = t_new - v
            return t - t_new

        V = self.voronoi_relevant_vectors()
        t = vector(t)
        p = 0
        while not (ZZ(2 ** p) * voronoi_cell).contains(t):
            p += 1
        t_new = t
        i = p
        while i >= 1:
            V_scaled = [v * (2 ** (i - 1)) for v in V]
            t_new = t_new - CVPP_2V(t_new, V_scaled, ZZ(2 ** (i - 1)) * voronoi_cell)
            i -= 1
        return t - t_new

    def __init__(self, A, elt=None, check=True):
        """
        TESTS::

            sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
            sage: A(QQ(4))
            Traceback (most recent call last):
            ...
            TypeError: elt should be a vector, a matrix, or an element of the base field

            sage: B = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [-1,0]])])
            sage: elt = B(Matrix([[1,1], [-1,1]])); elt
            e0 + e1
            sage: TestSuite(elt).run()
            sage: B(Matrix([[0,1], [1,0]]))
            Traceback (most recent call last):
            ...
            ValueError: matrix does not define an element of the algebra
        """
        AlgebraElement.__init__(self, A)
        k = A.base_ring()
        n = A.degree()
        if elt is None:
            self._vector = vector(k, n)
            self._matrix = Matrix(k, n)
        else:
            if isinstance(elt, int):
                elt = Integer(elt)
            elif isinstance(elt, list):
                elt = vector(elt)
            if A == elt.parent():
                self._vector = elt._vector.base_extend(k)
                self._matrix = elt._matrix.base_extend(k)
            elif k.has_coerce_map_from(elt.parent()):
                e = k(elt)
                if e == 0:
                    self._vector = vector(k, n)
                    self._matrix = Matrix(k, n)
                elif A.is_unitary():
                    self._vector = A._one * e
                    self._matrix = Matrix.identity(k, n) * e
                else:
                    raise TypeError("algebra is not unitary")
            elif is_Vector(elt):
                self._vector = elt.base_extend(k)
                self._matrix = Matrix(k, sum([elt[i] * A.table()[i] for i in xrange(n)]))
            elif is_Matrix(elt):
                if not A.is_unitary():
                    raise TypeError("algebra is not unitary")
                self._vector = A._one * elt
                if not check or sum([self._vector[i]*A.table()[i] for i in xrange(n)]) == elt:
                    self._matrix = elt
                else:
                    raise ValueError("matrix does not define an element of the algebra")
            else:
                raise TypeError("elt should be a vector, a matrix, " +
                                "or an element of the base field")

def _normalize(x):
    r"""
    make x of length 1
    """
    from sage.misc.functional import norm
    x=vector(RR, x)
    if norm(x) == 0:
        return x
    return vector(x/norm(x))

    def has_rational_point(self, point = False,
                           algorithm = 'default', read_cache = True):
        r"""
        Returns True if and only if the conic ``self``
        has a point over its base field `B`.

        If ``point`` is True, then returns a second output, which is
        a rational point if one exists.
        
        Points are cached whenever they are found. Cached information
        is used if and only if ``read_cache`` is True.
        
        EXAMPLES:

            sage: Conic(RR, [1, 1, 1]).has_rational_point()
            False
            sage: Conic(CC, [1, 1, 1]).has_rational_point()
            True
            
            sage: Conic(RR, [1, 2, -3]).has_rational_point(point = True)
            (True, (1.73205080756888 : 0.000000000000000 : 1.00000000000000))
        """
        if read_cache:
            if self._rational_point is not None:
                if point:
                    return True, self._rational_point
                else:
                    return True
        B = self.base_ring()
        if is_ComplexField(B):
            if point:
                [_,_,_,d,e,f] = self._coefficients
                if d == 0:
                    return True, self.point([0,1,0])
                return True, self.point([0, ((e**2-4*d*f).sqrt()-e)/(2*d), 1],
                                        check = False)
            return True
        if is_RealField(B):
            D, T = self.diagonal_matrix()
            [a, b, c] = [D[0,0], D[1,1], D[2,2]]
            if a == 0:
                ret = True, self.point(T*vector([1,0,0]), check = False)
            elif a*c <= 0:
                ret = True, self.point(T*vector([(-c/a).sqrt(),0,1]),
                                       check = False)
            elif b == 0:
                ret = True, self.point(T*vector([0,1,0]), check = False)
            elif b*c <= 0:
                ret = True, self.point(T*vector([0,(-c/b).sqrt(),0,1]),
                                       check = False)
            else:
                ret = False, None
            if point:
                return ret
            return ret[0]
        raise NotImplementedError, "has_rational_point not implemented for " \
                                   "conics over base field %s" % B 

def _stereo_arc(x,y, xy=None,  north=(1,0,0), right=(0,1,0), translation=-1, **kwds):
    from sage.misc.functional import n
    x=vector(x)
    y=vector(y)
    sx=n(_stereo_coordinates(x, north=north, right=right, translation=translation))
    sy=n(_stereo_coordinates(y, north=north, right=right, translation=translation))
    if xy == None:
        xy=x+y
    sxy=n(_stereo_coordinates(xy, north=north, right=right, translation=translation))
    return _arc(sx,sy,sxy,**kwds)