def test_stoch_eig_high_prec():
    n = 1e-100
    with mp.workdps(100):
        P = mp.matrix([[1-3*(mp.exp(n)-1), 3*(mp.exp(n)-1)],
                       [mp.exp(n)-1      , 1-(mp.exp(n)-1)]])

    run_stoch_eig(P, verbose=VERBOSE)

def test_gth_solve_high_prec():
    n = 1e-100
    with mp.workdps(100):
        P = mp.matrix([[-3*(mp.exp(n)-1), 3*(mp.exp(n)-1)],
                       [mp.exp(n)-1     , -(mp.exp(n)-1) ]])

    run_gth_solve(P, verbose=VERBOSE)

def test_g():
    # Test data for the g_k. See DLMF 5.11.4.
    with mp.workdps(30):
        g = [mp.mpf(1), mp.mpf(1)/12, mp.mpf(1)/288,
             -mp.mpf(139)/51840, -mp.mpf(571)/2488320,
             mp.mpf(163879)/209018880, mp.mpf(5246819)/75246796800]
        mp_assert_allclose(compute_g(7), g)

	def runTest(self):
		import math
		from .math import cos as pcos, sin as psin
		from .types import polynomial, kronecker_monomial, double, real
		self.assertEqual(math.cos(3),pcos(3))
		self.assertEqual(math.cos(3.1234),pcos(3.1234))
		self.assertEqual(math.sin(3),psin(3))
		self.assertEqual(math.sin(3.1234),psin(3.1234))
		pt = polynomial(double,kronecker_monomial())()
		self.assertEqual(math.cos(3),pcos(pt(3)))
		self.assertEqual(math.cos(-2.456),pcos(pt(2.456)))
		self.assertEqual(math.sin(3),psin(pt(3)))
		self.assertEqual(math.sin(-2.456),-psin(pt(2.456)))
		self.assertRaises(TypeError,lambda : pcos(""))
		self.assertRaises(TypeError,lambda : psin(""))
		try:
			from mpmath import mpf, workdps
			from mpmath import cos as mpcos, sin as mpsin
			pt = polynomial(real,kronecker_monomial())()
			self.assertEqual(mpcos(mpf("1.2345")),pcos(mpf("1.2345")))
			self.assertEqual(mpcos(mpf("3")),pcos(pt(mpf("3"))))
			self.assertEqual(mpcos(mpf("-2.456")),pcos(pt(mpf("-2.456"))))
			self.assertEqual(mpsin(mpf("1.2345")),psin(mpf("1.2345")))
			self.assertEqual(mpsin(mpf("3")),psin(pt(mpf("3"))))
			self.assertEqual(mpsin(mpf("-2.456")),psin(pt(mpf("-2.456"))))
			with workdps(500):
				self.assertEqual(mpcos(mpf("1.2345")),pcos(mpf("1.2345")))
				self.assertEqual(mpcos(mpf("3")),pcos(pt(mpf("3"))))
				self.assertEqual(mpcos(mpf("-2.456")),pcos(pt(mpf("-2.456"))))
				self.assertEqual(mpsin(mpf("1.2345")),psin(mpf("1.2345")))
				self.assertEqual(mpsin(mpf("3")),psin(pt(mpf("3"))))
				self.assertEqual(mpsin(mpf("-2.456")),psin(pt(mpf("-2.456"))))
		except ImportError:
			pass

def field_isomorphism_pslq(a, b):
    """Construct field isomorphism using PSLQ algorithm."""
    if not all(_.domain.is_RationalField and _.ext.is_real for _ in (a, b)):
        raise NotImplementedError("PSLQ doesn't support complex coefficients")

    f = a.minpoly
    g = b.minpoly.replace(f.gen)
    m = b.minpoly.degree()

    for n in mpmath.libmp.libintmath.giant_steps(32, 256):  # pragma: no branch
        with mpmath.workdps(n):
            A = lambdify((), a.ext, "mpmath")()
            B = lambdify((), b.ext, "mpmath")()
            basis = [B**i for i in range(m)] + [A]
            coeffs = mpmath.pslq(basis, maxcoeff=int(1e10), maxsteps=1000)

        if coeffs is None:
            break

        coeffs = [QQ(c, coeffs[-1]) for c in coeffs[:-1]]
        while not coeffs[-1]:
            coeffs.pop()
        coeffs.reverse()

        h = Poly(coeffs, f.gen, domain='QQ')

        if f.compose(h).rem(g).is_zero or f.compose(-h).rem(g).is_zero:
            return [-c for c in coeffs]

def _noncentral_chi_cdf(x, df, nc, dps=None):
    if dps is None:
        dps = mpmath.mp.dps
    x, df, nc = mpmath.mpf(x), mpmath.mpf(df), mpmath.mpf(nc)
    with mpmath.workdps(dps):
        res = mpmath.quad(lambda t: _noncentral_chi_pdf(t, df, nc), [0, x])
        return res

def taylor_series_at_1(N):
    coeffs = []
    with mpmath.workdps(100):
        coeffs.append(-mpmath.euler)
        for n in range(2, N + 1):
            coeffs.append((-1)**n*mpmath.zeta(n)/n)
    return coeffs

def gammaincc(a, x, dps=50, maxterms=10**8):
    """Compute gammaincc exactly like mpmath does but allow for more
    terms in hypercomb. See

    mpmath/functions/expintegrals.py#L187

    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a = a, x
        
        if mp.isint(z):
            try:
                # mpmath has a fast integer path
                return mpf2float(mp.gammainc(z, a=a, regularized=True))
            except mp.libmp.NoConvergence:
                pass
        nega = mp.fneg(a, exact=True)
        G = [z]
        # Use 2F0 series when possible; fall back to lower gamma representation
        try:
            def h(z):
                r = z-1
                return [([mp.exp(nega), a], [1, r], [], G, [1, -r], [], 1/nega)]
            return mpf2float(mp.hypercomb(h, [z], force_series=True))
        except mp.libmp.NoConvergence:
            def h(z):
                T1 = [], [1, z-1], [z], G, [], [], 0
                T2 = [-mp.exp(nega), a, z], [1, z, -1], [], G, [1], [1+z], a
                return T1, T2
            return mpf2float(mp.hypercomb(h, [z], maxterms=maxterms))

def test_alpha():
    # Test data for the alpha_k. See DLMF 8.12.14.
    with mp.workdps(30):
        alpha = [mp.mpf(0), mp.mpf(1), mp.mpf(1)/3, mp.mpf(1)/36,
                 -mp.mpf(1)/270, mp.mpf(1)/4320, mp.mpf(1)/17010,
                 -mp.mpf(139)/5443200, mp.mpf(1)/204120]
        mp_assert_allclose(compute_alpha(9), alpha)

def test_d():
    # Compare the d_{k, n} to the results in appendix F of [1].
    #
    # Sources
    # -------
    # [1] DiDonato and Morris, Computation of the Incomplete Gamma
    #     Function Ratios and their Inverse, ACM Transactions on
    #     Mathematical Software, 1986.

    with mp.workdps(50):
        dataset = [(0, 0, -mp.mpf('0.333333333333333333333333333333')),
                   (0, 12, mp.mpf('0.102618097842403080425739573227e-7')),
                   (1, 0, -mp.mpf('0.185185185185185185185185185185e-2')),
                   (1, 12, mp.mpf('0.119516285997781473243076536700e-7')),
                   (2, 0, mp.mpf('0.413359788359788359788359788360e-2')),
                   (2, 12, -mp.mpf('0.140925299108675210532930244154e-7')),
                   (3, 0, mp.mpf('0.649434156378600823045267489712e-3')),
                   (3, 12, -mp.mpf('0.191111684859736540606728140873e-7')),
                   (4, 0, -mp.mpf('0.861888290916711698604702719929e-3')),
                   (4, 12, mp.mpf('0.288658297427087836297341274604e-7')),
                   (5, 0, -mp.mpf('0.336798553366358150308767592718e-3')),
                   (5, 12, mp.mpf('0.482409670378941807563762631739e-7')),
                   (6, 0, mp.mpf('0.531307936463992223165748542978e-3')),
                   (6, 12, -mp.mpf('0.882860074633048352505085243179e-7')),
                   (7, 0, mp.mpf('0.344367606892377671254279625109e-3')),
                   (7, 12, -mp.mpf('0.175629733590604619378669693914e-6')),
                   (8, 0, -mp.mpf('0.652623918595309418922034919727e-3')),
                   (8, 12, mp.mpf('0.377358774161109793380344937299e-6')),
                   (9, 0, -mp.mpf('0.596761290192746250124390067179e-3')),
                   (9, 12, mp.mpf('0.870823417786464116761231237189e-6'))]
        d = compute_d(10, 13)
        res = [d[k][n] for k, n, std in dataset]
        std = map(lambda x: x[2], dataset)
        mp_assert_allclose(res, std)

def logReg_ObjectiveFunction(X,Y,W,u,s,coef_transfer,b):
    w=np.array([W[1:]]).T
    w_0=W[0]
    A=np.dot(X,w)+float(w_0)
    B=-np.multiply(Y,A)
    with mp.workdps(30):
        NLL=sum(log_one_plus_exp(B))+coef_transfer*0.5*GPriorWeightRegTerm(W,u,s)#+(1/float(np.shape(Y)[0]))*np.dot(b.T,np.array([W]).T)[0,0]
    return NLL

def zetac_series(N):
    coeffs = []
    with mpmath.workdps(100):
        coeffs.append(-1.5)
        for n in range(1, N):
            coeff = mpmath.diff(mpmath.zeta, 0, n)/mpmath.factorial(n)
            coeffs.append(coeff)
    return coeffs

def _student_t_cdf(df, t, dps=None):
    if dps is None:
        dps = mpmath.mp.dps
    with mpmath.workdps(dps):
        df, t = mpmath.mpf(df), mpmath.mpf(t)
        fac = mpmath.hyp2f1(0.5, 0.5*(df + 1), 1.5, -t**2/df)
        fac *= t*mpmath.gamma(0.5*(df + 1))
        fac /= mpmath.sqrt(mpmath.pi*df)*mpmath.gamma(0.5*df)
        return 0.5 + fac

 def idmap(self, *args):
     if self.spfunc_first:
         res = self.spfunc(*args)
         if np.isnan(res):
             return np.nan
         args = list(args)
         args[self.index] = res
         with mpmath.workdps(self.dps):
             res = self.mpfunc(*tuple(args))
             # Imaginary parts are spurious
             res = mpf2float(res.real)
     else:
         with mpmath.workdps(self.dps):
             res = self.mpfunc(*args)
             res = mpf2float(res.real)
         args = list(args)
         args[self.index] = res
         res = self.spfunc(*tuple(args))
     return res

def main():
    print(__doc__)
    with mpmath.workdps(50):
        p, q = lambertw_pade()
        p, q = p[::-1], q[::-1]
        print("p = {}".format(p))
        print("q = {}".format(q))

    x, y = np.linspace(-1.5, 1.5, 75), np.linspace(-1.5, 1.5, 75)
    x, y = np.meshgrid(x, y)
    z = x + 1j*y
    lambertw_std = []
    for z0 in z.flatten():
        lambertw_std.append(complex(mpmath.lambertw(z0)))
    lambertw_std = np.array(lambertw_std).reshape(x.shape)

    fig, axes = plt.subplots(nrows=3, ncols=1)
    # Compare Pade approximation to true result
    p = np.array([float(p0) for p0 in p])
    q = np.array([float(q0) for q0 in q])
    pade_approx = np.polyval(p, z)/np.polyval(q, z)
    pade_err = abs(pade_approx - lambertw_std)
    axes[0].pcolormesh(x, y, pade_err)
    # Compare two terms of asymptotic series to true result
    asy_approx = np.log(z) - np.log(np.log(z))
    asy_err = abs(asy_approx - lambertw_std)
    axes[1].pcolormesh(x, y, asy_err)
    # Compare two terms of the series around the branch point to the
    # true result
    p = np.sqrt(2*(np.exp(1)*z + 1))
    series_approx = -1 + p - p**2/3
    series_err = abs(series_approx - lambertw_std)
    im = axes[2].pcolormesh(x, y, series_err)

    fig.colorbar(im, ax=axes.ravel().tolist())
    plt.show()

    fig, ax = plt.subplots(nrows=1, ncols=1)
    pade_better = pade_err < asy_err
    im = ax.pcolormesh(x, y, pade_better)
    t = np.linspace(-0.3, 0.3)
    ax.plot(-2.5*abs(t) - 0.2, t, 'r')
    fig.colorbar(im, ax=ax)
    plt.show()

def multinomial(n,m,multiplePrecision=False,decimalPlaces=40):
    '''
    n = int
    m = list of integers summing to n
    '''
    mybinom = binom
    if multiplePrecision: #use multiple-precision arithmetics
        with mp.workdps(decimalPlaces):
            if sum(m) != n:
                return mp.mpf('0.0')
            else:
                return np.prod([mp.mpf(mybinom(sum(m[i:]),m[i])) for i in xrange(len(m)-1)])
    else: #use floating-point arithmetics
        if sum(m) != n:
            return 0.0
        else:
    #        if len(m)==1: m = list(m)+[0] #else the reduce statement does not work
    #        return reduce(lambda x,y:x*y,[mybinom(sum(m[i:]),m[i]) for i in xrange(len(m)-1)])
            return np.prod([mybinom(sum(m[i:]),m[i]) for i in xrange(len(m)-1)])

def main():
    print(__doc__)
    K = 25
    N = 25
    with mp.workdps(50):
        d = compute_d(K, N)
    fn = os.path.join(os.path.dirname(__file__), '..', 'cephes', 'igam.h')
    with open(fn + '.new', 'w') as f:
        f.write(header.format(K, N))
        for k, row in enumerate(d):
            row = map(lambda x: mp.nstr(x, 17, min_fixed=0, max_fixed=0), row)
            f.write('{')
            f.write(", ".join(row))
            if k < K - 1:
                f.write('},\n')
            else:
                f.write('}};\n')
        f.write(footer)
    os.rename(fn + '.new', fn)

	def runTest(self):
		import math
		from .math import cos as pcos, sin as psin
		from .types import polynomial, k_monomial, double, real
		# NOTE: according to
		# https://docs.python.org/3/library/math.html
		# the CPython math functions are wrappers around the corresponding
		# C functions. The results should thus be the same.
		self.assertEqual(math.cos(3.),pcos(3.))
		self.assertEqual(math.cos(3.1234),pcos(3.1234))
		self.assertEqual(math.sin(3.),psin(3.))
		self.assertEqual(math.sin(3.1234),psin(3.1234))
		pt = polynomial(double,k_monomial)()
		self.assertEqual(math.cos(3),pcos(pt(3)))
		self.assertEqual(math.cos(2.456),pcos(pt(2.456)))
		self.assertEqual(math.sin(3),psin(pt(3)))
		self.assertEqual(math.sin(-2.456),psin(pt(-2.456)))
		self.assertRaises(TypeError,lambda : pcos(""))
		self.assertRaises(TypeError,lambda : psin(""))
		try:
			from mpmath import mpf, workdps
			from mpmath import cos as mpcos, sin as mpsin
			pt = polynomial(real,k_monomial)()
			self.assertEqual(mpcos(mpf("1.2345")),pcos(mpf("1.2345")))
			self.assertEqual(mpcos(mpf("3")),pcos(pt(mpf("3"))))
			self.assertEqual(mpcos(mpf("-2.456")),pcos(pt(mpf("-2.456"))))
			self.assertEqual(mpsin(mpf("1.2345")),psin(mpf("1.2345")))
			self.assertEqual(mpsin(mpf("3")),psin(pt(mpf("3"))))
			self.assertEqual(mpsin(mpf("-2.456")),psin(pt(mpf("-2.456"))))
			with workdps(500):
				self.assertEqual(mpcos(mpf("1.2345")),pcos(mpf("1.2345")))
				self.assertEqual(mpcos(mpf("3")),pcos(pt(mpf("3"))))
				self.assertEqual(mpcos(mpf("-2.456")),pcos(pt(mpf("-2.456"))))
				self.assertEqual(mpsin(mpf("1.2345")),psin(mpf("1.2345")))
				self.assertEqual(mpsin(mpf("3")),psin(pt(mpf("3"))))
				self.assertEqual(mpsin(mpf("-2.456")),psin(pt(mpf("-2.456"))))
		except ImportError:
			pass
		self.binomialTest()
		self.sincosTest()
		self.evaluateTest()
		self.subsTest()
		self.invertTest()

def gammainc(a, x, dps=50, maxterms=10**8):
    """Compute gammainc exactly like mpmath does but allow for more
    summands in hypercomb. See

    mpmath/functions/expintegrals.py#L134
    
    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a, b = mp.mpf(a), mp.mpf(x), mp.mpf(x)
        G = [z]
        negb = mp.fneg(b, exact=True)

        def h(z):
            T1 = [mp.exp(negb), b, z], [1, z, -1], [], G, [1], [1+z], b
            return (T1,)

        res = mp.hypercomb(h, [z], maxterms=maxterms)
        return mpf2float(res)

 def test_log(self):
     with mp.workdps(30):
         logcoeffs = mp.taylor(lambda x: mp.log(1 + x), 0, 10)
         expcoeffs = mp.taylor(lambda x: mp.exp(x) - 1, 0, 10)
         invlogcoeffs = lagrange_inversion(logcoeffs)
         mp_assert_allclose(invlogcoeffs, expcoeffs)