def string2number(s):
    # a start to replace p2sage (used for the paramters in the FE)

    strs = str(s).replace(' ','')
    try:
        if 'e' in strs:
            # check for e(m/n) := exp(2*pi*i*m/n), used by Dirichlet characters, for example
            r = re.match('^\$?e\\\\left\(\\\\frac\{(?P<num>\d+)\}\{(?P<den>\d+)\}\\\\right\)\$?$',strs)
            if not r:
                r = re.match('^e\((?P<num>\d+)/(?P<den>\d+)\)$',strs)
            if r:
                q = Rational(r.groupdict()['num'])/Rational(r.groupdict()['den'])
                return CDF(exp(2*pi*I*q))
        if 'I' in strs:
            return CDF(strs)
        elif '/' in strs:
            return Rational(strs)
        elif strs=='0.5':  # Temporary fix because 0.5 in db for EC
            return Rational('1/2')
        elif '.' in strs:
            return float(strs)
        else:
            return Integer(strs)
    except:
        return s

def number_of_coefficients_needed(Q, kappa_fe, lambda_fe, max_t):
    # TODO: This doesn't work. Trouble when computing t0
    # We completely mimic what lcalc does when it decides whether
    # to print a warning.

    DIGITS = 14    # These are names of lcalc parameters, and we are
    DIGITS2 = 2    # mimicking them.

    logger.debug("Start NOC")
    theta = sum(kappa_fe)
    c = DIGITS2 * log(10.0)
    a = len(kappa_fe)

    c1 = 0.0
    for j in range(a):
        logger.debug("In loop NOC")
        t0 = kappa_fe[j] * max_t + complex(lambda_fe[j]).imag()
        logger.debug("In loop 2 NOC")
        if abs(t0) < 2 * c / (math.pi * a):
            logger.debug("In loop 3_1 NOC")
            c1 += kappa_fe[j] * pi / 2.0
        else:
            c1 += kappa_fe[j] * abs(c / (t0 * a))
            logger.debug("In loop 3_2 NOC")

    return int(round(Q * exp(log(2.3 * DIGITS * theta / c1) * theta) + 10))

def main4():
    print("MAIN 4:")
    zeta = CDF(exp(Integer(2) * pi * I / Integer(5)))
    print(zeta)
    v = [legendre_symbol(n, Integer(5)) * zeta**(Integer(2) * n)
         for n in range(Integer(1), Integer(5))]
    S = sum([point(tuple(z), pointsize=Integer(100)) for z in v])
    G = point(tuple(sum(v)), pointsize=Integer(100), rgbcolor='red')
    (S + G).save(filename="gauss_sum.png")
    print()

        def laplace(cls, self, parameters, variable, x='x', s='t'):
            r"""
            Returns the Laplace transform of self with respect to the variable
            var.

            INPUT:

            -  ``x`` - variable of self

            -  ``s`` - variable of Laplace transform.

            We assume that a piecewise function is 0 outside of its domain and
            that the left-most endpoint of the domain is 0.

            EXAMPLES::

                sage: x, s, w = var('x, s, w')
                sage: f = piecewise([[(0,1),1],[[1,2], 1-x]])
                sage: f.laplace(x, s)
                -e^(-s)/s + (s + 1)*e^(-2*s)/s^2 + 1/s - e^(-s)/s^2
                sage: f.laplace(x, w)
                -e^(-w)/w + (w + 1)*e^(-2*w)/w^2 + 1/w - e^(-w)/w^2

            ::

                sage: y, t = var('y, t')
                sage: f = piecewise([[[1,2], 1-y]])
                sage: f.laplace(y, t)
                (t + 1)*e^(-2*t)/t^2 - e^(-t)/t^2

            ::

                sage: s = var('s')
                sage: t = var('t')
                sage: f1(t) = -t
                sage: f2(t) = 2
                sage: f = piecewise([[[0,1],f1],[(1,infinity),f2]])
                sage: f.laplace(t,s)
                (s + 1)*e^(-s)/s^2 + 2*e^(-s)/s - 1/s^2
            """
            from sage.all import assume, exp, forget
            x = SR.var(x)
            s = SR.var(s)
            assume(s>0)
            result = 0
            for domain, f in parameters:
                for interval in domain:
                    a = interval.lower()
                    b = interval.upper()
                    result += (SR(f)*exp(-s*x)).integral(x,a,b)
            forget(s>0)
            return result

    def from_conjugacy_class_index_to_polynomial(self, index):
        """ A function converting from a conjugacy class index (starting at 1) to the local Euler polynomial.
            Saves a sequence of processed polynomials, obtained from the local factors table, so it can reuse computations from prime to prime
            This sequence is indexed by conjugacy class indices (starting at 1, filled with dummy first) and gives the corresponding polynomials in the form
            [coeff_deg_0, coeff_deg_1, ...], where coeff_deg_i is in ComplexField(). This could be changed later, or made parametrizable
        """
        try:
            return self._from_conjugacy_class_index_to_polynomial_fn(index)
        except AttributeError:
            local_factors = self.local_factors_table()
            field = ComplexField()
            root_of_unity = exp((field.gen()) * 2 * field.pi() / int(self.character_field()))
            local_factor_processed_pols = [0]   # dummy to account for the shift in indices
            for pol in local_factors:
                local_factor_processed_pols.append(
                    process_polynomial_over_algebraic_integer(pol, field, root_of_unity))

            def tmp(conjugacy_class_index_start_1):
                return local_factor_processed_pols[conjugacy_class_index_start_1]
            self._from_conjugacy_class_index_to_polynomial_fn = tmp
            return self._from_conjugacy_class_index_to_polynomial_fn(index)

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

def draw_fundamental_domain(N,group='Gamma0',model="H",axes=None,filename=None,**kwds):
    r""" Draw fundamental domain
    INPUT:
        - ''model'' -- (default ''H'')
        - ''H'' -- Upper halfplane
        - ''D'' -- Disk model
        - ''filename''-- filename to print to
        - ''**kwds''-- additional arguments to matplotlib 
        - ''axes''  -- set geometry of output
        =[x0,x1,y0,y1] -- restrict figure to [x0,x1]x[y0,y1]

    EXAMPLES::

        sage: G=MySubgroup(Gamma0(3))
        sage: G.draw_fundamental_domain()

    """
    G=eval(group+'('+str(N)+')')
    #print G
    name ="$"+latex(G)+"$"
    ## need a "nice" set of coset representatives to draw a connected fundamental domain. Only implemented for Gamma_0(N)
    coset_reps = nice_coset_reps(G)
    if(model=="D"):
        g=draw_funddom_d(coset_reps,format,I)
    else:
        g=draw_funddom(coset_reps,format)
    if(axes<>None):
        [x0,x1,y0,y1]=axes
    elif(model=="D"):
        x0=-1 ; x1=1 ; y0=-1.1 ; y1=1 
    else:
        # find the width of the fundamental domain
        w=0  #self.cusp_width(Cusp(Infinity))
        wmin=0 ; wmax=1
        max_x = RR(0.55)
        rho = CC( exp(2*pi*I/3))
        for V in coset_reps:
            ## we also compare the real parts of where rho and infinity are mapped
            r1 = (V.acton(rho)).real()
            if(V[1,0]<>0):
                inf1 = RR(V[0,0] / V[1,0])
            else:
                inf1 = 0
            if(V[1 ,0 ]==0  and V[0 ,0 ]==1 ):
                if(V[0 ,1 ]>wmax):
                    wmax=V[0 ,1 ]
                if(V[0 ,1 ]<wmin):
                    wmin=V[0 ,1 ]
            if( max(r1,inf1) > max_x):
                max_x = max(r1,inf1)
            #print "wmin,wmax=",wmin,wmax
            #x0=-1; x1=1; y0=-0.2; y1=1.5
        x0=RR(-max_x) ; x1=RR(max_x) ; y0=RR(-0.15) ; y1=RR(1.5) 
    ## Draw the axes manually (since  can't figure out how to do it automatically)
    ax = line([[x0,0.0],[x1,0.0]],color='black')
        #ax = ax + line([[0.0,y0],[0.0,y1]],color='black')
        ## ticks
    ax = ax + line([[-0.5,-0.01],[-0.5,0.01]],color='black')
    ax = ax + line([[0.5,-0.01],[0.5,0.01]],color='black')
    g = g + ax
    if model=="H":
        t = text(name, (0, -0.1), fontsize=16, color='black')
    else:
        t = text(name, (0, -1.1), fontsize=16, color='black')
        g = g + t
        g.set_aspect_ratio(1)
        g.set_axes_range(x0,x1,y0,y1)
        g.axes(False)
        if(filename<>None):
            fig = g.matplotlib()
            fig.set_canvas(FigureCanvasAgg(fig))
            axes = fig.get_axes()[0]
            axes.minorticks_off()
            axes.set_yticks([])
            fig.savefig(filename,**kwds)
        else:
            return g

 def value(self, z, embedding=0):
     if self.prec == 0:
         return 0
     else:
         q = exp(2*CC.pi()*CC(0,1)*z)
         return sum(self.coefficient_embedding(n,embedding)*q**n for n in range(self.prec))

    def dimension(self,k,ignore=False, debug = 0):
        if k < 2 and not ignore:
            raise NotImplementedError("k has to >= 2")
        s = self._signature
        if not (2*k in ZZ):
            raise ValueError("k has to be integral or half-integral")
        if (2*k+s)%4 != 0 and not ignore:
            raise NotImplementedError("2k has to be congruent to -signature mod 4")
        if self._v2.has_key(0):
            v2 = self._v2[0]
        else:
            v2 = 1

        if self._g != None:
            if not self._aniso_formula:
                vals = self._g.values()
            #else:
                #print "using aniso_formula"
            M = self._g
        else:
            vals = self._M.values()
            M = self._M
            
        prec = ceil(max(log(M.order(),2),52)+1)+17
        #print prec
        RR = RealField(prec)
        CC = ComplexField(prec)
        d = self._d
        m = self._m
        if debug > 0: print d,m
            
        if self._alpha3 == None:
            if self._aniso_formula:
                self._alpha4 = 1
                self._alpha3 = -sum([BB(a)*mm for a,mm in self._v2.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 += Integer(d) - Integer(1) - self._g.beta_formula()
                #print self._alpha3, self._g.a5prime_formula()
                self._alpha3 = self._alpha3/RR(2)
            else:
                self._alpha3 = sum([(1-a)*mm for a,mm in self._v2.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 += sum([(1-a)*mm for a,mm in vals.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 = self._alpha3 / Integer(2)
                self._alpha4 = 1/Integer(2)*(vals[0]+v2) # the codimension of SkL in MkL
        alpha3 = self._alpha3
        alpha4 = self._alpha4
        if debug > 0: print alpha3, alpha4
        g1=M.char_invariant(1)
        g1=CC(g1[0]*g1[1])
        #print g1
        g2=M.char_invariant(2)
        g2=RR(real(g2[0]*g2[1]))
        if debug > 0: print g2, g2.parent()
        g3=M.char_invariant(-3)
        g3=CC(g3[0]*g3[1])
        if debug > 0: print RR(d) / RR(4), sqrt(RR(m)) / RR(4), CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8)))
        alpha1 = RR(d) / RR(4) - (sqrt(RR(m)) / RR(4)  * CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8))) * g2)
        if debug > 0: print alpha1
        alpha2 = RR(d) / RR(3) + sqrt(RR(m)) / (3 * sqrt(RR(3))) * real(exp(CC(2 * CC.pi() * CC(0,1) * (4 * k + 3 * s - 10) / 24)) * (g1+g3))
        if debug > 0: print alpha1, alpha2, g1, g2, g3, d, k, s
        dim = real(d + (d * k / Integer(12)) - alpha1 - alpha2 - alpha3)
        if debug > 0:
            print "dimension:", dim
        if abs(dim-round(dim)) > 1e-6:
            raise RuntimeError("Error ({0}) too large in dimension formula for {1} and k={2}".format(abs(dim-round(dim)), self._M if self._M is not None else self._g, k))
        dimr = dim
        dim = Integer(round(dim))
        if k >=2 and dim < 0:
            raise RuntimeError("Negative dimension (= {0}, {1})!".format(dim, dimr))
        return dim

 def compute_cm_values_numeric(self,digits=12,insert_in_db=True):
     r"""
     Compute CM-values numerically.
     """
     if isinstance(self._cm_values,dict) and self._cm_values  <> {}:
         return self._cm_values
      # the points we want are i and rho. More can be added later...
     bits = ceil(int(digits) * int(4))
     CF = ComplexField(bits)
     RF = ComplexField(bits)
     eps = RF(10 ** - (digits + 1))
     if(self._verbose > 1):
         wmf_logger.debug("eps={0}".format(eps))
     K = self.base_ring()
     # recall that
     degree = self.degree()
     cm_vals = dict()
     rho = CyclotomicField(3).gen()
     zi = CyclotomicField(4).gen()
     points = [rho, zi]
     maxprec = 1000  # max size of q-expansion
     minprec = 10  # max size of q-expansion
     for tau in points:
         q = CF(exp(2 * pi * I * tau))
         fexp = dict()
         cm_vals[tau] = dict()
         if(tau == I and self.level() == -1):
             # cv=    #"Exact(soon...)" #_cohen_exact_formula(k)
             for h in range(degree):
                 cm_vals[tau][h] = cv
             continue
         if K.absolute_degree()==1:
             v1 = CF(0)
             v2 = CF(1)
             try:
                 for prec in range(minprec, maxprec, 10):
                     if(self._verbose > 1):
                         wmf_logger.debug("prec={0}".format(prec))
                     v2 = self.as_factor().q_eigenform(prec).truncate(prec)(q)
                     err = abs(v2 - v1)
                     if(self._verbose > 1):
                         wmf_logger.debug("err={0}".format(err))
                     if(err < eps):
                         raise StopIteration()
                     v1 = v2
                 cm_vals[tau][0] = None
             except StopIteration:
                 cm_vals[tau][0] = v2
         else:
             v1 = dict()
             v2 = dict()
             err = dict()
             for h in range(degree):
                 v1[h] = 1
                 v2[h] = 0
             try:
                 for prec in range(minprec, maxprec, 10):
                     if(self._verbose > 1):
                         wmf_logger.debug("prec={0}".format(prec))
                     c = self.coefficients(range(prec),insert_in_db=insert_in_db)
                     for h in range(degree):
                         fexp[h] = list()
                         v2[h] = 0
                         for n in range(prec):
                             cn = c[n]
                             if hasattr(cn, 'complex_embeddings'):
                                 cc = cn.complex_embeddings(CF.prec())[h]
                             else:
                                 cc = CF(cn)
                             v2[h] = v2[h] + cc * q ** n
                         err[h] = abs(v2[h] - v1[h])
                         if(self._verbose > 1):
                             wmf_logger.debug("v1[{0}]={1}".format(h,v1[h]))
                             wmf_logger.debug("v2[{0}]={1}".format(h,v2[h]))
                             wmf_logger.debug("err[{0}]={2}".format(h,err[h]))
                         if(max(err.values()) < eps):
                             raise StopIteration()
                         v1[h] = v2[h]
             except StopIteration:
                 pass
             for h in range(degree):
                 if(err[h] < eps):
                     cm_vals[tau][h] = v2[h]
                 else:
                     cm_vals[tau][h] = None
     self._cm_values = cm_vals

def gaussum(n, N, prec=53):
    CC = ComplexField(prec)
    return sum(CC(exp(2 * CC.pi() * CC(0, 1) * n * m ** 2 / N)) for m in range(N))

    def dimension(self,k,ignore=False, debug = 0):
        if k < 2 and not ignore:
            raise NotImplementedError("k has to >= 2")
        s = self._signature
        if not (2*k in ZZ):
            raise ValueError("k has to be integral or half-integral")
        if (2*k+s)%2 != 0:
            return 0
        m = self._m
        n2 = self._n2
        if self._v2.has_key(0):
            v2 = self._v2[0]
        else:
            v2 = 1

        if self._g != None:
            if not self._aniso_formula:
                vals = self._g.values()
            #else:
                #print "using aniso_formula"
            M = self._g
        else:
            vals = self._M.values()
            M = self._M

        if (2*k+s)%4 == 0:
            d = Integer(1)/Integer(2)*(m+n2) # |dimension of the Weil representation on even functions|
            self._d = d
            self._alpha4 = 1/Integer(2)*(vals[0]+v2) # the codimension of SkL in MkL
        else:
            d = Integer(1)/Integer(2)*(m-n2) # |dimension of the Weil representation on odd functions|
            self._d = d
            self._alpha4 = 1/Integer(2)*(vals[0]-v2) # the codimension of SkL in MkL
            
        prec = ceil(max(log(M.order(),2),52)+1)+17
        #print prec
        RR = RealField(prec)
        CC = ComplexField(prec)
        
        if debug > 0: print "d, m = {0}, {1}".format(d,m)
        eps = exp( 2 * CC.pi() * CC(0,1) * (s + 2*k) / Integer(4) )
        eps = round(real(eps))        
        if self._alpha3 is None or self._last_eps != eps:
            self._last_eps = eps
            if self._aniso_formula:
                self._alpha4 = 1
                self._alpha3 = -sum([BB(a)*mm for a,mm in self._v2.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 += Integer(d) - Integer(1) - self._g.beta_formula()
                #print self._alpha3, self._g.a5prime_formula()
                self._alpha3 = self._alpha3/RR(2)
            else:
                self._alpha3 = eps*sum([(1-a)*mm for a,mm in self._v2.iteritems() if a != 0])
                if debug>0: print "alpha3t = ", self._alpha3
                self._alpha3 += sum([(1-a)*mm for a,mm in vals.iteritems() if a != 0])
                #print self._alpha3
                self._alpha3 = self._alpha3 / Integer(2)
        alpha3 = self._alpha3
        alpha4 = self._alpha4
        if debug > 0: print alpha3, alpha4
        g1=M.char_invariant(1)
        g1=CC(g1[0]*g1[1])
        #print g1
        g2=M.char_invariant(2)
        g2=CC(g2[0]*g2[1])
        if debug > 0: print g2, g2.parent()
        g3=M.char_invariant(-3)
        g3=CC(g3[0]*g3[1])
        if debug > 0: print "eps = {0}".format(eps)
        if debug > 0: print "d/4 = {0}, m/4 = {1}, e^(2pi i (2k+s)/8) = {2}".format(RR(d) / RR(4), sqrt(RR(m)) / RR(4), CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8))))
        if eps == 1:
            g2_2 = real(g2)
        else:
            g2_2 = imag(g2)*CC(0,1)
        alpha1 = RR(d) / RR(4) - sqrt(RR(m)) / RR(4)  * CC(exp(2 * CC.pi() * CC(0,1) * (2 * k + s) / Integer(8)) * g2_2)
        if debug > 0: print alpha1
        alpha2 = RR(d) / RR(3) + sqrt(RR(m)) / (3 * sqrt(RR(3))) * real(exp(CC(2 * CC.pi() * CC(0,1) * (4 * k + 3 * s - 10) / 24)) * (g1 + eps*g3))
        if debug > 0: print "alpha1 = {0}, alpha2 = {1}, alpha3 = {2}, g1 = {3}, g2 = {4}, g3 = {5}, d = {6}, k = {7}, s = {8}".format(alpha1, alpha2, alpha3, g1, g2, g3, d, k, s)
        dim = real(d + (d * k / Integer(12)) - alpha1 - alpha2 - alpha3)
        if debug > 0:
            print "dimension:", dim
        if abs(dim-round(dim)) > 1e-6:
            raise RuntimeError("Error ({0}) too large in dimension formula for {1} and k={2}".format(abs(dim-round(dim)), self._M if self._M is not None else self._g, k))
        dimr = dim
        dim = Integer(round(dim))
        if k >=2 and dim < 0:
            raise RuntimeError("Negative dimension (= {0}, {1})!".format(dim, dimr))
        return dim

def draw_fundamental_domain(N, group="Gamma0", model="H", axes=None, filename=None, **kwds):
    r""" Draw fundamental domain
    INPUT:
    - ''model'' -- (default ''H'')
    = ''H'' -- Upper halfplane
    = ''D'' -- Disk model
    - ''filename''-- filename to print to
    - ''**kwds''-- additional arguments to matplotlib
    - ''axes''  -- set geometry of output
    =[x0,x1,y0,y1] -- restrict figure to [x0,x1]x[y0,y1]

    EXAMPLES::

    sage: G=MySubgroup(Gamma0(3))
    sage: G.draw_fundamental_domain()

    """
    if group.strip() == "Gamma0":
        G = Gamma0(N)
    elif group.strip() == "Gamma1":
        G = Gamma1(N)
    elif group.strip() == "Gamma":
        G = Gamma(N)
    else:
        raise ValueError('group must be one of: "Gamma0", "Gamma1", "Gamma"')
    s = "$" + latex(G) + "$"
    s = s.replace("mbox", "mathrm")
    s = s.replace("Bold", "mathbb")
    name = s
    # name ="$\mbox{SL}_{2}(\mathbb{Z})$"
    # need a "nice" set of coset representatives to draw a connected
    # fundamental domain. Only implemented for Gamma_0(N)
    coset_reps = nice_coset_reps(G)
    # if(group=='Gamma0'):
    # else:
    # coset_reps = list(G.coset_reps())
    from matplotlib.backends.backend_agg import FigureCanvasAgg

    if model == "D":
        g = _draw_funddom_d(coset_reps, format, I)
    else:
        g = _draw_funddom(coset_reps, format)
    if axes is not None:
        [x0, x1, y0, y1] = axes
    elif model == "D":
        x0 = -1
        x1 = 1
        y0 = -1.1
        y1 = 1
    else:
        # find the width of the fundamental domain
        # w = 0  # self.cusp_width(Cusp(Infinity)) FIXME: w is never used
        wmin = 0
        wmax = 1
        max_x = RR(0.55)
        rho = CC(exp(2 * pi * I / 3))
        for V in coset_reps:
            ## we also compare the real parts of where rho and infinity are mapped
            r1 = (V.acton(rho)).real()
            if V[1, 0] != 0:
                inf1 = RR(V[0, 0] / V[1, 0])
            else:
                inf1 = 0
            if V[1, 0] == 0 and V[0, 0] == 1:
                if V[0, 1] > wmax:
                    wmax = V[0, 1]
                if V[0, 1] < wmin:
                    wmin = V[0, 1]
            if max(r1, inf1) > max_x:
                max_x = max(r1, inf1)
            logging.debug("wmin,wmax=%s,%s" % (wmin, wmax))
        # x0=-1; x1=1; y0=-0.2; y1=1.5
        x0 = RR(-max_x)
        x1 = RR(max_x)
        y0 = RR(-0.15)
        y1 = RR(1.5)
        ## Draw the axes manually (since  can't figure out how to do it automatically)
    ax = line([[x0, 0.0], [x1, 0.0]], color="black")
    # ax = ax + line([[0.0,y0],[0.0,y1]],color='black')
    ## ticks
    ax = ax + line([[-0.5, -0.01], [-0.5, 0.01]], color="black")
    ax = ax + line([[0.5, -0.01], [0.5, 0.01]], color="black")
    g = g + ax
    if model == "H":
        t = text(name, (0, -0.1), fontsize=16, color="black")
        t = t + text("$ -\\frac{1}{2} $", (-0.5, -0.1), fontsize=12, color="black")
        t = t + text("$ \\frac{1}{2} $", (0.5, -0.1), fontsize=12, color="black")
    else:
        t = text(name, (0, -1.1), fontsize=16, color="black")
    g = g + t
    g.set_aspect_ratio(1)
    g.set_axes_range(x0, x1, y0, y1)
    g.axes(False)

    if filename is not None:
        fig = g.matplotlib()
        fig.set_canvas(FigureCanvasAgg(fig))
        axes = fig.get_axes()[0]
        axes.minorticks_off()
        axes.set_yticks([])
#.........这里部分代码省略.........

    def _hyperbolic_arc_d(self, z0, z3, first=False):
        """
        Function to construct Bezier path as an approximation to
        the hyperbolic arc between the complex numbers z0 and z3 in the
        hyperbolic plane.
        """

        w0 = self._cayley_transform(z0)
        w3 = self._cayley_transform(z3)
        if self._verbose>0:
            print "in plane z0,z3=",z0,z3
            print "in disc: ",w0,w3
        npts = 10
        if z0 == infinity or z0==CC(infinity):
            zm = [z3 + I*(j+0.5) for j in range(npts)]
            wm = [self._cayley_transform(x) for x in zm]
            pts = [w3]
            pts.extend(wm)
            pts.append(w0)
            opt = self._options
            opt['fill']=False
            self._graphics.add_primitive(BezierPath([[(x.real(),x.imag()) for x in pts ]],opt))
            return 
        if z3 == infinity  or z3==CC(infinity):
            zm = [z0 + I*(j+0.5) for j in range(npts)]
            wm = [self._cayley_transform(x) for x in zm]
            pts = [w0]
            pts.extend(wm)
            pts.append(w3)
            opt = self._options
            opt['fill']=False
            self._graphics.add_primitive(BezierPath([[(x.real(),x.imag()) for x in pts ]],opt))
            #self._graphics.add_primitive(Line([w1.real(),w2.real()],[w1.imag(),w2.imag()],self._options))
            #self.path.append([(0,w0.imag()),CC(0,y), (0,w3.imag())])
            return
        x0=z0.real(); y0 = z0.imag()
        x3=z3.real(); y3 = z3.imag()
        if y0 == 0 and y3 == 0:
            p = (z0.real()+z3.real())/2
            r = abs(z0-p)
            zm = CC(p, r)
            self._hyperbolic_arc_d(z0, zm, first)
            self._hyperbolic_arc_d(zm, z3)
            return
        else:
            p = ((x0+x3)*(x3-x0)+(y0+y3)*(y3-y0))/(2*(x3-x0)) 
            r = sqrt((p - x0)**2 + y0**2)  # radius of the circle in H
            zm = ((z0+z3)/2-p)/abs((z0+z3)/2-p)*r+p  # midpoint (at least approximately) of geodesic between z0 and z3
            t0 = CC(z0 - p).argument()
            t3 = CC(z3 - p).argument()
            if x0 <= x3:
                zm = [p + r*exp(I*(t0+t*(t3-t0)/npts)) for t in range(npts+1)]
                #print "zm=",zm
            else:
                zm = [p + r*exp(I*(t0+t*(t3-t0)/npts)) for t in range(npts+1)]            
                #print "zm=",zm
            wm = [self._cayley_transform(x) for x in zm]
            opt = self._options
            opt['fill']=False
            w0 = self._cayley_transform(z0)
            w3 = self._cayley_transform(z3)
            pts = [w0]
            pts.extend(wm)
            pts.append(w3)
            self._graphics.add_primitive(BezierPath([[(x.real(),x.imag()) for x in pts ]],opt))
            return 
            #print "z0_test=",(p+r*exp(t0*I))
                #print "z3_test=",(p+r*exp(t3*I))
#t = (8*zm-4*(z0+z3)).imag()/3/(z3-z0).real()
            # I have no idea what these points should represent....
            #z1 = z0 + t*CC(z0.imag(), (p-z0.real()))
            #z2 = z3 - t*CC(z3.imag(), (p-z3.real()))
            wm = [self._cayley_transform(x) for x in zm]
            pp = self._cayley_transform(CC(p,0))
            w1 = self._cayley_transform(z0)
            w2 = self._cayley_transform(z3)
            c = self._cayley_transform(CC(p,0)) # center of circle on the unit disk.
            if self._verbose>0:
                print "p,r=",p,r     
                print "zm=",zm
                #print "t=",t
                #print "tt=",(8*zm-4*(z0+z3)).imag()/3/(z3-z0).real()

                print "C(c)=",pp
                print "C(zm)=",wm
                print "C(z0)=",w1
                print "C(z3)=",w2
                #print "z2=",z2
            
            r = abs(w1-c) # radius
            t0 = CC(w0 - pp).argument()
            t3 = CC(w3 - pp).argument()
            t = abs(t0-t3)
        if self._verbose>0:
            print "adding a:rc ",zm.real(),zm.imag(),r,r,t,t0,t3
        self._graphics.add_primitive(Line([w1.real(),w2.real(),wm.real()],[w1.imag(),w2.imag(),wm.imag()],{'thickness':2,'alpha':1,
                                                                                       'rgbcolor':'blue',
                                                                                                         'legend_label':""}))
        self._graphics.add_primitive(Point([w1.real(),w2.real(),wm.real()],[w1.imag(),w2.imag(),wm.imag()],{'size':10,'alpha':1,
                                                                                                            'faceted':True,
#.........这里部分代码省略.........