def nu1_mu_portrait(n):
    """ returns an encoded scatter plot of the nth roots of unity in the complex plane """
    if n == 1:
        return db.gps_sato_tate.lookup('1.2.1.2.1a').get('trace_histogram')
    if n <= 120:
        plot =  sum([line2d([(-2*cos(2*pi*m/n),-2*sin(2*pi*m/n)),(2*cos(2*pi*m/n),2*sin(2*pi*m/n))],thickness=3) for m in range(n)]) + circle((0,0),0.1,rgbcolor=(0,0,0),fill=True)
    else:
        plot = circle((0,0),2,fill=True)
    plot.xmin(-2); plot.xmax(2); plot.ymin(-2); plot.ymax(2)
    plot.set_aspect_ratio(4.0/3.0)
    plot.axes(False)
    return encode_plot(plot)

def su2_mu_portrait(n):
    """ returns an encoded line plot of SU(2) x mu(n) in the complex plane """
    if n == 1:
        return db.gps_sato_tate.lookup('1.2.3.1.1a').get('trace_histogram')
    if n <= 120:
        plot =  sum([line2d([(-2*cos(2*pi*m/n),-2*sin(2*pi*m/n)),(2*cos(2*pi*m/n),2*sin(2*pi*m/n))],thickness=3) for m in range(n)])
    else:
        plot = circle((0,0),2,fill=True)
    plot.xmin(-2); plot.xmax(2); plot.ymin(-2); plot.ymax(2)
    plot.set_aspect_ratio(4.0/3.0)
    plot.axes(False)
    return encode_plot(plot)

        def fourier_series_partial_sum(cls, self, parameters, variable, N, L):
            r"""
            Returns the partial sum

            .. math::

               f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],

            as a string.

            EXAMPLE::

                sage: f(x) = x^2
                sage: f = piecewise([[(-1,1),f]])
                sage: f.fourier_series_partial_sum(3,1)
                cos(2*pi*x)/pi^2 - 4*cos(pi*x)/pi^2 + 1/3
                sage: f1(x) = -1
                sage: f2(x) = 2
                sage: f = piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
                sage: f.fourier_series_partial_sum(3,pi)
                -3*cos(x)/pi - 3*sin(2*x)/pi + 3*sin(x)/pi - 1/4
            """
            from sage.all import pi, sin, cos, srange
            x = self.default_variable()
            a0 = self.fourier_series_cosine_coefficient(0,L)
            result = a0/2 + sum([(self.fourier_series_cosine_coefficient(n,L)*cos(n*pi*x/L) +
                                  self.fourier_series_sine_coefficient(n,L)*sin(n*pi*x/L))
                                 for n in srange(1,N)])
            return SR(result).expand()

        def fourier_series_sine_coefficient(cls, self, parameters, variable, n, L):
            r"""
            Returns the n-th Fourier series coefficient of
            `\sin(n\pi x/L)`, `b_n`.

            INPUT:


            -  ``self`` - the function f(x), defined over -L x L

            -  ``n`` - an integer n0

            -  ``L`` - (the period)/2


            OUTPUT:
            `b_n = \frac{1}{L}\int_{-L}^L f(x)\sin(n\pi x/L)dx`

            EXAMPLES::

                sage: f(x) = x^2
                sage: f = piecewise([[(-1,1),f]])
                sage: f.fourier_series_sine_coefficient(2,1)  # L=1, n=2
                0
            """
            from sage.all import sin, pi
            x = SR.var('x')
            result = 0
            for domain, f in parameters:
                for interval in domain:
                    a = interval.lower()
                    b = interval.upper()
                    result += (f*sin(pi*x*n/L)/L).integrate(x, a, b)
            return SR(result).simplify_trig()

def get_matrix(matrix, mesh, argyris=False,
               convection_field=(sg.cos(sg.pi/3), sg.sin(sg.pi/3)),
               forcing_function=1, stabilization_coeff=0.0):
    """
    Build a finite element matrix.
    """
    if argyris:
        builder = {'mass': cfem.cf_ap_build_mass,
                   'stiffness': cfem.cf_ap_build_stiffness,
                   'betaplane': cfem.cf_ap_build_betaplane,
                   'biharmonic': cfem.cf_ap_build_biharmonic}[matrix]
        order = 5

    else:
        builder = {'mass': cfem.cf_build_mass,
                   'convection': cfem.cf_build_convection,
                   'stiffness': cfem.cf_build_stiffness,
                   'hessian': cfem.cf_build_hessian}[matrix]
        order = mesh.order
    convection, _, ref_data = _get_lookups(convection_field, forcing_function,
                                           stabilization_coeff, order,
                                           argyris=argyris)
    cmesh, elements = _get_cmesh(mesh)

    length = elements.shape[0]*elements.shape[1]**2
    rows = np.zeros(length, dtype=np.int32) + -99999
    columns = np.zeros(length, dtype=np.int32) + -99999
    values = np.zeros(length, dtype=np.double) + -99999
    triplet_matrix = CTripletMatrix(length=length,
                                    rows=rows.ctypes.data_as(PINT),
                                    columns=columns.ctypes.data_as(PINT),
                                    values=values.ctypes.data_as(PDOUBLE))

    builder(cmesh, ref_data, convection, triplet_matrix)
    return sparse.coo_matrix((values, (rows, columns))).tocsc()

def get_load_vector(mesh, time, argyris=False,
                    convection_field=(sg.cos(sg.pi/3), sg.sin(sg.pi/3)),
                    forcing_function=1, stabilization_coeff=0.0):
    """
    Form a load vector from symbolic expressions for the convection
    field and forcing function. The symbolic expressions are cached.
    """
    if argyris:
        order = 5
    else:
        order = mesh.order
    convection, forcing, ref_data = _get_lookups(
        convection_field, forcing_function, stabilization_coeff, order,
        argyris=argyris)
    cmesh, elements = _get_cmesh(mesh)

    load_vector = np.zeros(mesh.nodes.shape[0])
    cvector = CVector(length=mesh.nodes.shape[0],
                      values=load_vector.ctypes.data_as(PDOUBLE))

    if argyris:
        builder = cfem.cf_ap_build_load
    else:
        builder = cfem.cf_build_load
    builder(cmesh, ref_data, convection, forcing, time, cvector)

    return load_vector

def get_linearization(matrix, mesh, solution, argyris=False,
                      convection_field=(sg.cos(sg.pi/3), sg.sin(sg.pi/3)),
                      forcing_function=1, stabilization_coeff=0.0):
    if not argyris:
        raise NotImplementedError
    order = 5
    builder = {
        0: cfem.cf_ap_build_jacobian_0,
        1: cfem.cf_ap_build_jacobian_1,
    }[matrix]
    convection, _, ref_data = _get_lookups(convection_field, forcing_function,
                                           stabilization_coeff, order,
                                           argyris=True)
    cmesh, elements = _get_cmesh(mesh)

    length = elements.shape[0]*elements.shape[1]**2
    rows = np.zeros(length, dtype=np.int32) + -99999
    columns = np.zeros(length, dtype=np.int32) + -99999
    values = np.zeros(length, dtype=np.double) + -99999
    triplet_matrix = CTripletMatrix(length=length,
                                    rows=rows.ctypes.data_as(PINT),
                                    columns=columns.ctypes.data_as(PINT),
                                    values=values.ctypes.data_as(PDOUBLE))

    assert len(solution) == mesh.nodes.shape[0]
    cvector = CVector(length=mesh.nodes.shape[0],
                      values=solution.ctypes.data_as(PDOUBLE))
    builder(cmesh, ref_data, convection, cvector, triplet_matrix)
    answer = sparse.coo_matrix((values, (rows, columns))).tocsc()
    return answer

def mu_portrait(n):
    """ returns an encoded scatter plot of the nth roots of unity in the complex plane """
    if n <= 120:
        plot =  list_plot([(cos(2*pi*m/n),sin(2*pi*m/n)) for m in range(n)],pointsize=30+60/n, axes=False)
    else:
        plot = circle((0,0),1,thickness=3)
    plot.xmin(-1); plot.xmax(1); plot.ymin(-1); plot.ymax(1)
    plot.set_aspect_ratio(4.0/3.0)
    return encode_plot(plot)

def plotit(k):
    k = int(k[0])
    # FIXME there could be a filename collission
    fn = tempfile.mktemp(suffix=".png")
    x = var('x')
    p = plot(sin(k * x))
    p.save(filename=fn)
    data = file(fn).read()
    os.remove(fn)
    return data

 def _sage_(self):
     import sage.all as sage
     return sage.sin(self.args[0]._sage_())

        def fourier_series_partial_sum(self, parameters, variable, N,
                                       L=None):
            r"""
            Returns the partial sum up to a given order of the Fourier series
            of the periodic function `f` extending the piecewise-defined
            function ``self``.

            The Fourier partial sum of order `N` is defined as

            .. MATH::

                S_{N}(x) = \frac{a_0}{2} + \sum_{n=1}^{N} \left[
                      a_n\cos\left(\frac{n\pi x}{L}\right)
                    + b_n\sin\left(\frac{n\pi x}{L}\right)\right],

            where `L` is the half-period of `f` and the `a_n`'s and `b_n`'s
            are respectively the cosine coefficients and sine coefficients
            of the Fourier series of `f` (cf.
            :meth:`fourier_series_cosine_coefficient` and
            :meth:`fourier_series_sine_coefficient`).

            INPUT:

            - ``N`` -- a positive integer; the order of the partial sum

            - ``L`` -- (default: ``None``) the half-period of `f`; if none
              is provided, `L` is assumed to be the half-width of the domain
              of ``self``

            OUTPUT:

            - the partial sum `S_{N}(x)`, as a symbolic expression

            EXAMPLES:

            A square wave function of period 2::

                sage: f = piecewise([((-1,0), -1), ((0,1), 1)])
                sage: f.fourier_series_partial_sum(5)
                4/5*sin(5*pi*x)/pi + 4/3*sin(3*pi*x)/pi + 4*sin(pi*x)/pi

            If the domain of the piecewise-defined function encompasses
            more than one period, the half-period must be passed as the
            second argument; for instance::

                sage: f2 = piecewise([((-1,0), -1), ((0,1), 1),
                ....:                 ((1,2), -1), ((2,3), 1)])
                sage: bool(f2.restriction((-1,1)) == f)  # f2 extends f on (-1,3)
                True
                sage: f2.fourier_series_partial_sum(5, 1)  # half-period = 1
                4/5*sin(5*pi*x)/pi + 4/3*sin(3*pi*x)/pi + 4*sin(pi*x)/pi
                sage: bool(f2.fourier_series_partial_sum(5, 1) ==
                ....:      f.fourier_series_partial_sum(5))
                True

            The default half-period is 2, so that skipping the second
            argument yields a different result::

                sage: f2.fourier_series_partial_sum(5)  # half-period = 2
                4*sin(pi*x)/pi

            An example of partial sum involving both cosine and sine terms::

                sage: f = piecewise([((-1,0), 0), ((0,1/2), 2*x),
                ....:                ((1/2,1), 2*(1-x))])
                sage: f.fourier_series_partial_sum(5)
                -2*cos(2*pi*x)/pi^2 + 4/25*sin(5*pi*x)/pi^2
                 - 4/9*sin(3*pi*x)/pi^2 + 4*sin(pi*x)/pi^2 + 1/4

            """
            from sage.all import pi, sin, cos, srange
            if not L:
                L = (self.domain().sup() - self.domain().inf()) / 2
            x = self.default_variable()
            a0 = self.fourier_series_cosine_coefficient(0, L)
            result = a0/2 + sum([(self.fourier_series_cosine_coefficient(n, L)*cos(n*pi*x/L) +
                                  self.fourier_series_sine_coefficient(n, L)*sin(n*pi*x/L))
                                 for n in srange(1, N+1)])
            return SR(result).expand()

        def fourier_series_sine_coefficient(self, parameters, variable,
                                            n, L=None):
            r"""
            Return the `n`-th sine coefficient of the Fourier series of
            the periodic function `f` extending the piecewise-defined
            function ``self``.

            Given an integer `n\geq 0`, the `n`-th sine coefficient of
            the Fourier series of `f` is defined by

            .. MATH::

                b_n = \frac{1}{L}\int_{-L}^L
                        f(x)\sin\left(\frac{n\pi x}{L}\right) dx,

            where `L` is the half-period of `f`. The number `b_n` is
            the coefficient of `\sin(n\pi x/L)` in the Fourier
            series of `f` (cf. :meth:`fourier_series_partial_sum`).

            INPUT:

            - ``n`` -- a non-negative integer

            - ``L`` -- (default: ``None``) the half-period of `f`; if none
              is provided, `L` is assumed to be the half-width of the domain
              of ``self``

            OUTPUT:

            - the Fourier coefficient `b_n`, as defined above

            EXAMPLES:

            A square wave function of period 2::

                sage: f = piecewise([((-1,0), -1), ((0,1), 1)])
                sage: f.fourier_series_sine_coefficient(1)
                4/pi
                sage: f.fourier_series_sine_coefficient(2)
                0
                sage: f.fourier_series_sine_coefficient(3)
                4/3/pi

            If the domain of the piecewise-defined function encompasses
            more than one period, the half-period must be passed as the
            second argument; for instance::

                sage: f2 = piecewise([((-1,0), -1), ((0,1), 1),
                ....:                 ((1,2), -1), ((2,3), 1)])
                sage: bool(f2.restriction((-1,1)) == f)  # f2 extends f on (-1,3)
                True
                sage: f2.fourier_series_sine_coefficient(1, 1)  # half-period = 1
                4/pi
                sage: f2.fourier_series_sine_coefficient(3, 1)  # half-period = 1
                4/3/pi

            The default half-period is 2 and one has::

                sage: f2.fourier_series_sine_coefficient(1)  # half-period = 2
                0
                sage: f2.fourier_series_sine_coefficient(3)  # half-period = 2
                0

            The Fourier coefficients obtained from ``f`` are actually
            recovered for `n=2` and `n=6` respectively::

                sage: f2.fourier_series_sine_coefficient(2)
                4/pi
                sage: f2.fourier_series_sine_coefficient(6)
                4/3/pi

            """
            from sage.all import sin, pi
            L0 = (self.domain().sup() - self.domain().inf()) / 2
            if not L:
                L = L0
            else:
                m = L0 / L
                if not (m.is_integer() and m > 0):
                    raise ValueError("the width of the domain of " +
                                     "{} is not a multiple ".format(self) +
                                     "of the given period")
            x = SR.var('x')
            result = 0
            for domain, f in parameters:
                for interval in domain:
                    a = interval.lower()
                    b = interval.upper()
                    result += (f*sin(pi*x*n/L)).integrate(x, a, b)
            return SR(result/L0).simplify_trig()

def test_sage_conversions():
    try:
        import sage.all as sage
    except ImportError:
        return

    x, y = sage.SR.var('x y')
    x1, y1 = symbols('x, y')

    # Symbol
    assert x1._sage_() == x
    assert x1 == sympify(x)

    # Integer
    assert Integer(12)._sage_() == sage.Integer(12)
    assert Integer(12) == sympify(sage.Integer(12))

    # Rational
    assert (Integer(1) / 2)._sage_() == sage.Integer(1) / 2
    assert Integer(1) / 2 == sympify(sage.Integer(1) / 2)

    # Operators
    assert x1 + y == x1 + y1
    assert x1 * y == x1 * y1
    assert x1 ** y == x1 ** y1
    assert x1 - y == x1 - y1
    assert x1 / y == x1 / y1

    assert x + y1 == x + y
    assert x * y1 == x * y
    # Doesn't work in Sage 6.1.1ubuntu2
    # assert x ** y1 == x ** y
    assert x - y1 == x - y
    assert x / y1 == x / y

    # Conversions
    assert (x1 + y1)._sage_() == x + y
    assert (x1 * y1)._sage_() == x * y
    assert (x1 ** y1)._sage_() == x ** y
    assert (x1 - y1)._sage_() == x - y
    assert (x1 / y1)._sage_() == x / y

    assert x1 + y1 == sympify(x + y)
    assert x1 * y1 == sympify(x * y)
    assert x1 ** y1 == sympify(x ** y)
    assert x1 - y1 == sympify(x - y)
    assert x1 / y1 == sympify(x / y)

    # Functions
    assert sin(x1) == sin(x)
    assert sin(x1)._sage_() == sage.sin(x)
    assert sin(x1) == sympify(sage.sin(x))

    assert cos(x1) == cos(x)
    assert cos(x1)._sage_() == sage.cos(x)
    assert cos(x1) == sympify(sage.cos(x))

    assert function_symbol('f', x1, y1)._sage_() == sage.function('f', x, y)
    assert function_symbol('f', 2 * x1, x1 + y1).diff(x1)._sage_() == sage.function('f', 2 * x, x + y).diff(x)

    # For the following test, sage needs to be modified
    # assert sage.sin(x) == sage.sin(x1)

    # Constants
    assert pi._sage_() == sage.pi
    assert E._sage_() == sage.e
    assert I._sage_() == sage.I

    assert pi == sympify(sage.pi)
    assert E == sympify(sage.e)

    # SympyConverter does not support converting the following
    # assert I == sympify(sage.I)

    # Matrix
    assert DenseMatrix(1, 2, [x1, y1])._sage_() == sage.matrix([[x, y]])

    # SympyConverter does not support converting the following
    # assert DenseMatrix(1, 2, [x1, y1]) == sympify(sage.matrix([[x, y]]))

    # Sage Number
    a = sage.Mod(2, 7)
    b = PyNumber(a, sage_module)

    a = a + 8
    b = b + 8
    assert isinstance(b, PyNumber)
    assert b._sage_() == a

    a = a + x
    b = b + x
    assert isinstance(b, Add)
    assert b._sage_() == a

    # Sage Function
    e = x1 + wrap_sage_function(sage.log_gamma(x))
    assert str(e) == "x + log_gamma(x)"
    assert isinstance(e, Add)
    assert e + wrap_sage_function(sage.log_gamma(x)) == x1 + 2*wrap_sage_function(sage.log_gamma(x))

#.........这里部分代码省略.........

def test_sage_conversions():
    try:
        import sage.all as sage
    except ImportError:
        return

    x, y = sage.SR.var('x y')
    x1, y1 = symbols('x, y')

    # Symbol
    assert x1._sage_() == x
    assert x1 == sympify(x)

    # Integer
    assert Integer(12)._sage_() == sage.Integer(12)
    assert Integer(12) == sympify(sage.Integer(12))

    # Rational
    assert (Integer(1) / 2)._sage_() == sage.Integer(1) / 2
    assert Integer(1) / 2 == sympify(sage.Integer(1) / 2)

    # Operators
    assert x1 + y == x1 + y1
    assert x1 * y == x1 * y1
    assert x1 ** y == x1 ** y1
    assert x1 - y == x1 - y1
    assert x1 / y == x1 / y1

    assert x + y1 == x + y
    assert x * y1 == x * y
    # Doesn't work in Sage 6.1.1ubuntu2
    # assert x ** y1 == x ** y
    assert x - y1 == x - y
    assert x / y1 == x / y

    # Conversions
    assert (x1 + y1)._sage_() == x + y
    assert (x1 * y1)._sage_() == x * y
    assert (x1 ** y1)._sage_() == x ** y
    assert (x1 - y1)._sage_() == x - y
    assert (x1 / y1)._sage_() == x / y

    assert x1 + y1 == sympify(x + y)
    assert x1 * y1 == sympify(x * y)
    assert x1 ** y1 == sympify(x ** y)
    assert x1 - y1 == sympify(x - y)
    assert x1 / y1 == sympify(x / y)

    # Functions
    assert sin(x1) == sin(x)
    assert sin(x1)._sage_() == sage.sin(x)
    assert sin(x1) == sympify(sage.sin(x))

    assert cos(x1) == cos(x)
    assert cos(x1)._sage_() == sage.cos(x)
    assert cos(x1) == sympify(sage.cos(x))

    assert function_symbol('f', x1, y1)._sage_() == sage.function('f', x, y)
    assert function_symbol('f', 2 * x1, x1 + y1).diff(x1)._sage_() == sage.function('f', 2 * x, x + y).diff(x)

    # For the following test, sage needs to be modified
    # assert sage.sin(x) == sage.sin(x1)

    # Constants
    assert pi._sage_() == sage.pi
    assert E._sage_() == sage.e
    assert I._sage_() == sage.I

    assert pi == sympify(sage.pi)
    assert E == sympify(sage.e)

    # SympyConverter does not support converting the following
    # assert I == sympify(sage.I)

    # Matrix
    assert DenseMatrix(1, 2, [x1, y1])._sage_() == sage.matrix([[x, y]])