def _display_expression(self, keys, user_substitutes={}):
        """Helper function for human friendly display of the symbolic components."""
        # Create some pretty maths symbols for the display.
        sigma, alpha, nu, omega, l, variance = sym.var('\sigma, \alpha, \nu, \omega, \ell, \sigma^2')
        substitutes = {'scale': sigma, 'shape': alpha, 'lengthscale': l, 'variance': variance}
        substitutes.update(user_substitutes)

        function_substitutes = {normcdfln : lambda arg : sym.log(normcdf(arg)),
                                logisticln : lambda arg : -sym.log(1+sym.exp(-arg)),
                                logistic : lambda arg : 1/(1+sym.exp(-arg)),
                                erfcx : lambda arg : erfc(arg)/sym.exp(arg*arg),
                                gammaln : lambda arg : sym.log(sym.gamma(arg))}
        expr = getFromDict(self.expressions, keys)
        for var_name, sub in self.variable_sort(self.expressions['update_cache'], reverse=True):
            for var in self.variables['cache']:
                if var_name == var.name:
                    expr = expr.subs(var, sub)
                    break
        for var_name, sub in self.variable_sort(self.expressions['parameters_changed'], reverse=True):
            for var in self.variables['sub']:
                if var_name == var.name:
                    expr = expr.subs(var, sub)
                    break

        for var_name, sub in self.variable_sort(substitutes, reverse=True):
            for var in self.variables['theta']:
                if var_name == var.name:
                    expr = expr.subs(var, sub)
                    break
        for m, r in function_substitutes.iteritems():
            expr = expr.replace(m, r)#normcdfln, lambda arg : sym.log(normcdf(arg)))
        return expr.simplify()

def test_issue_1572_1364_1368():
    assert solve((sqrt(x**2 - 1) - 2)) in ([sqrt(5), -sqrt(5)],
                                           [-sqrt(5), sqrt(5)])
    assert solve((2**exp(y**2/x) + 2)/(x**2 + 15), y) == (
        [-sqrt(x)*sqrt(log((log(2) + I*pi)/log(2))),
          sqrt(x)*sqrt(log((log(2) + I*pi)/log(2)))]
          )

    C1, C2 = symbols('C1 C2')
    f = Function('f')
    assert solve(C1 + C2/x**2 - exp(-f(x)), f(x)) == [log(x**2/(C1*x**2 + C2))]
    a = symbols('a')
    E = S.Exp1
    assert solve(1 - log(a + 4*x**2), x) in (
                                        [-sqrt(-a + E)/2, sqrt(-a + E)/2],
                                        [sqrt(-a + E)/2, -sqrt(-a + E)/2]
                                        )
    assert solve(log(a**(-3) - x**2)/a, x) in (
                            [-sqrt(-1 + a**(-3)), sqrt(-1 + a**(-3))],
                            [sqrt(-1 + a**(-3)), -sqrt(-1 + a**(-3))],)
    assert solve(1 - log(a + 4*x**2), x) in (
                                             [-sqrt(-a + E)/2, sqrt(-a + E)/2],
                                             [sqrt(-a + E)/2, -sqrt(-a + E)/2],)
    assert solve((a**2 + 1) * (sin(a*x) + cos(a*x)), x) == [-pi/(4*a), 3*pi/(4*a)]
    assert solve(3 - (sinh(a*x) + cosh(a*x)), x) == [2*atanh(S.Half)/a]
    assert solve(3-(sinh(a*x) + cosh(a*x)**2), x) == \
             [
             2*atanh(-1 + sqrt(2))/a,
             2*atanh(S(1)/2 + sqrt(5)/2)/a,
             2*atanh(-sqrt(2) - 1)/a,
             2*atanh(-sqrt(5)/2 + S(1)/2)/a
             ]
    assert solve(atan(x) - 1) == [tan(1)]

def test_acceleration():
    e = (1 + 1/n)**n
    assert round(richardson(e, n, 10, 20).evalf(), 10) == round(E.evalf(), 10)

    A = Sum(Integer(-1)**(k + 1) / k, (k, 1, n))
    assert round(shanks(A, n, 25).evalf(), 4) == round(log(2).evalf(), 4)
    assert round(shanks(A, n, 25, 5).evalf(), 10) == round(log(2).evalf(), 10)

def test_manualintegrate_exponentials():
    assert manualintegrate(exp(2*x), x) == exp(2*x) / 2
    assert manualintegrate(2**x, x) == (2 ** x) / log(2)

    assert manualintegrate(1 / x, x) == log(x)
    assert manualintegrate(1 / (2*x + 3), x) == log(2*x + 3) / 2
    assert manualintegrate(log(x)**2 / x, x) == log(x)**3 / 3

def test_ci():
    m1 = exp_polar(I*pi)
    m1_ = exp_polar(-I*pi)
    pI = exp_polar(I*pi/2)
    mI = exp_polar(-I*pi/2)

    assert Ci(m1*x) == Ci(x) + I*pi
    assert Ci(m1_*x) == Ci(x) - I*pi
    assert Ci(pI*x) == Chi(x) + I*pi/2
    assert Ci(mI*x) == Chi(x) - I*pi/2
    assert Chi(m1*x) == Chi(x) + I*pi
    assert Chi(m1_*x) == Chi(x) - I*pi
    assert Chi(pI*x) == Ci(x) + I*pi/2
    assert Chi(mI*x) == Ci(x) - I*pi/2
    assert Ci(exp_polar(2*I*pi)*x) == Ci(x) + 2*I*pi
    assert Chi(exp_polar(-2*I*pi)*x) == Chi(x) - 2*I*pi
    assert Chi(exp_polar(2*I*pi)*x) == Chi(x) + 2*I*pi
    assert Ci(exp_polar(-2*I*pi)*x) == Ci(x) - 2*I*pi

    assert mytd(Ci(x), cos(x)/x, x)
    assert mytd(Chi(x), cosh(x)/x, x)

    assert mytn(Ci(x), Ci(x).rewrite(Ei),
                Ei(x*exp_polar(-I*pi/2))/2 + Ei(x*exp_polar(I*pi/2))/2, x)
    assert mytn(Chi(x), Chi(x).rewrite(Ei),
                Ei(x)/2 + Ei(x*exp_polar(I*pi))/2 - I*pi/2, x)

    assert tn_arg(Ci)
    assert tn_arg(Chi)

    from sympy import O, EulerGamma, log, limit
    assert Ci(x).nseries(x, n=4) == EulerGamma + log(x) - x**2/4 + x**4/96 + O(x**5)
    assert Chi(x).nseries(x, n=4) == EulerGamma + log(x) + x**2/4 + x**4/96 + O(x**5)
    assert limit(log(x) - Ci(2*x), x, 0) == -log(2) - EulerGamma

def test_messy():
    from sympy import (laplace_transform, Si, Shi, Chi, atan, Piecewise,
                       acoth, E1, besselj, acosh, asin, And, re,
                       fourier_transform, sqrt)
    assert laplace_transform(Si(x), x, s) == ((-atan(s) + pi/2)/s, 0, True)

    assert laplace_transform(Shi(x), x, s) == (acoth(s)/s, 1, True)

    # where should the logs be simplified?
    assert laplace_transform(Chi(x), x, s) == \
        ((log(s**(-2)) - log((s**2 - 1)/s**2))/(2*s), 1, True)

    # TODO maybe simplify the inequalities?
    assert laplace_transform(besselj(a, x), x, s)[1:] == \
        (0, And(S(0) < re(a/2) + S(1)/2, S(0) < re(a/2) + 1))

    # NOTE s < 0 can be done, but argument reduction is not good enough yet
    assert fourier_transform(besselj(1, x)/x, x, s, noconds=False) == \
        (Piecewise((0, 4*abs(pi**2*s**2) > 1),
                   (2*sqrt(-4*pi**2*s**2 + 1), True)), s > 0)
    # TODO FT(besselj(0,x)) - conditions are messy (but for acceptable reasons)
    #                       - folding could be better

    assert integrate(E1(x)*besselj(0, x), (x, 0, oo), meijerg=True) == \
        log(1 + sqrt(2))
    assert integrate(E1(x)*besselj(1, x), (x, 0, oo), meijerg=True) == \
        log(S(1)/2 + sqrt(2)/2)

    assert integrate(1/x/sqrt(1 - x**2), x, meijerg=True) == \
        Piecewise((-acosh(1/x), 1 < abs(x**(-2))), (I*asin(1/x), True))

def test_fraction():
    x, y, z = map(Symbol, 'xyz')
    A = Symbol('A', commutative=False)

    assert fraction(Rational(1, 2)) == (1, 2)

    assert fraction(x) == (x, 1)
    assert fraction(1/x) == (1, x)
    assert fraction(x/y) == (x, y)
    assert fraction(x/2) == (x, 2)

    assert fraction(x*y/z) == (x*y, z)
    assert fraction(x/(y*z)) == (x, y*z)

    assert fraction(1/y**2) == (1, y**2)
    assert fraction(x/y**2) == (x, y**2)

    assert fraction((x**2 + 1)/y) == (x**2 + 1, y)
    assert fraction(x*(y + 1)/y**7) == (x*(y + 1), y**7)

    assert fraction(exp(-x), exact=True) == (exp(-x), 1)
    assert fraction((1/(x + y))/2, exact=True) == (1, Mul(2,(x + y), evaluate=False))

    assert fraction(x*A/y) == (x*A, y)
    assert fraction(x*A**-1/y) == (x*A**-1, y)

    n = symbols('n', negative=True)
    assert fraction(exp(n)) == (1, exp(-n))
    assert fraction(exp(-n)) == (exp(-n), 1)

    p = symbols('p', positive=True)
    assert fraction(exp(-p)*log(p), exact=True) == (exp(-p)*log(p), 1)

def test_residue_reduce():
    a = Poly(2*t**2 - t - x**2, t)
    d = Poly(t**3 - x**2*t, t)
    DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1/x, t)], 'Tfuncs': [log]})
    assert residue_reduce(a, d, DE, z, invert=False) == \
        ([(Poly(z**2 - S(1)/4, z), Poly((1 + 3*x*z - 6*z**2 -
        2*x**2 + 4*x**2*z**2)*t - x*z + x**2 + 2*x**2*z**2 - 2*z*x**3, t))], False)
    assert residue_reduce(a, d, DE, z, invert=True) == \
        ([(Poly(z**2 - S(1)/4, z), Poly(t + 2*x*z, t))], False)
    assert residue_reduce(Poly(-2/x, t), Poly(t**2 - 1, t,), DE, z, invert=False) == \
        ([(Poly(z**2 - 1, z), Poly(-2*z*t/x - 2/x, t))], True)
    ans = residue_reduce(Poly(-2/x, t), Poly(t**2 - 1, t), DE, z, invert=True)
    assert ans == ([(Poly(z**2 - 1, z), Poly(t + z, t))], True)
    assert residue_reduce_to_basic(ans[0], DE, z) == -log(-1 + log(x)) + log(1 + log(x))

    DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(-t**2 - t/x - (1 - nu**2/x**2), t)]})
    # TODO: Skip or make faster
    assert residue_reduce(Poly((-2*nu**2 - x**4)/(2*x**2)*t - (1 + x**2)/x, t),
    Poly(t**2 + 1 + x**2/2, t), DE, z) == \
        ([(Poly(z + S(1)/2, z, domain='QQ'), Poly(t**2 + 1 + x**2/2, t, domain='EX'))], True)
    DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1 + t**2, t)]})
    assert residue_reduce(Poly(-2*x*t + 1 - x**2, t),
    Poly(t**2 + 2*x*t + 1 + x**2, t), DE, z) == \
        ([(Poly(z**2 + S(1)/4, z), Poly(t + x + 2*z, t))], True)
    DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t, t)]})
    assert residue_reduce(Poly(t, t), Poly(t + sqrt(2), t), DE, z) == \
        ([(Poly(z - 1, z), Poly(t + sqrt(2), t))], True)

def test_DifferentialExtension_all_attrs():
    # Test 'unimportant' attributes
    DE = DifferentialExtension(exp(x)*log(x), x, dummy=False, handle_first='exp')
    assert DE.f == exp(x)*log(x)
    assert DE.newf == t0*t1
    assert DE.x == x
    assert DE.cases == ['base', 'exp', 'primitive']
    assert DE.case == 'primitive'

    assert DE.level == -1
    assert DE.t == t1 == DE.T[DE.level]
    assert DE.d == Poly(1/x, t1) == DE.D[DE.level]
    raises(ValueError, lambda: DE.increment_level())
    DE.decrement_level()
    assert DE.level == -2
    assert DE.t == t0 == DE.T[DE.level]
    assert DE.d == Poly(t0, t0) == DE.D[DE.level]
    assert DE.case == 'exp'
    DE.decrement_level()
    assert DE.level == -3
    assert DE.t == x == DE.T[DE.level] == DE.x
    assert DE.d == Poly(1, x) == DE.D[DE.level]
    assert DE.case == 'base'
    raises(ValueError, lambda: DE.decrement_level())
    DE.increment_level()
    DE.increment_level()
    assert DE.level == -1
    assert DE.t == t1 == DE.T[DE.level]
    assert DE.d == Poly(1/x, t1) == DE.D[DE.level]
    assert DE.case == 'primitive'

def test_gruntz_evaluation_slow():
    sskip()
    assert gruntz((exp(exp(-x/(1+exp(-x))))*exp(-x/(1+exp(-x/(1+exp(-x)))))
                   *exp(exp(-x+exp(-x/(1+exp(-x))))))
                  / (exp(-x/(1+exp(-x))))**2 - exp(x) + x, x, oo) == 2
    assert gruntz(exp(exp(exp(x)/(1-1/x)))
                  - exp(exp(exp(x)/(1-1/x-log(x)**(-log(x))))), x, oo) == -oo

def test_ei():
    pos = Symbol('p', positive=True)
    neg = Symbol('n', negative=True)
    assert Ei(-pos) == Ei(polar_lift(-1)*pos) - I*pi
    assert Ei(neg) == Ei(polar_lift(neg)) - I*pi
    assert tn_branch(Ei)
    assert mytd(Ei(x), exp(x)/x, x)
    assert mytn(Ei(x), Ei(x).rewrite(uppergamma),
                -uppergamma(0, x*polar_lift(-1)) - I*pi, x)
    assert mytn(Ei(x), Ei(x).rewrite(expint),
                -expint(1, x*polar_lift(-1)) - I*pi, x)
    assert Ei(x).rewrite(expint).rewrite(Ei) == Ei(x)
    assert Ei(x*exp_polar(2*I*pi)) == Ei(x) + 2*I*pi
    assert Ei(x*exp_polar(-2*I*pi)) == Ei(x) - 2*I*pi

    assert mytn(Ei(x), Ei(x).rewrite(Shi), Chi(x) + Shi(x), x)
    assert mytn(Ei(x*polar_lift(I)), Ei(x*polar_lift(I)).rewrite(Si),
                Ci(x) + I*Si(x) + I*pi/2, x)

    assert Ei(log(x)).rewrite(li) == li(x)
    assert Ei(2*log(x)).rewrite(li) == li(x**2)

    assert gruntz(Ei(x+exp(-x))*exp(-x)*x, x, oo) == 1

    assert Ei(x).series(x) == EulerGamma + log(x) + x + x**2/4 + \
        x**3/18 + x**4/96 + x**5/600 + O(x**6)

 def ideal_mixing_energy(self, phase, symbols, param_search):
     #pylint: disable=W0613
     """
     Returns the ideal mixing energy in symbolic form.
     """
     # Normalize site ratios
     site_ratio_normalization = self._site_ratio_normalization(phase)
     site_ratios = phase.sublattices
     site_ratios = [c/site_ratio_normalization for c in site_ratios]
     ideal_mixing_term = S.Zero
     for subl_index, sublattice in enumerate(phase.constituents):
         active_comps = set(sublattice).intersection(self.components)
         if len(active_comps) == 1:
             continue # no mixing if only one species in sublattice
         ratio = site_ratios[subl_index]
         for comp in active_comps:
             sitefrac = \
                 v.SiteFraction(phase.name, subl_index, comp)
             # -35.8413*x term is there to keep derivative continuous
             mixing_term = Piecewise((sitefrac * log(sitefrac), \
                 sitefrac > 1e-12), ((3.0*log(10.0)/(2.5e11))+ \
                 (sitefrac-1e-12)*(1.0-log(1e-12))+(5e-11) * \
                 (sitefrac-1e-12)**2, True))
             ideal_mixing_term += (mixing_term*ratio)
     ideal_mixing_term *= (v.R * v.T)
     return ideal_mixing_term

def test_heurisch_trigonometric():
    assert heurisch(sin(x), x) == -cos(x)
    assert heurisch(pi*sin(x) + 1, x) == x - pi*cos(x)

    assert heurisch(cos(x), x) == sin(x)
    assert heurisch(tan(x), x) in [
        log(1 + tan(x)**2)/2,
        log(tan(x) + I) + I*x,
        log(tan(x) - I) - I*x,
    ]

    assert heurisch(sin(x)*sin(y), x) == -cos(x)*sin(y)
    assert heurisch(sin(x)*sin(y), y) == -cos(y)*sin(x)

    # gives sin(x) in answer when run via setup.py and cos(x) when run via py.test
    assert heurisch(sin(x)*cos(x), x) in [sin(x)**2 / 2, -cos(x)**2 / 2]
    assert heurisch(cos(x)/sin(x), x) == log(sin(x))

    assert heurisch(x*sin(7*x), x) == sin(7*x) / 49 - x*cos(7*x) / 7
    assert heurisch(1/pi/4 * x**2*cos(x), x) == 1/pi/4*(x**2*sin(x) -
                    2*sin(x) + 2*x*cos(x))

    assert heurisch(acos(x/4) * asin(x/4), x) == 2*x - (sqrt(16 - x**2))*asin(x/4) \
        + (sqrt(16 - x**2))*acos(x/4) + x*asin(x/4)*acos(x/4)

    assert heurisch(sin(x)/(cos(x)**2+1), x) == -atan(cos(x)) #fixes issue 13723
    assert heurisch(1/(cos(x)+2), x) == 2*sqrt(3)*atan(sqrt(3)*tan(x/2)/3)/3
    assert heurisch(2*sin(x)*cos(x)/(sin(x)**4 + 1), x) == atan(sqrt(2)*sin(x)
        - 1) - atan(sqrt(2)*sin(x) + 1)

    assert heurisch(1/cosh(x), x) == 2*atan(tanh(x/2))

def test_to_expr():
    x = symbols('x')
    R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx')
    p = HolonomicFunction(Dx - 1, x, 0, [1]).to_expr()
    q = exp(x)
    assert p == q
    p = HolonomicFunction(Dx**2 + 1, x, 0, [1, 0]).to_expr()
    q = cos(x)
    assert p == q
    p = HolonomicFunction(Dx**2 - 1, x, 0, [1, 0]).to_expr()
    q = cosh(x)
    assert p == q
    p = HolonomicFunction(2 + (4*x - 1)*Dx + \
        (x**2 - x)*Dx**2, x, 0, [1, 2]).to_expr().expand()
    q = 1/(x**2 - 2*x + 1)
    assert p == q
    p = expr_to_holonomic(sin(x)**2/x).integrate((x, 0, x)).to_expr()
    q = (sin(x)**2/x).integrate((x, 0, x))
    assert p == q
    C_0, C_1, C_2, C_3 = symbols('C_0, C_1, C_2, C_3')
    p = expr_to_holonomic(log(1+x**2)).to_expr()
    q = C_2*log(x**2 + 1)
    assert p == q
    p = expr_to_holonomic(log(1+x**2)).diff().to_expr()
    q = C_0*x/(x**2 + 1)
    assert p == q
    p = expr_to_holonomic(erf(x) + x).to_expr()
    q = 3*C_3*x - 3*sqrt(pi)*C_3*erf(x)/2 + x + 2*x/sqrt(pi)
    assert p == q
    p = expr_to_holonomic(sqrt(x), x0=1).to_expr()
    assert p == sqrt(x)
    assert expr_to_holonomic(sqrt(x)).to_expr() == sqrt(x)
    p = expr_to_holonomic(sqrt(1 + x**2)).to_expr()
    assert p == sqrt(1+x**2)
    p = expr_to_holonomic((2*x**2 + 1)**(S(2)/3)).to_expr()
    assert p == (2*x**2 + 1)**(S(2)/3)
    p = expr_to_holonomic(sqrt(-x**2+2*x)).to_expr()
    assert p == sqrt(x)*sqrt(-x + 2)
    p = expr_to_holonomic((-2*x**3+7*x)**(S(2)/3)).to_expr()
    q = x**(S(2)/3)*(-2*x**2 + 7)**(S(2)/3)
    assert p == q
    p = from_hyper(hyper((-2, -3), (S(1)/2, ), x))
    s = hyperexpand(hyper((-2, -3), (S(1)/2, ), x))
    D_0 = Symbol('D_0')
    C_0 = Symbol('C_0')
    assert (p.to_expr().subs({C_0:1, D_0:0}) - s).simplify() == 0
    p.y0 = {0: [1], S(1)/2: [0]}
    assert p.to_expr() == s
    assert expr_to_holonomic(x**5).to_expr() == x**5
    assert expr_to_holonomic(2*x**3-3*x**2).to_expr().expand() == \
        2*x**3-3*x**2
    a = symbols("a")
    p = (expr_to_holonomic(1.4*x)*expr_to_holonomic(a*x, x)).to_expr()
    q = 1.4*a*x**2
    assert p == q
    p = (expr_to_holonomic(1.4*x)+expr_to_holonomic(a*x, x)).to_expr()
    q = x*(a + 1.4)
    assert p == q
    p = (expr_to_holonomic(1.4*x)+expr_to_holonomic(x)).to_expr()
    assert p == 2.4*x

def test_asinh():
    x, y = symbols('x,y')
    assert asinh(x) == asinh(x)
    assert asinh(-x) == -asinh(x)
    assert asinh(nan) == nan
    assert asinh( 0) == 0
    assert asinh(+1) == log(sqrt(2) + 1)

    assert asinh(-1) == log(sqrt(2) - 1)
    assert asinh(I) == pi*I/2
    assert asinh(-I) == -pi*I/2
    assert asinh(I/2) == pi*I/6
    assert asinh(-I/2) == -pi*I/6

    assert asinh(oo) == oo
    assert asinh(-oo) == -oo

    assert asinh(I*oo) == oo
    assert asinh(-I *oo) == -oo

    assert asinh(zoo) == zoo

    assert asinh(I *(sqrt(3) - 1)/(2**(S(3)/2))) == pi*I/12
    assert asinh(-I *(sqrt(3) - 1)/(2**(S(3)/2))) == -pi*I/12

    assert asinh(I*(sqrt(5) - 1)/4) == pi*I/10
    assert asinh(-I*(sqrt(5) - 1)/4) == -pi*I/10

    assert asinh(I*(sqrt(5) + 1)/4) == 3*pi*I/10
    assert asinh(-I*(sqrt(5) + 1)/4) == -3*pi*I/10

def test_seriesbug2c():
    w = Symbol("w", real=True)
    #more complicated case, but sin(x)~x, so the result is the same as in (1)
    e=(sin(2*w)/w)**(1+w)
    assert e.nseries(w,0,1) == 2 + O(w)
    assert e.nseries(w,0,3) == 2-Rational(4,3)*w**2+w**2*log(2)**2+2*w*log(2)+O(w**3, w)
    assert e.nseries(w,0,2).subs(w,0) == 2

def test_acosh():
    # TODO please write more tests  -- see issue 3751
    # From http://functions.wolfram.com/ElementaryFunctions/ArcCosh/03/01/
    # at specific points
    assert acosh(1) == 0
    assert acosh(-1) == pi*I
    assert acosh(0) == I*pi/2
    assert acosh(Rational(1, 2)) == I*pi/3
    assert acosh(Rational(-1, 2)) == 2*pi*I/3

    assert acosh(zoo) == oo

    assert acosh(I) == log(I*(1 + sqrt(2)))
    assert acosh(-I) == log(-I*(1 + sqrt(2)))
    assert acosh((sqrt(3) - 1)/(2*sqrt(2))) == 5*pi*I/12
    assert acosh(-(sqrt(3) - 1)/(2*sqrt(2))) == 7*pi*I/12
    assert acosh(sqrt(2)/2) == I*pi/4
    assert acosh(-sqrt(2)/2) == 3*I*pi/4
    assert acosh(sqrt(3)/2) == I*pi/6
    assert acosh(-sqrt(3)/2) == 5*I*pi/6
    assert acosh(sqrt(2 + sqrt(2))/2) == I*pi/8
    assert acosh(-sqrt(2 + sqrt(2))/2) == 7*I*pi/8
    assert acosh(sqrt(2 - sqrt(2))/2) == 3*I*pi/8
    assert acosh(-sqrt(2 - sqrt(2))/2) == 5*I*pi/8
    assert acosh((1 + sqrt(3))/(2*sqrt(2))) == I*pi/12
    assert acosh(-(1 + sqrt(3))/(2*sqrt(2))) == 11*I*pi/12
    assert acosh((sqrt(5) + 1)/4) == I*pi/5
    assert acosh(-(sqrt(5) + 1)/4) == 4*I*pi/5

def test_ratsimp():
    f, g = 1/x + 1/y, (x + y)/(x*y)

    assert f != g and ratsimp(f) == g

    f, g = 1/(1 + 1/x), 1 - 1/(x + 1)

    assert f != g and ratsimp(f) == g

    f, g = x/(x + y) + y/(x + y), 1

    assert f != g and ratsimp(f) == g

    f, g = -x - y - y**2/(x + y) + x**2/(x + y), -2*y

    assert f != g and ratsimp(f) == g

    f = (a*c*x*y + a*c*z - b*d*x*y - b*d*z - b*t*x*y - b*t*x - b*t*z + e*x)/(x*y + z)
    G = [a*c - b*d - b*t + (-b*t*x + e*x)/(x*y + z),
         a*c - b*d - b*t - ( b*t*x - e*x)/(x*y + z)]

    assert f != g and ratsimp(f) in G

    A = sqrt(pi)

    B = log(erf(x) - 1)
    C = log(erf(x) + 1)

    D = 8 - 8*erf(x)

    f = A*B/D - A*C/D + A*C*erf(x)/D - A*B*erf(x)/D + 2*A/D

    assert ratsimp(f) == A*B/8 - A*C/8 - A/(4*erf(x) - 4)

def test_derivative_by_array():
    from sympy.abc import a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z

    bexpr = x*y**2*exp(z)*log(t)
    sexpr = sin(bexpr)
    cexpr = cos(bexpr)

    a = Array([sexpr])

    assert derive_by_array(sexpr, t) == x*y**2*exp(z)*cos(x*y**2*exp(z)*log(t))/t
    assert derive_by_array(sexpr, [x, y, z]) == Array([bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr, bexpr*cexpr])
    assert derive_by_array(a, [x, y, z]) == Array([[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr], [bexpr*cexpr]])

    assert derive_by_array(sexpr, [[x, y], [z, t]]) == Array([[bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr], [bexpr*cexpr, bexpr/log(t)/t*cexpr]])
    assert derive_by_array(a, [[x, y], [z, t]]) == Array([[[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr]], [[bexpr*cexpr], [bexpr/log(t)/t*cexpr]]])
    assert derive_by_array([[x, y], [z, t]], [x, y]) == Array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]])
    assert derive_by_array([[x, y], [z, t]], [[x, y], [z, t]]) == Array([[[[1, 0], [0, 0]], [[0, 1], [0, 0]]],
                                                                         [[[0, 0], [1, 0]], [[0, 0], [0, 1]]]])

    assert diff(sexpr, t) == x*y**2*exp(z)*cos(x*y**2*exp(z)*log(t))/t
    assert diff(sexpr, Array([x, y, z])) == Array([bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr, bexpr*cexpr])
    assert diff(a, Array([x, y, z])) == Array([[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr], [bexpr*cexpr]])

    assert diff(sexpr, Array([[x, y], [z, t]])) == Array([[bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr], [bexpr*cexpr, bexpr/log(t)/t*cexpr]])
    assert diff(a, Array([[x, y], [z, t]])) == Array([[[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr]], [[bexpr*cexpr], [bexpr/log(t)/t*cexpr]]])
    assert diff(Array([[x, y], [z, t]]), Array([x, y])) == Array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]])
    assert diff(Array([[x, y], [z, t]]), Array([[x, y], [z, t]])) == Array([[[[1, 0], [0, 0]], [[0, 1], [0, 0]]],
                                                                         [[[0, 0], [1, 0]], [[0, 0], [0, 1]]]])

def test_trigsimp1():
    x, y = symbols('x,y')

    assert trigsimp(1 - sin(x)**2) == cos(x)**2
    assert trigsimp(1 - cos(x)**2) == sin(x)**2
    assert trigsimp(sin(x)**2 + cos(x)**2) == 1
    assert trigsimp(1 + tan(x)**2) == 1/cos(x)**2
    assert trigsimp(1/cos(x)**2 - 1) == tan(x)**2
    assert trigsimp(1/cos(x)**2 - tan(x)**2) == 1
    assert trigsimp(1 + cot(x)**2) == 1/sin(x)**2
    assert trigsimp(1/sin(x)**2 - 1) == cot(x)**2
    assert trigsimp(1/sin(x)**2 - cot(x)**2) == 1

    assert trigsimp(5*cos(x)**2 + 5*sin(x)**2) == 5
    assert trigsimp(5*cos(x/2)**2 + 2*sin(x/2)**2) in \
                [2 + 3*cos(x/2)**2, 5 - 3*sin(x/2)**2]

    assert trigsimp(sin(x)/cos(x)) == tan(x)
    assert trigsimp(2*tan(x)*cos(x)) == 2*sin(x)
    assert trigsimp(cot(x)**3*sin(x)**3) == cos(x)**3
    assert trigsimp(y*tan(x)**2/sin(x)**2) == y/cos(x)**2
    assert trigsimp(cot(x)/cos(x)) == 1/sin(x)

    assert trigsimp(cos(0.12345)**2 + sin(0.12345)**2) == 1
    e = 2*sin(x)**2 + 2*cos(x)**2
    assert trigsimp(log(e), deep=True) == log(2)