def point3d(v, size=5, **kwds):
    """
    Plot a point or list of points in 3d space.
    
    INPUT:
    
    
    -  ``v`` - a point or list of points
    
    -  ``size`` - (default: 5) size of the point (or
       points)
    
    -  ``color`` - a word that describes a color
    
    -  ``rgbcolor`` - (r,g,b) with r, g, b between 0 and 1
       that describes a color
    
    -  ``opacity`` - (default: 1) if less than 1 then is
       transparent
    
    
    EXAMPLES::
    
        sage: sum([point3d((i,i^2,i^3), size=5) for i in range(10)])

    We check to make sure this works with vectors::

        sage: pl = point3d([vector(ZZ,(1, 0, 0)), vector(ZZ,(0, 1, 0)), (-1, -1, 0)])
        sage: p = point(vector((2,3,4)))
        sage: print p
        Graphics3d Object
        

    We check to make sure the options work::

        sage: point3d((4,3,2),size=20,color='red',opacity=.5)
        
    """
    if (
        (isinstance(v, (list, tuple)) or is_Vector(v))
        and len(v) == 3
        and not (isinstance(v[0], (list, tuple)) or is_Vector(v[0]))
    ):
        return Point(v, size, **kwds)
    else:
        A = sum([Point(z, size, **kwds) for z in v])
        A._set_extra_kwds(kwds)
        return A

    def eval(self, Vobj):
        r"""
        Evaluates the left hand side `A\vec{x}+b` on the given
        vertex/ray/line.

        NOTES:

          * Evaluating on a vertex returns `A\vec{x}+b`
          * Evaluating on a ray returns `A\vec{r}`. Only the sign or
            whether it is zero is meaningful.
          * Evaluating on a line returns `A\vec{l}`. Only whether it
            is zero or not is meaningful.

        EXAMPLES::

            sage: triangle = Polyhedron(vertices=[[1,0],[0,1],[-1,-1]])
            sage: ineq = next(triangle.inequality_generator())
            sage: ineq
            An inequality (2, -1) x + 1 >= 0
            sage: [ ineq.eval(v) for v in triangle.vertex_generator() ]
            [0, 0, 3]
            sage: [ ineq * v for v in triangle.vertex_generator() ]
            [0, 0, 3]

        If you pass a vector, it is assumed to be the coordinate vector of a point::

            sage: ineq.eval( vector(ZZ, [3,2]) )
            5
        """
        if is_Vector(Vobj):
            return self.A() * Vobj + self.b()
        return Vobj.evaluated_on(self)

 def pyobject(self, ex, obj):
     from mathics.core import expression
     from mathics.core.expression import Number
     
     if obj is None:
         return expression.Symbol('Null')
     elif isinstance(obj, (list, tuple)) or is_Vector(obj):
         return expression.Expression('List', *(from_sage(item, self.subs) for item in obj))
     elif isinstance(obj, Constant):
         return expression.Symbol(obj._conversions.get('mathematica', obj._name))
     elif is_Integer(obj):
         return expression.Integer(str(obj))
     elif isinstance(obj, sage.Rational):
         rational = expression.Rational(str(obj))
         if rational.value.denom() == 1:
             return expression.Integer(rational.value.numer())
         else:
             return rational
     elif isinstance(obj, sage.RealDoubleElement) or is_RealNumber(obj):
         return expression.Real(str(obj))
     elif is_ComplexNumber(obj):
         real = Number.from_string(str(obj.real())).value
         imag = Number.from_string(str(obj.imag())).value
         return expression.Complex(real, imag)
     elif isinstance(obj, NumberFieldElement_quadratic):
         # TODO: this need not be a complex number, but we assume so!
         real = Number.from_string(str(obj.real())).value
         imag = Number.from_string(str(obj.imag())).value
         return expression.Complex(real, imag)
     else:
         return expression.from_python(obj)

    def point(self, v, check=True):
        r"""
        Constructs a point on ``self`` corresponding to the input ``v``.

        If ``check`` is True, then checks if ``v`` defines a valid
        point on ``self``.

        If no rational point on ``self`` is known yet, then also caches the point
        for use by ``self.rational_point()`` and ``self.parametrization()``.

        EXAMPLES ::

            sage: c = Conic([1, -1, 1])
            sage: c.point([15, 17, 8])
            (15/8 : 17/8 : 1)
            sage: c.rational_point()
            (15/8 : 17/8 : 1)
            sage: d = Conic([1, -1, 1])
            sage: d.rational_point()
            (1 : 1 : 0)
        """
        if is_Vector(v):
            v = Sequence(v)
        p = ProjectiveCurve_generic.point(self, v, check=check)
        if self._rational_point is None:
            self._rational_point = p
        return p

    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 jacobian(functions, variables):
    """
    Return the Jacobian matrix, which is the matrix of partial
    derivatives in which the i,j entry of the Jacobian matrix is the
    partial derivative diff(functions[i], variables[j]).

    EXAMPLES::

        sage: x,y = var('x,y')
        sage: g=x^2-2*x*y
        sage: jacobian(g, (x,y))
        [2*x - 2*y      -2*x]

    The Jacobian of the Jacobian should give us the "second derivative", which is the Hessian matrix::

        sage: jacobian(jacobian(g, (x,y)), (x,y))
        [ 2 -2]
        [-2  0]
        sage: g.hessian()
        [ 2 -2]
        [-2  0]

        sage: f=(x^3*sin(y), cos(x)*sin(y), exp(x))
        sage: jacobian(f, (x,y))
        [  3*x^2*sin(y)     x^3*cos(y)]
        [-sin(x)*sin(y)  cos(x)*cos(y)]
        [           e^x              0]
        sage: jacobian(f, (y,x))
        [    x^3*cos(y)   3*x^2*sin(y)]
        [ cos(x)*cos(y) -sin(x)*sin(y)]
        [             0            e^x]

    """
    if is_Matrix(functions) and (functions.nrows()==1 or functions.ncols()==1):
        functions = functions.list()
    elif not (isinstance(functions, (tuple, list)) or is_Vector(functions)):
        functions = [functions]

    if not isinstance(variables, (tuple, list)) and not is_Vector(variables):
        variables = [variables]

    return matrix([[diff(f, v) for v in variables] for f in functions])

#.........这里部分代码省略.........
       tuple is returned of the form (range_min, range_max,
       range_step) such that ``srange(range_min, range_max,
       range_step, include_endpoint=True)`` gives the correct points
       for evaluation.

    EXAMPLES::

        sage: x,y,z=var('x,y,z')
        sage: f(x,y)=x+y-z
        sage: g(x,y)=x+y
        sage: h(y)=-y
        sage: sage.plot.misc.setup_for_eval_on_grid(f, [(0, 2),(1,3),(-4,1)], plot_points=5)
        (<sage.ext...>, [(0.0, 2.0, 0.5), (1.0, 3.0, 0.5), (-4.0, 1.0, 1.25)])
        sage: sage.plot.misc.setup_for_eval_on_grid([g,h], [(0, 2),(-1,1)], plot_points=5)
        ((<sage.ext...>, <sage.ext...>), [(0.0, 2.0, 0.5), (-1.0, 1.0, 0.5)])
        sage: sage.plot.misc.setup_for_eval_on_grid([sin,cos], [(-1,1)], plot_points=9)
        ((<sage.ext...>, <sage.ext...>), [(-1.0, 1.0, 0.25)])
        sage: sage.plot.misc.setup_for_eval_on_grid([lambda x: x^2,cos], [(-1,1)], plot_points=9)
        ((<function <lambda> ...>, <sage.ext...>), [(-1.0, 1.0, 0.25)])
        sage: sage.plot.misc.setup_for_eval_on_grid([x+y], [(x,-1,1),(y,-2,2)])
        ((<sage.ext...>,), [(-1.0, 1.0, 2.0), (-2.0, 2.0, 4.0)])
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(x,-1,1),(y,-1,1)], plot_points=[4,9])
        (<sage.ext...>, [(-1.0, 1.0, 0.6666666666666666), (-1.0, 1.0, 0.25)])
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(x,-1,1),(y,-1,1)], plot_points=[4,9,10])
        Traceback (most recent call last):
        ...
        ValueError: plot_points must be either an integer or a list of integers, one for each range
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(1,-1),(y,-1,1)], plot_points=[4,9,10])
        Traceback (most recent call last):
        ...
        ValueError: Some variable ranges specify variables while others do not
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(1,-1),(-1,1)], plot_points=5)
        doctest:...: DeprecationWarning:
        Unnamed ranges for more than one variable is deprecated and
        will be removed from a future release of Sage; you can used
        named ranges instead, like (x,0,2)
        See http://trac.sagemath.org/7008 for details.
        (<sage.ext...>, [(1.0, -1.0, 0.5), (-1.0, 1.0, 0.5)])
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(y,1,-1),(x,-1,1)], plot_points=5)
        (<sage.ext...>, [(1.0, -1.0, 0.5), (-1.0, 1.0, 0.5)])
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(x,1,-1),(x,-1,1)], plot_points=5)
        Traceback (most recent call last):
        ...
        ValueError: range variables should be distinct, but there are duplicates
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(x,1,1),(y,-1,1)])
        Traceback (most recent call last):
        ...
        ValueError: plot start point and end point must be different
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(x,1,-1),(y,-1,1)], return_vars=True)
        (<sage.ext...>, [(1.0, -1.0, 2.0), (-1.0, 1.0, 2.0)], [x, y])
        sage: sage.plot.misc.setup_for_eval_on_grid(x+y, [(y,1,-1),(x,-1,1)], return_vars=True)
        (<sage.ext...>, [(1.0, -1.0, 2.0), (-1.0, 1.0, 2.0)], [y, x])
    """
    if max(map(len, ranges)) != min(map(len, ranges)):
        raise ValueError("Some variable ranges specify variables while others do not")

    if len(ranges[0])==3:
        vars = [r[0] for r in ranges]
        ranges = [r[1:] for r in ranges]
        if len(set(vars))<len(vars):
            raise ValueError("range variables should be distinct, but there are duplicates")
    else:
        vars, free_vars = unify_arguments(funcs)
        if len(free_vars)>1:
            from sage.misc.superseded import deprecation
            deprecation(7008, "Unnamed ranges for more than one variable is deprecated and will be removed from a future release of Sage; you can used named ranges instead, like (x,0,2)")

    # pad the variables if we don't have enough
    nargs = len(ranges)
    if len(vars)<nargs:
        vars += ('_',)*(nargs-len(vars))

    ranges = [[float(z) for z in r] for r in ranges]

    if plot_points is None:
        plot_points=2

    if not isinstance(plot_points, (list, tuple)):
        plot_points = [plot_points]*len(ranges)
    elif len(plot_points)!=nargs:
        raise ValueError("plot_points must be either an integer or a list of integers, one for each range")

    plot_points = [int(p) if p>=2 else 2 for p in plot_points]
    range_steps = [abs(range[1] - range[0])/(p-1) for range, p in zip(ranges, plot_points)]
    if min(range_steps) == float(0):
        raise ValueError("plot start point and end point must be different")

    options={}
    if nargs==1:
        options['expect_one_var']=True

    if is_Vector(funcs):
        funcs = list(funcs)

    #TODO: raise an error if there is a function/method in funcs that takes more values than we have ranges

    if return_vars:
        return fast_float(funcs, *vars,**options), [tuple(range+[range_step]) for range,range_step in zip(ranges, range_steps)], vars
    else:
        return fast_float(funcs, *vars,**options), [tuple(range+[range_step]) for range,range_step in zip(ranges, range_steps)]

    def __call__(self, v):
        """
        Evaluate this quadratic form Q on a vector or matrix of elements
        coercible to the base ring of the quadratic form.  If a vector
        is given then the output will be the ring element Q(`v`), but if a
        matrix is given then the output will be the quadratic form Q'
        which in matrix notation is given by:

        .. math::
                Q' = v^t * Q * v.


        EXAMPLES::

            ## Evaluate a quadratic form at a vector:
            ## --------------------------------------
            sage: Q = QuadraticForm(QQ, 3, range(6))
            sage: Q
            Quadratic form in 3 variables over Rational Field with coefficients:
            [ 0 1 2 ]
            [ * 3 4 ]
            [ * * 5 ]
            sage: Q([1,2,3])
            89
            sage: Q([1,0,0])
            0
            sage: Q([1,1,1])
            15

    ::

            ## Evaluate a quadratic form using a column matrix:
            ## ------------------------------------------------
            sage: Q = QuadraticForm(QQ, 2, range(1,4))
            sage: A = Matrix(ZZ,2,2,[-1,0,0,1])
            sage: Q(A)
            Quadratic form in 2 variables over Rational Field with coefficients:
            [ 1 -2 ]
            [ * 3 ]
            sage: Q([1,0])
            1
            sage: type(Q([1,0]))
            <type 'sage.rings.rational.Rational'>
            sage: Q = QuadraticForm(QQ, 2, range(1,4))
            sage: Q(matrix(2, [1,0]))
            Quadratic form in 1 variables over Rational Field with coefficients:
            [ 1 ]

        ::

            ## Simple 2x2 change of variables:
            ## -------------------------------
            sage: Q = QuadraticForm(ZZ, 2, [1,0,1])
            sage: Q
            Quadratic form in 2 variables over Integer Ring with coefficients:
            [ 1 0 ]
            [ * 1 ]
            sage: M = Matrix(ZZ, 2, 2, [1,1,0,1])
            sage: M
            [1 1]
            [0 1]
            sage: Q(M)
            Quadratic form in 2 variables over Integer Ring with coefficients:
            [ 1 2 ]
            [ * 2 ]

        ::

            ## Some more tests:
            ## ----------------
            sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
            sage: Q([1,2,3])
            14
            sage: v = vector([1,2,3])
            sage: Q(v)
            14
            sage: t = tuple([1,2,3])
            sage: Q(v)
            14
            sage: M = Matrix(ZZ, 3, [1,2,3])
            sage: Q(M)
            Quadratic form in 1 variables over Integer Ring with coefficients:
            [ 14 ]

        """
        ## If we are passed a matrix A, return the quadratic form Q(A(x))
        ## (In matrix notation: A^t * Q * A)
        n = self.dim()

        if is_Matrix(v):
            ## Check that v has the correct number of rows
            if v.nrows() != n:
                raise TypeError("Oops!  The matrix must have " + str(n) + " rows. =(")

            ## Create the new quadratic form
            m = v.ncols()
            Q2 = QuadraticForm(self.base_ring(), m)
            return QFEvaluateMatrix(self, v, Q2)

        elif (is_Vector(v) or isinstance(v, (list, tuple))):
#.........这里部分代码省略.........

#.........这里部分代码省略.........
        f_x = cos(u) * cos(v) * G1 * G2
        f_y = cos(u) * sin(v) * G1 * G2
        f_z = sin(u) * G1
        sphinx_plot(parametric_plot3d([f_x, f_y, f_z], (u,-pi/2,pi/2), (v,0,2*pi), plot_points=[50,50], frame=False, color="green"))

    A cool self-intersecting surface (Eppener surface?)::

        sage: f_x = u - u^3/3 + u*v^2
        sage: f_y = v - v^3/3 + v*u^2
        sage: f_z = u^2 - v^2
        sage: parametric_plot3d([f_x, f_y, f_z], (u,-25,25), (v,-25,25), plot_points=[50,50], frame=False, color="green")
        Graphics3d Object

    .. PLOT::

        u, v = var('u,v')
        f_x = u - u**3/3 + u*v**2
        f_y = v - v**3/3 + v*u**2
        f_z = u**2 - v**2
        sphinx_plot(parametric_plot3d([f_x, f_y, f_z], (u,-25,25), (v,-25,25), plot_points=[50,50], frame=False, color="green"))

    The breather surface
    (http://en.wikipedia.org/wiki/Breather_surface)::

        sage: K = sqrt(0.84)
        sage: G = (0.4*((K*cosh(0.4*u))^2 + (0.4*sin(K*v))^2))
        sage: f_x = (2*K*cosh(0.4*u)*(-(K*cos(v)*cos(K*v)) - sin(v)*sin(K*v)))/G
        sage: f_y = (2*K*cosh(0.4*u)*(-(K*sin(v)*cos(K*v)) + cos(v)*sin(K*v)))/G
        sage: f_z = -u + (2*0.84*cosh(0.4*u)*sinh(0.4*u))/G
        sage: parametric_plot3d([f_x, f_y, f_z], (u,-13.2,13.2), (v,-37.4,37.4), plot_points=[90,90], frame=False, color="green")
        Graphics3d Object

    .. PLOT::

        u, v = var('u,v')
        K = sqrt(0.84)
        G = (0.4*((K*cosh(0.4*u))**2 + (0.4*sin(K*v))**2))
        f_x = (2*K*cosh(0.4*u)*(-(K*cos(v)*cos(K*v)) - sin(v)*sin(K*v)))/G
        f_y = (2*K*cosh(0.4*u)*(-(K*sin(v)*cos(K*v)) + cos(v)*sin(K*v)))/G
        f_z = -u + (2*0.84*cosh(0.4*u)*sinh(0.4*u))/G
        sphinx_plot(parametric_plot3d([f_x, f_y, f_z], (u,-13.2,13.2), (v,-37.4,37.4), plot_points=[90,90], frame=False, color="green"))

    TESTS::

        sage: u, v = var('u,v')
        sage: plot3d(u^2-v^2, (u,-1,1), (u,-1,1))
        Traceback (most recent call last):
        ...
        ValueError: range variables should be distinct, but there are duplicates


    From :trac:`2858`::

        sage: parametric_plot3d((u,-u,v), (u,-10,10),(v,-10,10))
        Graphics3d Object
        sage: f(u)=u; g(v)=v^2; parametric_plot3d((g,f,f), (-10,10),(-10,10))
        Graphics3d Object

    From :trac:`5368`::

        sage: x, y = var('x,y')
        sage: plot3d(x*y^2 - sin(x), (x,-1,1), (y,-1,1))
        Graphics3d Object

    """
    # TODO:
    #   * Surface -- behavior of functions not defined everywhere -- see note above
    #   * Iterative refinement

    # color_function -- (default: "automatic") how to determine the color of curves and surfaces
    # color_function_scaling -- (default: True) whether to scale the input to color_function
    # exclusions -- (default: "automatic") u points or (u,v) conditions to exclude.
    #         (E.g., exclusions could be a function e = lambda u, v: False if u < v else True
    # exclusions_style -- (default: None) what to draw at excluded points
    # max_recursion -- (default: "automatic") maximum number of recursive subdivisions,
    #                   when ...
    # mesh -- (default: "automatic") how many mesh divisions in each direction to draw
    # mesh_functions -- (default: "automatic") how to determine the placement of mesh divisions
    # mesh_shading -- (default: None) how to shade regions between mesh divisions
    # plot_range -- (default: "automatic") range of values to include

    if is_Vector(f):
        f = tuple(f)

    if isinstance(f, (list, tuple)) and len(f) > 0 and isinstance(f[0], (list, tuple)):
        return sum([parametric_plot3d(v, urange, vrange, plot_points=plot_points, **kwds) for v in f])

    if not isinstance(f, (tuple, list)) or len(f) != 3:
        raise ValueError("f must be a list, tuple, or vector of length 3")

    if vrange is None:
        if plot_points == "automatic":
            plot_points = 75
        G = _parametric_plot3d_curve(f, urange, plot_points=plot_points, **kwds)
    else:
        if plot_points == "automatic":
            plot_points = [40, 40]
        G = _parametric_plot3d_surface(f, urange, vrange, plot_points=plot_points, boundary_style=boundary_style, **kwds)
    G._set_extra_kwds(kwds)
    return G

#.........这里部分代码省略.........
        sage: fx = (sinh(v)*cos(3*u))/(1+cosh(u)*cosh(v))
        sage: fy = (sinh(v)*sin(3*u))/(1+cosh(u)*cosh(v))
        sage: fz = (cosh(v)*sinh(u))/(1+cosh(u)*cosh(v))
        sage: parametric_plot3d([fx, fy, fz], (u, -pi, pi), (v, -pi, pi), plot_points = [50,50], frame=False, color="red")
    
    A Helicoid (lines through a helix,
    http://en.wikipedia.org/wiki/Helix)::
    
        sage: u, v = var('u,v')
        sage: fx = sinh(v)*sin(u)
        sage: fy = -sinh(v)*cos(u)
        sage: fz = 3*u
        sage: parametric_plot3d([fx, fy, fz], (u, -pi, pi), (v, -pi, pi), plot_points = [50,50], frame=False, color="red")
    
    Kuen's surface
    (http://www.math.umd.edu/research/bianchi/Gifccsurfs/ccsurfs.html)::
    
        sage: fx = (2*(cos(u) + u*sin(u))*sin(v))/(1+ u^2*sin(v)^2)
        sage: fy = (2*(sin(u) - u*cos(u))*sin(v))/(1+ u^2*sin(v)^2)
        sage: fz = log(tan(1/2 *v)) + (2*cos(v))/(1+ u^2*sin(v)^2)
        sage: parametric_plot3d([fx, fy, fz], (u, 0, 2*pi), (v, 0.01, pi-0.01), plot_points = [50,50], frame=False, color="green")
    
    A 5-pointed star::
    
        sage: fx = cos(u)*cos(v)*(abs(cos(1*u/4))^0.5 + abs(sin(1*u/4))^0.5)^(-1/0.3)*(abs(cos(5*v/4))^1.7 + abs(sin(5*v/4))^1.7)^(-1/0.1)
        sage: fy = cos(u)*sin(v)*(abs(cos(1*u/4))^0.5 + abs(sin(1*u/4))^0.5)^(-1/0.3)*(abs(cos(5*v/4))^1.7 + abs(sin(5*v/4))^1.7)^(-1/0.1)
        sage: fz = sin(u)*(abs(cos(1*u/4))^0.5 + abs(sin(1*u/4))^0.5)^(-1/0.3)
        sage: parametric_plot3d([fx, fy, fz], (u, -pi/2, pi/2), (v, 0, 2*pi), plot_points = [50,50], frame=False, color="green")
    
    A cool self-intersecting surface (Eppener surface?)::
    
        sage: fx = u - u^3/3 + u*v^2
        sage: fy = v - v^3/3 + v*u^2
        sage: fz = u^2 - v^2
        sage: parametric_plot3d([fx, fy, fz], (u, -25, 25), (v, -25, 25), plot_points = [50,50], frame=False, color="green")
    
    The breather surface
    (http://en.wikipedia.org/wiki/Breather_surface)::
    
        sage: fx = (2*sqrt(0.84)*cosh(0.4*u)*(-(sqrt(0.84)*cos(v)*cos(sqrt(0.84)*v)) - sin(v)*sin(sqrt(0.84)*v)))/(0.4*((sqrt(0.84)*cosh(0.4*u))^2 + (0.4*sin(sqrt(0.84)*v))^2))
        sage: fy = (2*sqrt(0.84)*cosh(0.4*u)*(-(sqrt(0.84)*sin(v)*cos(sqrt(0.84)*v)) + cos(v)*sin(sqrt(0.84)*v)))/(0.4*((sqrt(0.84)*cosh(0.4*u))^2 + (0.4*sin(sqrt(0.84)*v))^2))
        sage: fz = -u + (2*0.84*cosh(0.4*u)*sinh(0.4*u))/(0.4*((sqrt(0.84)*cosh(0.4*u))^2 + (0.4*sin(sqrt(0.84)*v))^2))
        sage: parametric_plot3d([fx, fy, fz], (u, -13.2, 13.2), (v, -37.4, 37.4), plot_points = [90,90], frame=False, color="green")
    
    TESTS::
    
        sage: u, v = var('u,v')
        sage: plot3d(u^2-v^2, (u, -1, 1), (u, -1, 1))
        Traceback (most recent call last):
        ...
        ValueError: range variables should be distinct, but there are duplicates


    From Trac #2858::
    
        sage: parametric_plot3d((u,-u,v), (u,-10,10),(v,-10,10))
        sage: f(u)=u; g(v)=v^2; parametric_plot3d((g,f,f), (-10,10),(-10,10))

    From Trac #5368::

        sage: x, y = var('x,y')
        sage: plot3d(x*y^2 - sin(x), (x,-1,1), (y,-1,1))

    """
    # TODO:
    #   * Surface -- behavior of functions not defined everywhere -- see note above
    #   * Iterative refinement

    
    # color_function -- (default: "automatic") how to determine the color of curves and surfaces
    # color_function_scaling -- (default: True) whether to scale the input to color_function
    # exclusions -- (default: "automatic") u points or (u,v) conditions to exclude.
    #         (E.g., exclusions could be a function e = lambda u, v: False if u < v else True
    # exclusions_style -- (default: None) what to draw at excluded points
    # max_recursion -- (default: "automatic") maximum number of recursive subdivisions,
    #                   when ...
    # mesh -- (default: "automatic") how many mesh divisions in each direction to draw
    # mesh_functions -- (default: "automatic") how to determine the placement of mesh divisions
    # mesh_shading -- (default: None) how to shade regions between mesh divisions
    # plot_range -- (default: "automatic") range of values to include

    if is_Vector(f):
        f = tuple(f)

    if isinstance(f, (list,tuple)) and len(f) > 0 and isinstance(f[0], (list,tuple)):
        return sum([parametric_plot3d(v, urange, vrange, plot_points=plot_points, **kwds) for v in f])
            
    if not isinstance(f, (tuple, list)) or len(f) != 3:
        raise ValueError, "f must be a list, tuple, or vector of length 3"

    if vrange is None:
        if plot_points == "automatic":
            plot_points = 75
        G = _parametric_plot3d_curve(f, urange, plot_points=plot_points, **kwds)
    else:
        if plot_points == "automatic":
            plot_points = [40,40]
        G = _parametric_plot3d_surface(f, urange, vrange, plot_points=plot_points, boundary_style=boundary_style, **kwds)
    G._set_extra_kwds(kwds)
    return G