def Pprime(self,z):
		# A+S 18.9.
		from mpmath import ellipfun, sqrt, cos, sin, mpc, mpf
		Delta = self.Delta
		e1, e2, e3 = self.__roots
		if self.__ng3:
			z = mpc(0,1) * z
		if Delta > 0:
			zs = sqrt(e1 - e3) * z
			m = (e2 - e3) / (e1 - e3)
			retval = -2 * sqrt((e1 - e3)**3) * ellipfun('cn',u=zs,m=m) * ellipfun('dn',u=zs,m=m) / (ellipfun('sn',u=zs,m=m)**3)
		elif Delta < 0:
			H2 = (sqrt((e2 - e3) * (e2 - e1))).real
			assert(H2 > 0)
			m = mpf(1) / mpf(2) - 3 * e2 / (4 * H2)
			zp = 2 * z * sqrt(H2)
			retval = -4 * sqrt(H2**3) * ellipfun('sn',u=zp,m=m) * ellipfun('dn',u=zp,m=m) / ((1 - ellipfun('cn',u=zp,m=m))**2)
		else:
			g2, g3 = self.__invariants
			if g2 == 0 and g3 == 0:
				retval = -2 / (z**3)
			else:
				c = e1 / 2
				A = sqrt(3 * c)
				retval = -6 * c * A * cos(A * z) / (sin(A * z))**3
		if self.__ng3:
			return mpc(0,-1) * retval
		else:
			return retval

 def __cubic_roots(self,a,c,d):
     from mpmath import mpf, mpc, sqrt, cbrt
     assert(all([isinstance(_,mpf) for _ in [a,c,d]]))
     Delta = -4 * a * c*c*c - 27 * a*a * d*d
     self.__Delta = Delta
     # NOTE: this was the original function used for root finding.
     # proots, err = polyroots([a,0,c,d],error=True,maxsteps=5000000)
     # Computation of the cubic roots.
     u_list = [mpf(1),mpc(-1,sqrt(3))/2,mpc(-1,-sqrt(3))/2]
     Delta0 = -3 * a * c
     Delta1 = 27 * a * a * d
     # Here we have two choices for the value from sqrt, positive or negative. Since C
     # is used as a denominator below, we pick the choice with the greatest absolute value.
     # http://en.wikipedia.org/wiki/Cubic_function
     # This should handle the case g2 = 0 gracefully.
     C1 = cbrt((Delta1 + sqrt(Delta1 * Delta1 - 4 * Delta0 * Delta0 * Delta0)) / 2)
     C2 = cbrt((Delta1 - sqrt(Delta1 * Delta1 - 4 * Delta0 * Delta0 * Delta0)) / 2)
     if abs(C1) > abs(C2):
         C = C1
     else:
         C = C2
     proots = [(-1 / (3 * a)) * (u * C + Delta0 / (u * C)) for u in u_list]
     # NOTE: we ignore any residual imaginary part that we know must come from numerical artefacts.
     if Delta < 0:
         # Sort the roots following the convention: complex with negative imaginary, real, complex with positive imaginary.
         # Then assign the roots following the P convention (e2 real root, e1 complex with positive imaginary).
         e3,e2,e1 = sorted(proots,key = lambda c: c.imag)
     else:
         # The convention in this case is to sort in descending order.
         e1,e2,e3 = sorted([_.real for _ in proots],reverse = True)
     return e1,e2,e3

def e_ratio(a,b,e,x):
  # Get S
  bt2 = mp.beta(a,b-1.0)             # Beta function
  bix = mp.betainc(a,b+1.0,0.0,e)    # Incomplete Beta function
  hf = mp.hyp2f1(1.0,a,a+b-1.0,-1.0) # 2F1, Gauss' hypergeometric function
  hfre = mp.re(hf)
  Sval = bix - x*bt2*hfre
  # Get U
  c1 = mp.mpc(1.0 + a)
  c2 = mp.mpc(-b)
  c3 = mp.mpc(1.0)
  c4 = mp.mpc(2.0 + a)
  Uval = mp.appellf1(c1,c2,c3,c4,e,-e)
  Ure = mp.re(Uval)
  # Get P & Q
  Pval = mp.hyp2f1(a+1.0,1.0-b,a+2.0,e) # 2F1, Gauss' hypergeometric function
  Pre = mp.re(Pval)
  Qval = mp.hyp2f1(a+1.0,2.0-b,a+2.0,e) # 2F1, Gauss' hypergeometric function
  Qre = mp.re(Qval)
  # Get T
  Tval = ( (e**(1.0+a)) / (1.0+a) )*( 3.0*Pre + 2.0*Qre - Ure )
  Tval = Tval + 4.0*Sval
  # Get Rval (ratio)
  Rval = 0.25*(1.0-e*e)*( (1.0-e)**(1.0-b) )*( e**(1.0-a) )*Tval
  return Rval

	def __compute_periods(self):
		# A+S 18.9.
		from mpmath import sqrt, ellipk, mpc, pi, mpf
		Delta = self.Delta
		e1, e2, e3 = self.__roots
		if Delta > 0:
			m = (e2 - e3) / (e1 - e3)
			Km = ellipk(m)
			Kpm = ellipk(1 - m)
			om = Km / sqrt(e1 - e3)
			omp = mpc(0,1) * om * Kpm / Km
		elif Delta < 0:
			# NOTE: the expression in the sqrt has to be real and positive, as e1 and e3 are
			# complex conjugate and e2 is real.
			H2 = (sqrt((e2 - e3) * (e2 - e1))).real
			assert(H2 > 0)
			m = mpf(1) / mpf(2) - 3 * e2 / (4 * H2)
			Km = ellipk(m)
			Kpm = ellipk(1 - m)
			om2 = Km / sqrt(H2)
			om2p = mpc(0,1) * Kpm * om2 / Km
			om = (om2 - om2p) / 2
			omp = (om2 + om2p) / 2
		else:
			g2, g3 = self.__invariants
			if g2 == 0 and g3 == 0:
				om = mpf('+inf')
				omp = mpc(0,'+inf')
			else:
				# NOTE: here there is no need for the dichotomy on the sign of g3 because
				# we are already working in a regime in which g3 >= 0 by definition.
				c = e1 / 2
				om = 1 / sqrt(12 * c) * pi()
				omp = mpc(0,'+inf')
		return 2 * om, 2 * omp

	def __compute_roots(self,a,c,d):
		from mpmath import mpf, mpc, sqrt, cbrt
		assert(all([isinstance(_,mpf) for _ in [a,c,d]]))
		Delta = self.__Delta
		# NOTE: this was the original function used for root finding.
		# proots, err = polyroots([a,0,c,d],error=True,maxsteps=5000000)
		# Computation of the cubic roots.
		# TODO special casing.
		u_list = [mpf(1),mpc(-1,sqrt(3))/2,mpc(-1,-sqrt(3))/2]
		Delta0 = -3 * a * c
		Delta1 = 27 * a * a * d
		C1 = cbrt((Delta1 + sqrt(Delta1 * Delta1 - 4 * Delta0 * Delta0 * Delta0)) / 2)
		C2 = cbrt((Delta1 - sqrt(Delta1 * Delta1 - 4 * Delta0 * Delta0 * Delta0)) / 2)
		if abs(C1) > abs(C2):
			C = C1
		else:
			C = C2
		proots = [(-1 / (3 * a)) * (u * C + Delta0 / (u * C)) for u in u_list]
		# NOTE: we ignore any residual imaginary part that we know must come from numerical artefacts.
		if Delta < 0:
			# Sort the roots following the convention: complex with negative imaginary, real, complex with positive imaginary.
			# Then assign the roots following the P convention (e2 real root, e1 complex with positive imaginary).
			e3,e2,e1 = sorted(proots,key = lambda c: c.imag)
		else:
			# The convention in this case is to sort in descending order.
			e1,e2,e3 = sorted([_.real for _ in proots],reverse = True)
		return e1,e2,e3

 def test_mpArray_init_(self):
     n = random.randint(1,5000)
     #x = [random.random() for i in range(n)]
     #self.assertRaises(TypeError, mpArray.mpArray, x)
     
     y = [mpmath.rand() for i in range(n)]
     arr = mpArray.mpArray(y)
     self.assertTrue(numpy.all(arr==y))
     
     z = [mpmath.mpc(mpmath.rand(),mpmath.rand()) for i in range(n)]
     
     z_real = [zz.real for zz in z]
     z_imag = [zz.imag for zz in z]
     arr = mpArray.mpArray(z)
     self.assertTrue(numpy.all(arr==z))
     self.assertTrue(numpy.all(arr.real() == z_real))
     self.assertTrue(numpy.all(arr.imag() == z_imag))
     
     start = time.clock()
     x=[mpmath.mpc('0','0')]*n
     t1 = time.clock() - start
     start = time.clock()
     x=[mpmath.mpc('0','0') for i in range(n)]
     t2 = time.clock() -start
     if self.verbose:
         print("[]*",n,t1,"vs.[for i in range(n)]",n,t2)

def dedekind(tau, floatpre):
    """
    Algorithm 22 (Dedekind eta)
    Input : tau in the upper half-plane, k in N
    Output : eta(tau)
    """
    a = 2 * mpmath.pi / mpmath.mpf(24)
    b = mpmath.exp(mpmath.mpc(0, a))
    p = 1
    m = 0
    while m <= 0.999:
        n = nearest_integer(tau.real)
        if n != 0:
            tau -= n
            p *= b ** n
        m = tau.real * tau.real + tau.imag * tau.imag
        if m <= 0.999:
            ro = mpmath.sqrt(mpmath.power(tau, -1) * 1j)
            if ro.real < 0:
                ro = -ro
            p = p * ro
            tau = (-p.real + p.imag * 1j) / m
    q1 = mpmath.exp(a * tau * 1j)
    q = q1 ** 24
    s = 1
    qs = mpmath.mpc(1, 0)
    qn = 1
    des = mpmath.mpf(10) ** (-floatpre)
    while abs(qs) > des:
        t = -q * qn * qn * qs
        qn = qn * q
        qs = qn * t
        s += t + qs
    return p * q1 * s

def QScomplex(val):
    if QSMODE == MODE_NORM:
        return complex(val)
    else:
        if type(val) is str or type(val) is unicode:
            if 'nan' in val:
                return mpmath.mpc(real='nan',imag='nan')
            real = None
            imag = None
            delim = None
            if '+' in val[1:]:
                delim = '+'
            elif '-' in val[1:]:
                delim = '-'
            if delim is None:
                if 'j' in val:
                    imag = val.replace('j','')
                else:
                    real = val
            else:
                index = val[1:].find(delim) + 1
                real = val[:index]
                imag = val[index:].replace('j','')
            return mpmath.mpc(real=real,imag=imag)
        else:
            return mpmath.mpc(val)

def sph_h2n_exact(n, z):
    """Return the value of h^{(2)}_n computed using the exact formula.

    The expression used is http://dlmf.nist.gov/10.49.E7 .

    """
    zm = mpmathify(z)
    s = sum(mpc(0,-1)**(k-n-1)*_a(k, n)/zm**(k+1) for k in xrange(n+1))
    return exp(mpc(0,-1)*zm)*s

def sph_i2n_exact(n, z):
    """Return the value of i^{(2)}_n computed using the exact formula.

    The expression used is http://dlmf.nist.gov/10.49.E10 .

    """
    zm = mpmathify(z)
    s1 = sum(mpc(-1,0)**k * _a(k, n)/zm**(k+1) for k in xrange(n+1))
    s2 = sum(_a(k, n)/zm**(k+1) for k in xrange(n+1))
    return exp(zm)/2 * s1 + mpc(-1,0)**n*exp(-zm)/2 * s2

def BSLaplace(S,K,T,t,r,sig,N,phi): 
        """Solving the Black Scholes PDE in the Laplace domain"""
        x   = ln(S/K)     
        r = mpf(r);sig = mpf(sig);T = mpf(T);t=mpf(t)
        S = mpf(S);K = mpf(K);x=mpf(x)
        mu  = r - 0.5*(sig**2)
       
        tau = T - t   
        c1 = mpf('0.5017')
        c2 = mpf('0.6407')
        c3 = mpf('0.6122')
        c4 = mpc('0','0.2645')        
        
        ans = 0.0
        h = 2*pi/N
        h = mpf(h)
        for k in range(N/2): # Use symmetry
            theta = -pi + (k+0.5)*h
            z     =  N/tau*(c1*theta/tan(c2*theta) - c3 + c4*theta)
            dz    =  N/tau*(-c1*c2*theta/(sin(c2*theta)**2) + c1/tan(c2*theta)+c4)
            eps1  =  (-mu + sqrt(mu**2 + 2*(sig**2)*(z+r)))/(sig**2)
            eps2  =  (-mu - sqrt(mu**2 + 2*(sig**2)*(z+r)))/(sig**2)
            b1    =  1/(eps1-eps2)*(eps2/(z+r) + (1 - eps2)/z)
            b2    =  1/(eps1-eps2)*(eps1/(z+r) + (1 - eps1)/z)
            ans  +=  exp(z*tau)*bs(x,b1,b2,eps1,eps2,z,r,phi)*dz
            val = (K*(h/(2j*pi)*ans)).real
           
            
        return 2*val

	def Pinv(self,P):
		from mpmath import ellipf, sqrt, asin, acos, mpc, mpf
		Delta = self.Delta
		e1, e2, e3 = self.__roots
		if self.__ng3:
			P = -P
		if Delta > 0:
			m = (e2 - e3) / (e1 - e3)
			retval = (1 / sqrt(e1 - e3)) * ellipf(asin(sqrt((e1 - e3)/(P - e3))),m=m)
		elif Delta < 0:
			H2 = (sqrt((e2 - e3) * (e2 - e1))).real
			assert(H2 > 0)
			m = mpf(1) / mpf(2) - 3 * e2 / (4 * H2)
			retval = 1 / (2 * sqrt(H2)) * ellipf(acos((e2-P+H2)/(e2-P-H2)),m=m)
		else:
			g2, g3 = self.__invariants
			if g2 == 0 and g3 == 0:
				retval = 1 / sqrt(P)
			else:
				c = e1 / 2
				retval = (1 / sqrt(3 * c)) * asin(sqrt((3 * c)/(P + c)))
		if self.__ng3:
			retval /= mpc(0,1)
		alpha, beta, _, _ = self.reduce_to_fpp(retval)
		T1, T2 = self.periods
		return T1 * alpha + T2 * beta

 def test_mpArray_toarray(self):
     n = random.randint(1,5000)
     input = numpy.random.random(n) + 1.j*numpy.random.random()
     y = [mpmath.mpc(x.real, x.imag) for x in input]
     arr = mpArray.mpArray(y)
     out = arr.toarray()
     self.assertTrue(numpy.all(out == input))

 def test_Matrix01(self):
     a = mpmath.mpc(mpmath.pi,"0")
     b = mpmath.mpc("0","1")
     m = mpArray.mpArray([[a,b],[b,a]])
     mat = mpArray.mpMatrix(m)
     evals, evecs = mat.eigen(False, 'qeispack', verbose=True)
     
     self.assertTrue(mpmath.fabs(mpmath.mpc(mpmath.pi,"1")  - evals[0]) < 1e-32)
     self.assertTrue(mpmath.fabs(mpmath.mpc(mpmath.pi,"-1")  - evals[1]) < 1e-32)
     self.assertTrue(numpy.all([mpmath.fabs(mpmath.mpf(1)/mpmath.sqrt("2") - mpmath.fabs(x)) < 1e-32 for x in evecs[0]]))
     self.assertTrue(numpy.all([mpmath.fabs(mpmath.mpf(1)/mpmath.sqrt("2") - mpmath.fabs(x)) < 1e-32 for x in evecs[1]]))
     
     
     self.assertTrue(mpmath.fabs(mpmath.fsum(evecs[0]*evecs[0].conj())) - mpmath.mpf(1) < 1e-32)
     self.assertTrue(numpy.all([vec.abs2() -mpmath.mpf(1) < 1e-32 for vec in evecs]))
     self.assertTrue(mpmath.fabs(evecs[0].inner(evecs[1])) < 1e-32)

 def __new__(cls, input=[], dtype='object'):
     if isinstance(input, int):
         data = [mpmath.mpc('0','0')]*input
         obj = numpy.asarray(data, dtype='object').view(cls)
     else:
         obj = numpy.asarray(input, dtype=dtype).view(cls)
     return obj

	def P(self,z):
		# A+S 18.9.
		from mpmath import sqrt, mpc, sin, ellipfun, mpf
		Delta = self.Delta
		e1, e2, e3 = self.__roots
		if self.__ng3:
			z = mpc(0,1) * z
		if Delta > 0:
			zs = sqrt(e1 - e3) * z
			m = (e2 - e3) / (e1 - e3)
			retval = e3 + (e1 - e3) / ellipfun('sn',u=zs,m=m)**2
		elif Delta < 0:
			H2 = (sqrt((e2 - e3) * (e2 - e1))).real
			assert(H2 > 0)
			m = mpf(1) / mpf(2) - 3 * e2 / (4 * H2)
			zp = 2 * z * sqrt(H2)
			retval = e2 + H2 * (1 + ellipfun('cn',u=zp,m=m)) / (1 - ellipfun('cn',u=zp,m=m))
		else:
			g2, g3 = self.__invariants
			if g2 == 0 and g3 == 0:
				retval = 1 / (z**2)
			else:
				c = e1 / 2
				retval = -c + 3 * c / (sin(sqrt(3 * c) * z))**2
		if self.__ng3:
			return -retval
		else:
			return retval

def zp(x):
    """
    plasma dispersion function                                
    using complementary error function in mpmath library.                       
                                                                                
    """
    return -mp.sqrt(mp.pi) * mp.exp(-x**2) * mp.erfi(x) + mpc(0, 1) * mp.sqrt(mp.pi) * mp.exp(-x**2)

def find_Y_and_M(G,R,ndigs=12,Yset=None,Mset=None):
    r"""
    Compute a good value of M and Y for Maass forms on G

    INPUT:

    - ''G'' -- group
    - ''R'' -- real
    - ''ndigs'' -- integer (number of desired digits of precision)
    - ''Yset'' -- real (default None) if set we return M corr. to this Y
    - ''Mset'' -- integer (default None) if set we return Y corr. to this M

    OUTPUT:

    - [Y,M] -- good values of Y (real) and M (integer)

    EXAMPLES::

    

    TODO:
    Better and more effective bound
    """
    l=G._level
    if(Mset <> None):
        # then we get Y corr. to this M
        Y0=mpmath.sqrt(3.0)/mpmath.mpf(2*l)
        
    if(Yset==None):
        Y0=mpmath.sqrt(3.0)/mpmath.mpf(2*l)
        Y=mpmath.mpf(0.95*Y0)
    else:
        Y=Yset
    #print "Y=",Y,"Yset=",Yset
    IR=mpmath.mpc(0,R)
    eps=mpmath.mpf(10 **-ndigs)
    twopiY=mpmath.pi()*Y*mpmath.mpf(2)

    M0=get_M_for_maass(R,Y,eps) 
    if(M0<10):
        M0=10
    ## Do this in low precision
    dold=mpmath.mp.dps
    #print "Start M=",M0
    #print "dold=",dold
    #mpmath.mp.dps=100
    try:
        for n in range(M0,10000 ):
            X=mpmath.pi()*Y*mpmath.mpf(2 *n)
            #print "X,IR=",X,IR
            test=mpmath.fp.besselk(IR,X)
            if(abs(test)<eps):
                raise StopIteration()
    except StopIteration:
        M=n
    else:
        M=n
        raise Exception,"Error: Did not get small enough error:=M=%s gave err=%s" % (M,test)
    mpmath.mp.dps=dold
    return [Y,M]

def testing_kbes(Rt,Xt):
    [R0,R1,NR]=Rt
    [X0,X1,NX]=Xt
    NRr=mpmath.mpf(NR)
    NXr=mpmath.mpf(NX)
    for j in range(1,NR):
        rj=mpmath.mpf(j)
        R=R0+R1*rj/NRr
        iR=mpmath.mpc(0,R)
        for k in range(1,NX):
            rk=mpmath.mpf(k)
            x=X0+X1*rk/NXr
            print "r,x=",R,x
            if(x>R):
                print "kbes_asymp="
                timeit( "kbes_asymp(R,x)",repeat=1)
            else:
                print "kbes_rec="
                timeit( "kbes_rec(R,x)",repeat=1)
            print "mpmath.besselk="
            timeit("mpmath.besselk(iR,x)",repeat=1)
            

            #print "t1(",R,x,")=",t1
            #print "t2(",R,x,")=",t2
            if(R<15.0):
                if(x<0.3 *R):
                    print "Case 1"
                elif(x<=max(10.0 +1.2*R,2 *R)):
                    print "Case 2"
            elif(R>20  and x>4 *R):
                print "Case 3"
            else:
                print "Case 4"

def getNthPadovanNumber( arg ):
    n = fadd( real( arg ), 4 )

    a = root( fsub( fdiv( 27, 2 ), fdiv( fmul( 3, sqrt( 69 ) ), 2 ) ), 3 )
    b = root( fdiv( fadd( 9, sqrt( 69 ) ), 2 ), 3 )
    c = fadd( 1, fmul( mpc( 0, 1 ), sqrt( 3 ) ) )
    d = fsub( 1, fmul( mpc( 0, 1 ), sqrt( 3 ) ) )
    e = power( 3, fdiv( 2, 3 ) )

    r = fadd( fdiv( a, 3 ), fdiv( b, e ) )
    s = fsub( fmul( fdiv( d, -6 ), a ), fdiv( fmul( c, b ), fmul( 2, e ) ) )
    t = fsub( fmul( fdiv( c, -6 ), a ), fdiv( fmul( d, b ), fmul( 2, e ) ) )

    return nint( re( fsum( [ fdiv( power( r, n ), fadd( fmul( 2, r ), 3 ) ),
                             fdiv( power( s, n ), fadd( fmul( 2, s ), 3 ) ),
                             fdiv( power( t, n ), fadd( fmul( 2, t ), 3 ) ) ] ) ) )