def test_mpmath_lambda():
    mpmath.mp.dps = 50
    sin02 = mpmath.mpf("0.19866933079506121545941262711838975037020672954020")
    f = lambdify(x, sin(x), "mpmath")
    prec = 1e-49  # mpmath precision is around 50 decimal places
    assert -prec < f(mpmath.mpf("0.2")) - sin02 < prec
    raises(TypeError, lambda: f(x))

def _optimize_single(opt_params):
    '''mapped function for optimization of a single function in a single domain
    
    opt_params is a tuple (expression, subs, bound, use_basinhopping, use_corners)
    returns the interval
    '''

    (e, subs, bound, use_basinhopping, use_corners) = opt_params
    assert not (use_basinhopping and use_corners)

    rv = None
    e = _simplify_eq(e)

    if bound is None:
        rv = _eval_eq_direct(e, subs)
    else:
        under_approx = None

        if use_corners:
            under_approx = _eval_at_corners(e, subs)
        elif use_basinhopping:
            under_approx = _basinhopping(e, subs)
        else:
            under_approx = _eval_at_middle(e, subs)

        # under_approx param should be an array of size 1, since it gets updated
        rv = _eval_eq_bounded(e, subs, bound, [under_approx])

    # convert rv from interval to tuple (for pickling)
    return [float(mpf(rv.a)), float(mpf(rv.b))]

def daubechis(N):
	# make polynomial
	q_y = [sm.binomial(N-1+k,k) for k in reversed(range(N))]

	# get polynomial roots y[k]
	y = sm.mp.polyroots(q_y, maxsteps=200, extraprec=64)

	z = []
	for yk in y:
		# subustitute y = -1/4z + 1/2 - 1/4/z to factor f(y) = y - y[k]
		f = [sm.mpf('-1/4'), sm.mpf('1/2') - yk, sm.mpf('-1/4')]

		# get polynomial roots z[k]
		z += sm.mp.polyroots(f)

	# make polynomial using the roots within unit circle
	h0z = sm.sqrt('2')
	for zk in z:
		if sm.fabs(zk) < 1:
			h0z *= sympy.sympify('(z-zk)/(1-zk)').subs('zk',zk)

	# adapt vanising moments
	hz = (sympy.sympify('(1+z)/2')**N*h0z).expand()

	# get scaling coefficients
	return [sympy.re(hz.coeff('z',k)) for k in reversed(range(N*2))]

def uniform_nodes(n):
  """Return uniform nodes."""

  r = set()
  for i in range(n):
    r.add(mpmath.mpf('-1.0') + i * mpmath.mpf('2.0')/(n-1))

  return map_to_zero_one(r)

def mp_real(x):
    if type(x) == mpc:
        if x.imag == mpf(0):
            return mpf(x.real)
        else:
            raise ValueError
    else:
        return mpf(x)

def gauss_lobatto(n):
  """Return Gauss-Lobatto nodes.

  Gauss-Lobatto nodes are roots of P'_{n-1}(x) and -1, 1.
  """

  x = sympy.var('x')
  p = legendre_poly(n-1).diff(x)
  r = find_roots(p)
  r = [mpmath.mpf('-1.0'), mpmath.mpf('1.0')] + r
  return sorted(r)

def test_mod():
    assert mpf(234) % 1 == 0
    assert mpf(-3) % 256 == 253
    assert mpf(0.25) % 23490.5 == 0.25
    assert mpf(0.25) % -23490.5 == -23490.25
    assert mpf(-0.25) % 23490.5 == 23490.25
    assert mpf(-0.25) % -23490.5 == -0.25
    # Check that these cases are handled efficiently
    assert mpf('1e10000000000') % 1 == 0
    assert mpf('1.23e-1000000000') % 1 == mpf('1.23e-1000000000')
    # test __rmod__
    assert 3 % mpf('1.75') == 1.25

def gauss_lobatto_nodes(n):
  """Return Gauss-Lobatto nodes.

  Before mapping to [0,1], the nodes are: roots of :math:`P'_{n-1}(x)`
  and -1, 1.
  """

  x = sympy.var('x')
  p = legendre_poly(n-1).diff(x)
  r = find_roots(p)
  r = [mpmath.mpf('-1.0'), mpmath.mpf('1.0')] + r

  return map_to_zero_one(r)

def as_mpmath(x, prec, options):
    x = sympify(x)
    if isinstance(x, C.Zero):
        return mpf(0)
    if isinstance(x, C.Infinity):
        return mpf("inf")
    if isinstance(x, C.NegativeInfinity):
        return mpf("-inf")
    # XXX
    re, im, _, _ = evalf(x, prec, options)
    if im:
        return mpc(re or fzero, im)
    return mpf(re)

def mc_compute_stationary(P, precision=None, tol=None):
    n = P.shape[0]

    if precision is None:
        # Compute eigenvalues and eigenvectors
        eigvals, eigvecs = la.eig(P, left=True, right=False)

        # Find the index for unit eigenvalues
        index = np.where(abs(eigvals - 1.) < 1e-12)[0]

        # Pull out the eigenvectors that correspond to unit eig-vals
        uniteigvecs = eigvecs[:, index]

        stationary_dists = uniteigvecs/np.sum(uniteigvecs, axis=0)

    else:
        # Create a list to store eigvals
        stationary_dists_list = []
        if tol is None:
            # If tolerance isn't specified then use 2*precision
            tol = mp.mpf(2 * 10**(-precision + 1))

        with mp.workdps(precision):
            eigvals, eigvecs = mp.eig(mp.matrix(P), left=True, right=False)

            for ind, el in enumerate(eigvals):
                if el>=(mp.mpf(1)-mp.mpf(tol)) and el<=(mp.mpf(1)+mp.mpf(tol)):
                    stationary_dists_list.append(eigvecs[ind, :])

            stationary_dists = np.asarray(stationary_dists_list).T
            stationary_dists = (stationary_dists/sum(stationary_dists)).astype(np.float)


    # Check to make sure all of the elements of invar_dist are positive
    if np.any(stationary_dists < -1e-16):
        warn("Elements of your invariant distribution were negative; " +
              "Re-trying with additional precision")

        if precision is None:
            stationary_dists = mc_compute_stationary(P, precision=18, tol=tol)

        elif precision is not None:
            raise ValueError("Elements of your stationary distribution were" +
                             "negative.  Try computing with higher precision")

    # Since we will be accessing the columns of this matrix, we
    # might consider adding .astype(np.float, order='F') to make it
    # column major at beginning
    return stationary_dists.squeeze()

def test_mpmath_lambda():
    dps = mpmath.mp.dps
    mpmath.mp.dps = 50
    try:
        sin02 = mpmath.mpf("0.19866933079506121545941262711838975037020672954020")
        f = lambdify(x, sin(x), "mpmath")
        prec = 1e-49 # mpmath precision is around 50 decimal places
        assert -prec < f(mpmath.mpf("0.2")) - sin02 < prec
        try:
            f(x) # if this succeeds, it can't be a mpmath function
            assert False
        except TypeError:
            pass
    finally:
        mpmath.mp.dps = dps

 def f(i, j):
     gamDict = self.__genGammaDict(j, gDict)
     gamDict[-1] = mpf('0')
     def g(n):
         return gamDict[n]
     _g = vectorize(g)
     return _g(i)

def test_3250():
    from sympy.mpmath import mpf
    assert str(Float(mpf((1,22,2,22)), '')) == '-88.000'
    assert Float('23.e3', '')._prec == 10
    assert Float('23e3', '')._prec == 20
    assert Float('23000', '')._prec == 20
    assert Float('-23000', '')._prec == 20

    def log(self):
        """
        Logaritmo de un intervalo: 'self.log()'

        NOTA: Si el intervalo contiene al 0, pero no es estrictamente negativo,
        se calcula el logaritmo de la intersecci\'on  del intervalo con el dominio
        natural del logaritmo, i.e., [0,+inf].
        """
        
        if 0 in self:

            domainNatural = Intervalo( 0, mpf('inf') )
            intervalRestricted = self.intersection( domainNatural )

            txt_warning = "\nWARNING: Interval {} contains 0 or negative numbers.\n".format(self)
            print txt_warning
            
            txt_warning = "Restricting to the intersection "\
                "with the natural domain of log(x), i.e. {}\n".format(intervalRestricted)
            print txt_warning
            
            return Intervalo( mp.log(intervalRestricted.lo), mp.log(intervalRestricted.hi) )

        elif 0 > self.hi:
            txt_error = 'Interval {} < 0\nlog(x) cannot be computed '\
                'for negative numbers.'.format(self)
            raise ValueError( txt_error )

        else:
            return Intervalo( mp.log(self.lo), mp.log(self.hi) )

 def __genGammaDict(self, j, gDict):
     if self.isZeroT:
         return dict((k, self.scaledVolt[j]**k /gDict[k]) for k in self.wijnTerms)
     else:
         fg = fgamma(self.scaledVolt[j])
         return dict((k, fgamma(self.scaledVolt[j] + mpf(str(k))) / \
                             (fg*gDict[k])) for k in self.wijnTerms)

def test_math_lambda():
    mpmath.mp.dps = 50
    sin02 = mpmath.mpf("0.19866933079506121545941262711838975037020672954020")
    f = lambdify(x, sin(x), "math")
    prec = 1e-15
    assert -prec < f(0.2) - sin02 < prec
    raises(ValueError,"f(x)") # if this succeeds, it can't be a python math function

 def __init__(self, parameters, maxParameter, 
              g, V, prefac, gtot,
              scaledVolt, T, Vq):
     self.Vq = Vq
     self.parameters = parameters
     self.V = V
     self.T = T
     self.isZeroT = mp.almosteq(self.T, mpf(0))
     self.scaledVolt = scaledVolt
     self.g = g
     self.maxParameter = maxParameter
     self.prefac = prefac
     self.gtot = gtot
     ### input check ###            
     if len(self.maxParameter.shape) == 1:
         self.maxParameter = self.maxParameter[newaxis,:]
     if len(self.prefac.shape) == 1:
         self.prefac = self.prefac[newaxis,:]
     if self.g.shape != self.parameters.shape or \
                 self.V.shape != self.scaledVolt.shape or \
                 self.maxParameter.shape[1] != self.parameters.shape[0]:
         
         raise ValueError
     if len(self.g.shape) == 1:
         self.g = self.g[:, newaxis]
         self.parameters = self.parameters[:, newaxis]
     self.prefac = self.prefac * np.power(self.maxParameter, 
                               -self.scaledVolt[...,newaxis])

 def genScaledVoltage(self):
     if self.isZeroT:
         self.Vq = self.Q*ELEC*self.V / HBAR
         self.scaledVolt = -1j*self.Vq
     else:
         self.Vq = self.Q*ELEC*self.V/(2*pi*BLTZMN*self.T)        
         self.scaledVolt = fmpc(self.gtot /mpf(2), - self.Vq )

 def _f(i, j):
     ldaDICT = dict((k, self.lda[k][j])  for k in self.wijnTerms)
     ldaDICT[-1] = mpf(0)
     def g(n): 
         return ldaDICT[n]
     _g = np.vectorize(g)
     return _g(i)

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision. """
        _prec, mp.prec = mp.prec, prec

        try:
            func = lambdify(self.poly.gen, self.expr)

            interval = self._get_interval()
            refined =  False

            while True:
                if self.is_real:
                    x0 = mpf(str(interval.center))
                else:
                    x0 = mpc(*map(str, interval.center))

                try:
                    root = findroot(func, x0)
                except ValueError:
                    interval = interval.refine()
                    refined = True
                    continue
                else:
                    if refined:
                        self._set_interval(interval)

                    break
        finally:
            mp.prec = _prec

        return Float._new(root.real._mpf_, prec) + I*Float._new(root.imag._mpf_, prec)