def getNthOctagonalTriangularNumber( n ):
    sign = power( -1, real( n ) )

    return nint( floor( fdiv( fmul( fsub( 7, fprod( [ 2, sqrt( 6 ), sign ] ) ),
                                    power( fadd( sqrt( 3 ), sqrt( 2 ) ),
                                           fsub( fmul( 4, real_int( n ) ), 2 ) ) ),
                              96 ) ) )

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 dfdb (y,x,b,c):

    global FT,XO,XT,XF

    ft=FT

    v=x[0]
    i=y[0]

    iss=IS=b[0]
    n=N=b[1]
    ikf=IKF=b[2]
    isr=ISR=b[3]
    nr=NR=b[4]
    vj=VJ=b[5]
    m=M=b[6]
    rs=RS=b[7]

    #ещё раз: ВСЕ константы должны быть в перспективе либо глобальными переменными, как FT, либо составляющими вектора c, но он кривой, потому лучше уж глобальными.

    return mpm.matrix ([[iss*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*(-mpm.exp((-i*rs + v)/(ft*n)) + XO)*(mpm.exp((-i*rs + v)/(ft*n)) - XO)/(XT*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO))) + mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*(mpm.exp((-i*rs + v)/(ft*n)) - XO),
                         iss**XT*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*(-i*rs + v)*(mpm.exp((-i*rs + v)/(ft*n)) - XO)*mpm.exp((-i*rs + v)/(ft*n))/(XT*ft*n**XT*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO))) - iss*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*(-i*rs + v)*mpm.exp((-i*rs + v)/(ft*n))/(ft*n**XT),
                        iss*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO))*(-ikf/(XT*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO))**XT) + XO/(XT*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO))))*(mpm.exp((-i*rs + v)/(ft*n)) - XO)/ikf,
                         ((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*(mpm.exp((-i*rs + v)/(ft*nr)) - XO),
                         -isr*(-i*rs + v)*((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*mpm.exp((-i*rs + v)/(ft*nr))/(ft*nr**XT),
                        isr*m*(XO - (-i*rs + v)/vj)*(-i*rs + v)*((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*(mpm.exp((-i*rs + v)/(ft*nr)) - XO)/(vj**XT*((XO - (-i*rs + v)/vj)**XT + XF)),
                         isr*((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*(mpm.exp((-i*rs + v)/(ft*nr)) - XO)*mpm.log((XO - (-i*rs + v)/vj)**XT + XF)/XT,
                         i*iss**XT*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*(mpm.exp((-i*rs + v)/(ft*n)) - XO)*mpm.exp((-i*rs + v)/(ft*n))/(XT*ft*n*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO))) - i*iss*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*mpm.exp((-i*rs + v)/(ft*n))/(ft*n)+ \
                         i*isr*m*(XO - (-i*rs + v)/vj)*((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*(mpm.exp((-i*rs + v)/(ft*nr)) - XO)/(vj*((XO - (-i*rs + v)/vj)**XT + XF)) - i*isr*((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*mpm.exp((-i*rs + v)/(ft*nr))/(ft*nr)
                        ]])

def getNthNonagonalOctagonalNumber( n ):
    sqrt6 = sqrt( 6 )
    sqrt7 = sqrt( 7 )

    return nint( floor( fdiv( fmul( fsub( fmul( 11, sqrt7 ), fmul( 9, sqrt6 ) ),
                                    power( fadd( sqrt6, sqrt7 ), fsub( fmul( 8, real_int( n ) ), 5 ) ) ),
                              672 ) ) )

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 gamma_shot(self, wrel, l, n, t, tc):
        """
        Calculate 
        1, transfer gain
        2, electron shot noise
        
        """
        if wrel > 1.0 and wrel < 1.2:
            mp.mp.dps = 80
        else:
            mp.mp.dps = 40
        wc = wrel * mp.sqrt(1+n)
        za = self.za_l(wc, l, n, t, tc)
        mp.mp.dps = 15
        zr = self.zr_mp(wc, l, tc)

        # below calculating shot noise:
        ldc = self.ant_len/l
        nc  =permittivity * boltzmann * tc/ ldc**2 / echarge**2
        vtc = np.sqrt(2 * boltzmann * tc/ emass)
        ne = nc * vtc * (1 + n * mp.sqrt(t)) * 2 * np.pi * self.ant_rad * self.ant_len / np.sqrt(4 * np.pi)
        ###################
        ## a: coefficient in shot noise. see Issautier et al. 1999
        ###################
        a = 1 + echarge * 3.6 / boltzmann/tc
        shot_noise = 2 * a * echarge**2 * mp.fabs(za)**2 * ne        
        return [mp.fabs((zr+za)/zr)**2, shot_noise]

def resolver_equacao(a, b, c):
    u'''
    Resolve uma equação do segundo grau.
    '''
    a, b, c = mpf(a), mpf(b), mpf(c)
    d = b ** 2 - 4 * a * c
    return ( (-b+sqrt(d))/(2*a), (-b-sqrt(d))/(2*a) )

 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 fermihalf(x, sgn):
    """ Series approximation to the F_{1/2}(x) or F_{-1/2}(x) 
        Fermi-Dirac integral """

    f = lambda k: mp.sqrt(x ** 2 + np.pi ** 2 * (2 * k - 1) ** 2)

    # if x < -100:
    #    return 0.0
    if x < -9 or True:
        if sgn > 0:
            return mp.exp(x)
        else:
            return mp.exp(x)

    if sgn > 0:  # F_{1/2}(x)
        a = np.array((1.0 / 770751818298, -1.0 / 3574503105, -13.0 / 184757992,
                      85.0 / 3603084, 3923.0 / 220484, 74141.0 / 8289, -5990294.0 / 7995))
        g = lambda k: mp.sqrt(f(k) - x)

    else:  # F_{-1/2}(x)
        a = np.array((-1.0 / 128458636383, -1.0 / 714900621, -1.0 / 3553038,
                      27.0 / 381503, 3923.0 / 110242, 8220.0 / 919))
        g = lambda k: -0.5 * mp.sqrt(f(k) - x) / f(k)

    F = np.polyval(a, x) + 2 * np.sqrt(2 * np.pi) * sum(map(g, range(1, 21)))
    return F

    def gamma_shot(self, wrel, l, n, t, k, tc):
        """
        Calculate 
        1, transfer gain
        2, electron shot noise
        
        """
        if wrel > 1.0 and wrel < 1.1:
            mp.mp.dps = 40
        else:
            mp.mp.dps = 20
        wc = wrel * mp.sqrt(1+n)
        za_val = self.za(wrel, l, n, t, k, tc)
        mp.mp.dps = 15
        zr_val = self.zr(wc, l, tc)

        # below calculating shot noise:
        ldc = self.ant_len/l
        nc  =permittivity * boltzmann * tc/ ldc**2 / echarge**2
        vtc = np.sqrt(2 * boltzmann * tc/ emass)
        ne = nc * vtc * (1 + n * mp.sqrt(t)) * 2 * np.pi * self.ant_rad * self.ant_len / np.sqrt(4 * np.pi)
        ###################################
        ## a: coefficient in shot noise. ##
        ###################################
        scpot= 4
        a = 1 + echarge * scpot / boltzmann/tc
        shot_noise = 2 * a * echarge**2 * mp.fabs(za_val)**2 * ne        
        return [mp.fabs((zr_val+za_val)/zr_val)**2, shot_noise]

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 __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 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 mobility(self, z=1000, E=0, T=300, pn=None):
     if pn is None:
         Eg = self.band_gap(T, symbolic=False, electron_volts=False)
         # print Eg, self.__to_numeric(-Eg/(k*T)), mp.exp(self.__to_numeric(-Eg/(k*T)))
         pn = self.Nc(T, symbolic=False) * self.Nv(T, symbolic=False) * mp.exp(
             self.__to_numeric(-Eg / (k * T))) * 1e-12
         # print pn
     N = 0
     for dopant in self.dopants:
         N += dopant.concentration(z)
     N *= 1e-6
     # print N
     mobility = {'mobility_e': {'mu_L': 0, 'mu_I': 0, 'mu_ccs': 0, 'mu_tot': 0},
                 'mobility_h': {'mu_L': 0, 'mu_I': 0, 'mu_ccs': 0, 'mu_tot': 0}}
     for key in mobility.keys():
         mu_L = self.reference[key]['mu_L0'] * (T / 300.0) ** (-self.reference[key]['alpha'])
         mu_I = (self.reference[key]['A'] * (T ** (3 / 2)) / N) / (
         mp.log(1 + self.reference[key]['B'] * (T ** 2) / N) - self.reference[key]['B'] * (T ** 2) / (
         self.reference[key]['B'] * (T ** 2) + N))
         try:
             mu_ccs = (2e17 * (T ** (3 / 2)) / mp.sqrt(pn)) / (mp.log(1 + 8.28e8 * (T ** 2) * (pn ** (-1 / 3))))
             X = mp.sqrt(6 * mu_L * (mu_I + mu_ccs) / (mu_I * mu_ccs))
         except:
             mu_ccs = np.nan
             X = 0
         # print X
         mu_tot = mu_L * (1.025 / (1 + ((X / 1.68) ** (1.43))) - 0.025)
         Field_coeff = (1 + (mu_tot * E * 1e-2 / self.reference[key]['v_s']) ** self.reference[key]['beta']) ** (
         -1 / self.reference[key]['beta'])
         mobility[key]['mu_L'] = mu_L * 1e-4
         mobility[key]['mu_I'] = mu_I * 1e-4
         mobility[key]['mu_ccs'] = mu_ccs * 1e-4
         mobility[key]['mu_tot'] = mu_tot * 1e-4 * Field_coeff
     return mobility

	def __compute_t_r(self,n_lobes,lobe_idx,H_in,Hr,d_eval,p4roots,lead_cf):
		from pyranha import math
		from mpmath import asin, sqrt, ellipf, mpf
		assert(n_lobes == 1 or n_lobes == 2)
		C = -lead_cf
		assert(C > 0)
		# First determine if we are integrating in the upper or lower plane.
		# NOTE: we don't care about eps, we are only interested in the sign.
		if (self.__F1 * math.sin(pt('h_\\ast'))).trim().evaluate(d_eval) > 0:
			sign = 1
		else:
			sign = -1
		if n_lobes == 2:
			assert(lobe_idx == 0 or lobe_idx == 1)
			r0,r1,r2,r3 = p4roots
			# k is the same in both cases.
			k = sqrt(((r3-r2)*(r1-r0))/((r3-r1)*(r2-r0)))
			if lobe_idx == 0:
				assert(Hr == r0)
				phi = asin(sqrt(((r3-r1)*(H_in-r0))/((r1-r0)*(r3-H_in))))
			else:
				assert(Hr == r2)
				phi = asin(sqrt(((r3-r1)*(H_in-r2))/((H_in-r1)*(r3-r2))))
			return -sign * mpf(2) / sqrt(C * (r3 - r1) * (r2 - r0)) * ellipf(phi,k**2)
		else:
			# TODO: single lobe case.
			assert(False)
			pass

def eta(lam):
    """Function from DLMF 8.12.1 shifted to be centered at 0."""
    if lam > 0:
        return mp.sqrt(2*(lam - mp.log(lam + 1)))
    elif lam < 0:
        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
    else:
        return 0

def fp(ne):
    """ 
    Return the plasma frequency. 
    ~ 8.98 * sqrt(ne)

    """
    
    return mp.sqrt(echarge**2/emass/permittivity)/2/mp.pi * mp.sqrt(ne)

 def bimax_integrand(self, z, wc, l, n, t):
     """
     Integrand of electron-noise integral.
     
     """
     return f1(wc*l/z/mp.sqrt(2)) * z * \
         (mp.exp(-z**2) + n/mp.sqrt(t)*mp.exp(-z**2 / t)) / \
         (mp.fabs(BiMax.d_l(z, wc, n, t))**2 * wc**2)

def findCenteredPolygonalNumber( n, k ):
    if real_int( k ) < 3:
        raise ValueError( 'the number of sides of the polygon cannot be less than 3,' )

    s = fdiv( k, 2 )

    return nint( fdiv( fadd( sqrt( s ),
                       sqrt( fsum( [ fmul( 4, real( n ) ), s, -4 ] ) ) ), fmul( 2, sqrt( s ) ) ) )