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

        EXAMPLES::

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

        Infinite geodesics cannot be bisected::

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

    def approximate(self, x, parent=None):
        r"""
        Approximate using de Bruijn's formula

        .. math::

             \rho(x) \sim \frac{exp(-x \xi + Ei(\xi))}{\sqrt{2\pi x}\xi}

        which is asymptotically equal to Dickman's function, and is much
        faster to compute.

        REFERENCES:

        - N. De Bruijn, "The Asymptotic behavior of a function
          occurring in the theory of primes." J. Indian Math Soc. v 15.
          (1951)

        EXAMPLES::

            sage: dickman_rho.approximate(10)
            2.41739196365564e-11
            sage: dickman_rho(10)
            2.77017183772596e-11
            sage: dickman_rho.approximate(1000)
            4.32938809066403e-3464
        """
        log, exp, sqrt, pi = math.log, math.exp, math.sqrt, math.pi
        x = float(x)
        xi = log(x)
        y = (exp(xi)-1.0)/xi - x
        while abs(y) > 1e-12:
            dydxi = (exp(xi)*(xi-1.0) + 1.0)/(xi*xi)
            xi -= y/dydxi
            y = (exp(xi)-1.0)/xi - x
        return (-x*xi + RR(xi).eint()).exp() / (sqrt(2*pi*x)*xi)

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

            sage: gamma_inc(2.,0)
            1.00000000000000
            sage: gamma_inc(2,0)
            1
            sage: gamma_inc(1/2,2)
            -(erf(sqrt(2)) - 1)*sqrt(pi)
            sage: gamma_inc(1/2,1)
            -(erf(1) - 1)*sqrt(pi)
            sage: gamma_inc(1/2,0)
            sqrt(pi)
            sage: gamma_inc(x,0)
            gamma(x)
            sage: gamma_inc(1,2)
            e^(-2)
            sage: gamma_inc(0,2)
            -Ei(-2)
        """
        if not isinstance(x, Expression) and not isinstance(y, Expression) and \
               (is_inexact(x) or is_inexact(y)):
            x, y = coercion_model.canonical_coercion(x, y)
            return self._evalf_(x, y, parent(x))

        if y == 0:
            return gamma(x)
        if x == 1:
            return exp(-y)
        if x == 0:
            return -Ei(-y)
        if x == Rational(1)/2: #only for x>0
            return sqrt(pi)*(1-erf(sqrt(y)))
        return None

    def partition_function(self, beta, epsilon):
        r"""
        Return the partition function of ``self``.

        The partition function of a 6 vertex model is defined by:

        .. MATH::

            Z = \sum_{\nu} e^{-\beta E(\nu)}

        where we sum over all configurations and `E` is the energy function.
        The constant `\beta` is known as the *inverse temperature* and is
        equal to `1 / k_B T` where `k_B` is Boltzmann's constant and `T` is
        the system's temperature.

        INPUT:

        - ``beta`` -- the inverse temperature constant `\beta`
        - ``epsilon`` -- the energy constants, see
          :meth:`~sage.combinat.six_vertex_model.SixVertexConfiguration.energy()`

        EXAMPLES::

            sage: M = SixVertexModel(3, boundary_conditions='ice')
            sage: M.partition_function(2, [1,2,1,2,1,2])
            e^(-24) + 2*e^(-28) + e^(-30) + 2*e^(-32) + e^(-36)

        REFERENCES:

        :wikipedia:`Partition_function_(statistical_mechanics)`
        """
        from sage.functions.log import exp
        return sum(exp(-beta * nu.energy(epsilon)) for nu in self)

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

            sage: exp_integral_e(1.0, x)
            exp_integral_e(1.00000000000000, x)
            sage: exp_integral_e(x, 1.0)
            exp_integral_e(x, 1.00000000000000)
            sage: exp_integral_e(3, 0)
            1/2

        TESTS:

        Check that Python ints work (:trac:`14766`)::

            sage: exp_integral_e(int(3), 0)
            1/2
        """
        z_zero = False
        # special case: z == 0 and n > 1
        if isinstance(z, Expression):
            if z.is_trivial_zero():
                z_zero = True # for later
                if n > 1:
                    return 1/(n-1)
        else:
            if not z:
                z_zero = True
                if n > 1:
                    return 1/(n-1)

        # special case: n == 0
        if isinstance(n, Expression):
            if n.is_trivial_zero():
                if z_zero:
                    return None
                else:
                    return exp(-z)/z
        else:
            if not n:
                if z_zero:
                    return None
                else:
                    return exp(-z)/z

        return None # leaves the expression unevaluated

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

            sage: exp_integral_e(1.0, x)
            exp_integral_e(1.00000000000000, x)
            sage: exp_integral_e(x, 1.0)
            exp_integral_e(x, 1.00000000000000)
            sage: exp_integral_e(1.0, 1.0)
            0.219383934395520

        """
        if not isinstance(n, Expression) and not isinstance(z, Expression) and \
               (is_inexact(n) or is_inexact(z)):
            coercion_model = sage.structure.element.get_coercion_model()
            n, z = coercion_model.canonical_coercion(n, z)
            return self._evalf_(n, z, parent(n))

        z_zero = False
        # special case: z == 0 and n > 1
        if isinstance(z, Expression):
            if z.is_trivial_zero():
                z_zero = True # for later
                if n > 1:
                    return 1/(n-1)
        else:
            if not z:
                z_zero = True
                if n > 1:
                    return 1/(n-1)

        # special case: n == 0
        if isinstance(n, Expression):
            if n.is_trivial_zero():
                if z_zero:
                    return None
                else:
                    return exp(-z)/z
        else:
            if not n:
                if z_zero:
                    return None
                else:
                    return exp(-z)/z

        return None # leaves the expression unevaluated

    def xseries(self, all_conjugates=True):
        r"""Returns the corresponding x-series.

        Parameters
        ----------
        all_conjugates : bool
            (default: True) If ``True``, returns all conjugates
            x-representations of this Puiseux t-series. If ``False``,
            only returns one representative.

        Returns
        -------
        list
            List of PuiseuxXSeries representations of this PuiseuxTSeries.

        """
        # obtain relevant rings:
        #   o R = parent ring of curve
        #   o L = parent ring of T-series
        #   o S = temporary polynomial ring over base ring of T-series
        #   o P = Puiseux series ring
        L = self.ypart.parent()
        t = L.gen()
        S = L.base_ring()['z']
        z = S.gen()

        R = self.f.parent()
        x,y = R.gens()
        P = PuiseuxSeriesRing(L.base_ring(), str(x))
        x = P.gen()

        # given x = alpha + lambda*t^e solve for t. this involves finding an
        # e-th root of either (1/lambda) or of lambda, depending on e's sign
        e = self.ramification_index
        lamb = S(self.xcoefficient)
        order = self.order
        if e > 0:
            phi = lamb*z**e - 1
        else:
            phi = z**abs(e) - lamb
        mu = phi.roots(QQbar, multiplicities=False)[0]

        if all_conjugates:
            conjugates = [mu*exp(2*pi*I*k/abs(e)) for k in range(abs(e))]
        else:
            conjugates = [mu]
        map(lambda x: x.exactify(), conjugates)

        # determine the resulting x-series
        xseries = []
        for c in conjugates:
            t = self.ypart.parent().gen()
            fconj = self.ypart(c*t)
            p = P(fconj(x**(QQ(1)/e)))
            p = p.add_bigoh(QQ(order+1)/abs(e))
            xseries.append(p)
        return xseries

def circle_image(A,B):
    G = Graphics()
    G += circle((0,0), 1 , color = 'grey')
    from collections import defaultdict
    tmp = defaultdict(int)
    for a in A:
        for j in range(a):
            if gcd(j,a) == 1:
                rational = Rational(j)/Rational(a)
                tmp[(rational.numerator(),rational.denominator())] += 1
    
    for b in B:
        for j in range(b):
            if gcd(j,b) == 1:
                rational = Rational(j)/Rational(b)
                tmp[(rational.numerator(),rational.denominator())] -= 1
    C = ComplexField()
    for val in tmp:
        if tmp[val] > 0:
            G += text(str(tmp[val]),exp(C(2*3.14159*I*val[0]/val[1])), fontsize = 30, axes = False, color = "green")
        if tmp[val] < 0:
            G += text(str(abs(tmp[val])),exp(C(2*3.14159*I*val[0]/val[1])), fontsize = 30, axes = False, color = "red")
    return G

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

            sage: bessel_K(1,0)
            bessel_K(1, 0)
            sage: bessel_K(1.0, 0.0)
            +infinity
            sage: bessel_K(-1, 1).n(128)
            0.60190723019723457473754000153561733926
        """
        # special identity
        if n == Integer(1) / Integer(2) and x > 0:
            return sqrt(pi / 2) * exp(-x) * x ** (-Integer(1) / Integer(2))

    def _derivative_(self, x, diff_param=None):
        """
        EXAMPLES::

            sage: Ei(x).diff(x)
            e^x/x
            sage: Ei(x).diff(x).subs(x=1)
            e
            sage: Ei(x^2).diff(x)
            2*e^(x^2)/x
            sage: f = function('f')
            sage: Ei(f(x)).diff(x)
            e^f(x)*D[0](f)(x)/f(x)
        """
        return exp(x)/x

    def __getattr__(self, name):
        """
        EXAMPLE::

            sage: from sage.crypto.lwe import DiscreteGaussianSamplerRejection
            sage: DiscreteGaussianSamplerRejection(3.0).foo
            Traceback (most recent call last):
            ...
            AttributeError: 'DiscreteGaussianSamplerRejection' object has no attribute 'foo'
        """
        if name == "rho":
            # we delay the creation of rho until we actually need it
            R = RealField(self.precision)
            self.rho = [round(self.max_precs * exp((-(R(x) / R(self.stddev))**2)/R(2))) for x in range(0,self.upper_bound)]
            self.rho[0] = self.rho[0] / 2
            return self.rho
        else:
            raise AttributeError("'%s' object has no attribute '%s'"%(self.__class__.__name__, name))

    def _derivative_(self, x, diff_param=None):
        """
        Derivative of erf function

        EXAMPLES::

            sage: erf(x).diff(x)
            2*e^(-x^2)/sqrt(pi)

        TESTS::

        Check if #8568 is fixed::

            sage: var('c,x')
            (c, x)
            sage: derivative(erf(c*x),x)
            2*c*e^(-c^2*x^2)/sqrt(pi)
            sage: erf(c*x).diff(x)._maxima_init_()
            '((%pi)^(-1/2))*(c)*(exp(((c)^(2))*((x)^(2))*(-1)))*(2)'
        """
        return 2*exp(-x**2)/sqrt(pi)

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

            sage: bessel_K(1,0)
            bessel_K(1, 0)
            sage: bessel_K(1.0, 0.0)
            +infinity
            sage: bessel_K(-1, 1).n(128)
            0.60190723019723457473754000153561733926
        """
        if (not isinstance(n, Expression) and not isinstance(x, Expression) and
                (is_inexact(n) or is_inexact(x))):
            coercion_model = get_coercion_model()
            n, x = coercion_model.canonical_coercion(n, x)
            return self._evalf_(n, x, parent(n))

        # special identity
        if n == Integer(1) / Integer(2) and x > 0:
            return sqrt(pi / 2) * exp(-x) * x ** (-Integer(1) / Integer(2))

        return None  # leaves the expression unevaluated

    def midpoint(self):  # UHP
        r"""
        Return the (hyperbolic) midpoint of ``self`` if it exists.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: g = UHP.random_geodesic()
            sage: m = g.midpoint()
            sage: d1 = UHP.dist(m, g.start())
            sage: d2 = UHP.dist(m, g.end())
            sage: bool(abs(d1 - d2) < 10**-9)
            True

        Infinite geodesics have no midpoint::

            sage: UHP.get_geodesic(0, 2).midpoint()
            Traceback (most recent call last):
            ...
            ValueError: the length must be finite
        """
        if self.length() == infinity:
            raise ValueError("the length must be finite")

        start = self._start.coordinates()
        end = self._end.coordinates()
        d = self._model._dist_points(start, end) / 2
        S = self.complete()._to_std_geod(start)
        T = matrix([[exp(d), 0], [0, 1]])
        M = S.inverse() * T * S
        if ((real(start - end) < EPSILON)
                or (abs(real(start - end)) < EPSILON
                    and imag(start - end) < EPSILON)):
            end_p = start
        else:
            end_p = end
        return self._model.get_point(mobius_transform(M, end_p))

def subexpressions_list(f, pars=None):
    """
    Construct the lists with the intermediate steps on the evaluation of the
    function.

    INPUT:

    - ``f`` -- a symbolic function of several components.

    - ``pars`` -- a list of the parameters that appear in the function
      this should be the symbolic constants that appear in f but are not
      arguments.

    OUTPUT:

    - a list of the intermediate subexpressions that appear in the evaluation
      of f.

    - a list with the operations used to construct each of the subexpressions.
      each element of this list is a tuple, formed by a string describing the
      operation made, and the operands.

    For the trigonometric functions, some extra expressions will be added.
    These extra expressions will be used later to compute their derivatives.


    EXAMPLES::

        sage: from sage.interfaces.tides import subexpressions_list
        sage: var('x,y')
        (x, y)
        sage: f(x,y) = [x^2+y, cos(x)/log(y)]
        sage: subexpressions_list(f)
        ([x^2, x^2 + y, sin(x), cos(x), log(y), cos(x)/log(y)],
        [('mul', x, x),
        ('add', y, x^2),
        ('sin', x),
        ('cos', x),
        ('log', y),
        ('div', log(y), cos(x))])

    ::

        sage: f(a)=[cos(a), arctan(a)]
        sage: from sage.interfaces.tides import subexpressions_list
        sage: subexpressions_list(f)
        ([sin(a), cos(a), a^2, a^2 + 1, arctan(a)],
        [('sin', a), ('cos', a), ('mul', a, a), ('add', 1, a^2), ('atan', a)])

    ::

        sage: from sage.interfaces.tides import subexpressions_list
        sage: var('s,b,r')
        (s, b, r)
        sage: f(t,x,y,z)= [s*(y-x),x*(r-z)-y,x*y-b*z]
        sage: subexpressions_list(f,[s,b,r])
        ([-y,
        x - y,
        s*(x - y),
        -s*(x - y),
        -z,
        r - z,
        (r - z)*x,
        -y,
        (r - z)*x - y,
        x*y,
        b*z,
        -b*z,
        x*y - b*z],
        [('mul', -1, y),
        ('add', -y, x),
        ('mul', x - y, s),
        ('mul', -1, s*(x - y)),
        ('mul', -1, z),
        ('add', -z, r),
        ('mul', x, r - z),
        ('mul', -1, y),
        ('add', -y, (r - z)*x),
        ('mul', y, x),
        ('mul', z, b),
        ('mul', -1, b*z),
        ('add', -b*z, x*y)])

    ::

        sage: var('x, y')
        (x, y)
        sage: f(x,y)=[exp(x^2+sin(y))]
        sage: from sage.interfaces.tides import *
        sage: subexpressions_list(f)
        ([x^2, sin(y), cos(y), x^2 + sin(y), e^(x^2 + sin(y))],
        [('mul', x, x),
        ('sin', y),
        ('cos', y),
        ('add', sin(y), x^2),
        ('exp', x^2 + sin(y))])


    """
    from sage.functions.trig import sin, cos, arcsin, arctan, arccos
#.........这里部分代码省略.........

    def __init__(self, B, sigma=1, c=None, precision=None):
        r"""
        Construct a discrete Gaussian sampler over the lattice `Λ(B)`
        with parameter ``sigma`` and center `c`.

        INPUT:

        - ``B`` -- a basis for the lattice, one of the following:

          - an integer matrix,
          - an object with a ``matrix()`` method, e.g. ``ZZ^n``, or
          - an object where ``matrix(B)`` succeeds, e.g. a list of vectors.

        - ``sigma`` -- Gaussian parameter `σ>0`.
        - ``c`` -- center `c`, any vector in `\ZZ^n` is supported, but `c ∈ Λ(B)` is faster.
        - ``precision`` -- bit precision `≥ 53`.

        EXAMPLE::

            sage: from sage.stats.distributions.discrete_gaussian_lattice import DiscreteGaussianDistributionLatticeSampler
            sage: n = 2; sigma = 3.0; m = 5000
            sage: D = DiscreteGaussianDistributionLatticeSampler(ZZ^n, sigma)
            sage: f = D.f
            sage: c = D._normalisation_factor_zz(); c
            56.2162803067524

            sage: l = [D() for _ in xrange(m)]
            sage: v = vector(ZZ, n, (-3,-3))
            sage: l.count(v), ZZ(round(m*f(v)/c))
            (39, 33)

            sage: target = vector(ZZ, n, (0,0))
            sage: l.count(target), ZZ(round(m*f(target)/c))
            (116, 89)

            sage: from sage.stats.distributions.discrete_gaussian_lattice import DiscreteGaussianDistributionLatticeSampler
            sage: qf = QuadraticForm(matrix(3, [2, 1, 1,  1, 2, 1,  1, 1, 2]))
            sage: D = DiscreteGaussianDistributionLatticeSampler(qf, 3.0); D
            Discrete Gaussian sampler with σ = 3.000000, c=(0, 0, 0) over lattice with basis
            <BLANKLINE>
            [2 1 1]
            [1 2 1]
            [1 1 2]
            sage: D()
            (0, 1, -1)
        """
        precision = DiscreteGaussianDistributionLatticeSampler.compute_precision(precision, sigma)

        self._RR = RealField(precision)
        self._sigma = self._RR(sigma)

        try:
            B = matrix(B)
        except ValueError:
            pass

        try:
            B = B.matrix()
        except AttributeError:
            pass

        self.B = B
        self._G = B.gram_schmidt()[0]

        try:
            c = vector(ZZ, B.ncols(), c)
        except TypeError:
            try:
                c = vector(QQ, B.ncols(), c)
            except TypeError:
                c = vector(RR, B.ncols(), c)

        self._c = c

        self.f = lambda x: exp(-(vector(ZZ, B.ncols(), x)-c).norm()**2/(2*self._sigma**2))

        # deal with trivial case first, it is common
        if self._G == 1 and self._c == 0:
            self._c_in_lattice = True
            D = DiscreteGaussianDistributionIntegerSampler(sigma=sigma)
            self.D = tuple([D for _ in range(self.B.nrows())])
            self.VS = FreeModule(ZZ, B.nrows())
            return

        w = B.solve_left(c)
        if w in ZZ**B.nrows():
            self._c_in_lattice = True
            D = []
            for i in range(self.B.nrows()):
                sigma_ = self._sigma/self._G[i].norm()
                D.append( DiscreteGaussianDistributionIntegerSampler(sigma=sigma_) )
            self.D = tuple(D)
            self.VS = FreeModule(ZZ, B.nrows())
        else:
            self._c_in_lattice = False

    def Stirling(var, precision=None, skip_constant_factor=False):
        r"""
        Return Stirling's approximation formula for factorials.

        INPUT:

        - ``var`` -- a string for the variable name.

        - ``precision`` -- (default: ``None``) an integer `\ge 3`. If ``None``, then
          the default precision of the asymptotic ring is used.

        - ``skip_constant_factor`` -- (default: ``False``) a
          boolean. If set, then the constant factor `\sqrt{2\pi}` is left out.
          As a consequence, the coefficient ring of the output changes
          from ``Symbolic Constants Subring`` (if ``False``) to
          ``Rational Field`` (if ``True``).

        OUTPUT:

        An asymptotic expansion.

        EXAMPLES::

            sage: asymptotic_expansions.Stirling('n', precision=5)
            sqrt(2)*sqrt(pi)*e^(n*log(n))*(e^n)^(-1)*n^(1/2) +
            1/12*sqrt(2)*sqrt(pi)*e^(n*log(n))*(e^n)^(-1)*n^(-1/2) +
            1/288*sqrt(2)*sqrt(pi)*e^(n*log(n))*(e^n)^(-1)*n^(-3/2) +
            O(e^(n*log(n))*(e^n)^(-1)*n^(-5/2))
            sage: _.parent()
            Asymptotic Ring <(e^(n*log(n)))^QQ * (e^n)^QQ * n^QQ * log(n)^QQ>
            over Symbolic Constants Subring

        .. SEEALSO::

            :meth:`log_Stirling`,
            :meth:`~sage.rings.asymptotic.asymptotic_ring.AsymptoticExpansion.factorial`.

        TESTS::

            sage: expansion = asymptotic_expansions.Stirling('n', precision=5)
            sage: n = expansion.parent().gen()
            sage: expansion.compare_with_values(n, lambda x: x.factorial(), [5, 10, 20])  # rel tol 1e-6
            [(5, 0.00675841118?), (10, 0.0067589306?), (20, 0.006744925?)]
            sage: asymptotic_expansions.Stirling('n', precision=5,
            ....:                                skip_constant_factor=True)
            e^(n*log(n))*(e^n)^(-1)*n^(1/2) +
            1/12*e^(n*log(n))*(e^n)^(-1)*n^(-1/2) +
            1/288*e^(n*log(n))*(e^n)^(-1)*n^(-3/2) +
            O(e^(n*log(n))*(e^n)^(-1)*n^(-5/2))
            sage: _.parent()
            Asymptotic Ring <(e^(n*log(n)))^QQ * (e^n)^QQ * n^QQ * log(n)^QQ>
            over Rational Field
            sage: asymptotic_expansions.Stirling('m', precision=4)
            sqrt(2)*sqrt(pi)*e^(m*log(m))*(e^m)^(-1)*m^(1/2) +
            O(e^(m*log(m))*(e^m)^(-1)*m^(-1/2))
            sage: asymptotic_expansions.Stirling('m', precision=3)
            O(e^(m*log(m))*(e^m)^(-1)*m^(1/2))
            sage: asymptotic_expansions.Stirling('m', precision=2)
            Traceback (most recent call last):
            ...
            ValueError: precision must be at least 3
        """
        if precision < 3:
            raise ValueError("precision must be at least 3")
        log_Stirling = AsymptoticExpansionGenerators.log_Stirling(
            var, precision=precision, skip_constant_summand=True)

        P = log_Stirling.parent().change_parameter(
            growth_group='(e^({n}*log({n})))^QQ * (e^{n})^QQ * {n}^QQ * log({n})^QQ'.format(n=var))
        from sage.functions.log import exp
        result = exp(P(log_Stirling))

        if not skip_constant_factor:
            from sage.symbolic.ring import SR
            SCR = SR.subring(no_variables=True)
            result *= (2*SCR('pi')).sqrt()

        return result

#.........这里部分代码省略.........
            dList += [identity_matrix(n).column(i) for i in range(n)]

    elif directions['select'] == 'oct':

        if n != 2:
            raise NotImplementedError('Directions select octagon not implemented for n other than 2. Try box.')

        theta = [i*pi/4 for i in range(8)] # octagon
        dList = [vector(RR,[cos(t), sin(t)]) for t in theta]

    elif directions['select'] == 'random':

        order = directions['order'] if 'order' in directions else 12

        if n == 2:
            theta = [random.uniform(0, 2*pi.n(digits=5)) for i in range(order)]
            dList = [vector(RR,[cos(theta[i]), sin(theta[i])]) for i in range(order)]
        else:
            raise NotImplementedError('Directions select random not implemented for n greater than 2. Try box.')

    elif directions['select'] == 'custom':

        dList = directions['dList']

    else:

        raise TypeError('Template directions not understood.')


    # transform directions to numpy array, and get number of directions
    dArray = np.array(dList)
    k = len(dArray)

    global Phi_tau, expX0, alpha_tau_B

    if got_homogeneous: # dx/dx = Ax

        # #######################################################
        # Compute first element of the approximating sequence   #
        # #######################################################

        # compute matrix exponential exp(A*tau)
        Phi_tau = expm(np.multiply(A, tau))

        # compute exp(tau*A)X0
        expX0 = Phi_tau * X0

        # compute the bloating factor
        Ainfty = A.norm(Infinity)
        RX0 = radius(X0)

        unitBall = BoxInfty(center = zero_vector(n), radius = 1, base_ring = base_ring)
        alpha_tau = (exp(tau*Ainfty) - 1 - tau*Ainfty)*(RX0)
        alpha_tau_B = (alpha_tau*np.identity(n)) * unitBall

        # now we have that:
        # Omega0 = X0.convex_hull(expX0.Minkowski_sum(alpha_tau_B))

        # compute the first element of the approximating sequence, Omega_0
        #if first_element_evaluation == 'exact':
        #    Omega0 = X0.convex_hull(expX0.Minkowski_sum(alpha_tau_B))

        #elif first_element_evaluation == 'approximate': # NOT TESTED!!!
            #Omega0_A = dArray
        #    Omega0_b = np.zeros(k)

    def _closed_form(hyp):
        a, b, z = hyp.operands()
        a, b = a.operands(), b.operands()
        p, q = len(a), len(b)

        if z == 0:
            return Integer(1)
        if p == q == 0:
            return exp(z)
        if p == 1 and q == 0:
            return (1 - z) ** (-a[0])

        if p == 0 and q == 1:
            # TODO: make this require only linear time
            def _0f1(b, z):
                F12 = cosh(2 * sqrt(z))
                F32 = sinh(2 * sqrt(z)) / (2 * sqrt(z))
                if 2 * b == 1:
                    return F12
                if 2 * b == 3:
                    return F32
                if 2 * b > 3:
                    return ((b - 2) * (b - 1) / z * (_0f1(b - 2, z) -
                            _0f1(b - 1, z)))
                if 2 * b < 1:
                    return (_0f1(b + 1, z) + z / (b * (b + 1)) *
                            _0f1(b + 2, z))
                raise ValueError
            # Can evaluate 0F1 in terms of elementary functions when
            # the parameter is a half-integer
            if 2 * b[0] in ZZ and b[0] not in ZZ:
                return _0f1(b[0], z)

        # Confluent hypergeometric function
        if p == 1 and q == 1:
            aa, bb = a[0], b[0]
            if aa * 2 == 1 and bb * 2 == 3:
                t = sqrt(-z)
                return sqrt(pi) / 2 * erf(t) / t
            if a == 1 and b == 2:
                return (exp(z) - 1) / z
            n, m = aa, bb
            if n in ZZ and m in ZZ and m > 0 and n > 0:
                rf = rising_factorial
                if m <= n:
                    return (exp(z) * sum(rf(m - n, k) * (-z) ** k /
                            factorial(k) / rf(m, k) for k in
                            xrange(n - m + 1)))
                else:
                    T = sum(rf(n - m + 1, k) * z ** k /
                            (factorial(k) * rf(2 - m, k)) for k in
                            xrange(m - n))
                    U = sum(rf(1 - n, k) * (-z) ** k /
                            (factorial(k) * rf(2 - m, k)) for k in
                            xrange(n))
                    return (factorial(m - 2) * rf(1 - m, n) *
                            z ** (1 - m) / factorial(n - 1) *
                            (T - exp(z) * U))

        if p == 2 and q == 1:
            R12 = QQ('1/2')
            R32 = QQ('3/2')

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

                if (a, b, c) == (R12, 1, 2):
                    return (2 - 2 * sqrt(1 - z)) / z
                if (a, b, c) == (1, 1, 2):
                    return -log(1 - z) / z
                if (a, b, c) == (1, R32, R12):
                    return (1 + z) / (1 - z) ** 2
                if (a, b, c) == (1, R32, 2):
#.........这里部分代码省略.........

 def f(t):
     return exp(I * ((cos(t) + I * sin(t)) *
                     arc_radius + arc_center)) * radius