def test_airyai():
    z = Symbol('z', real=False)
    t = Symbol('t', negative=True)
    p = Symbol('p', positive=True)

    assert isinstance(airyai(z), airyai)

    assert airyai(0) == 3**(S(1)/3)/(3*gamma(S(2)/3))
    assert airyai(oo) == 0
    assert airyai(-oo) == 0

    assert diff(airyai(z), z) == airyaiprime(z)

    assert series(airyai(z), z, 0, 3) == (
        3**(S(5)/6)*gamma(S(1)/3)/(6*pi) - 3**(S(1)/6)*z*gamma(S(2)/3)/(2*pi) + O(z**3))

    assert airyai(z).rewrite(hyper) == (
        -3**(S(2)/3)*z*hyper((), (S(4)/3,), z**S(3)/9)/(3*gamma(S(1)/3)) +
         3**(S(1)/3)*hyper((), (S(2)/3,), z**S(3)/9)/(3*gamma(S(2)/3)))

    assert isinstance(airyai(z).rewrite(besselj), airyai)
    assert airyai(t).rewrite(besselj) == (
        sqrt(-t)*(besselj(-S(1)/3, 2*(-t)**(S(3)/2)/3) +
                  besselj(S(1)/3, 2*(-t)**(S(3)/2)/3))/3)
    assert airyai(z).rewrite(besseli) == (
        -z*besseli(S(1)/3, 2*z**(S(3)/2)/3)/(3*(z**(S(3)/2))**(S(1)/3)) +
         (z**(S(3)/2))**(S(1)/3)*besseli(-S(1)/3, 2*z**(S(3)/2)/3)/3)
    assert airyai(p).rewrite(besseli) == (
        sqrt(p)*(besseli(-S(1)/3, 2*p**(S(3)/2)/3) -
                 besseli(S(1)/3, 2*p**(S(3)/2)/3))/3)

    assert expand_func(airyai(2*(3*z**5)**(S(1)/3))) == (
        -sqrt(3)*(-1 + (z**5)**(S(1)/3)/z**(S(5)/3))*airybi(2*3**(S(1)/3)*z**(S(5)/3))/6 +
         (1 + (z**5)**(S(1)/3)/z**(S(5)/3))*airyai(2*3**(S(1)/3)*z**(S(5)/3))/2)

def test_from_sympy():
    x = symbols("x")
    R, Dx = DifferentialOperators(QQ.old_poly_ring(x), "Dx")
    p = from_sympy((sin(x) / x) ** 2)
    q = HolonomicFunction(
        8 * x + (4 * x ** 2 + 6) * Dx + 6 * x * Dx ** 2 + x ** 2 * Dx ** 3,
        x,
        1,
        [sin(1) ** 2, -2 * sin(1) ** 2 + 2 * sin(1) * cos(1), -8 * sin(1) * cos(1) + 2 * cos(1) ** 2 + 4 * sin(1) ** 2],
    )
    assert p == q
    p = from_sympy(1 / (1 + x ** 2) ** 2)
    q = HolonomicFunction(4 * x + (x ** 2 + 1) * Dx, x, 0, 1)
    assert p == q
    p = from_sympy(exp(x) * sin(x) + x * log(1 + x))
    q = HolonomicFunction(
        (4 * x ** 3 + 20 * x ** 2 + 40 * x + 36)
        + (-4 * x ** 4 - 20 * x ** 3 - 40 * x ** 2 - 36 * x) * Dx
        + (4 * x ** 5 + 12 * x ** 4 + 14 * x ** 3 + 16 * x ** 2 + 20 * x - 8) * Dx ** 2
        + (-4 * x ** 5 - 10 * x ** 4 - 4 * x ** 3 + 4 * x ** 2 - 2 * x + 8) * Dx ** 3
        + (2 * x ** 5 + 4 * x ** 4 - 2 * x ** 3 - 7 * x ** 2 + 2 * x + 5) * Dx ** 4,
        x,
        0,
        [0, 1, 4, -1],
    )
    assert p == q
    p = from_sympy(x * exp(x) + cos(x) + 1)
    q = HolonomicFunction(
        (-x - 3) * Dx + (x + 2) * Dx ** 2 + (-x - 3) * Dx ** 3 + (x + 2) * Dx ** 4, x, 0, [2, 1, 1, 3]
    )
    assert p == q
    assert (x * exp(x) + cos(x) + 1).series(n=10) == p.series(n=10)
    p = from_sympy(log(1 + x) ** 2 + 1)
    q = HolonomicFunction(Dx + (3 * x + 3) * Dx ** 2 + (x ** 2 + 2 * x + 1) * Dx ** 3, x, 0, [1, 0, 2])
    assert p == q
    p = from_sympy(erf(x) ** 2 + x)
    q = HolonomicFunction(
        (32 * x ** 4 - 8 * x ** 2 + 8) * Dx ** 2 + (24 * x ** 3 - 2 * x) * Dx ** 3 + (4 * x ** 2 + 1) * Dx ** 4,
        x,
        0,
        [0, 1, 8 / pi, 0],
    )
    assert p == q
    p = from_sympy(cosh(x) * x)
    q = HolonomicFunction((-x ** 2 + 2) - 2 * x * Dx + x ** 2 * Dx ** 2, x, 0, [0, 1])
    assert p == q
    p = from_sympy(besselj(2, x))
    q = HolonomicFunction((x ** 2 - 4) + x * Dx + x ** 2 * Dx ** 2, x, 0, [0, 0])
    assert p == q
    p = from_sympy(besselj(0, x) + exp(x))
    q = HolonomicFunction(
        (-2 * x ** 2 - x + 1)
        + (2 * x ** 2 - x - 3) * Dx
        + (-2 * x ** 2 + x + 2) * Dx ** 2
        + (2 * x ** 2 + x) * Dx ** 3,
        x,
        0,
        [2, 1, 1 / 2],
    )
    assert p == q

def test_airybiprime():
    z = Symbol('z', real=False)
    t = Symbol('t', negative=True)
    p = Symbol('p', positive=True)

    assert isinstance(airybiprime(z), airybiprime)

    assert airybiprime(0) == 3**(S(1)/6)/gamma(S(1)/3)
    assert airybiprime(oo) == oo
    assert airybiprime(-oo) == 0

    assert diff(airybiprime(z), z) == z*airybi(z)

    assert series(airybiprime(z), z, 0, 3) == (
        3**(S(1)/6)/gamma(S(1)/3) + 3**(S(5)/6)*z**2/(6*gamma(S(2)/3)) + O(z**3))

    assert airybiprime(z).rewrite(hyper) == (
        3**(S(5)/6)*z**2*hyper((), (S(5)/3,), z**S(3)/9)/(6*gamma(S(2)/3)) +
        3**(S(1)/6)*hyper((), (S(1)/3,), z**S(3)/9)/gamma(S(1)/3))

    assert isinstance(airybiprime(z).rewrite(besselj), airybiprime)
    assert airyai(t).rewrite(besselj) == (
        sqrt(-t)*(besselj(-S(1)/3, 2*(-t)**(S(3)/2)/3) +
                  besselj(S(1)/3, 2*(-t)**(S(3)/2)/3))/3)
    assert airybiprime(z).rewrite(besseli) == (
        sqrt(3)*(z**2*besseli(S(2)/3, 2*z**(S(3)/2)/3)/(z**(S(3)/2))**(S(2)/3) +
                 (z**(S(3)/2))**(S(2)/3)*besseli(-S(2)/3, 2*z**(S(3)/2)/3))/3)
    assert airybiprime(p).rewrite(besseli) == (
        sqrt(3)*p*(besseli(-S(2)/3, 2*p**(S(3)/2)/3) + besseli(S(2)/3, 2*p**(S(3)/2)/3))/3)

    assert expand_func(airybiprime(2*(3*z**5)**(S(1)/3))) == (
        sqrt(3)*(z**(S(5)/3)/(z**5)**(S(1)/3) - 1)*airyaiprime(2*3**(S(1)/3)*z**(S(5)/3))/2 +
        (z**(S(5)/3)/(z**5)**(S(1)/3) + 1)*airybiprime(2*3**(S(1)/3)*z**(S(5)/3))/2)

def test_to_hyper():
    x = symbols('x')
    R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
    p = HolonomicFunction(Dx - 2, x, 0, [3]).to_hyper()
    q = 3 * hyper([], [], 2*x)
    assert p == q
    p = hyperexpand(HolonomicFunction((1 + x) * Dx - 3, x, 0, [2]).to_hyper()).expand()
    q = 2*x**3 + 6*x**2 + 6*x + 2
    assert p == q
    p = HolonomicFunction((1 + x)*Dx**2 + Dx, x, 0, [0, 1]).to_hyper()
    q = -x**2*hyper((2, 2, 1), (2, 3), -x)/2 + x
    assert p == q
    p = HolonomicFunction(2*x*Dx + Dx**2, x, 0, [0, 2/sqrt(pi)]).to_hyper()
    q = 2*x*hyper((1/2,), (3/2,), -x**2)/sqrt(pi)
    assert p == q
    p = hyperexpand(HolonomicFunction(2*x*Dx + Dx**2, x, 0, [1, -2/sqrt(pi)]).to_hyper())
    q = erfc(x)
    assert p.rewrite(erfc) == q
    p =  hyperexpand(HolonomicFunction((x**2 - 1) + x*Dx + x**2*Dx**2,
        x, 0, [0, S(1)/2]).to_hyper())
    q = besselj(1, x)
    assert p == q
    p = hyperexpand(HolonomicFunction(x*Dx**2 + Dx + x, x, 0, [1, 0]).to_hyper())
    q = besselj(0, x)
    assert p == q

def test_diff():
    assert besselj(n, z).diff(z) == besselj(n - 1, z)/2 - besselj(n + 1, z)/2
    assert bessely(n, z).diff(z) == bessely(n - 1, z)/2 - bessely(n + 1, z)/2
    assert besseli(n, z).diff(z) == besseli(n - 1, z)/2 + besseli(n + 1, z)/2
    assert besselk(n, z).diff(z) == -besselk(n - 1, z)/2 - besselk(n + 1, z)/2
    assert hankel1(n, z).diff(z) == hankel1(n - 1, z)/2 - hankel1(n + 1, z)/2
    assert hankel2(n, z).diff(z) == hankel2(n - 1, z)/2 - hankel2(n + 1, z)/2

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_expr_to_holonomic():
    x = symbols('x')
    R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
    p = expr_to_holonomic((sin(x)/x)**2)
    q = HolonomicFunction(8*x + (4*x**2 + 6)*Dx + 6*x*Dx**2 + x**2*Dx**3, x, 0, \
        [1, 0, -2/3])
    assert p == q
    p = expr_to_holonomic(1/(1+x**2)**2)
    q = HolonomicFunction(4*x + (x**2 + 1)*Dx, x, 0, 1)
    assert p == q
    p = expr_to_holonomic(exp(x)*sin(x)+x*log(1+x))
    q = HolonomicFunction((2*x**3 + 10*x**2 + 20*x + 18) + (-2*x**4 - 10*x**3 - 20*x**2 \
        - 18*x)*Dx + (2*x**5 + 6*x**4 + 7*x**3 + 8*x**2 + 10*x - 4)*Dx**2 + \
        (-2*x**5 - 5*x**4 - 2*x**3 + 2*x**2 - x + 4)*Dx**3 + (x**5 + 2*x**4 - x**3 - \
        7*x**2/2 + x + 5/2)*Dx**4, x, 0, [0, 1, 4, -1])
    assert p == q
    p = expr_to_holonomic(x*exp(x)+cos(x)+1)
    q = HolonomicFunction((-x - 3)*Dx + (x + 2)*Dx**2 + (-x - 3)*Dx**3 + (x + 2)*Dx**4, x, \
        0, [2, 1, 1, 3])
    assert p == q
    assert (x*exp(x)+cos(x)+1).series(n=10) == p.series(n=10)
    p = expr_to_holonomic(log(1 + x)**2 + 1)
    q = HolonomicFunction(Dx + (3*x + 3)*Dx**2 + (x**2 + 2*x + 1)*Dx**3, x, 0, [1, 0, 2])
    assert p == q
    p = expr_to_holonomic(erf(x)**2 + x)
    q = HolonomicFunction((8*x**4 - 2*x**2 + 2)*Dx**2 + (6*x**3 - x/2)*Dx**3 + \
        (x**2+ 1/4)*Dx**4, x, 0, [0, 1, 8/pi, 0])
    assert p == q
    p = expr_to_holonomic(cosh(x)*x)
    q = HolonomicFunction((-x**2 + 2) -2*x*Dx + x**2*Dx**2, x, 0, [0, 1])
    assert p == q
    p = expr_to_holonomic(besselj(2, x))
    q = HolonomicFunction((x**2 - 4) + x*Dx + x**2*Dx**2, x, 0, [0, 0])
    assert p == q
    p = expr_to_holonomic(besselj(0, x) + exp(x))
    q = HolonomicFunction((-x**2 - x/2 + 1/2) + (x**2 - x/2 - 3/2)*Dx + (-x**2 + x/2 + 1)*Dx**2 +\
        (x**2 + x/2)*Dx**3, x, 0, [2, 1, 1/2])
    assert p == q
    p = expr_to_holonomic(sin(x)**2/x)
    q = HolonomicFunction(4 + 4*x*Dx + 3*Dx**2 + x*Dx**3, x, 0, [0, 1, 0])
    assert p == q
    p = expr_to_holonomic(sin(x)**2/x, x0=2)
    q = HolonomicFunction((4) + (4*x)*Dx + (3)*Dx**2 + (x)*Dx**3, x, 2, [sin(2)**2/2,
        sin(2)*cos(2) - sin(2)**2/4, -3*sin(2)**2/4 + cos(2)**2 - sin(2)*cos(2)])
    assert p == q
    p = expr_to_holonomic(log(x)/2 - Ci(2*x)/2 + Ci(2)/2)
    q = HolonomicFunction(4*Dx + 4*x*Dx**2 + 3*Dx**3 + x*Dx**4, x, 0, \
        [-log(2)/2 - EulerGamma/2 + Ci(2)/2, 0, 1, 0])
    assert p == q
    p = p.to_expr()
    q = log(x)/2 - Ci(2*x)/2 + Ci(2)/2
    assert p == q
    p = expr_to_holonomic(x**(S(1)/2), x0=1)
    q = HolonomicFunction(x*Dx - 1/2, x, 1, 1)
    assert p == q
    p = expr_to_holonomic(sqrt(1 + x**2))
    q = HolonomicFunction((-x) + (x**2 + 1)*Dx, x, 0, 1)
    assert p == q

def test_rewrite():
    from sympy import polar_lift, exp, I
    assert besselj(n, z).rewrite(jn) == sqrt(2*z/pi)*jn(n - S(1)/2, z)
    assert bessely(n, z).rewrite(yn) == sqrt(2*z/pi)*yn(n - S(1)/2, z)
    assert besseli(n, z).rewrite(besselj) == \
        exp(-I*n*pi/2)*besselj(n, polar_lift(I)*z)
    assert besselj(n, z).rewrite(besseli) == \
        exp(I*n*pi/2)*besseli(n, polar_lift(-I)*z)
    nu = randcplx()
    assert tn(besselj(nu, z), besselj(nu, z).rewrite(besseli), z)
    assert tn(besseli(nu, z), besseli(nu, z).rewrite(besselj), z)

def test_laplace_transform():
    LT = laplace_transform
    a, b, c, = symbols('a b c', positive=True)
    t = symbols('t')
    w = Symbol("w")
    f = Function("f")

    # Test unevaluated form
    assert laplace_transform(f(t), t, w) == LaplaceTransform(f(t), t, w)
    assert inverse_laplace_transform(f(w), w, t, plane=0) == InverseLaplaceTransform(f(w), w, t, 0)

    # test a bug
    spos = symbols('s', positive=True)
    assert LT(exp(t), t, spos)[:2] == (1/(spos - 1), True)

    # basic tests from wikipedia

    assert LT((t-a)**b*exp(-c*(t-a))*Heaviside(t-a), t, s) \
           == ((s + c)**(-b - 1)*exp(-a*s)*gamma(b + 1), -c, True)
    assert LT(t**a, t, s) == (s**(-a - 1)*gamma(a + 1), 0, True)
    assert LT(Heaviside(t), t, s) == (1/s, 0, True)
    assert LT(Heaviside(t - a), t, s) == (exp(-a*s)/s, 0, True)
    assert LT(1 - exp(-a*t), t, s) == (a/(s*(a + s)), 0, True)

    assert LT((exp(2*t)-1)*exp(-b - t)*Heaviside(t)/2, t, s, noconds=True) \
           == exp(-b)/(s**2 - 1)

    assert LT(exp(t), t, s)[:2] == (1/(s-1), 1)
    assert LT(exp(2*t), t, s)[:2] == (1/(s-2), 2)
    assert LT(exp(a*t), t, s)[:2] == (1/(s-a), a)

    assert LT(log(t/a), t, s) == ((log(a) + log(s) + EulerGamma)/(-s), 0, True)

    assert LT(erf(t), t, s) == ((-erf(s/2) + 1)*exp(s**2/4)/s, 0, True)

    assert LT(sin(a*t), t, s) == (a/(a**2 + s**2), 0, True)
    assert LT(cos(a*t), t, s) == (s/(a**2 + s**2), 0, True)
    # TODO would be nice to have these come out better
    assert LT(exp(-a*t)*sin(b*t), t, s) == (1/b/(1 + (a + s)**2/b**2), -a, True)
    assert LT(exp(-a*t)*cos(b*t), t, s) == \
           (1/(s + a)/(1 + b**2/(a + s)**2), -a, True)
    # TODO sinh, cosh have delicate cancellation

    assert LT(besselj(0, t), t, s) == (1/sqrt(1 + s**2), 0, True)
    assert LT(besselj(1, t), t, s) == (1 - 1/sqrt(1 + 1/s**2), 0, True)
    # TODO general order works, but is a *mess*
    # TODO besseli also works, but is an even greater mess

    # test a bug in conditions processing
    # TODO the auxiliary condition should be recognised/simplified
    assert LT(exp(t)*cos(t), t, s)[:-1] in [
        ((s - 1)/(s**2 - 2*s + 2), -oo),
        ((s - 1)/((s - 1)**2 + 1), -oo),
        ]

def test_besselsimp():
    from sympy import besselj, besseli, besselk, bessely, jn, yn, exp_polar, cosh
    assert besselsimp(exp(-I*pi*y/2)*besseli(y, z*exp_polar(I*pi/2))) == \
           besselj(y, z)
    assert besselsimp(exp(-I*pi*a/2)*besseli(a, 2*sqrt(x)*exp_polar(I*pi/2))) == \
           besselj(a, 2*sqrt(x))
    assert besselsimp(sqrt(2)*sqrt(pi)*x**(S(1)/4)*exp(I*pi/4)*exp(-I*pi*a/2) * \
                      besseli(-S(1)/2, sqrt(x)*exp_polar(I*pi/2)) * \
                      besseli(a, sqrt(x)*exp_polar(I*pi/2))/2) == \
           besselj(a, sqrt(x)) * cos(sqrt(x))
    assert besselsimp(besseli(S(-1)/2, z)) == sqrt(2)*cosh(z)/(sqrt(pi)*sqrt(z))
    assert besselsimp(besseli(a, z*exp_polar(-I*pi/2))) == exp(-I*pi*a/2)*besselj(a, z)

def test_branching():
    from sympy import exp_polar, polar_lift, Symbol, I, exp
    assert besselj(polar_lift(k), x) == besselj(k, x)
    assert besseli(polar_lift(k), x) == besseli(k, x)

    n = Symbol('n', integer=True)
    assert besselj(n, exp_polar(2*pi*I)*x) == besselj(n, x)
    assert besselj(n, polar_lift(x)) == besselj(n, x)
    assert besseli(n, exp_polar(2*pi*I)*x) == besseli(n, x)
    assert besseli(n, polar_lift(x)) == besseli(n, x)

    def tn(func, s):
        from random import uniform
        c = uniform(1, 5)
        expr = func(s, c*exp_polar(I*pi)) - func(s, c*exp_polar(-I*pi))
        eps = 1e-15
        expr2 = func(s + eps, -c + eps*I) - func(s + eps, -c - eps*I)
        return abs(expr.n() - expr2.n()).n() < 1e-10

    nu = Symbol('nu')
    assert besselj(nu, exp_polar(2*pi*I)*x) == exp(2*pi*I*nu)*besselj(nu, x)
    assert besseli(nu, exp_polar(2*pi*I)*x) == exp(2*pi*I*nu)*besseli(nu, x)
    assert tn(besselj, 2)
    assert tn(besselj, pi)
    assert tn(besselj, I)
    assert tn(besseli, 2)
    assert tn(besseli, pi)
    assert tn(besseli, I)

def test_rewrite():
    from sympy import polar_lift, exp, I

    assert besselj(n, z).rewrite(jn) == sqrt(2*z/pi)*jn(n - S(1)/2, z)
    assert bessely(n, z).rewrite(yn) == sqrt(2*z/pi)*yn(n - S(1)/2, z)
    assert besseli(n, z).rewrite(besselj) == \
        exp(-I*n*pi/2)*besselj(n, polar_lift(I)*z)
    assert besselj(n, z).rewrite(besseli) == \
        exp(I*n*pi/2)*besseli(n, polar_lift(-I)*z)

    nu = randcplx()

    assert tn(besselj(nu, z), besselj(nu, z).rewrite(besseli), z)
    assert tn(besselj(nu, z), besselj(nu, z).rewrite(bessely), z)

    assert tn(besseli(nu, z), besseli(nu, z).rewrite(besselj), z)
    assert tn(besseli(nu, z), besseli(nu, z).rewrite(bessely), z)

    assert tn(bessely(nu, z), bessely(nu, z).rewrite(besselj), z)
    assert tn(bessely(nu, z), bessely(nu, z).rewrite(besseli), z)

    assert tn(besselk(nu, z), besselk(nu, z).rewrite(besselj), z)
    assert tn(besselk(nu, z), besselk(nu, z).rewrite(besseli), z)
    assert tn(besselk(nu, z), besselk(nu, z).rewrite(bessely), z)

    # check that a rewrite was triggered, when the order is set to a generic
    # symbol 'nu'
    assert yn(nu, z) != yn(nu, z).rewrite(jn)
    assert hn1(nu, z) != hn1(nu, z).rewrite(jn)
    assert hn2(nu, z) != hn2(nu, z).rewrite(jn)
    assert jn(nu, z) != jn(nu, z).rewrite(yn)
    assert hn1(nu, z) != hn1(nu, z).rewrite(yn)
    assert hn2(nu, z) != hn2(nu, z).rewrite(yn)

    # rewriting spherical bessel functions (SBFs) w.r.t. besselj, bessely is
    # not allowed if a generic symbol 'nu' is used as the order of the SBFs
    # to avoid inconsistencies (the order of bessel[jy] is allowed to be
    # complex-valued, whereas SBFs are defined only for integer orders)
    order = nu
    for f in (besselj, bessely):
        assert hn1(order, z) == hn1(order, z).rewrite(f)
        assert hn2(order, z) == hn2(order, z).rewrite(f)

    assert jn(order, z).rewrite(besselj) == sqrt(2)*sqrt(pi)*sqrt(1/z)*besselj(order + S(1)/2, z)/2
    assert jn(order, z).rewrite(bessely) == (-1)**nu*sqrt(2)*sqrt(pi)*sqrt(1/z)*bessely(-order - S(1)/2, z)/2

    # for integral orders rewriting SBFs w.r.t bessel[jy] is allowed
    N = Symbol('n', integer=True)
    ri = randint(-11, 10)
    for order in (ri, N):
        for f in (besselj, bessely):
            assert yn(order, z) != yn(order, z).rewrite(f)
            assert jn(order, z) != jn(order, z).rewrite(f)
            assert hn1(order, z) != hn1(order, z).rewrite(f)
            assert hn2(order, z) != hn2(order, z).rewrite(f)

    for func, refunc in product((yn, jn, hn1, hn2),
                                (jn, yn, besselj, bessely)):
        assert tn(func(ri, z), func(ri, z).rewrite(refunc), z)

def test_hyper():
    for x in sorted(exparg):
        test("erf", x, N(sp.erf(x)))
    for x in sorted(exparg):
        test("erfc", x, N(sp.erfc(x)))

    gamarg = FiniteSet(*(x+S(1)/12 for x in exparg))
    betarg = ProductSet(gamarg, gamarg)
    for x in sorted(gamarg):
        test("lgamma", x, N(sp.log(abs(sp.gamma(x)))))
    for x in sorted(gamarg):
        test("gamma", x, N(sp.gamma(x)))
    for x, y in sorted(betarg, key=lambda (x, y): (y, x)):
        test("beta", x, y, N(sp.beta(x, y)))

    pgamarg = FiniteSet(S(1)/12, S(1)/3, S(3)/2, 5)
    pgamargp = ProductSet(gamarg & Interval(0, oo, True), pgamarg)
    for a, x in sorted(pgamargp):
        test("pgamma", a, x, N(sp.lowergamma(a, x)))
    for a, x in sorted(pgamargp):
        test("pgammac", a, x, N(sp.uppergamma(a, x)))
    for a, x in sorted(pgamargp):
        test("pgammar", a, x, N(sp.lowergamma(a, x)/sp.gamma(a)))
    for a, x in sorted(pgamargp):
        test("pgammarc", a, x, N(sp.uppergamma(a, x)/sp.gamma(a)))
    for a, x in sorted(pgamargp):
        test("ipgammarc", a, N(sp.uppergamma(a, x)/sp.gamma(a)), x)

    pbetargp = [(a, b, x) for a, b, x in ProductSet(betarg, pgamarg)
                if a > 0 and b > 0 and x < 1]
    pbetargp.sort(key=lambda (a, b, x): (b, a, x))
    for a, b, x in pbetargp:
        test("pbeta", a, b, x, mp.betainc(mpf(a), mpf(b), x2=mpf(x)))
    for a, b, x in pbetargp:
        test("pbetar", a, b, x, mp.betainc(mpf(a), mpf(b), x2=mpf(x),
                                           regularized=True))
    for a, b, x in pbetargp:
        test("ipbetar", a, b, mp.betainc(mpf(a), mpf(b), x2=mpf(x),
                                         regularized=True), x)

    for x in sorted(posarg):
        test("j0", x, N(sp.besselj(0, x)))
    for x in sorted(posarg):
        test("j1", x, N(sp.besselj(1, x)))
    for x in sorted(posarg-FiniteSet(0)):
        test("y0", x, N(sp.bessely(0, x)))
    for x in sorted(posarg-FiniteSet(0)):
        test("y1", x, N(sp.bessely(1, x)))

def test_inverse_laplace_transform():
    from sympy import sinh, cosh, besselj, besseli, simplify, factor_terms

    ILT = inverse_laplace_transform
    a, b, c, = symbols("a b c", positive=True, finite=True)
    t = symbols("t")

    def simp_hyp(expr):
        return factor_terms(expand_mul(expr)).rewrite(sin)

    # just test inverses of all of the above
    assert ILT(1 / s, s, t) == Heaviside(t)
    assert ILT(1 / s ** 2, s, t) == t * Heaviside(t)
    assert ILT(1 / s ** 5, s, t) == t ** 4 * Heaviside(t) / 24
    assert ILT(exp(-a * s) / s, s, t) == Heaviside(t - a)
    assert ILT(exp(-a * s) / (s + b), s, t) == exp(b * (a - t)) * Heaviside(-a + t)
    assert ILT(a / (s ** 2 + a ** 2), s, t) == sin(a * t) * Heaviside(t)
    assert ILT(s / (s ** 2 + a ** 2), s, t) == cos(a * t) * Heaviside(t)
    # TODO is there a way around simp_hyp?
    assert simp_hyp(ILT(a / (s ** 2 - a ** 2), s, t)) == sinh(a * t) * Heaviside(t)
    assert simp_hyp(ILT(s / (s ** 2 - a ** 2), s, t)) == cosh(a * t) * Heaviside(t)
    assert ILT(a / ((s + b) ** 2 + a ** 2), s, t) == exp(-b * t) * sin(a * t) * Heaviside(t)
    assert ILT((s + b) / ((s + b) ** 2 + a ** 2), s, t) == exp(-b * t) * cos(a * t) * Heaviside(t)
    # TODO sinh/cosh shifted come out a mess. also delayed trig is a mess
    # TODO should this simplify further?
    assert ILT(exp(-a * s) / s ** b, s, t) == (t - a) ** (b - 1) * Heaviside(t - a) / gamma(b)

    assert ILT(exp(-a * s) / sqrt(1 + s ** 2), s, t) == Heaviside(t - a) * besselj(
        0, a - t
    )  # note: besselj(0, x) is even

    # XXX ILT turns these branch factor into trig functions ...
    assert simplify(
        ILT(a ** b * (s + sqrt(s ** 2 - a ** 2)) ** (-b) / sqrt(s ** 2 - a ** 2), s, t).rewrite(exp)
    ) == Heaviside(t) * besseli(b, a * t)
    assert ILT(a ** b * (s + sqrt(s ** 2 + a ** 2)) ** (-b) / sqrt(s ** 2 + a ** 2), s, t).rewrite(exp) == Heaviside(
        t
    ) * besselj(b, a * t)

    assert ILT(1 / (s * sqrt(s + 1)), s, t) == Heaviside(t) * erf(sqrt(t))
    # TODO can we make erf(t) work?

    assert ILT(1 / (s ** 2 * (s ** 2 + 1)), s, t) == (t - sin(t)) * Heaviside(t)

    assert ILT((s * eye(2) - Matrix([[1, 0], [0, 2]])).inv(), s, t) == Matrix(
        [[exp(t) * Heaviside(t), 0], [0, exp(2 * t) * Heaviside(t)]]
    )

def test_bessel_rand():
    assert td(besselj(randcplx(), z), z)
    assert td(bessely(randcplx(), z), z)
    assert td(besseli(randcplx(), z), z)
    assert td(besselk(randcplx(), z), z)
    assert td(hankel1(randcplx(), z), z)
    assert td(hankel2(randcplx(), z), z)
    assert td(jn(randcplx(), z), z)
    assert td(yn(randcplx(), z), z)

def test_inverse_laplace_transform():
    from sympy import (expand, sinh, cosh, besselj, besseli, exp_polar,
                       unpolarify, simplify)
    ILT = inverse_laplace_transform
    a, b, c, = symbols('a b c', positive=True)
    t = symbols('t')

    def simp_hyp(expr):
        return expand(expand(expr).rewrite(sin))

    # just test inverses of all of the above
    assert ILT(1/s, s, t) == Heaviside(t)
    assert ILT(1/s**2, s, t) == t*Heaviside(t)
    assert ILT(1/s**5, s, t) == t**4*Heaviside(t)/24
    assert ILT(exp(-a*s)/s, s, t) == Heaviside(t - a)
    assert ILT(exp(-a*s)/(s + b), s, t) == exp(b*(a - t))*Heaviside(-a + t)
    assert ILT(a/(s**2 + a**2), s, t) == sin(a*t)*Heaviside(t)
    assert ILT(s/(s**2 + a**2), s, t) == cos(a*t)*Heaviside(t)
    # TODO is there a way around simp_hyp?
    assert simp_hyp(ILT(a/(s**2 - a**2), s, t)) == sinh(a*t)*Heaviside(t)
    assert simp_hyp(ILT(s/(s**2 - a**2), s, t)) == cosh(a*t)*Heaviside(t)
    assert ILT(a/((s + b)**2 + a**2), s, t) == exp(-b*t)*sin(a*t)*Heaviside(t)
    assert ILT(
        (s + b)/((s + b)**2 + a**2), s, t) == exp(-b*t)*cos(a*t)*Heaviside(t)
    # TODO sinh/cosh shifted come out a mess. also delayed trig is a mess
    # TODO should this simplify further?
    assert ILT(exp(-a*s)/s**b, s, t) == \
        (t - a)**(b - 1)*Heaviside(t - a)/gamma(b)

    assert ILT(exp(-a*s)/sqrt(1 + s**2), s, t) == \
        Heaviside(t - a)*besselj(0, a - t)  # note: besselj(0, x) is even

    # XXX ILT turns these branch factor into trig functions ...
    assert simplify(ILT(a**b*(s + sqrt(s**2 - a**2))**(-b)/sqrt(s**2 - a**2),
                    s, t).rewrite(exp)) == \
        Heaviside(t)*besseli(b, a*t)
    assert ILT(a**b*(s + sqrt(s**2 + a**2))**(-b)/sqrt(s**2 + a**2),
               s, t).rewrite(exp) == \
        Heaviside(t)*besselj(b, a*t)

    assert ILT(1/(s*sqrt(s + 1)), s, t) == Heaviside(t)*erf(sqrt(t))

def test_inversion():
    from sympy import piecewise_fold, besselj, sqrt, I, sin, cos, Heaviside
    def inv(f): return piecewise_fold(meijerint_inversion(f, s, t))
    assert inv(1/(s**2 + 1)) == sin(t)*Heaviside(t)
    assert inv(s/(s**2 + 1)) == cos(t)*Heaviside(t)
    assert inv(exp(-s)/s) == Heaviside(t - 1)
    assert inv(1/sqrt(1 + s**2)) == besselj(0, t)*Heaviside(t)

    # Test some antcedents checking.
    assert meijerint_inversion(sqrt(s)/sqrt(1 + s**2), s, t) is None
    assert inv(exp(s**2)) is None
    assert meijerint_inversion(exp(-s**2), s, t) is None

def test_laplace_transform():
    LT = laplace_transform
    a, b, c, = symbols('a b c', positive=True)
    t = symbols('t')

    # basic tests from wikipedia

    assert LT((t-a)**b*exp(-c*(t-a))*Heaviside(t-a), t, s) \
           == ((s + c)**(-b - 1)*exp(-a*s)*gamma(b + 1), -c, True)
    assert LT(t**a, t, s) == (s**(-a - 1)*gamma(a + 1), 0, True)
    assert LT(Heaviside(t), t, s) == (1/s, 0, True)
    assert LT(Heaviside(t - a), t, s) == (exp(-a*s)/s, 0, True)
    assert LT(1 - exp(-a*t), t, s) == (a/(s*(a + s)), 0, True)

    assert LT((exp(2*t)-1)*exp(-b - t)*Heaviside(t)/2, t, s, noconds=True) \
           == exp(-b)/(s**2 - 1)

    assert LT(exp(t), t, s)[0:2] == (1/(s-1), 1)
    assert LT(exp(2*t), t, s)[0:2] == (1/(s-2), 2)
    assert LT(exp(a*t), t, s)[0:2] == (1/(s-a), a)

    lt = LT(log(t/a), t, s)
    assert lt[1:] == (0, True)
    # TODO hyperexpand is not clever enough to recognise this on its own
    assert hyperexpand(lt[0], allow_hyper=True) == (-log(a*s) - EulerGamma)/s

    assert LT(erf(t), t, s) == ((-erf(s/2) + 1)*exp(s**2/4)/s, 0, True)

    assert LT(sin(a*t), t, s) == (a/(a**2 + s**2), 0, True)
    assert LT(cos(a*t), t, s) == (s/(a**2 + s**2), 0, True)
    # TODO would be nice to have these come out better
    assert LT(exp(-a*t)*sin(b*t), t, s) == (1/b/(1 + (a + s)**2/b**2), -a, True)
    assert LT(exp(-a*t)*cos(b*t), t, s) == \
           (sqrt((a + s)**2)/b**2/(1 + (a + s)**2/b**2), -a, True)
    # TODO sinh, cosh have delicate cancellation

    # TODO conditions are a mess
    assert LT(besselj(0, t), t, s, noconds=True) == 1/sqrt(1 + s**2)
    assert LT(besselj(1, t), t, s, noconds=True) == 1 - 1/sqrt(1 + 1/s**2)

def test_bessel_eval():
    from sympy import I, Symbol
    n, m, k = Symbol('n', integer=True), Symbol('m'), Symbol('k', integer=True, zero=False)

    for f in [besselj, besseli]:
        assert f(0, 0) == S.One
        assert f(2.1, 0) == S.Zero
        assert f(-3, 0) == S.Zero
        assert f(-10.2, 0) == S.ComplexInfinity
        assert f(1 + 3*I, 0) == S.Zero
        assert f(-3 + I, 0) == S.ComplexInfinity
        assert f(-2*I, 0) == S.NaN
        assert f(n, 0) != S.One and f(n, 0) != S.Zero
        assert f(m, 0) != S.One and f(m, 0) != S.Zero
        assert f(k, 0) == S.Zero

    assert bessely(0, 0) == S.NegativeInfinity
    assert besselk(0, 0) == S.Infinity
    for f in [bessely, besselk]:
        assert f(1 + I, 0) == S.ComplexInfinity
        assert f(I, 0) == S.NaN

    for f in [besselj, bessely]:
        assert f(m, S.Infinity) == S.Zero
        assert f(m, S.NegativeInfinity) == S.Zero

    for f in [besseli, besselk]:
        assert f(m, I*S.Infinity) == S.Zero
        assert f(m, I*S.NegativeInfinity) == S.Zero

    for f in [besseli, besselk]:
        assert f(-4, z) == f(4, z)
        assert f(-3, z) == f(3, z)
        assert f(-n, z) == f(n, z)
        assert f(-m, z) != f(m, z)

    for f in [besselj, bessely]:
        assert f(-4, z) == f(4, z)
        assert f(-3, z) == -f(3, z)
        assert f(-n, z) == (-1)**n*f(n, z)
        assert f(-m, z) != (-1)**m*f(m, z)

    for f in [besselj, besseli]:
        assert f(m, -z) == (-z)**m*z**(-m)*f(m, z)

    assert besseli(2, -z) == besseli(2, z)
    assert besseli(3, -z) == -besseli(3, z)

    assert besselj(0, -z) == besselj(0, z)
    assert besselj(1, -z) == -besselj(1, z)

    assert besseli(0, I*z) == besselj(0, z)
    assert besseli(1, I*z) == I*besselj(1, z)
    assert besselj(3, I*z) == -I*besseli(3, z)

def test_laplace_transform():
    LT = laplace_transform
    a, b, c, = symbols('a b c', positive=True)
    t = symbols('t')

    # test a bug
    spos = symbols('s', positive=True)
    assert LT(exp(t), t, spos)[:2] == (1/(spos - 1), True)

    # basic tests from wikipedia

    assert LT((t-a)**b*exp(-c*(t-a))*Heaviside(t-a), t, s) \
           == ((s + c)**(-b - 1)*exp(-a*s)*gamma(b + 1), -c, True)
    assert LT(t**a, t, s) == (s**(-a - 1)*gamma(a + 1), 0, True)
    assert LT(Heaviside(t), t, s) == (1/s, 0, True)
    assert LT(Heaviside(t - a), t, s) == (exp(-a*s)/s, 0, True)
    assert LT(1 - exp(-a*t), t, s) == (a/(s*(a + s)), 0, True)

    assert LT((exp(2*t)-1)*exp(-b - t)*Heaviside(t)/2, t, s, noconds=True) \
           == exp(-b)/(s**2 - 1)

    assert LT(exp(t), t, s)[0:2] == (1/(s-1), 1)
    assert LT(exp(2*t), t, s)[0:2] == (1/(s-2), 2)
    assert LT(exp(a*t), t, s)[0:2] == (1/(s-a), a)

    assert LT(log(t/a), t, s) == ((log(a) + log(s) + EulerGamma)/(-s), 0, True)

    assert LT(erf(t), t, s) == ((-erf(s/2) + 1)*exp(s**2/4)/s, 0, True)

    assert LT(sin(a*t), t, s) == (a/(a**2 + s**2), 0, True)
    assert LT(cos(a*t), t, s) == (s/(a**2 + s**2), 0, True)
    # TODO would be nice to have these come out better
    assert LT(exp(-a*t)*sin(b*t), t, s) == (1/b/(1 + (a + s)**2/b**2), -a, True)
    assert LT(exp(-a*t)*cos(b*t), t, s) == \
           (1/(s + a)/(1 + b**2/(a + s)**2), -a, True)
    # TODO sinh, cosh have delicate cancellation

    assert LT(besselj(0, t), t, s) == (1/sqrt(1 + s**2), 0, True)
    assert LT(besselj(1, t), t, s) == (1 - 1/sqrt(1 + 1/s**2), 0, True)