Python mpmath.mpc函数代码示例

	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