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_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_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_hyperexpand_bases():
    assert (
        hyperexpand(hyper([2], [a], z))
        == a + z ** (-a + 1) * (-a ** 2 + 3 * a + z * (a - 1) - 2) * exp(z) * lowergamma(a - 1, z) - 1
    )
    # TODO [a+1, a-S.Half], [2*a]
    assert hyperexpand(hyper([1, 2], [3], z)) == -2 / z - 2 * log(exp_polar(-I * pi) * z + 1) / z ** 2
    assert hyperexpand(hyper([S.Half, 2], [S(3) / 2], z)) == -1 / (2 * z - 2) + log((sqrt(z) + 1) / (-sqrt(z) + 1)) / (
        4 * sqrt(z)
    )
    assert hyperexpand(hyper([S(1) / 2, S(1) / 2], [S(5) / 2], z)) == (-3 * z + 3) / 4 / (z * sqrt(-z + 1)) + (
        6 * z - 3
    ) * asin(sqrt(z)) / (4 * z ** (S(3) / 2))
    assert hyperexpand(hyper([1, 2], [S(3) / 2], z)) == -1 / (2 * z - 2) - asin(sqrt(z)) / (
        sqrt(z) * (2 * z - 2) * sqrt(-z + 1)
    )
    assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == sqrt(z) * (6 * z / 7 - S(6) / 5) * atanh(
        sqrt(z)
    ) + (-30 * z ** 2 + 32 * z - 6) / 35 / z - 6 * log(-z + 1) / (35 * z ** 2)
    assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == -4 * log(sqrt(-z + 1) / 2 + S(1) / 2) / z
    # TODO hyperexpand(hyper([a], [2*a + 1], z))
    # TODO [S.Half, a], [S(3)/2, a+1]
    assert hyperexpand(hyper([2], [b, 1], z)) == z ** (-b / 2 + S(1) / 2) * besseli(b - 1, 2 * sqrt(z)) * gamma(
        b
    ) + z ** (-b / 2 + 1) * besseli(b, 2 * sqrt(z)) * gamma(b)

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_bessel_eval():
    from sympy import I
    assert besselj(-4, z) == besselj(4, z)
    assert besselj(-3, z) == -besselj(3, z)

    assert bessely(-2, z) == bessely(2, z)
    assert bessely(-1, z) == -bessely(1, 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_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_bessel():
    from sympy import (besselj, Heaviside, besseli, polar_lift, exp_polar,
                       powdenest)
    assert simplify(integrate(besselj(a, z)*besselj(b, z)/z, (z, 0, oo),
                     meijerg=True, conds='none')) == \
           2*sin(pi*a/2 - pi*b/2)/(pi*(a - b)*(a + b))
    assert simplify(integrate(besselj(a, z)*besselj(a, z)/z, (z, 0, oo),
                     meijerg=True, conds='none')) == 1/(2*a)

    # TODO more orthogonality integrals

    # TODO there is some improvement possible here:
    #  - the result can be simplified to besselj(y, z))
    assert powdenest(simplify(integrate(sin(z*x)*(x**2-1)**(-(y+S(1)/2)),
                              (x, 1, oo), meijerg=True, conds='none')
                              *2/((z/2)**y*sqrt(pi)*gamma(S(1)/2-y))),
                     polar=True) == \
           exp(-I*pi*y/2)*besseli(y, z*exp_polar(I*pi/2))

    # Werner Rosenheinrich
    # SOME INDEFINITE INTEGRALS OF BESSEL FUNCTIONS

    assert integrate(x*besselj(0, x), x, meijerg=True) == x*besselj(1, x)
    assert integrate(x*besseli(0, x), x, meijerg=True) == x*besseli(1, x)
    # TODO can do higher powers, but come out as high order ... should they be
    #      reduced to order 0, 1?
    assert integrate(besselj(1, x), x, meijerg=True) == -besselj(0, x)
    assert integrate(besselj(1, x)**2/x, x, meijerg=True) == \
           -(besselj(0, x)**2 + besselj(1, x)**2)/2
    # TODO more besseli when tables are extended or recursive mellin works
    assert integrate(besselj(0, x)**2/x**2, x, meijerg=True) == \
           -2*x*besselj(0, x)**2 - 2*x*besselj(1, x)**2 \
           + 2*besselj(0, x)*besselj(1, x) - besselj(0, x)**2/x
    assert integrate(besselj(0, x)*besselj(1, x), x, meijerg=True) == \
           -besselj(0, x)**2/2
    assert integrate(x**2*besselj(0, x)*besselj(1, x), x, meijerg=True) == \
           x**2*besselj(1, x)**2/2
    assert integrate(besselj(0, x)*besselj(1, x)/x, x, meijerg=True) == \
           (x*besselj(0, x)**2 + x*besselj(1, x)**2 - \
            besselj(0, x)*besselj(1, x))
    # TODO how does besselj(0, a*x)*besselj(0, b*x) work?
    # TODO how does besselj(0, x)**2*besselj(1, x)**2 work?
    # TODO sin(x)*besselj(0, x) etc come out a mess
    # TODO can x*log(x)*besselj(0, x) be done?
    # TODO how does besselj(1, x)*besselj(0, x+a) work?
    # TODO more indefinite integrals when struve functions etc are implemented

    # test a substitution
    assert integrate(besselj(1, x**2)*x, x, meijerg=True) == \
        -besselj(0, x**2)/2

def test_chi_noncentral():
    k = Symbol("k", integer=True)
    l = Symbol("l")

    X = ChiNoncentral("x", k, l)
    assert density(X)(x) == (x**k*l*(x*l)**(-k/2)*
                          exp(-x**2/2 - l**2/2)*besseli(k/2 - 1, x*l))

def test_meijerg_eval():
    from sympy import besseli, exp_polar
    from sympy.abc import l
    a = randcplx()
    arg = x*exp_polar(k*pi*I)
    expr1 = pi*meijerg([[], [(a + 1)/2]], [[a/2], [-a/2, (a + 1)/2]], arg**2/4)
    expr2 = besseli(a, arg)

    # Test that the two expressions agree for all arguments.
    for x_ in [0.5, 1.5]:
        for k_ in [0.0, 0.1, 0.3, 0.5, 0.8, 1, 5.751, 15.3]:
            assert abs((expr1 - expr2).n(subs={x: x_, k: k_})) < 1e-10
            assert abs((expr1 - expr2).n(subs={x: x_, k: -k_})) < 1e-10

    # Test continuity independently
    eps = 1e-13
    expr2 = expr1.subs(k, l)
    for x_ in [0.5, 1.5]:
        for k_ in [0.5, S(1)/3, 0.25, 0.75, S(2)/3, 1.0, 1.5]:
            assert abs((expr1 - expr2).n(
                       subs={x: x_, k: k_ + eps, l: k_ - eps})) < 1e-10
            assert abs((expr1 - expr2).n(
                       subs={x: x_, k: -k_ + eps, l: -k_ - eps})) < 1e-10

    expr = (meijerg(((0.5,), ()), ((0.5, 0, 0.5), ()), exp_polar(-I*pi)/4)
            + meijerg(((0.5,), ()), ((0.5, 0, 0.5), ()), exp_polar(I*pi)/4)) \
        /(2*sqrt(pi))
    assert (expr - pi/exp(1)).n(chop=True) == 0

def test_bessel():
    from sympy import besselj, besseli

    assert simplify(integrate(besselj(a, z) * besselj(b, z) / z, (z, 0, oo), meijerg=True, conds="none")) == 2 * sin(
        pi * (a / 2 - b / 2)
    ) / (pi * (a - b) * (a + b))
    assert simplify(integrate(besselj(a, z) * besselj(a, z) / z, (z, 0, oo), meijerg=True, conds="none")) == 1 / (2 * a)

    # TODO more orthogonality integrals

    assert simplify(
        integrate(sin(z * x) * (x ** 2 - 1) ** (-(y + S(1) / 2)), (x, 1, oo), meijerg=True, conds="none")
        * 2
        / ((z / 2) ** y * sqrt(pi) * gamma(S(1) / 2 - y))
    ) == besselj(y, z)

    # Werner Rosenheinrich
    # SOME INDEFINITE INTEGRALS OF BESSEL FUNCTIONS

    assert integrate(x * besselj(0, x), x, meijerg=True) == x * besselj(1, x)
    assert integrate(x * besseli(0, x), x, meijerg=True) == x * besseli(1, x)
    # TODO can do higher powers, but come out as high order ... should they be
    #      reduced to order 0, 1?
    assert integrate(besselj(1, x), x, meijerg=True) == -besselj(0, x)
    assert integrate(besselj(1, x) ** 2 / x, x, meijerg=True) == -(besselj(0, x) ** 2 + besselj(1, x) ** 2) / 2
    # TODO more besseli when tables are extended or recursive mellin works
    assert (
        integrate(besselj(0, x) ** 2 / x ** 2, x, meijerg=True)
        == -2 * x * besselj(0, x) ** 2
        - 2 * x * besselj(1, x) ** 2
        + 2 * besselj(0, x) * besselj(1, x)
        - besselj(0, x) ** 2 / x
    )
    assert integrate(besselj(0, x) * besselj(1, x), x, meijerg=True) == -besselj(0, x) ** 2 / 2
    assert integrate(x ** 2 * besselj(0, x) * besselj(1, x), x, meijerg=True) == x ** 2 * besselj(1, x) ** 2 / 2
    assert integrate(besselj(0, x) * besselj(1, x) / x, x, meijerg=True) == (
        x * besselj(0, x) ** 2 + x * besselj(1, x) ** 2 - besselj(0, x) * besselj(1, x)
    )
    # TODO how does besselj(0, a*x)*besselj(0, b*x) work?
    # TODO how does besselj(0, x)**2*besselj(1, x)**2 work?
    # TODO sin(x)*besselj(0, x) etc come out a mess
    # TODO can x*log(x)*besselj(0, x) be done?
    # TODO how does besselj(1, x)*besselj(0, x+a) work?
    # TODO more indefinite integrals when struve functions etc are implemented

    # test a substitution
    assert integrate(besselj(1, x ** 2) * x, x, meijerg=True) == -besselj(0, x ** 2) / 2

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_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_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_hyperexpand_bases():
    assert hyperexpand(hyper([2], [a], z)) == \
  a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)*lowergamma(a - 1, z) - 1
    # TODO [a+1, a-S.Half], [2*a]
    assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
    assert hyperexpand(hyper([S.Half, 2], [S(3)/2], z)) == \
      -1/(2*z - 2) + log((z**(S(1)/2) + 1)/(-z**(S(1)/2) + 1))/(4*z**(S(1)/2))
    assert hyperexpand(hyper([S(1)/2, S(1)/2], [S(5)/2], z)) == \
               (-3*z + 3)/(4*z*(-z + 1)**(S(1)/2)) \
               + (6*z - 3)*asin(z**(S(1)/2))/(4*z**(S(3)/2))
    assert hyperexpand(hyper([1, 2], [S(3)/2], z)) == -1/(2*z - 2) \
            - asin(z**(S(1)/2))/(z**(S(1)/2)*(2*z - 2)*(-z + 1)**(S(1)/2))
    assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == \
             z**(S(1)/2)*(6*z/7 - S(6)/5)*atanh(z**(S(1)/2)) \
           + (-30*z**2 + 32*z - 6)/(35*z) - 6*log(-z + 1)/(35*z**2)
    assert hyperexpand(hyper([1+S.Half, 1, 1], [2, 2], z)) == \
           -4*log((-z + 1)**(S(1)/2)/2 + S(1)/2)/z
    # TODO hyperexpand(hyper([a], [2*a + 1], z))
    # TODO [S.Half, a], [S(3)/2, a+1]
    assert hyperexpand(hyper([2], [b, 1], z)) == \
             z**(-b/2 + S(1)/2)*besseli(b - 1, 2*z**(S(1)/2))*gamma(b) \
           + z**(-b/2 + 1)*besseli(b, 2*z**(S(1)/2))*gamma(b)

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_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_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 hopkins(u):
    return sm.besseli(1,x) * sm.exp(-(x/2)**2 / L) 

def test_mellin_transform_bessel():
    from sympy import Max, Min, hyper, meijerg
    MT = mellin_transform

    # 8.4.19
    assert MT(besselj(a, 2*sqrt(x)), x, s) == \
           (gamma(a/2 + s)/gamma(a/2 - s + 1), (-re(a)/2, S(3)/4), True)
    assert MT(sin(sqrt(x))*besselj(a, sqrt(x)), x, s) == \
           (2**a*gamma(S(1)/2 - 2*s)*gamma((a+1)/2 + s) \
                / (gamma(1 - s- a/2)*gamma(1 + a - 2*s)),
            (-(re(a) + 1)/2, S(1)/4), True)
    # TODO why does this 2**(a+2)/4 not cancel?
    assert MT(cos(sqrt(x))*besselj(a, sqrt(x)), x, s) == \
           (2**(a+2)*gamma(a/2 + s)*gamma(S(1)/2 - 2*s)
                / (gamma(S(1)/2 - s - a/2)*gamma(a - 2*s + 1)) / 4,
            (-re(a)/2, S(1)/4), True)
    assert MT(besselj(a, sqrt(x))**2, x, s) == \
           (gamma(a + s)*gamma(S(1)/2 - s)
                / (sqrt(pi)*gamma(1 - s)*gamma(1 + a - s)),
            (-re(a), S(1)/2), True)
    assert MT(besselj(a, sqrt(x))*besselj(-a, sqrt(x)), x, s) == \
           (gamma(s)*gamma(S(1)/2 - s)
                / (sqrt(pi)*gamma(1 - a - s)*gamma(1 + a - s)),
            (0, S(1)/2), True)
    # NOTE: prudnikov gives the strip below as (1/2 - re(a), 1). As far as
    #       I can see this is wrong (since besselj(z) ~ 1/sqrt(z) for z large)
    assert MT(besselj(a - 1, sqrt(x))*besselj(a, sqrt(x)), x, s) == \
           (gamma(1-s)*gamma(a + s - S(1)/2)
                / (sqrt(pi)*gamma(S(3)/2 - s)*gamma(a - s + S(1)/2)),
            (S(1)/2 - re(a), S(1)/2), True)
    assert MT(besselj(a, sqrt(x))*besselj(b, sqrt(x)), x, s) == \
           (2**(2*s)*gamma(1 - 2*s)*gamma((a+b)/2 + s)
                / (gamma(1 - s + (b-a)/2)*gamma(1 - s + (a-b)/2)
                   *gamma( 1 - s + (a+b)/2)),
            (-(re(a) + re(b))/2, S(1)/2), True)
    assert MT(besselj(a, sqrt(x))**2 + besselj(-a, sqrt(x))**2, x, s)[1:] == \
           ((Max(re(a), -re(a)), S(1)/2), True)

    # Section 8.4.20
    assert MT(bessely(a, 2*sqrt(x)), x, s) == \
           (-cos(pi*a/2 - pi*s)*gamma(s - a/2)*gamma(s + a/2)/pi,
            (Max(-re(a)/2, re(a)/2), S(3)/4), True)
    assert MT(sin(sqrt(x))*bessely(a, sqrt(x)), x, s) == \
           (-2**(2*s)*sin(pi*a/2 - pi*s)*gamma(S(1)/2 - 2*s)
                * gamma((1-a)/2 + s)*gamma((1+a)/2 + s)
                / (sqrt(pi)*gamma(1 - s - a/2)*gamma(1 - s + a/2)),
            (Max(-(re(a) + 1)/2, (re(a) - 1)/2), S(1)/4), True)
    assert MT(cos(sqrt(x))*bessely(a, sqrt(x)), x, s) == \
           (-2**(2*s)*cos(pi*a/2 - pi*s)*gamma(s - a/2)*gamma(s + a/2)*gamma(S(1)/2 - 2*s)
                / (sqrt(pi)*gamma(S(1)/2 - s - a/2)*gamma(S(1)/2 - s + a/2)),
            (Max(-re(a)/2, re(a)/2), S(1)/4), True)
    assert MT(besselj(a, sqrt(x))*bessely(a, sqrt(x)), x, s) == \
           (-cos(pi*s)*gamma(s)*gamma(a + s)*gamma(S(1)/2 - s)
                / (pi**S('3/2')*gamma(1 + a - s)),
            (Max(-re(a), 0), S(1)/2), True)
    assert MT(besselj(a, sqrt(x))*bessely(b, sqrt(x)), x, s) == \
           (-2**(2*s)*cos(pi*a/2 - pi*b/2 + pi*s)*gamma(1 - 2*s)
                * gamma(a/2 - b/2 + s)*gamma(a/2 + b/2 + s)
                / (pi*gamma(a/2 - b/2 - s + 1)*gamma(a/2 + b/2 - s + 1)),
            (Max((-re(a) + re(b))/2, (-re(a) - re(b))/2), S(1)/2), True)
    # NOTE bessely(a, sqrt(x))**2 and bessely(a, sqrt(x))*bessely(b, sqrt(x))
    # are a mess (no matter what way you look at it ...)
    assert MT(bessely(a, sqrt(x))**2, x, s)[1:] == \
            ((Max(-re(a), 0, re(a)), S(1)/2), True)

    # Section 8.4.22
    # TODO we can't do any of these (delicate cancellation)

    # Section 8.4.23
    assert MT(besselk(a, 2*sqrt(x)), x, s) == \
           (gamma(s - a/2)*gamma(s + a/2)/2, (Max(-re(a)/2, re(a)/2), oo), True)
    assert MT(besselj(a, 2*sqrt(2*sqrt(x)))*besselk(a, 2*sqrt(2*sqrt(x))), x, s) == \
           (4**(-s)*gamma(2*s)*gamma(a/2 + s)/gamma(a/2 - s + 1)/2,
            (Max(-re(a)/2, 0), oo), True)
    # TODO bessely(a, x)*besselk(a, x) is a mess
    assert MT(besseli(a, sqrt(x))*besselk(a, sqrt(x)), x, s) == \
           (gamma(s)*gamma(a + s)*gamma(-s + S(1)/2)/(2*sqrt(pi)*gamma(a - s + 1)),
            (Max(-re(a), 0), S(1)/2), True)
    assert MT(besseli(b, sqrt(x))*besselk(a, sqrt(x)), x, s) == \
           (2**(2*s - 1)*gamma(-2*s + 1)*gamma(-a/2 + b/2 + s)*gamma(a/2 + b/2 + s) \
               /(gamma(-a/2 + b/2 - s + 1)*gamma(a/2 + b/2 - s + 1)),
            (Max(-re(a)/2 - re(b)/2, re(a)/2 - re(b)/2), S(1)/2), True)
    # TODO products of besselk are a mess

    # TODO this can be simplified considerably (although I have no idea how)
    mt = MT(exp(-x/2)*besselk(a, x/2), x, s)
    assert not mt[0].has(meijerg, hyper)
    assert mt[1:] == ((Max(-re(a), re(a)), oo), True)